2. 中南大学资源与安全工程学院，湖南 长沙 410083；
3. 商丘师范学院 环境与规划学院，河南 商丘476100
2. School of Resources and Safety Engineering, Central South University， Changsha 410083,China;
3. Environment and Planning School, Shangqiu Normal University，Shangqiu 476100,China
1 引 言
时空去相干和大气折射延迟是限制DInSAR地面沉降监测精度的主要因素，为了削弱大气扰动对差分干涉相位质量的影响，文献[1]提出永久散射体雷达干涉测量技术(permanent scatterers InSAR ,PS-InSAR)，较好地解决了大气扰动问题，随后，该技术得到迅速发展，但由于要求数据量较多(通常20幅以上)，且数据处理范围通常要求在5km2以内[2, 3, 4, 5, 6, 7, 8, 9, 10]，因此，该技术目前很难被广泛应用。文献[11—19]采用基于相干点目标的DInSAR技术进行削弱大气扰动和时空去相干的研究与试验分析，得到的试验区平均沉降图能有效揭示地面形变的空间展布，得到较可靠的监测结果。可见，采用相干点目标分析方法能降低大气延迟影响、提高信噪比，但这些研究没有对相干点目标影像对的质量进行合理评价，如果其中部分影像质量太差，相干性太低，必将会影响到最后的沉降监测结果；文献[20]利用6幅ENVISAT ASAR影像，应用干涉图叠加法对珠江三角洲的地面沉降进行研究，发现因过去20年的城市化发展导致的广州、佛山、东莞等地的地面沉降现象非常明显，但其研究是在非单一公共主影像的基础上进行干涉图叠加，会减弱干涉图形变的线性关系，增大干涉图叠加方法(Stacking)的假定条件要求(其一，单副雷达影像中大气延迟相位为随机分布，各差分干涉图中大气延迟相位影响也为随机分布；其二，研究区的地表形变为线性特征)[11, 17]，同时，没有考虑干涉图质量权重问题，影响了监测结果的精度。基于此，本文提出一种基于单一公用主影像且考虑干涉图质量权重的干涉图叠加方法，解决以上技术存在问题，提高监测精度，获得较好的地面沉降监测结果。
2 公用主影像干涉图加权叠加方法 2.1 基本原理及流程公用主影像干涉图加权叠加方法是选取最优公用主影像进行配准，将配准到同一主影像几何空间的差分干涉解缠相位图进行加权叠加，利用大气延迟相位的随机分布特性和地表形变信号近似线性表现的特点，提高叠加结果的信噪比。该方法假设在干涉图中，大气扰动的误差相位是随机的，而形变为线性形变[13, 18, 20]，并根据相干性来定权。基于这种假设，将多幅干涉图对应的解缠相位叠加起来，所得的形变相位信息对应着累加时间基线内加权的变形量；叠加后的大气延迟相位，却不是单幅干涉图中大气相位误差随干涉图数量倍数增长的结果，而只是干涉图数量的平方根倍增长的结果。这样，叠加相位图中形变信息和大气延迟项之间的信噪比就得到提高。
公用主影像干涉图叠加方法数据处理的基本流程包括：主影像的优化选取[21, 22, 23]、辐射定标等数据预处理；影像配准与重采样到主影像空间；高相干区域的选取；干涉图生成；去平地、地形生成差分干涉图；相位解缠；根据相干图质量定权，干涉图解缠相位叠加；地理编码、形变图生成等。具体流程见图 1。
|
| 图 1 公用主影像干涉图加权叠加方法处理流程 Fig. 1 Flow of weighted stacking based on common master image |
从图 1中可以看出，公用主影像干涉图加权叠加方法基本包含了常规DInSAR的所有关键步骤。但要进行有效的干涉图加权叠加，首先，基于同一公用主影像的干涉图叠加方法，应进行公用主影像的优化选取，选择大气扰动等误差影响最小的影像作为公用主影像，以便最大程度减少大气延迟和其他噪声的影响。对于公用主影像的优化选取，既要使时间基线最佳(尽量位于时间跨度中间) ，又要顾及干涉对的有效空间基线和Doppler 质心频率、大气延迟影响因子以及季节因素等对干涉图相关性的显著影响，寻找综合影响最小的影像为公用主影像[21, 22, 23]。其次，要根据干涉图相干性质量(本文根据各影像相干点的数量成正比来定权，数量最多的相干图叠加权重qji设定为1)进行定权，如果存在有与其他干涉图相差明显较大的影像(即相干图相干点的数量≤1/2最高相干图相干点数时)，将其权重设为0，不让它参加干涉图叠加计算与后续的数据处理。
2.2 干涉图加权叠加方法的函数模型所有干涉图相位解缠叠加后，高相干点总的形变相位为所有干涉图解缠形变相位的加权和，且与时间基线加权累加量成正比。则高相干点的平均形变相位变化速率Videf可表示为
式中，N为干涉图的数目; Φij与qji分别为第i个像元第j幅差分干涉图的解缠形变相位与其叠加权重;ΔTj为第j幅差分干涉图主辅影像之间的时间间隔(时间基线);λ为雷达波的波长。
由于单幅雷达影像中的大气延迟相位分布是随机的，各差分干涉图中的大气延迟相位影响也随机分布，根据误差传播定律，则所有干涉图叠加后，高相干点受到的大气延迟影响可表示为
式中，Φji为第j幅差分干涉像对对应的大气相位延迟；ΔΦi为i像元在N个干涉图叠加后的大气延迟影响总和。
则干涉图叠加后，高相干点大气延迟对线性形变速率的影响可表示为
由式(1)和式(3)可得，干涉图叠加像元中点目标i的平均速率
可表示为
从式(4)可以看出，叠加后平均形变速率中，形变速率的信噪比得到了提高。因此，干涉图叠加像元中平均形变速率的标准偏差为
式中，σVi为像元点i的平均形变速率标准偏差；Vji为像元点i在第j幅差分干涉图中的相位变化速率。这个相位变化速率的标准偏差可以反映出该研究区域地表形变的线性变化特征，标准偏差越小，形变线性特征明显，监测精度越高。
干涉图加权叠加最有效的方法是对相干性较高的点相位进行叠加，在选择高相干像元时，首先计算各雷达影像对应的相干系数，然后可用函数模型(6)[4, 13]进行高相干目标点的提取
式中，m、n表示选择的移动窗口方位向、距离向的像元数；γi为某像元i在移动窗口中相干系数的均值，是通过包含该像元的移动窗口的所有像素复数信息来进行估算的，当γi高于给定阈值γc的点选定为干涉图相位解缠的高相干目标点，将所有满足该模型的点组合成高质量相干点集。 3 公用主影像加权叠加方法在沉降监测中的应用 3.1 研究区域与数据来源
自20世纪初以来，天津市区及近郊一直存在不同程度的沉降，某些地段沉降量已经超过3m，但1986年开始对过度开采地下水采取控沉措施后，市区地段沉降速度明显减小，可市区以外的近郊及县区，某些地段沉降不但没有减弱，相反沉降量增大，沉降速率有进一步加大的趋势[24]。
考虑到公用主影像干涉图加权叠加方法的优越性，选择天津市主城区及近郊作为试验区，其中心经纬度为(39.11288°N，117.06963°E)，面积约24km×28km，利用2003—2006年间获取的11景ENVISAT ASAR数据(Track2447，Frame22817)进行试验。
首先从这11幅SAR影像中择优选择1幅SAR影像作为公用主影像，可利用综合函数模型[3, 4]式(7)来进行确定
式中，Tc、B⊥c、FDCc分别表示时间基线、空间垂直基线和多普勒质心频率差异的极限基线；Tmk、B⊥mk、FDCmk分别表示主影像m与第k幅辅影像对应的时间基线、空间垂直基线和多普勒质心频率差异；α、β、ω为对应的指数因子。建立该模型的目的是通过求取函数综合相关系数的最大值来确定公用主影像，即γ均m为最大值时，对应的主影像为最优公用主影像。根据式(7)，成功选取2004-12-10的14529影像作为主影像，其他影像作为辅影像并重采样到主影像空间，组成10个干涉像对，其相关参数见表 1。
| 主-辅 影像 | 辅影像 成像时间 | 时间 间隔/d | 垂直 基线/m |
| 1 | 14529-10521 | 2004-03-05 | -8×35 | -70.99 |
| 2 | 14529-11523 | 2004-05-14 | -6×35 | -351.17 |
| 3 | 14529-12024 | 2004-06-18 | -5×35 | 94.67 |
| 4 | 14529-12525 | 2004-07-23 | -4×35 | 360.23 |
| 5 | 14529-13026 | 2004-08-27 | -3×35 | 112.76 |
| 6 | 14529-16032 | 2005-03-25 | 3×35 | -678.37 |
| 7 | 14529-17034 | 2005-06-03 | 5×35 | -214.71 |
| 8 | 14529-18537 | 2005-09-16 | 8×35 | 95.55 |
| 9 | 14529-20040 | 2005-12-30 | 11×35 | 166.96 |
| 10 | 14529-21042 | 2006-03-10 | 13×35 | 49.48 |
在雷达差分干涉测量中，通常利用美国航天飞机测图任务SRTM 3″分辨率(90m)DEM来去除地形相位。天津地处华北平原上，地势比较平坦，地形起伏在10m以内，具体见图 2，图中白色方框为所选雷达影像覆盖范围。该试验区地形的相位贡献不大[11]。
|
| 图 2 天津地表高程变化图 Fig. 2 Tendency chart of topographic elevation in Tianjin area |
根据DInSAR的基本原理，以14529作为干涉公用主影像，将其他10幅辅影像与14529影像进行粗配准和精配准，并重采样到主影像14529几何空间。根据式(6)高相干点的选取模型，选定0.2作为相干系数阈值，则可得各相干图见图 3。从图 3中可以看出，14529-10521、14529-12024、14529-13026、14529-17034、14529-18537和14529-21042相干性较好，相干点多且相干系数高，相干性基本一致，它们对应的有效空间基线较短，几何去相干较小；而14529-16032相干性最差，它对应的有效空间基线最大，为-678.37m，是其他有效空间基线的近2倍以上，造成的几何去相干太严重，高相干点不到其他相干图的一半。如果将其直接与其他相干性好的相干图进行叠加解算，必将会影响最终的监测精度。因此，合理设置相位叠加时各影像的权重尤为重要，权重的设置，可根据相干图 3中各影像相干点的数量成正比来定权，数量最多的相干图叠加权重qij设定为1，当相干点的数量明显低于其他相干图时(经试验证明，相干图相干点的数量≤1/2最高相干图相干点数，即qij≤1/2时)，需将此相干图叠加权重qij定权为零，可见14529-16032的叠加权重应设为0，否则会降低最终沉降监测结果。
|
| 图 3 相干图集 Fig. 3 Coherent images |
再将配准的干涉图与SRTM3 DEM进行差分干涉，得到对应的差分干涉图，见图 4。
|
| 图 4 差分干涉图集 Fig. 4 Differential interferograms |
从图 4中可以明显地看出，14529-16032的差分干涉图具有图明显的误差。主要是由于其干涉图相干点分布稀少，解缠相位不连续，很难形成正确的干涉影像。
高相干目标在空间分布上是离散的、间断的，采用Delaunay法则能将所有离散相干点用若干个没有重叠的三角形连接起来[20]，保证每一个点至少是一个三角形的顶点，假定每一个三角形的重心是一个相位解缠的节点，沿该三角形3条边做3个顶点的相位闭合线积分，如果任意两顶点间相位差的绝对值不大于π，则其积分值为0，若积分值为-2π或+2π，则表示该节点有一个负留数或正留数，计算每一个三角形对应的留数，为了防止解缠误差扩散，将正负留数用弧段成对地连接起来，在解缠时，避免积分路径穿过这些连有正负留数的弧段[18]，再利用最小费用流法[6]来进行相位解缠，可得各解缠相位图 5。根据前面定权的理论，计算各影像的权重，其中14529-16032干涉对解缠相位的权重为0。通过以上10幅相位解缠图叠加解算后，得到的年平均沉降速率见图 6，从图 6中可以看出，天津市中心城区的沉降趋势基本一致，大部分沉降速率在2cm/a左右，但北辰区和西南角存在明显的沉降中心，其最大沉降速率达8cm/a，可见天津郊区沉降速率要比市中心城区大，其主要原因是2003—2006年间天津工业向郊区发展和转移，导致郊区地下水的抽取严重，地面下沉加速，与官方监测[24](水准测量结果)一致。
|
| 图 5 解缠相位图集 Fig. 5 Unwrapping phase diagrams |
|
| 图 6 平均形变速率 Fig. 6 Mean rate of deformation |
根据式(5)，对解算的平均形变速率计算标准偏差，可得如下形变速率标准偏差图 7。由图 7可知，市区线性形变标准偏差较小且比较均匀，最大不超过0.5cm/a，表明该地区在2003—2006年内地表沉降速率的线性特征明显，监测结果可靠。
|
| 图 7 形变速率标准差 Fig. 7 Standard deviation of deformation rate |
天津沉降办公室提供了47个水准监测点2003年、2004年和2005年数据[24]，对获得的这些离散监测水准数据，求取2003—2006年度内平均沉降速率，基于Matlab软件，采用Kriging方法进行内插计算，其内插结果见图 8。
|
| 图 8 平均沉降速率Kriging内插结果 Fig. 8 Mean rate of subsidence by Kriging interpolation |
从图 8中，可以看出北辰区的水准监测内插结果与雷达干涉图加权叠加监测结果一致，其最大沉降速率都为8cm/a，沉降漏斗明显。表明公用主影像干涉图加权叠加方法与高精度水准监测精度相当。
3.4 与DInSAR及等权Stacking方法的对比根据DInSAR的基本原理与流程，对以上11景SAR影像数据进行优化组合试验，得到该地区最优监测结果见图 9。与公用主影像加权叠加方法图 6相比，DInSAR处理结果明显存在噪声的影响，且与水准监测数据最大沉降速率8cm/a相比，存在2 cm/a的差异。同时，进行了随机主影像等权Stacking方法，得到了最后形变监测结果如图 10所示，与本文公用主影像干涉图加权叠加方法监测结果图 6相比，其连续性较差且偏差稍大，且其沉降漏斗北辰区宜兴埠的最大形变监测值为7cm/a，与水准监测结果和新加权Stacking方法的监测结果有1cm/a的差异。基于以上试验分析，初步推论公用主影像加权叠加方法具有更优的地表形变监测能力。
|
| 图 9 DInSAR监测结果 Fig. 9 Monitoring results of conventional DInSAR |
|
| 图 10 传统Stacking方法监测结果 Fig. 10 Monitoring results of conventional Stacking |
在现有Stacking方法的基础上，提出一种削弱大气延迟影响的公用主影像干涉图加权叠加方法。该方法主要贡献：① 考虑时间基线、空间垂直基线、多普勒质心频率差异的综合影响，选取最优公用主影像，消弱主影像大气延迟影响；② 顾及相干图质量(即高相干点数量)进行定权，解决低相干影像图对监测结果的负面影响，通过加权叠加后，相位图的形变信息和大气噪声之间的信噪比得到明显提高。该方法要求干涉影像数不少于3幅，且研究区域的形变为缓慢形变，它能在数据量较少的情况下(与PS方法相比)，得到较好的监测结果。通过对天津市区2003—2006年11景ENVISAT ASAR影像进行试验，经比较分析，监测结果优于传统DInSAR和Stacking方法,其精度和高精度水准测量相当。可见，公用主影像干涉图加权叠加方法解决了PS-InSAR数据量要求较多的问题，削除传统Stacking方法中低相干影像对监测结果的影响，提高信噪比。
| [1] | FERRETTI A,PRATI C,ROEEA F. Permanent Scatterers InSAR Interefromety[J].IEEE Transactions on Geoscience and Remote Sensing, 2001,39(l):8-19. |
| [2] | COLESANTI C, FERRETTI A. SAR Monitoring of Progressive and Seasonal Ground Deformation Using the Permanent Scatterers Technique[J]. IEEE Transactions on Geoscience and Remote Sensing,2003,41(7): 1685-1701. |
| [3] | KAMPES B M. Displacement Parameter Estimation Using Permanent Scatterer Interferometry[D]. Delft:Delft University of Technology,2005. |
| [4] | HOOPER. Persistent Scatterer Radar Interferometry for Crustal Deformation Studies and Modeling of Volcanic Deformation[D]. Stanford :Stanford University,2006. |
| [5] | PERISSIN D, ROCCA F. High-Accuracy Urban DEM Using Permanent Scatterers[J]. IEEE Transactions on Geoscience and Remote Sensing, 2006:44(11): 3338-3347. |
| [6] | KAMPES B M. Radar Interferometry Persistent Scatterer Technique [M]. Berlin:Springer, 2006. |
| [7] | FERRETTI A, SAVIO. Submillimeter Accuracy of InSAR Time Series: Experimental Validation[J]. IEEE Transactions on Geoscience and Remote Sensing,2007, 45(5):1142-1153. |
| [8] | JUNG H, KIMB S, JUNGA H, et al. Satellite Observation of Coal Mining Subsidence by Persistent Scatterer Analysis[J]. Engineering Geology,2007: 92:1-13. |
| [9] | BELL J, AMELUNG F, FERRETTI A, et al. Permanent Scatterer InSAR Reveals Seasonal and Long-term Aquifer-system Response to Ground Water Pumping and Artificial Recharge[J]. Water Resources Research, 2008,44(2):402-407. |
| [10] | LONG Sichun,LI Tao,LIU Jingnan. Application Study of PS-DInSAR Technique Fusing Multi-metadata in Urban Ground Deformation Survey[C]//Proceedings of International Conference on Environmental Science and Information Application Technology.[S.l.]:ESIAT,2009: 326-330. |
| [11] | LONG Sichun. Advanced DInSAR Techniques Based on Coherent Targets and Its Application in Ground Subsidence Monitoring[D].Wuhan:Wuhan University , 2009. (龙四春.基于相干目标的DInSAR高级技术及其在地面沉降监测中的应用[D]武汉:武汉大学,2009.) |
| [12] | STROZZI T, WEGMVLLER U, WERNER C, et al. Measurement of Slow Uniform Surface Displacement with mm/year Accuracy [C]//Proceedings of IGARSS. Hawaii: IEEE,2000. |
| [13] | WERNER C, WEGMULLER U, WIESMANN A, et al. Interferometric Point Target Analysis with JERS-1 L-band SAR Data[C]//Proceedings of IGARSS.Muri: IEEE,2003. |
| [14] | BERARDINO P, CASU F, FOMARO G, et al. Small Baseline DInSAR Techniques for Earth Surface Deformation Analysis[C]//Proceedings of FRING 2003 Workshop.Franscati:[s.n.],2003. |
| [15] | BERARDINO P, CASU F. A Quantitative Analysis of the SBAS Algorithm Performance. Geoscience and Remote Sensing Symposium[C]//Proceedings of IGARSS.[S.l.]: IEEE , 2004. |
| [16] | WEGMVLLER U. GAMMA IPTA Processing Example Luxemburg,GAMMA Technical Report[R].Luxeburg:GAMMA,2005. |
| [17] | CASU F, MANZO M, LANARI R. A Quantitative Assessment of the SBAS Algorithm Performance for Surface Deformation Retrieval from DInSAR Data[J]. Remote Sensing of Environment,2006,102:195-210. |
| [18] | GE Daqing, WANG Yan, GUO Xiaofang,et al. Surface Deformation Monitoring with Multi-Baseline D-InSAR Based on Coherent Point Target[J]. Journal of Remote Sensing, 2007,11(4):574-580.(葛大庆,王艳,郭小方,等.基于相干点目标的多基线D-InSAR技术与地表形变监测[J].遥感学报,2007,11(4):574-580.) |
| [19] | ZHANG Yonghong,ZHANG Jixian,GONG Wenyu. Monitoring Urban Subsidence Based on SAR Interferometric Point Target Analysis[J]. Acta Geodaetica et Cartographica Sinica,2009,38(6):482-493.(张永红,张继贤,龚文瑜.基于SAR 干涉点目标分析技术的城市地表形变监测[J].测绘学报,2009,38(6):482-493.) |
| [20] | ZHAO Qing, LIN Hui, JIANG Liming. Ground Deformation Monitoring in Pearl River Delta Region with Stacking D-InSAR Technique[J]. Geoinformatics and Joint Conference on GIS and Built Enviroment, 2008,7145:1-9. |
| [21] | LONG Sichun,LIU Jingnan,LI Tao,et al.Method for Optimum Selection of Common Master Acquisition for PS-DInSAR Fusing GPS Data[J]. Journal of Tongji University (Natural Science),2010,38(3): 453-458. (龙四春,刘经南,李陶,等.融合GPS数据的PS-DInSAR公用主影像的优化选取[J].同济大学学报:自然科学版,2010,38(3):453-458.) |
| [22] | CHEN Qiang, DING Xiaoli, LIU Guoxiang. Optimum Selection of Common Master Acquisition for PS-DInSAR[J]. Acta Geodaetica et Cartographica Sinica, 2007,36 (4):395-399.(陈强,丁晓利,刘国祥. PS-DInSAR公共主影像的优化选取[J].测绘学报,2007年,36(4):395-399.) |
| [23] | ZHANG Hua,ZENG Qiming,LIU Yihua, et al. The Optimum Selection of Common Master Image for Series of Differential SAR Processing to Estimate Long and Slow Ground Deformation [C]//Proceedings of IGARSS. Seoul:IEEE, 2005:4586-4589. |
| [24] | Tianjin Office of Ground Subsidence Monitoring. Annual Report on Tianjin Ground Subsidence[R].Tianjin: Tianjin Water Resources Bureau,2006. (天津市控制地面沉降工作办公室.天津市地面沉降年报[R]. 天津:天津市水利局,2006.) |
| [25] | YU Xiaoxin, YANG Honglei,PENG Junhuan.A Modified Goldstein Algorithm for InSAR Interferogram Filtering[J].Geomatics and Information Science of Wuhan University,2011,36(9):1051-1054.(于晓歆,杨红磊,彭军还.一种改进的Goldstein InSAR干涉图滤波算法[J]. 武汉大学学报:信息科学版,2011,36(9):1051-1054.) |
| [26] | XIA Y, KAUFMANN H,GUO X F. Differential SAR Interferometry Using Corner Reflectors[C]//Proceedings of IGARSS.Toronto:IEEE, 2002:1243-1246. |
| [27] | XIA Y, KAUFMANN H, GUO X F. Landslide Monitoring in the Three Georges Area Using D-InSAR and Corner Reflectors[J]. Photogrammetric Engineering and Remote Sensing,2004,(70)10:1167-1172. |
| [28] | Xia Y,Bam. Earthquake: Surface Deformation Measurement Using Radar Interferometry[J]. Acta Seismologica Sinica, 2005,(18)4:451-459. |