利用雷达干涉时序分析方法研究地面沉降

吴文豪

(武汉大学卫星定位系统工程技术研究中心，湖北 武汉 430079)

利用雷达干涉时序分析方法研究地面沉降

吴文豪

(武汉大学卫星定位系统工程技术研究中心，湖北 武汉 430079)

天津地区因过度开采地下水导致地面沉降现象严重，对城市建设和发展产生了巨大的负面影响，因此科学分析其地面沉降机理及相关因素非常必要。本文以天津南郊为试验区，基于短基线集干涉技术实施地面沉降监测，并根据误差源的空间相关性采用干涉图的半变异值检测误差是否有效剔除，实现了短基线集处理结果的内部验证，减少了对水准结果检校的依赖。结果表明，相对于比较稳定的天津城区，天津南郊出现了大面积的沉降漏斗。最后结合水文地质数据分析了地面沉降区域与地下水系关联。

地面沉降；短基线；半变异函数

一、引 言

地面沉降是城市的主要地质灾害之一。随着我国城市化进程的加快，地面沉降规模扩大，危害加剧。为有效控制地面沉降的速度，减小其危害性，沉降监测是极为重要的一个环节，也是研究地面沉降机理与发展趋势必不可少的基础性工作。传统的地面沉降测量方法包括水准测量、基岩标测量和分层标测量。随着现代测绘技术的发展，布设GPS地面沉降观测站也成为沉降监测的数据来源，但这些技术都是获取地面沉降的点信息。合成孔径雷达干涉技术是近年来发展的新技术，它能够在广域范围内监测地面沉降。

差分雷达干涉技术对微小的地形变化敏感度极高，可以监测雷达视线方向厘米级的微小位移变化。然而进行长时间的地表微小形变监测中，时间和几何去相关，以及大气扰动能因素限制了该技术的应用。为实现长时间和高精度监测地表形变，许多学者提出了相应的方法来解决去相干问题。1999年Ferretti提出了一种“永久散射体”的时间序列处理方法[1]，随后Berardino提出了短基线集方法，研究低分辨率、大尺度上的形变[2]。与永久散射体方法相比，短基线集方法降低了几何和时间去相关效应。2003年Mora结合PS和短基线集方法的优点，提出了基于高相干点的形变分析方法[3]，其优点是采用少量的SAR图像也可以得到毫米级的结果。

与传统的沉降监测相比，差分干涉雷达广域监测地表变化有着不可比拟的优势。伴随着雷达干涉技术的发展，相关团队利用该技术进行广域范围内的城市沉降监测研究。国外陆续开展了采用雷达干涉技术对热点城市进行地面沉降的研究，这些城市包括洛杉矶、旧金山、拉斯维加斯、墨西哥的墨西哥城、挪威首都奥斯陆、意大利Naples和北部Bologna地区等。国内诸多学者也利用该技术开展了地面沉降监测实验，主要集中在苏州、上海、天津地区。如王超利用欧空局ERS1/2获取苏州市地面沉降测量值[4]；廖明生等探讨了相干点目标分析技术原理和方法，进行了上海城区的地表沉降监测实验[5]；张诗玉、范景辉等分别采用不同的雷达干涉时序分析方法获取天津地面沉降速率[6-7]；葛大庆以河北廊坊城区为实验区，运用了短基线集方法获取了形变时间序列和形变速率[8]；伍吉仓利用27景ENVISAT ASAR数据和18景ALOS PALSAR数据进行短基线集时间序列分析方法处理后得到上海磁悬浮列车线路周边地区的地面沉降的年沉降速率都超过了20 mm/a[9]。无论在国内还是在国外，雷达干涉时序分析方法已经成为研究城市沉降的一个技术热点[10]。

鉴于短基线集技术具有克服时间和空间失相关的优势，本文采用短基线方法获取天津南部郊区沉降结果，并根据差分干涉影像的空间相关性采用半变异函数检验结果的正确性，最后根据天津水文地质资料分析发现沉降漏斗空间分布与地下水系的分布存在相似性。

二、短基线集数据处理

1.短基线集理论

短基线集技术采用基线较短的干涉对组合进行干涉测量，不但继承了常规差分干涉的优点，并且避免了长空间基线的几何去相干问题和长时间基线的时间去相干问题，能够获取长时间缓慢地表形变的演变规律，提高了形变监测的时间分辨率。在短基线集中，高相干点像元依然是通过它们的相位属性来提取的[11]。差分相位干涉图中的两个高相干点(p，q)的观测模型可表示为

