基于SBAS的兰西县城区地面沉降监测研究

朱叶飞，陈火根，李向前，詹雅婷，苏一鸣

(1.江苏省地质调查研究院，南京 210018；2.地裂缝地质灾害重点实验室，南京 210018)

1 引 言

兰西县位于黑龙江省中西部，松嫩平原东部，隶属于黑龙江省绥化市。全县共辖18个乡、镇，总人口46.1万人。地方工业有雄厚的基础，形成了食品、纺织、机械、服装、医药、造纸、化工等十个行业为主体的工业格局，尤其是亚麻系列产品的开发已走到了全国的前列。但兰西县是否存地下水开采等原因导致的地面沉降尚不得而知，为了获取兰西县城区地表形变幅度和分布、采取合理措施控制地面沉降，本文采用InSAR技术对2012年兰西县城区的地面沉降的状况进行了监测。

合成孔径雷达差分干涉测量(DInSAR)是近年来发展起来的一种监测地表形变的手段[1-4]，能够对地表进行全天候、大范围、厘米甚至毫米级精度的监测。然而对于长时间缓慢地面沉降而言，由于SAR干涉像对的时间和空间基线距较大，导致影像时间、空间失相干严重[5-6]，利用传统DInSAR技术难以反演其形变序列。近年来国际上一些学者提出小基线集(SBAS)方法，并先后开展了基于SBAS的DInSAR监测地表沉降的应用和研究[7-12]，通过采用SBAS组合进行干涉测量，有效地减弱空间基线引起的失相干问题，提高了DInSAR技术形变计算的精度。SBAS方法利用奇异值分解(SVD)方法将多个小基线集联合起来求解，有效地解决不同SAR数据集之间空间基线过长造成的时间不连续问题，提高监测的时间分辨率。

本文引入SBAS方法，采用C波段的Radarsat-2数据对兰西县城区2012年2月～11月的地面沉降进行监测，获取该地区时间序列形变场，发现和定位沉降漏斗，对研究区沉降漏斗的时空状况进行分析，揭示沉降漏斗随时间演化的规律。

2 数据选取与计算方法

2.1 数据的选取

本文共收集从2012年2月～2012年11月覆盖研究区的11景Radarsat-2卫星宽模式的单视数据，各图像具体参数信息如表1所示。该地区的空间垂直基线和时间基线分布较为理想，时间基线最大间隔288天，最大的垂直基线为483m。本文要根据SBAS的原则计算出小基线集划分的情况，并估算出比较合理的小基线集组合方式(图1)，共选择了53个干涉组合，每个组合的垂直基线都小于400m。

表1 SAR图像获取信息表

图1 SBAS干涉图组合方式

2.2 基于相干点的SBAS计算方法

由于SBAS方法相干目标的提取是基于各个时期的单视数据，相干目标选择前需要对单视数据进行配准一样，本文选用2012年6月28日获取的SAR影像作为参考主影像，将所有影像以复相干系数作为测度配准至主影像，配准的精度控制在1/10个像元内。

在利用SBAS时序分析模型解算形变速率之前，首先要进行相干目标的筛选。目前相干目标像元的筛选方法主要有：基于像元强度的稳定性[13-14]、基于相位稳定性[15]、基于空间相干性[16]、点目标检测[17-18]等。像元强度稳定性检测算法是根据像元强度的稳定性来判断是否为相干点，一般要求SAR影像不少于30景。相位稳定性方法算法复杂，需采用大量模拟试验来确定阈值。结合本文数据量小的特点并考虑到算法的复杂性，采用空间相干性阈值和点目标检测的方法来选取相干目标点，这样既提高相干目标识别准确性，也增大相干目标的数量。本文共选取研究区初始相干点6365个参与运算，初选的相干点分布如图2所示。

图2 初选的相干点分布图

根据SBAS组合原则对收集到的11幅Radarsat-2图像进行组合，共选择了53个干涉对，利用研究区工作区DEM去除地形效应后得到差分干涉图。依据相干点的差分相位δφ可建立矩阵方程：

Dv=δφ

(1)

