利用永久散射体雷达干涉技术进行太原市地面沉降监测

王 琪

(1. 山西省煤炭地质物探测绘院，山西 晋中 030600; 2. 太原市基础地理数据中心，山西 太原 030009)

一、引 言

地面沉降是指由于自然因素或人为活动引发地壳表层松散土层压缩并导致地面标高降低的地质现象。近年来，差分干涉测量(DInSAR)技术已成为地面沉降监测的重要手段，具有全天候、大范围监测的优势，理论上可以达到厘米级的精度。特别是近几年发展的时间序列InSAR技术(time series InSAR)，如永久散射体雷达干涉测量PS-InSAR方法和短基线集雷达干涉测量SBAS InSAR方法，是在传统DInSAR基础上发展起来的更具潜力的新的监测方法，其克服了常规DInSAR中的时空失相干，并减弱了大气效应等因素的影响，大大提高了监测的精度，在城市等缓慢形变区域具有极大的应用潜力和优势。

本文将覆盖太原市2003-08—2010-09年间39景Envisat ASAR影像作为试验数据，采用Andy Hooper开发的StaMPS软件的PS方法对太原市的地面沉降进行试验研究，提取该时期的沉降速率图，研究该区域地表形变特征，并对地面沉降成因进行探讨分析。

二、StaMPS永久散射体雷达干涉技术(PS-InSAR)

A Ferretti，C Prati和F Rocca等于1999年共同提出并验证了永久散射体(PS-InSAR)方法，即在一组长时间序列SAR影像中，选取保持高相干性且较为稳定的点为PS点目标，其受时间空间去相干的影响小[1]。通过分析这些点目标的相位组成，利用各相位分量的时频特性，可以分离出大气相位，从而获得较为精确的地表形变信息。这种算法比较适合有大量人工建筑的城区，而在非城镇区域或地表形变速率不规则的地区选取的可靠PS点较少。Andy Hooper于2004年提出了一种新的PS算法，该算法采用幅度离散指数和相位空间相关性相结合的方法来识别PS点，不需要先验知识，可以在任何地形区域选取PS点，极大扩展了PS技术的应用范围。本文正是基于这个方法对太原市的地面沉降进行了试验。

假设N+1幅SAR影像，选择一幅各个方面(时间基线、空间基线、多普勒质心频率差等)最优的作为主影像，生成N幅干涉条纹图，引入外部参考DEM去除地形相位的影响,生成差分干涉图。首先基于Ferretti提出的幅度离散指数DA(amplitude dispersion index)提取PS候选点

式中，σA，μA为N幅干涉条纹图振幅的标准差和均值。设置DA阈值(0.4～0.42)初步选取PS候选点，对选出的PS候选点集利用相位的空间相关性分析其相位稳定性，据此筛选最佳PS点，基本模型如下：第i幅干涉图中的第x个PS候选点的差分干涉相位可以表示为

ψx,i=W{φD,x,i+φA,x,i+ΔφS,x,i+Δφθ,x,i+φN,x,i}

(2)

(3)

为此，可以设定用以筛选PS点的相关系数γx为

(5)

式中，N为干涉图的数量；γx即为判定PS点目标的指标。但仅仅如此还不够，StaMPS结合信噪比、γx和PS概率进行循环迭代来选取最终的PS点。

(6)

三、太原市地面沉降监测试验

1． 试验区概况和试验数据

太原位于山西省太原断陷盆地的北缘，属于汾渭盆地，是我国中西部地区地面沉降最为严重的城市之一。太原市从20世纪50年代末发现地面沉降，80年代至90年代为地面沉降快速发展阶段，至2004年形成了西张和城区两个沉降区，西张、万柏林、下元、吴家堡4个沉降漏斗中心。沉降区北起上兰镇，南至刘家堡乡郝村，西抵西镇，东达榆次西河堡村；南北长约39 km，东西宽约15 km，沉降面积达548 km2。90年代以来，沉降范围逐年向盆地边缘扩展，沉降漏斗面积逐年扩大，南部有向晋中盆地延伸趋势，已影响到该地区社会经济的可持续发展[3]。

试验选取时间跨度从2003年8月17日至2010年9月19日的39景欧空局30 m分辨率的Envisat ASAR单视复数影像(SLC)。ASAR的中心工作频率为C波段(频率为5.331 GHz)，成像中心入射角为23°，影像的覆盖范围如图1所示。

图1 ASAR数据覆盖范围

2. 差分干涉处理