式中，ψx，i是第i副干涉图中像元x的缠绕相位；N为干涉图的数量；x，i是空间相关部分的估计值，如地形沉降信息，大气扰动，轨道误差；是空间非相关部分视角误差，如地形误差；γ即为像元在时间维上的相干系数，它的x值越大，越有可能成为高相干点。在数据处理过程中，大气扰动、轨道误差及地表沉降都被认为在小区域内是强相关的。

2.数据处理

本文选取天津南郊作为试验区(图1中灰色框图区域)，采用该区域降轨ASAR影像22景，时间跨度为2007年5月至2009年6月。

剪切后的数据基于短基线集理论进行处理，具体处理流程如下：

(1)差分干涉图生成

1)配准：选取2008年3月14日的ASAR数据作为主图像，其他影像与主影像进行配准，对于多普勒中心频率差异超过100 Hz的影像对，配准前进行多普勒中心频率滤波。配准使主影像和辅影像具有统一的几何投影，是地面沉降干涉处理的关键环节，配准精度影响地面沉降的精度。

2)干涉对组合：一般认为时间间隔为35 d的干涉对沉降量很小，干涉图中只含有噪声，因此选取时间基线大于35 d且小于600 d，轨道基线小于450 m的影像对进行组合。干涉对组合如图2所示。

图1 研究区域与SAR影像覆盖范围

图2 短基线干涉对组合

3)差分干涉处理：首先将干涉对影像共轭相乘获得干涉相位，干涉相位由地形相位和形变相位组成；然后利用DEM数据去除地形相位信息，剩余的形变相位即为干涉条纹图。差分干涉图如图3所示。

(2)高相干点及其时间序列的获取

高相干点就是多景差分干涉图中保持稳定相干特征的目标点。通过监测相干点的相位变化来分析地表沉降，因此相干点提取是差分干涉技术的前提条件。获取高相干点的主要任务是完成部分误差的分离，评估干涉图相位在时间维的稳定性。具体步骤如下：

1)相干系数计算：理论上高相干点是基于相位信息定义的。为了减少数据运算量，首先根据影像强度采用振幅离差法获取高相干候选点。以某一个候选点为中心，指定半径圆形区域内使用自适应滤波器消除大气扰动和轨道误差，以及沉降信息，即x，i。地形误差产生的相位与基线呈线性关系，可采用最小二乘法进行估计，据此可以将不同的相位成分进行分离，然后根据式(2)计算γx值。基于已选择的高相干点和改正后的地形误差，重复进行式(2)的计算，直至计算γx值的均方根变化很小。

图3 短基线集差分干涉图

2)高相干点提取：由于试验数据量有限，时间维相干系数γx的直方图分布与大样本容量的直方图相比会出现偏差。为了消除这种偏差，文中选用模拟随机噪声的方法修正试验数据中的相干系数分布规律。然后设定γx大于0.7的像元为高相干点。

3)相位解缠：差分图中的相位是缠绕的，为了获取真实的相位信息，需要将缠绕的相位信息还原。相位解缠的基本假设是临近相位差分绝对值小于π，但是非空间相关相位噪声很可能导致相位差分出现绝对值大于π的情况，因此解缠前应尽量消除非空间相关误差。高相干点在空间维是稀疏点，而开源解缠软件snaphu只能对规则格网进行解缠。为解决这一问题，首先针对高相干点在空间维建立狄洛尼三角网，在时间维建立方格网；然后采用最邻近插值法将三角网重采样成为规则格网，这样就将三角网解缠转换为规则格网解缠；最后采用snaphu软件完成短基线集的三维解缠。

4)大气效应、轨道误差的消除。空间非相关噪声在沉降模型中可以直接视为随机噪声，但大气效应和轨道误差因其空间相关性会严重影响地面沉降计算结果。因此本文通过在时间维上的高通滤波和空间维上的低通滤波消除大气效应和轨道误差。与其他雷达干涉时间序列分析方法相比，该方法不需要提前预知地面沉降的时间序列模型，直接利用信号的空间相关性来估计像元相位稳定性[12]，如图4所示。

图4 短基线集方法获得的研究区域平均沉降速率图

3.结果验证

如式(1)所述，大气效应和轨道误差是雷达干涉的主要误差源，具有较强的空间相关性。因此分析差分干涉图的空间相关性可以判断大气扰动和轨道误差是否被有效剔除。半变异函数作为表示区域化变量空间相关性的重要函数，可以定量化描述差分干涉图相位的空间相关性[13－14]。因此本文采用半变异函数分析差分干涉图在空间相关性误差消除前后的差异。半变异函数假设某种属性的的空间变异既不是完全随机的，也不是完全确定的。一般认为空间变异可能包含3种影响因素：表征区域变量变异的空间相关因素；表征趋势的“漂移”或结构，以及随机误差[15]。半变异的数学表达式为