假设由N+1幅SAR图像得到M个干涉图，D是一个M×N矩阵，对第j行，位于主辅图像获取时间之间的列，D(j，k)=tk+1-tk，其他的D(j，k)=0；υ[N×1]为相干点的平均相位速度。

将SVD分解应用于矩阵方程(1)，就可以得到速度矢量v的最小范数解。地表形变最小范数意义上的最小二乘相位估计可表示为：

v=D+δφ

(2)

D+为D的伪逆矩阵：

(3)

式中ui为DDT的特征向量；vi为DTD的特征向量；λi为DDT的特征值。

从差分相位的组成出发，除了形变相位分量外，还有高程误差Δq的相位分量，建立方程组如下：

Dv+C·Δq=δφ

(4)

其中，C[MX1]是与基线距相关的系数矩阵，由此可以得到DEM误差。

由于大气相位在空间上表现为低频、在时间上表现为高频，因此对以上线性模型的残余相位在空间和时间上进行适当滤波，分离出大气相位和非线性形变相位。将非线性形变量与线性形变量相加，从而得到各个时间段累计的沉降量。

3 兰西县城区地面沉降分析

将研究区域的SAR数据采用SBAS方法进行时序分析得到兰西县城区2012年2月5日～2012年11月19日的高相干点的平均沉降速率(图3)。通过图4可以看出，2012年兰西县城区存在一定范围的地面沉降，主要呈漏斗形分布，另外在城区东侧存在局部小范围的沉降。沉降中心位于兰西县城区南侧5号点附近，最大地面沉降速率达到2.43cm/年，以此为中心地面沉降幅度向四周逐渐趋缓。笔者对研究区进行实地野外调查，发现沉降主要原因为地下水的开采导致的地面沉降，包括各类企业(如食品厂等)利用深水井进行地下水开采、开发区大量新建高层建筑基坑排水、城镇居民生活用地下水开采等，在局部地区表现为墙体开裂等沉降迹象。

图3 研究区地表形变散点图

在反演线性形变速率的基础上，通过空间滤波和时间滤波相结合的方式去除大气延迟相位的影响，反演得到非线性形变量。将线性形变量和非线性形变量结果叠加，得到研究区各时间段的累计形变量，选取散点图中沉降漏斗南北向剖面线上7个典型相干点，各相干点的地面沉降累计沉降序列如图4所示。从时间序列上看虽然不同点的沉降曲线有所差异，但兰西县城区沉降漏斗的地面沉降基本上是呈现线性的沉降。

图4 研究区典型高相干点累计沉降序列图

4 相邻轨道重复区的精度验证

图5 相邻轨道重复区地面沉降散点图

由于没有研究区同期的水准观测数据对InSAR结果进行验证，本文选取东侧相邻轨道的9景Radarsat-2数据，该幅图完全覆盖研究区范围。对相邻轨道的研究区进行SBAS方法的独立计算，并对重复研究区的InSAR计算结果进行相互对比验证。选取的9景SAR数据获取时间在2012年1月29日～2012年10月19日之间，依据SBAS组合原则共选择了28个干涉组合，计算得到的2012年线性形变速率如图5所示。在重复区域共选取研究区7个典型的同名高相干点，通过比较同名点的地表形变速率(表2)发现：两个相邻轨道独立计算的同名点2012年地表形变线性速率的相对误差都在0.4cm/年以内，误差的均方差为0.132cm/年。由于不同轨道图像获取时间的差异、轨道误差及大气延迟相位的差异，则计算的结果存在一定的差异，但上述精度验证表明相邻轨道重复区上的计算结果非常一致，表明本文的计算结果准确地反映了研究区地面沉降的幅度和分布状况。

表2 高相干同名点的相互验证(单位：cm/年)

5 结束语

本文引入SBAS方法监测兰西县城区地面沉降，突破InSAR观测周期以及空间和时间基线的限制，对DInSAR单次离散的观测进行时间序列分析，准确地获取研究区的地表形变速率图，发现了该研究区沉降漏斗中心速率达到2.43cm/年，并获得了该研究区沉降漏斗地表形变累积量的时序变化，发现该漏斗区地面沉降基本呈现线性沉降的规律。经实地调查发现研究区沉降的原因主要是地下水的过度开采。传统的以“点”为基础的研究区水准、GPS测量，由于空间采样率低，很难准确发现和定位沉降漏斗的中心、范围、大小和演化规律，而本文研究结果对预防兰西县城区地质灾害和环境破坏具有重要的意义。