首先利用Doris软件对39景影像经过配准、干涉、去平等步骤进行差分干涉处理，在去除地形相位时，引入了3弧秒分辨率的SRTM DEM数据，平面精度为90 m，高程精度约为10 m。影像间的空间基线越长，对应的地形相位误差越大。以2009年7月26日为主影像，39景影像的空间基线分布最大856 m，最小40 m，平均为197 m，可以看出基线分布比较均匀，且太原城区地形平坦，从而由地形引起的误差相位将非常小。经过差分干涉获得了38对干涉图，如图2所示。

图2 差分干涉图

3. PS干涉处理

StaMPS的PS处理分为以下几个步骤：①PS候选点初步选取；②估计PS候选点相位噪声；③根据相位噪声特性选取PS点；④进一步筛选出最终PS点；⑤对PS点的缠绕相位进行DEM误差校正；⑥三维相位解缠；⑦估计空间相关相位误差；⑧估计其他空间相关误差项[4]。

步骤①中，设置DA阈值为0.42，取较大的值可以选出更多的PS候选点，但同时存在较多的伪PS点，通过后续的步骤将不断进行筛选，最终将伪PS点去除。步骤②中，γx是通过估计空间相关相位项和空间不相关地形误差来定义的(式(5))，将像素采样到50 m格网空间从而使其满足空间相关条件，然后通过迭代方法估计空间相关相位项，设置迭代收敛条件为γx的变化值小于0.005。步骤③选取PS点的γx阈值是离差指数DA及伪PS点密度的函数，通过概率统计的方法设定，在处理中设置可接受的伪PS点密度为18/km2。步骤④首先剔除那些由于地理编码误差导致的伪PS点，以及受邻近PS点的影响而被认为是PS点的伪PS点，然后估计各PS点相位噪声，删除噪声比较低的点，经过这一步筛选之后剩余36 895个点，即为最终的PS点。第⑤步对上述的PS点的缠绕相位进行空间非相关地形误差校正。第⑥步利用对上述校正后的相位进行三维解缠，解缠的参考为以地理坐标(112.530 2°，37.959 5°)为中心，200 m范围内的29个PS点。可以看出，第1、3、4、9景影像空间垂直基线比较大，致使其基线失相关严重，从解缠的相位图中也可以看出这4景影像存在较大解缠误差，因此在解缠处理时通过“drop_ifg_index”参数将这4景影像排除，重新解缠并估计空间相关相位项。这样获得的解缠相位图如图3所示。利用解缠的相位观测值，反演了地表雷达视线向平均形变速率(MLV)，如图4所示；同时得到了形变速率的标准差，如图5所示。

图3 解缠相位图

图4 地表形变速率(mm/a) 图5 形变速率标准差(mm/a)

4. 地面沉降空间分布特征

根据太原市沉降监测资料及沉降历史特征(如图6所示)，PS-InSAR方法获得的地表形变(如图5所示)在空间分布上与以往的研究成果基本一致[5-6]。将PS-InSAR监测结果投影到Google影像上并标注了太原市沉降漏斗中心的分布(如图6所示)。由图6可以看出，在2003-08—2010-09时间段内太原市地面沉降主要分布在中部和南部地区。北部平原地区沉降速率较小，多数不超过-10 mm/a，部分地区存在反弹回升现象。东西两侧太行山、吕梁山山缘洪积扇冲积平原的沉降速率也较小，大多不超过-25 mm/a。总体来看，太原市地面沉降空间分布从下元沿着汾河往南呈喇叭状展开，沉降速率由北向南逐渐增大，研究区范围内最大沉降速率位于刘家堡乡附近，达到-45 mm/a；其次为孙家寨，达到-40 mm/a。地面沉降漏斗分布由北到南依次为：西张，大部分地区已经出现反弹现象，平均反弹速率为2 mm/a；万柏林，平均沉降速率为-21.5 mm/a；下元，平均沉降速率-32.4 mm/a；吴家堡，平均沉降速率为-25.9 mm/a；小店镇，平均沉降速率为-30.9 mm/a;孙家寨，平均沉降速率为-34.8 mm/a；刘家堡，平均沉降速率为-36.9 mm/a。从图5形变速率的标准差可以看出，最大的形变速率标准差发生在小店沉降中心，最大为2 mm/a，这可能是该地区存在大量耕地，由耕地的季节性翻耕导致地面散射特性变化引起的；其他地区则更小，这说明估计的形变结果是可靠的。

图6 地面沉降速率在Google地图上的分布(mm/a)