式中，γ(h)为已知点xi和xj的半变异；h为这两个点的空间距离；z为空间区域变量，在差分干涉图中z为相位值。如果差分干涉相位不存在空间相关，那么近距离范围内点之间的半变异很小；相反，较远距离的已知点的半变异变大。考虑到大气扰动的各向异性，相同的距离值可能出现多个半变异。因此在文中采用平均半变异和平均距离来绘制半变异图。为了简化计算半变异函数的数据运算量，又进一步对高相干点在空间维上作了降采样处理[16]。

考虑到地面沉降有着很强的空间相关性，本文选取时间间隔较近的4幅干涉图进行半变异分析，半变异图如图5所示。当差分干涉图因轨道和大气扰动的存在，解缠后相位值存在斜坡效应时，半变异表达式可变化为

式中，{E[z(xi)－z(xj)]}2称为偏项，如果相位值不是平稳随机分布，区域变量存在趋势项，则偏项值不为0，半变异值γ(h)呈现开口朝上的抛物线分布。因此图中深色曲线则预示了差分干涉图存在轨道误差和大气扰动，浅色曲线则表示偏项很小或者为0，大气扰动和轨道误差消除效果较好。所以利用半变异值可以对短基线干涉处理结果进行内部检验，验证误差项是否有效剔除，弥补因缺少GPS数据或者水准数据无法进行结果验证的不足。

图5 距离分组后的半变异图

三、结果分析

根据短基线集处理结果，天津东南郊区均出现明显的沉降现象，其中大寺镇和南河镇出现两个典型的沉降漏斗，并且与团泊洼水库组成一个更大的沉降漏斗。天津地区在气候上属于半干旱地区，水资源一直是制约该地区经济发展的重要因素，过度抽取地下水已经造成地下水位大幅度下降，出现大面积的地面沉降。近年来，尽管天津市区的地面沉降已经得到了基本控制，但是随着天津乡镇企业的迅速发展，郊县和农村地区的地下水开采量大幅度增长，天津市区的外围郊区出现了新的地面沉降。根据天津水务局所提供的资料，天津地区主要由堆积平原区组成，第四系覆盖层厚度巨大，含水层组分带比较复杂，埋深400 m范围自上而下可划分为5个含水层组：从上到下依次为第Ⅰ、Ⅱ、Ⅲ、Ⅳ含水层组和咸水含水层组，受地质条件影响，咸水含水层组一般位于第一、二层含水组中[17]。天津的主要采水层为第Ⅱ含水层组，以及第Ⅲ含水层组一部分，水位埋深在100 m左右，地表水很难进行有效补充。短基线集结果显示即使在团泊洼水库周边地区同样出现了大范围的沉降，这也证明了天津地下水的抽取主要是深层地下水，仅仅依靠土壤渗透，即使地表水充足，地表补给量远远小于地下开采量导致地面沉降，地面同样会出现沉降。

一般来说，地下水开采量、水位降深在时空分布上和地面沉降幅度上有明显的一致性，即开采量和水位降深越大的地区和层位，地面沉降量越大，反之亦然。但是短基线集结果与水务局提供的水文资料对比表明津南区含水组漏斗与地面沉降区域并不重合。水文资料显示第Ⅱ、Ⅲ含水组的漏斗中心在津南区的辛庄镇、咸水沽镇以及西青区的杨柳青和中北镇，然而这些区域地面沉降并不严重。并且第Ⅱ、Ⅲ含水组水位埋深一致处于增加的状态(见表1)，不存在地下水回灌的可能。我们推断含水组漏斗和地面沉降区域不一致现象的出现与沉降区主抽水层的厚度、固结度和土体压缩差异性有关。

天津地下水流动系统可分为多个地下水子系统[18]，图4中黑色虚线为海河干流冲海积地下水系统子区与子牙河古河道带地下水系统子区的分界线。本文研究的大寺镇、南河镇、团泊洼沉降区域以及整个东南郊区均属于海河干流冲海积地下水系统子区，而左侧的子牙河古河道带地下水系统子区并未出现明显的沉降。一般情况下，沉降区域应该是中心向外辐射分布，而该处沉降区域却是沿着地下水子系统分布的。而图中的沉降区域轮廓恰恰与地下水子系统的分界线是一致的，这说明地下水子系统之间的侧向补给与各个子系统内部补给还是有差异的，地面沉降区域往往与水文地质的构造有关。