在没有同期的水准测量数据的情况下，通过选用相邻轨道SAR数据进行InSAR独立计算并进行精度的对比验证，发现两者的计算结果非常一致，表明本文计算结果准确地反映了研究区的地表形变的状况。下一步工作是获取更多的SAR图像参与SBAS的运算以提高InSAR计算的精度，并开展外业水准测量工作与InSAR计算的结果进行对比。

参考文献：

[1] GABRIEL A K，GOLDSTEIN R M，ZEBKER H A.Mapping small elevation changes over large areas：Differential radar interferometry[J].Journal of Geophysical Research，1989，94(B7)：9183-9191.

[2] MASSONNET D，ROSSI M，CARMONA C，et al.The displacement field of the landers earthquake mapped by radar interferometry[J].Nature，1993，364(8)：138-142.

[3] DING X L，LIU G X，LI Z W，et al.Ground subsidence monitoring in Hong Kong with satellite SAR interferometry[J].Photogrammetry Engineering and Remote Sensing，2004，70(10)：1151-1156.

[4] HU J，LI Z W，DING X L，et al.Two-dimensional co-seismic surface displacements field of the chi-chi earthquake inferred from SAR image matching[J].Sensors，2008，8(10)：6484-6495.

[5] GATELLI F，GUAMIERI A M，PARIZZI F，et al.The wavenumber shift in SAR interfereometry[J].IEEE Transactions on Geosciences and Remote Sensing，1994，32(4)：855-864.

[6] ZEBKER H A，ROSEN P A，HENSLEY S.Atmospheric effects in interferometric synthetic aperture radar surface deformation and topographic maps[J].Journal of Geophysical Research，1997，102(B4)：7547-7563.

[7] BERARDINO P，FORNARO G，LANARI R，et al.A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms[J].IEEE Transcations on Geoscience and Remote Sensing，2002，40 (11)：2375-2383.

[8] LANARI R，MORA O，MANUNTA M，et al.A small baseline approach for investigating deformations on full resolution differential SAR interferograms[J].IEEE Transactions on Geoscience and Remote Sensing，2004，42(7)：1377-1386.

[9] 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.

[10] BERARDINO P，CASU F，FOMARO G，et al.Small baseline DIFSAR techniques for earth surface deformation analysis[C].Istituto Per Il Rilevamento Elettromagnetico Dell’Ambiente，IREA-National Research Council of Italy (CNR)，Napoli：2003.

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

[12] 尹宏杰，朱建军，李志伟，等.基于SBAS 的矿区形变监测研究[J].测绘学报，2011，40(1)：52-58.

[13] FERRETTI A，PRATI C，ROCCA F.Nonlinear subsidence rate estimation using permanent scatters in differential SAR interferometry[J].IEEE Transactions on Geoscience and Remote Sensing，2000，38(5)：2202-2212.

[14] KAMPES B M.Displacement parameter estimation using permanent scatterer interferometry[D].Delft：Delft University of Technology，2005.

[15] HOOPER A，ZEBKER H A，SEGALL P，et al.A new method for measuring feformation on volcanoes and other natural terrains using InSAR persistent scatterers[J].Geophysical Research Letters，2004，31 (23)：1-5.

[16] MORA O，MALLORQUI J J，BROQUETAS A.Linear and nonlinear terrain deformation maps from a reduced det of interferometric SAR images[J].IEEE Transactions on Geoscience and Remote Sensing，2003，41(10)：2243-2253.

[17] WNMER C，WEGMULLER U，STROZZI T，et al.Interferometric point target analysis for deformation mapping[A].IGARSS’03，Toulouse，France[C].21-25，July，2003，7：4362-4364.

[18] ANDREAS W，WEMER C，STROZZI T，et al.Combination of point and extended target based interferometric techniques[A].IGARSS’04，Anchorage，Alaska，USA[C].20-24，Sep，2004，2：989-991.