2011, Vol. 54
2. 香港理工大学 土地测量与地理资讯,香港 九龙
, DING Xiaoi-Li1,2, ZHU Jian-Jun1, WANG Chang-Cheng1, DING Wei1, YANG Ya-Fu1, WANG Yong-Zhe1
2. Dept. of Land Surveying & Geo-Informatics, the Hong Kong Polytechnic University, Hong Kong
永久散射体雷达干涉测量技术(PSInSAR)通过对一定数量的后向散射特性在长时间内保持稳定的点(PS点)进行相位分析，分离出各相位分量，从而提取出形变信息.由于其不受时间和空间失相干影响，近年来被学者们广泛应用于区域地表形变研究[1, 2].
PSInSAR 技术依据一定的原则将探测得出的所有PS点共同布设成PS网络，非常类似于传统大地测量中的水准控制网，网络中任意相邻两PS 点中的连线被称为一条PS基线，包含所有基线的PS网络被称为PS基线网络.通过建立相邻PS点间差分相位函数模型，可进一步求得相邻PS 点间的形变速率差及高程改正值差.在PS基线网络中，观测值为任一PS基线所连接的两相邻点间的干涉相位差，未知参数则为其对应的形变速率差Δv及高程改正值差ΔδH.在进行空间维解缠时，许多学者采用选取研究区域内一个假设稳定点(无形变)为起算点，通过一定的路径搜索原则将所有PS 点归算到该点上以实现空间维解缠[3, 4].但这种方法容易导致PS网络中闭合环出现未知参数增量之和，即形变速率差之和ΣΔv≠0和高程改正数差之和ΣΔδH≠0的情况.鉴于此，陈强等[5, 6]在2007年提出可通过对PS基线网络进行间接平差的方法实现空间维解缠，并在2009年应用此方法对香港地区形变进行监测，应用区域内的两个GPS控制点作为起算数据对PS网进行约束平差[7].但由于GPS 控制点分布的区域局限性，在对PS 基线网络进行间接平差法求解时很可能缺乏GPS数据.如果通过单纯选取稳定区域的PS点作为参考点，会受到天然PS点分布的局限，存在较大不确定性.
人工角反射器干涉测量技术(CRInSAR)通过在监测地区布设稳定可靠的CR 点，建立CR 点上的相位变化与形变量之间关系的模型以达到对离散点进行形变监测的目的.由于CR 点较强的后向散射系数使其理论上在SAR 图像上呈现亮点，且不受时间、空间失相关影响，适用于低相干区的长时间序列形变监测.CR 点安装灵活，可自由选取研究区域进行安装，且方便在影像上识别，可避免常规点目标识别时由地图坐标系投影至影像坐标系时引入的编码误差.同时，CR 点上的GPS接收机可准确测得其高程信息，避免差分干涉处理中由于外部DEM 精度限制及地理编码过程引入的误差.具有较强的优势[8, 9].鉴于此，本文提出一种将CR 点作为PS基线网络约束点进而进行约束平差的方法，实现了CRInSAR 与PSInSAR 技术的结合，可适用于研究区域内PS基线网络空间维解缠时无外部起算数据的情况，同时可增加基线网络间接平差模型的多余观测数.
2 CRInSAR与PSInSAR联合解算算法 2.1 CRInSAR技术基本原理及处理流程CRInSAR 技术进行形变监测时，需要预先在影像上提取出CR 点的行列号信息，并建立如下函数关系模型[10]：
|
(1) |
式中p表示CR 点号，m表示干涉对序号；ΔΦpm为第p个CR 点相对于参考CR 点的干涉相位差；Δkpm为整周模糊度差；vp为CR 点相对于参考点的形变速率；Tm为时间基线；Sp为系统误差参数；Δφatmm 是大气延迟引起的相位；ε 是残余相位和噪声；Δφtopop，m为地形引起的相位，可根据下式计算：
|
(2) |
其中λ 是雷达波长，Bm为干涉对的垂直基线长度，Rp为CR 点与卫星位置间距离，θ 为雷达入射角，ΔHp为CR 点与参考点的相对高程；
实际计算中，假设有N+1 景SAR 影像，有M+1个CR 点，选取其中一个CR 点为参考点，一幅SAR 影像为主影像，将其余N景影像与主影像配准并重采样，再进行干涉处理，可形成N组干涉对，计算结果为其余M个CR 点相对于参考点在时间序列上的形变值.当CR 点间距在1000 m 范围内时，大气相位影响Δφatmm 可忽略[7, 8].而λ、Bj、Ri、θ、ΔHij均为已知量，则在(1)式中未知参数仅为整周模糊度Δk、形变速率v和系统误差参数S.求解模型(1)式中未知参数的过程即为相位解缠，在这里我们可应用LAMDBA 方法实现，逐步分离出未知参数，进而得出CR 点上的最终的形变速率值v(其详细操作流程可参见文献[9~11])，但计算得出的CR 点形变速率需要严格评定精度，保证其精度的前提下才可作为后期PS基线网络空间维解缠网的约束.
2.2 PSInSAR技术基本原理与处理流程应用PSInSAR 算法进行形变解算的前期处理过程与CRInSAR 算法流程类似，同样需在时序SAR 影像中选取出一幅主影像，其余影像为从影像，将从影像分别配准并重采样到主影像上，与主影像进行干涉、去平地效应.与CRInSAR 算法不同之处是需要应用二轨法进行外部地形相位的去除，经过这些处理之后可得到一组时间序列差分干涉图.结合相干系数法与振幅离差法综合进行PS点目标的提取，从而在主影像上选取出一定数量的PS 候选点.对选取出的PS候选点进行筛选后，需要在研究区域内对PS点进行布网，将相邻PS点连接成一条PS基线，所有PS基线共同构成PS基线网络[5].
对于PS网络中任意一条PS基线(含PS点i和j)，对这两点的相位值作差，可建立如下模型[3~5]：
|
(3) |
上式中m表示干涉对序号.
|
(4) |
|
(5) |
|
(6) |
|
(7) |
分别表示第m幅干涉图中第j点相对于第i点的干涉相位差、相位整周模糊度差、沉降速率差以及高程改正差.βi，jm是相位高程转换系数，可依据(2)式计算；Tm表示时间基线；εi，jm则是随机误差，包括了大气延迟、噪声等因素引起的误差.
(3)式中，未知参数为Δki，jm，ΔδHi，j，Δvi，j，与CRInSAR 算法类似地可应用LAMDBA 方法对缠绕的差分干涉相位进行时间维相位解缠[12, 13]，求解出相邻两PS点间的形变速率差Δvi，j和高程修正值差ΔδHi，j.但此时得出的只是相邻点上的相对形变速率和高程改正结果，还需要进行空间维解缠求解出每个PS点上的形变速率值和高程改正值.
以求解形变速率值为例，在进行PS 基线网络空间维解缠时，由(6)式，可将每条PS 基线上相邻点间形变速率差Δvi，j看作观测值，绝对形变速率vi，vj为未知参数，依据测量平差中间接平差法结合最小二乘原则求解出最终未知参数的全局最优解.
此种算法的前提是有必要起算数据，否则基线网络间接平差网求解时法方程为秩亏[14].因此待求PS点中需要有已知形变速率点才可实现网络的空间维解缠.
2.3 CRInSAR与PSInSAR联合解算处理方法 2.3.1 基本原理由前两节可知，对PS 基线网络进行空间维解缠时，需要有起算数据点才可应用间接平差最小二乘法进行准确求解，因此本文将CR 点上求得的形变速率值及其通过GPS 接收机获取的高程值作为PS基线网络的约束数据，进行约束平差，最终求解出未知参数v和δH.这样可有效地将CRInSAR 与PSInSAR 两种技术结合，解决研究区域内无起算数据点的问题.
由式(6)，(7)，对于任一条基线边i，j，观测值为Δvi，j，ΔδHi，j，未知参数为vi、vj，δHi、δHj；假设部分PS基线网络分布如图 1 所示，可建立如下的间接平差观测模型[14]：
|
(8) |
|
(9) |
(8) 式中，vcr 为CR 点上的形变速率值，由2.1节中讨论可知，通过CRInSAR 处理流程可求解得出各CR 点相对于参考点的形变速率值[9~11]；(9)式中δHcr 为CR 点上高程改正值，可通过下式求解：
|
(10) |
Hcr 为CR 点上GPS 接收机测得的准确高程值，Hdem 为CR 点在外部DEM 中对应的高程值，可通过CR 点坐标信息在外部DEM 上进行提取；
由前述可知，vcr 和δHcr 为已知量，将其作为约束数据，应用最小二乘法可求解出PS网络中各PS点上的形变速率值及高程改正值.
|
图 1 PS基线网络模拟分布 Fig. 1 The simulate distribution of PS network |
由上述讨论，将CRInSAR 与PSInSAR 联合解算算法操作流程图表示如图 2.在实际操作中假设有N+1景SAR 影像，M+1个CR 点，对影像进行干涉处理，可生成N幅干涉图，结合外部DEM 数据去除N幅干涉图的地形相位后，可生成对应的N幅差分干涉图[15].根据2.1 节中介绍的CRInSAR算法，通过逐一单独提取干涉图(含地形相位)上各CR 点的相位信息，建立模型，应用LAMDBA 法进行相位解缠，可求解出各CR 点上的线性形变速率值vcr，同时计算出每个CR点上的高程改正值δHcr，以此作为后面PS网络空间维解缠的约束数据.
|
图 2 CRInSAR与PSInSAR联合解算算法流程图 Fig. 2 Flowchart of PSInSAR and CRInSAR federative calculation algorithm |
综合考虑相干系数和振幅离差指数信息，通过设定阈值可在主影像上选取一定数目的PS 候选点，并储存相应PS点位置信息.处理中，我们将PS点与CR 点位置信息分别存储.应用Delaunay三角网原则对选取出来的PS点与CR 点一起布网，并根据邻接矩阵模型存储各基线边的信息[7]，信息包括基线边上对应两PS 点的位置信息.布网时为尽量减弱大气延迟相位的影响，将基线边控制在1000m内[8, 16].
根据2.2节中介绍的方法，对基线网中所有基线边分别建立(3)式的模型.与CRInSAR 技术不同的是此时需要在差分干涉图中提取基线边上对应两PS点的相位信息.模型中未知参数为ΔδHi，j，Δvi，j及整周模糊度差Δki，j.为进一步分离缠绕相位，采用LAMDBA 法进行相位解缠，在原式基础上补充虚拟观测方程，解决方程观测个数不足的问题[12, 13].通过求解，可得出任一基线边上两PS点间的ΔδHi，j，Δvi，j，以此作为空间维解缠的观测值.
我们将求解得出的CR 点上的线性形变速率值vcr 作为间接平差形变速率网的约束数据，同理将CR 点上的高程改正值δHcr 作为PS高程改正网的约束数据.根据(6)，(7)式的模型，依据最小二乘原则建立法方程，可求解出PS网中各PS点上的形变速率值及高程改正值.
3 实验数据与处理方法 3.1 实验数据为验证上述算法流程，选取河南省境内约250km2范围为研究区域，所覆盖的中心大地坐标约为(113°8′E，34°21′N).研究区域位于郑州市西南方向，登封市与禹州市交界地带，由于距离城区较远，区域内居民区多以村落为主(如图 3所示).由图可见，区域覆盖了白沙水库、券门水库两大水库，区域南部地形以山体为主，且山体上植被覆盖稀少，多为裸露山体，而平坦地区多为河流及植被覆盖区.研究区域矿产资源丰富，矿山较多.沿研究区域内高速公路安装了三组共12个人工角反射器.实验中数据选取2007 年2月到2009年12月期间14 景PALSAR 影像.影像采样间隔为距离向3.17m，方位向2.34m.影像详细参数信息见表 1.
|
图 3 研究区域 Fig. 3 The study are |
|
表 1 SAR影像基本参数 Table 1 Basic parameters of the SAR images used in the study |
实验中选取09年8月9日的影像为主影像，采用JPL/NASA 的3″空间分辨率SRTM DEM 数据用于去除地形相位.首先将其余13景影像分别配准到主影像上并重采样，与主影像进行干涉，去平地效应、滤波、二轨法去地形相位处理后可生成对应的13幅时间序列干涉图和去除了地形相位的差分图.在处理过程中，为保证影像分辨率符合PS 点提取的要求，干涉过程中均采用单视处理.图 4 为2008年6月21日获取的SAR 影像与主影像处理生成的差分干涉图，从图中可明显看出一沉降漏斗(A 区).实验中沿高速公路安装了三组共12 个人工角反射器，最初目的是对高速公路进行沉降监测，其位置分布如图 5所示.人工角反射器为铝板材质的等腰直角四面体，直角边长为1.2m.安装过程中根据实地调查获取的先验信息预先将CR01，CR05 和CR12三个角反射器点安装在相对稳定的区域，分别作为三组人工角反射器的参考点(见图 5).在实验区域内沿高速公路同样进行了水准测量.水准路线同样分为三段，每段水准路线的起算点与CR01，CR05和CR12位置重合，以保证与CRInSAR 方法中的各组参考点一致，以方便对CR 点形变结果进行验证.人工角反射器在雷达影像上的识别综合考虑了强度与相关系数信息，由于实际施测中公路沿线施工导致了CR10号角反射器意外毁坏，最终在影像上识别出其余十一个人工角反射器点，其在主影像上的表现为图 6.
|
图 4 差分干涉图（2008-06-21 〜2009-08-09) Fig. 4 Differential interferogram(2008-06-21 〜2009-08-09) |
|
图 5 研究区域内CR点位置分布及三组CR点在PS基线网络中的网络连接图 Fig. 5 The distribution of CR points over the study area and the network connections of three teams of CR points in PS network |
|
图 6 CR点在主影像上的表现 Fig. 6 The CR points on the master image |
根据2.1 节和2.3 节中介绍的算法，应用CRInSAR 算法求解CR 点上的形变速率时，将各组CR 点分别与参考点进行连接，类似于构建PS基线网络的方式构建简单的CR 点网络建模求解.提取出各组CR 点相对于参考CR 点的干涉相位差，并应用CR 点上准确高程信息计算地形相位分量Δφtopoi，j 后，根据(1)式建立模型.本文应用LAMDBA方法以实现时间维的相位解缠，最终可求得各组CR 点相对于参考点的形变速率值.由于参考点的选取实际是结合了水准数据作为先验信息，因此实际的解算过程类似于大地测量水准网中的拟稳平差处理方法.根据CR 点的地图坐标，可提取出其在外部DEM 上的高程信息，结合CR 点上卫星接收机测得的绝对高程值，即可根据(9)式计算出各CR 点上的高程改正值.为保证CR 点上计算得出的形变速率结果的精度，需要对CRInSAR 技术处理的结果进行精度评定.内部精度评定采用传统测量平差理论中的精度评定方法[14].外部精度评定则将水准实测数据结果进行比较验证.在严格保证CR 点结果精度合理的前提下，才可将求解出来的CR 点上的形变速率及高程改正值结果作为PS基线网络空间维解缠的约束数据.
在进行PS点选取时，综合考虑相干系数，振幅离差指数和振幅信息[6]，选取在时间序列上相干系数均值高于0.5的点为初选点.由于选取出来的点多成聚堆分布，在此基础上应用振幅离差指数法精选，选取出振幅离差指数低于0.40 的点.经过两重筛选后所得的点仍有部分分布在水体中，因此采用振幅阈值法进一步筛选，并通过人工比对，删除一些误判点.通过这样的处理后，提取出研究区域内除CR 点外共13393 个PS 候选点.通过在Google地图上比对发现探测出的PS点多为山体上裸露的岩石、城区边界的建筑物、农村的石质房屋及公路附近的金属栅栏.将这些PS点与CR点共同依据Delaunay三角网原则布网，将三组CR 点与PS点之间的网络连接情况放大显示在图 5 右侧，根据邻接矩阵模型识别出各基线边[7]，并存储基线边上相应PS 点的位置信息.根据2.3节中操作流程，对PS候选点与CR 点构建的网络应用联合解算方法进行建模，依次进行时间维相位解缠、PS 基线网络空间约束平差，求得各PS 点上形变速率值及高程改正值的初始值.
为保证选取出来的PS点精度，在求出PS点参数初始值后需要进行残差迭代删除精度较差的PS点.具体操作时需要求出所有候选PS 点的相位残差wi，jm，删除残差值大于0.8的PS点后，重复前面的步骤，直到所有PS 候选点的相位残差值均小于0.8.保证了所有PS 点的残差值小于0.8 后，再将这些筛选出来的点进行布网.将迭代筛选后确定的PS点与CR 点共同布设好PS 基线网络后，需要计算所有基线边的时序相关系数均值，为确保网络中基线边的精度，我们删除了其中相关系数值小于0.3的基线边，这一过程也是迭代过程[17]，将最终确定出来的基线边进行存储，删除孤立的PS点后，最终在研究区域内筛选得出除CR 点外共11731 个PS点.此时对这些最终确定出来的精确的PS 点进行时间维解缠、基线网络布设及空间维解缠相关步骤后，可求得PS点上最终的参数值.
4 实验结果分析与验证 4.1 实验结果分析图 7为通过CRInSAR 算法获取的CR 点上年形变量(形变速率)结果，由图可见，实际计算得出的CR 点基本保持稳定，但由于高速公路附近矿区的持续开采，导致角反射点所在位置均发生沉降，其中CR04号点的年形变量达到了厘米级.实验过程中水准测量路线同样分为三段，各段路线的起始点均选取为三组角反射器的参考点位置(CR01，CR05和CR12)，各段路线中，有四个水准测点位置与CR04，CR06，CR07和CR11四个点位非常接近，因此本文将这四个水准测点选取为实验结果的验证点.水准实际施测时间分别为2008年12月22日，2009年6月30日和2009年11月9日，其中在2008年12月22日和2009年11月9日分别为SAR 影像获取时间(见表 1)，因此我们将这一时间段内的四个水准测点的形变结果作为CR 点上的外部验证数据，水准结果同样表示于图 7 中.从图中可以看出CR 点上计算所得的形变结果与水准测量的结果吻合较好.为进一步验证CR 点形变结果的精度，应用平差理论中的精度评定方法对12 个CR 点一年内的形变量进行了内部精度评定.由表 2可知，CR 点在一年内形变量的中误差均值(表中用MSE：MeanSquareError表示)为0.02mm，而与水准数据接近的四个CR 点的外部均方根误差(表中用RMS：RootMeanSquare表示)仅为±1.2 mm，可以将其作为后期PS网络的约束数据.
|
图 7 通过CRInSAR算法获取的CR点年形变量 Fig. 7 Annual deformation on CR points within CRInSAR algorithm |
|
|
表 2 CR点形变结果误差（mm) Table 2 Deformation error on CR points (mm) |
将CR 点上获取的参数值作为约束数据，依据2.3节中介绍的联合解算算法对整个研究区域内进行PS基线网络的空间维解缠.为进一步突出本文提出的联合解算算法的优势，我们同时对研究区域内进行了传统PSInSAR 算法实验，即在空间维解缠中不考虑CR 点信息，选取一个假设稳定点作为参考点得以实现.具体操作中，结合对研究区域进行的DInSAR 实验结果中稳定区域分布情况作为先验信息，选取右上角处C 点为参考点(图 8)，作为整个PS基线网络中空间维解缠的起算数据点.空间维解缠过程中同样采用间接平差方法以避免出现基线网络闭合环中参数增量和为零的情况(详见引言).将两种方法所获取的研究区域内PS点沉降速率结果表示在图 8(地图坐标系，以外部地形图为底图)中，从图中可以看出，两种算法的结果中，沉降明显区域分布位置总体类似，尤其是在A1 区内均有一明显沉降漏斗发育.但由于参考点选取的不确定性及参考点本身不稳定，导致两种算法的结果仍存在较大差异：(1)从整体效果来看，本文获取的沉降速率场总体表现为沉降，尤其是在南部约一半面积区域内颜色以红色和黄色为主，而传统方法获取的沉降速率场总体为抬升；(2)东南角处的沉降漏斗非常明显，沉降速率值达4cm/yr，而传统方法获取的结果中，东南角位置有一沉降漏斗发育，但沉降速率仅约为2cm/yr;(3)人工角反射器分布的高速公路附近区域沉降区势差异也较大(放大效果图显示于图 9中).根据本文结果，高速率公路上的点基本保持稳定，仅仅CR04号点所在位置存在约1cm 的年累积沉降，但在公路两侧由于有矿区分布导致2cm左右的年累积沉降(两侧黄色点分布位置)，由于处理过程中我们没有针对大气相位进行去除，只是在PS基线网络布设时限制了基线边距离，尽量减弱其影响，但由于大气延迟残留相位依然可将误差引入到最终形变速率结果中，因此在图 9a中高速公路西南方向，研究区域的西北角处可看到个别PS 点存在轻微抬升.而传统方法获取的结果，在高速公路上发生成片区域的明显抬升，仅仅在西北角位置处有轻微沉降.
|
图 8 本文算法获取地PS点形变速率结果与传统方法结果对比(斜距向） (a)本文算法结果；（b)传统算法结果. Fig. 8 Deformation velocities on PS points through the algorithm proposed in this paper (a) compared to that of traditional algorithm (b) |
|
图 9 角反射器附近PS点分布放大图 Fig. 9 Enlarged drawing of the PS distribution near the CRpoints |
为合理显示研究区域内煤矿的沉降情况，我们对研究区域内的煤矿位置和实际沉降情况进行了资料搜集，图 10为研究区域内一些主要煤矿的位置分布图，其底图为本文算法获取的PS点沉降速率结果.由图 9和图 10可见，在一年时间间隔内，部分沉降区可在图中明显看出沉降趋势.其位置分布总体上与差分干涉图中条纹明显位置吻合(A1区处沉降漏斗对应于图 4中A区).研究区域内沉降量总体上分布在-20~5mm/yr，存在缓慢沉降，对照煤矿在地图上位置可知，沉降较为明显的区域均为矿区分布密集区域.其中，最为明显的A 区为方山镇二矿所在位置.根据本文实验结果，其年累积沉降量在斜距方向约4cm，如果认为其均为垂直向沉降量的贡献，其垂直向的年累积沉降量为约5cm;在研究区域中部位置，由于祁沟村第三煤矿、庙庄煤矿及建生新煤矿的持续矿产资源开采，导致中部地区的PS点速率颜色分布以黄色和红色为主，其年累积沉降量约为3cm;与中部和东南部的明显沉降量相比，北部区域的沉降缓慢，基本分布在毫米级(颜色以绿色和蓝色为主)，但在人工角反射器分布的高速公路两侧依然存在两个沉降漏斗，其位置刚好为告城煤矿和朝阳沟煤矿的分布位置，由于持续的矿产开采，沉降漏斗已经开始发育.
|
图 10 研究区域煤矿位置分布图（以PS点沉降速率图为底图） Fig. 10 Distribution of collieries over the study area (with the deformation velocities image on PS points as background) |
图 11为根据本文算法得出的高程改正场.由图 11可知，研究区域内高程改正值基本分布在±10m内，与外部SRTM 高程精度为10m 吻合.但由于外部DEM 数据时效性不强，且去地形处理过程中容易引入配准、噪声等引起的误差，导致部分区域如白沙水库西南角、高速公路西南段的高程改正值达20m左右.根据本文算法将PS 点沉降速率量值概率分布统计如图 12中，由图 12a可知，研究区域内PS点的沉降速率值大部分集中在-20~5 mm/a，所有PS点上的沉降速率均值为10.192mm/yr，矿区分布区域的PS点由于矿产资源的开采活动导致沉降速率值较大，最大沉降速率值达-8cm/a，图 12b为高程改正值的统计结果，高程改正值总体分布在-15~10 m，所有PS 点上的高程改正值均值为3.6m，最大高程改正值为-22m.
|
图 11 研究区域PS点高程改正结果 Fig. 11 Height correction on PS points over the study area |
|
图 12 PS点沉降速率值及高程改正值统计 Fig. 12 Statistics of deformation velocities and height corrections on PS points |
为进一步证明联合算法获取的PS点沉降速率结果的准确程度，我们仍然利用4.1 节中采用的四个水准验证点作为PS 网沉降速率的外部验证数据，对本文算法和传统算法结果分别进行精度估计.具体操作时，采用100m 内平均值法进行验证[3]，将水准点位置附近100m 范围内的所有PS点沉降值转换为垂直向形变结果并取平均，以与水准验证点位置沉降值进行比较.由于应用外部水准测量方法对水准点进行施测的时间分别为2008 年12 月22日，2009年6月30日和2009年11月9日，其中在2008年12 月22 日和2009 年11 月9 日分别为SAR 影像获取时间(见表 1)，而2009 年6 月30 日时对应的形变量可通过插值方法计算得出.因此最终用于与水准点进行比较的PS点沉降量均值的时间段分别为2008 年12 月22 日至2009 年6 月30日，2008年12月22日至2009 年11 月9 日.图 13为两种算法获取PS点沉降量与这四个水准验证点上的沉降量对比结果.在2008 年12 月22 日到2009年6月30日期间，本文算法结果中PS点沉降量均值与水准验证点位置沉降量间误差均分布在3mm范围内，而传统算法结果中SZ-4 号水准验证点位置出现抬升；而在2008 年12 月22 日到2009年11月9日期间，传统算法结果中SZ-1与SZ-4号水准验证点位置均发生了较大量的抬升.将其具体的均方误差统计结果表示在表 3 中.从表中明显看出，加入了CR 点数据约束后，研究区域内的PS 点上的沉降结果估计精度明显提高，同时纠正了一些不合理的抬升点，特别是在2008 年12 月22 日至2009年11 月9 日时间段内，由于传统算法中参考点选取的不确定性，参考点位置处存在的沉降将误差引入了研究区域内PS点沉降结果中，误差较大.而加入CR 点约束后，PS点测量结果与各水准点的差异均方误差在两时间段内均值为±2.0mm，低于CR 点上获取的形变速率的精度±1.2 mm(表 2)，证明了应用CR 点速率作为约束数据的合理性.从验证的结果上明显体现出了加入CR 点约束后的优势.
|
图 13 约束前后PSInSAR结果与水准实测结果对比 Fig. 13 Results comparison between the PSInSAR algorithm and leveling |
|
|
表 3 约束前后均方误差结果比较（mm) Table 3 RMS comparison between the two algorithms (mm) |
本文提出一种CRInSAR 与PSInSAR 联合解算算法，有效实现了两种基于高相干目标的时间序列形变监测方法的结合.应用该方法具有以下优势：(1)CR 点安装灵活，可自由选取研究区域，可适用于研究区域内无外部精密实测数据点的情况，避免了人为选取参考点的不确定性；(2)CR 点可为研究区域空间维解缠提供数据约束，同时可增加间接平差模型的冗余观测数，使求解更加稳定；(3)CR 点具强后向散射特性，在SAR 影像上方便识别，可避免坐标换算到SAR 坐标系时引入的误差，提高求解精度；(4)CR 点上卫星接收机获取的准确高程信息可用于PS高程改正网的起算数据，在缺少时效性的外部DEM 数据时，该方法获取高程信息更可靠.
本文应用该方法对河南地区一矿区分布密集区域进行了形变监测，首先应用CRInSAR 方法获取了安装在研究区域内的三组角反射器上的形变速率值及高程改正值，应用外部水准数据进行精度验证，其形变速率精度可达±1.2mm，将CR 点上形变参数作为PS网络空间维解缠的约束数据，获取了该区域最终线性形变速率场和高程改正场，应用水准验证数据对加入CR 点约束前后的PS 网速率场进行外部精度对比分析，结果表明加入CR 点约束后，PS点测量结果与各水准点的差异均方误差在两时间段内均值为±2.0mm，相比于传统方法精度有明显提高，证明了应用CR 点速率作为约束数据的合理性，可修正由于参考点选取不当而出现的不合理抬升点，更加合理地反应研究区域内矿区的真实沉降情况.通过实地考查证实本文方法得出的研究区域内主要沉降区与实际矿区分布位置吻合，进一步证明了本文方法可有效地应用于研究区域内无外部起算数据的情况，适用性更广，精度更高.由于本文实验中只针对区域线性沉降进行探讨，在后续研究工作中可进一步分离大气相位，提取出研究区域内的非线性沉降场.
致谢感谢课题组胡俊、蒋弥、许文斌、尹宏杰、朱臖等同学开展的周期性野外GPS和水准测量，为本文提供了验证数据.本研究用的ALOSPALSAR 数据由日本空间局提供(AO-430).
| [1] | Ferretti A, Prati C, Rocca F. Permanent scatterers in SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing , 2001, 39(1): 8-20. DOI:10.1109/36.898661 |
| [2] | Ferretti A, Prati C, Rocca F. Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing , 2000, 38(5): 2202-2212. DOI:10.1109/36.868878 |
| [3] | 卢丽君. 基于时序SAR影像的地表形变监测方法及其应用. 武汉:武汉大学, 2008 Lu L J. Method for land deformation detection based on time series SAR images and its application (in Chinese). Wuhan: Wu Han University, 2008 |
| [4] | 龙江平. 基于PS-InSAR技术的New Orleans 地区的变形监测的研究. 长沙:中南大学信息物理工程学院, 2008 Long J P. Monitoring ground subsidence in New Orleans with Persistent Scatterers interferometry (in Chinese). Changsha: School of Info-Physics and Geomatics Engineering, Centrol South University, 2008 |
| [5] | 陈 强. 基于永久散射体雷达差分干涉探测区域地表形变的研究. 成都:西南交通大学, 2006 Chen Q. Detecting regional ground deformation by Differential SAR Interferometry based on Permanent Scatterers (in Chinese). Chengdu: South West Jiaotong University, 2006 http://cn.bing.com/academic/profile?id=2667732963&encoded=0&v=paper_preview&mkt=zh-cn |
| [6] | 陈强, 丁晓利, 刘国祥. 永久散射体雷达差分干涉应用于区域地表沉降探测. 地球物理学报 , 2007, 50(3): 737–743. Chen Q, Ding X L, Liu G X. Radar differential interferometry based on permanent scatterers and its application to detecting regional ground subsidence. Chinese J. Geophys. (in Chinese) , 2007, 50(3): 737-743. DOI:10.1002/cjg2.1088 |
| [7] | 陈强, 丁晓利, 刘国祥. 雷达干涉PS网络的基线识别与解算方法. 地球物理学报 , 2009, 52(9): 2230–2236. Chen Q, Ding X L, Liu G X. Baseline recognition and parameter estimation of persistent-scatterer network in radar interferometry. Chinese J. Geophys. (in Chinese) , 2009, 52(9): 2230-2236. |
| [8] | Burgmann R, Hilley G E, Ferretti A, et al. Resolving vertical tectonics in the San Francisco Bay area from permanent scatter InSAR and GPS analysis. Geology , 2006, 34(3): 221-224. DOI:10.1130/G22064.1 |
| [9] | Xia Y, Kaufmann H, Guo X F. Differential SAR interferometry using Corner Reflectors. IEEE 2002 International Geoscience and Remote Sensing Symposium , 2002, 2: 1243-1246. |
| [10] | Xia Y, Kaufmann H, Guo X F. Landslide monitoring in the Three Gorges area using D-InSAR and corner reflectors. Photogrammetric Engineering and Remote Sensing , 2004, 70(10): 1167-1172. DOI:10.14358/PERS.70.10.1167 |
| [11] | Xia Y. CR-based SAR-interferometry for landslide monitoring. IEEE 2008 International Geoscience and Remote Sensing Symposium , 2008, 2: 1239-1242. |
| [12] | Kampes B M. Displacement parameter estimation using permernent scatterer interferometry. Delft University of Technology , 2005. |
| [13] | Kampes B M, Hanssen R F. Ambiguity resolution for permanent scatterer interferometry. IEEE Transactions on Geoscience and Remote Sensing , 2004, 42(11): 2446-2453. DOI:10.1109/TGRS.2004.835222 |
| [14] | 张后苏. 误差理论与测量平差基础. 长沙: 中南工业大学出版社, 1995 . Zhang H S. Error Theory and Surveying Adjustment Basis (in Chinese). Changsha: Central South University Press, 1995 . |
| [15] | 廖明生, 林珲. 雷达干涉测量学:原理与信号处理基础. 北京: 测绘出版社, 2003 . Liao M S, Lin H. Synthetic Aperture Radar Intereferometry—Principle and Signal Processing (in Chinese). Beijing: Publishing House of Surveying and Mapping, 2003 . |
| [16] | 王艳, 廖明生, 李德仁. 利用长时间序列相干目标获取地面沉降场. 地球物理学报 , 2007, 50(2): 598–604. Wang Y, Liao M S, Li D R. Subsidence velocity retrieval from long term coherent targets in radar interferometric stacks. Chinese J. Geophys. (in Chinese) , 2007, 50(2): 598-604. |
| [17] | Liu G X, Sean M. B, Ding X L, et al. Estimating spatiotemporal ground deformation with improved persistent-Scatterer Radar Interferometry. IEEE Transactions on Geoscience and Remote Sensing , 2009, 47(9): 3209-3219. DOI:10.1109/TGRS.2009.2028797 |