与2000年之前相比，该时间段内太原的地面沉降无论是在空间分布还是在沉降快慢程度上都发生了变化。1956—2000年以前太原的地面沉降如图7所示，含西张、万柏林、下元、吴家堡4个沉降中心，沉降中心沿汾河呈串状分布，最大沉降中心位于吴家堡，南部的小店镇和孙家寨沉降漏斗则没有形成。在研究时间范围2003-08—2010-09内，太原地面沉降分布向南移动，出现了小店和孙家寨等新的沉降中心，发展成了一条小店至孙家寨的带状沉降区域。过去的西张沉降中心趋于稳定，部分地区出现了反弹回升现象；万柏林沉降速率大为减少，沉降面积也缩小很多；下元和吴家堡仍然为较明显的沉降中心，且两者有连成一片的趋势，但沉降速率与1989—2000年相比已经明显减缓。表1对比了太原市各沉降中心在1956—2000年与2003—2010年各时期内沉降速率的变化情况。总之，该时期地面累积沉降量由北往南逐渐减小，沉降速率逐渐增大。这说明随着2003年太原市采取关井压采及引黄入晋等一系列措施以来，原有地面沉降中心得到了一定程度的控制。

图7 太原不同历史时段累积地面沉降演化(太原市地面沉降监测报告，2004)

表1 不同时期太原地面沉降中心平均沉降速率对比 (mm/a)

注：1956—2000年沉降数据来自文献[7]。

5. 地面沉降时间演变特征

同时，在沉降中心(下元、万柏林、吴家堡、小店、孙家寨)提取PS点的地面沉降时间序列，通过简单的线性回归模型拟合出形变值随时间的变化，如图8—图11所示。小圆圈代表PS点，实线表示拟合的线性形变时间序列，两条虚线表示两个连续的PS点之间的相位差不能超过π，这是由于相位存在2π模糊度的原因。在这里所用数据为Envisat ASAR数据，波长5.6 cm，因此虚线代表±14 mm。从时间序列图可以看到，2007-04-11、2007-12-09、2008-01-13、2008-02-17这4景干涉图上PS点之间都超出了虚线，发现这4景干涉图都是在冬天季节，很有可能由于下雪覆盖的原因使得相干性下降，从而探测到的信号变得不可靠。

图8 下元沉降时间序列

图9 万柏林沉降时间序列

图10 吴家堡沉降时间序列

图11 小店镇沉降时间序列

四、结论与讨论

本文利用2003-08—2010-09的39景Envisat ASAR 数据基于StaMPS的PS-InSAR技术计算了该时期太原市的地面沉降，获得了沉降速率分布图。结果表明，StaMPS不仅在城区能够探测到较多可靠的PS点，并且在有植被的地区也能提取到PS点目标，但是在这些地区提取的形变结果并不是很可靠。通过与原有的沉降资料对比验证，证实该方法在城市地面沉降监测中的可靠性。但是，该方法仅能够探测到雷达视线向(LOS)的线性形变，未能对非线性形变进行估计，这对于StaMPS的推广应用将是不利的。由于缺少必要的数据，没有获得如(U,E,N)等方向的形变值。未来可联合GPS连续运行参考站(CORS)的数据进行解算，从而全面掌握更为精细的地表三维形变场。

参考文献：

[1] 廖明生，林晖.雷达干涉测量—原理与信号处理基础[M].北京：测绘出版社，2003.

[2] HOOPER A,SEGALL P,ZEBKER H. Persistent Scatterer Interferometric Synthetic Aperture Radar for Crustal Deformation Analysis, with Application to Volcán Alcedo,Galápagos[J].Journal of Geophysical Research,2007,112(B7).DOI:10.1029/2006JB004763.

[3] 孙自永，马腾，马军，等.太原市地层空间异质性对地面沉降分布的影响[J].岩土力学,2007,28(2):399-405.

[4] FERRETTI A,PRATI C,ROCCA F. Permanent Scatterers in SAR Interferometry[J].IEEE Transactions on Geoscience and Remote Sensing, 2001,39(1):8-21.

[5] 吴宏安，张永红，陈晓勇，等.基于小基线DInSAR的技术监测太原市2003～2009年地表形变场[J].地球物理学报，2011,54(3):673-680.

[6] 郭清海. 山西太原盆地孔隙地下水系统演化与相关环境问题成因分析[D].武汉：中国地质大学，2005.

[7] 王艳，廖明生，李德仁，等.利用长时间序列相干目标获取地面沉降场[J].地球物理学报，2007，50(2)：598-604.