四、结束语

本文采用短基线集方法监测天津南郊地表变化，结果表明包含大寺镇、南河镇及团泊洼库区的城区东南郊出现范围的沉降。为保证结果的准确性，基于差分干涉图相位的空间相关性，对差分干涉图进行半变异分析，验证了大气扰动和轨道误差的有效剔除。但是在进行半变异参数的计算时，降采样处理没有考虑数据的空间相关性，在沉降区域没有作掩膜处理，这有可能会影响到半变异函数的结果，需要作进一步的改进。最后结合水文地质资料对津南地区的地面沉降作了分析，发现地面沉降区域轮廓与地下水系子系统分布是相似的，本文推断这是地下水子系统之间的侧向补给障碍造成的。这一结论对治理地面沉降，合理开采地下水有着指导意义。

表1 含水组漏斗中心水位埋深m

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

[2] BERARDINO P，FORNARO G.A New Algorithm for Surface Deformation Monitoring Based on Small Baseline Differential SAR Interferograms[J].IEEE Transactions on Geoscience and Remote Sensing，2002，40：2375-2383.

[3] MORA O，MALLRQUI J J，BROQUETAS A.Linear and Nonlinear Terrain Deformation Maps from a Reduced Set of Interferometric SAR Images[J].IEEE Transactions on Geoscience and Remote Sensing，2003，40(10)：2243-2253.

[4] 王超，张红，刘智，等.星载合成孔径雷达干涉测量[M].北京：科学出版社，2002：198-205.

[5] 廖明生，卢丽君，王艳，等.基于点目标分析的InSAR技术检测地表微小形变的研究[J].城市地质，2006，1(2)：38-41.

[6] 张诗玉，李陶，夏耶.基于InSAR技术的城市地面沉降灾害监测研究[J].武汉大学学报：信息科学版，2008，33(8)：850-853.

[7] 范景辉，郭华东，郭小方.基于相干目标的干涉图叠加方法监测天津地区地面沉降[J].遥感学报，2008，12(1)：111-116.

[8] GE Daqing，WANG Yan，ZHANG Ling，et al.Mapping Urban Subsidence with TerraSAR-X Data by PSI Analysis[J].IGARSS，2010：3322-3325.

[9] WU Jicang，ZHANG Lina，LI Tao.Monitoring Ground Subsisdence in Shanghai Maglev Area Using PALSAR and ASAR Data[C]∥Proceeding of the Joint Interantional Symposium on Deformation Monitoring.Hongkong：[s. n.]，2011.

[10] 张永红，吴宏安，孙广通.时间序列InSAR技术中的形变模型研究[J].测绘学报，2012，41(6)：865-876.

[11] 张红，王超，吴涛.基于相干目标的DInSAR方法研究[M].北京：科学出版社，2009：30-39.

[12] SOUSA J J，HOOPER A，HASSEN R F.Comparative Study of Two Different PS_INSAR Approaches：DePSI VS StaMPS[C]∥Proceedings of Fringe 2009 Workshop，Frascati，Italy：[s.n.]，2009.

[13] BONCORI J P M，MOHR J J.Statistical Description of Tropospheric Delay for InSAR：Overview and a New Model[J].IGARSS，2007：4483-4486.

[14] HASSEN R F.Radar Interferometry：Data Interpretation and Error Analysis[M].[S.l.]：Kluwer Academic Publishers，2001：154-159.

[15] 陈健飞，张筱林.地理信息系统导论[M].北京：科学出版社，2010：341-346.

[16] LOHMAN R B，SIMONS M.Some Thoughts on the Use of InSAR Data to Constrain Models of Surface Deformation：Noise Structure and Data Downsampling[J].Geochemistry Geophysics Geosystems，2005，6(1)：1-12.

[17] 袁宏利，董民，贾国臣，等.天津平原水利工程地质环境概论[M].郑州：黄河水利出版社，2008：158-161.

[18] 董克刚，王威，于强，等.土水比指标在天津沉降区地下水资源管理中的应用[J].水资源保护，2009，25(6)：51-55.

Research on Ground Subsidence by InSAR Time Series

WU Wenhao

P237

B

0494-0911(2014)11-0011-05

2013-11-18

国家自然科学基金(41274048)

吴文豪(1987—)，男，河南驻马店人，博士生，主要从事雷达干涉测量的研究和应用。

吴文豪.利用雷达干涉时序分析方法研究地面沉降[J].测绘通报，2014(11)：11-15.

10.13474/j.cnki.11-2246.2014.0352