2011, Vol. 54
2. 东华理工大学地球科学与测绘工程学院, 江西抚州 344000;
3. U. S. Geological Survey Center for EROS, Sioux Falls, SD 57198, USA;
4. 太原市水务局, 太原 030002
2. School of Geosciences and Surveying and Mapping, East China Institute of Technology, Fuzhou Jiangxi 344000, China;
3. U. S. Geological Survey Center for EROS, Sioux Falls, SD 57198, USA;
4. Taiyuan Water Resources Bureau, Taiyuan 030002, China
地面沉降又称地面下沉或地陷.它是自然因素和人类工程经济活动因素作用下，由于地下松散地层固结压缩，导致区域性地面标高降低的一种地质现象[1].自1921年上海出现地面沉降以来，至今我国已有96个城市和地区发生了不同程度的地面沉降[2].太原市位于山西省太原断陷盆地的北缘，由于水资源短缺导致长期超量开采地下水，成为我国中西部地区地面沉降最为严重的城市之一.太原市地面沉降发现于20 世纪50 年代，至80 年代沉降加剧，至2003年已形成西张和城区两个沉降区，西张､万柏林､下元､吴家堡4 个沉降中心.沉降区南北长约39km, 东西宽约15km, 沉降面积达548km2.20世纪90年代以来，沉降范围逐年向盆地边缘扩展，沉降漏斗面积逐年扩大，南部有向晋中盆地延伸趋势[3， 4]，已影响到该地区社会经济的可持续发展.
为了采取合理的控制沉降措施，需要及时､有效地监测地面沉降.传统的监测方法主要有水准测量､基岩标和分层标测量.这些方法精度很高，但只能在比较小的范围内开展工作.对于大区域地面沉降监测可采用先进的全球定位系统(Global Positioning System, GPS)进行全方位的测量.但是由于GPS监测站的分布密度较稀疏，且成本昂贵，不适合大规模推广.随着雷达卫星遥感技术的快速发展，自20世纪90年代以来，利用星载合成孔径雷达差分干涉测量(Differential Interferometric SAR，DInSAR)方法进行大区域的地面沉降监测已取得可喜进展，可以精确地监测出地面沉降的变化[5].然而对于长时间缓慢地面沉降而言，由于SAR 干涉像对的时间和空间基线距较大，导致影像时间､空间失相干严重[6， 7]，利用传统DInSAR 技术难以反演其形变序列.为了克服这些困难，近年来，研究人员基于干涉点目标分析思想陆续提出了多种新的DInSAR 地表形变反演方法.其中代表性的有:永久散射体法(Permanent Scatterers, PS) [8~13]､小基线集法(Small Baseline Subsets, SBAS)[14~17]､相干目标法(Coherent Targets, CT)[18， 19].这几种方法都是在长时间尺度上，选择那些相位和幅度变化稳定的点作为干涉点目标，基于这些离散点目标反演地表形变信息.其中PS 法主要采用幅度离差指数阈值选取点目标，通过单一主影像组合干涉像对，它要求SAR 影像数目较多(通常超过25 景)，否则难以达到足够的形变反演精度，在提取非线性形变和消除大气相位影响方面PS 法给出很好的思路；小基线集法则将所有获得的SAR 数据组成若干个子集，子集内的SAR 影像基线距小，子集间的基线距大，它在反演线性形变时没有考虑大气影响；相干目标法则以相干系数阈值选取高相干点目标，通过多主影像组合干涉像对，只需少量SAR 影像就可以生成足够的干涉图，有利于形变反演，在提取非线性形变方面该方法计算复杂，通过分别计算低､高分辨率的非线性形变部分实现，容易引入新的算法误差.
针对上述方法的特点，本文拟采用小基线DInSAR 方法，它通过融合PS法和相干目标法的优点实现地面沉降反演，其中线性形变反演主要采用相干目标法，而非线性形变提取采用了类似PS 法的计算思想.基于该方法，我们反演了太原市2003~2009年长时间尺度的地表形变场，并利用研究区的水准实测数据对反演结果进行分析验证.
2 研究方法与数据准备 2.1 研究方法本文采用基于高相干点目标的小基线DInSAR技术，将地表形变反演分为两个部分，即线性形变反演和非线性形变反演.其主要内容包括以下几个方面:
(1) 干涉组合
假设在(t1，t2，…，tN)时间获取同一地区N幅SAR 图像，根据一定干涉条件组合(如时间基线距和垂直基线距不超过一定范围)，得到M对干涉组合有如下不等式关系:
|
(1) |
干涉组合不要求具有共同主图像，仅要求主辅图像都是按同一个时间顺序排列.这样处理使得少量SAR 图像也能组合较多的干涉图.对ENVISATASAR 数据来说，在干涉组合时，大基线距干涉像对要做距离向频域滤波，以消除距离向频谱偏移，增强干涉图相干性.
(2) 高相干点目标提取
点目标的后向散射特性几乎不随时间变化而变化，在长时间尺度上可以保持较高的相干性，因此可以将高相干点作为点目标.计算M幅干涉图中各像元的相干系数，并得到其平均相干系数，选择适当的相干系数阈值，将平均相干系数大于该阈值的像元认为是高相干点，如下式:
|
(2) |
其中ri表示单个像素在第i幅干涉图上的相干系数，thd_r是相干系数阈值.要注意，高阈值会舍弃许多有用的相干点，低的阈值又会带来许多噪声，因此选择合适的阈值是非常重要的.一般对于ERS数据或ASAR数据而言，当估计窗口大小为20×4或15×3时，阈值可以在0.25~0.35之间选择；但不是所有选择的点都是有效点，其中可能存在噪声点[18， 19].
(3) 线性形变反演
M幅干涉相位图经去平地和去地形后，可相应得到M幅差分相位图，各像素的相位模型可表示为
|
(3) |
式中，i为第i幅干涉图，φdef 为视线向的地表形变相位，包括线性形变和非线性形变，φerrortopo 为DEM引起的高程误差相位，φatm为大气影响相位，φnoise 为噪声相位，包括热噪声､时间和基线去相干噪声.其中地表形变相位φdef 和高程误差相位φerrortopo 可分别表示为:
|
(4) |
|
(5) |
式中，λ 为雷达载波波长，r为目标到雷达传感器斜距长，b为垂直基线，θ为雷达入射角，ε和v分别为高程误差和线性形变速率，T表示两次SAR 影像获取的时间基线.
在此基础上，通过Delaunay三角网连接所有的点目标.对于三角网上的任一条边，其两顶点(xm，ym)，(xn，yn)之间的相位差为:
|
(6) |
式中(xm，ym)､(xn，yn)为两顶点的位置坐标，Ti为第i幅干涉图的时间基线，β 为非线性形变相位项，α 为大气影响相位项，n为噪声相位项.
考虑到大气影响相位是一个低频信号，在空间上存在一个相关距离，一般约1~3km, 当三角网两点间距在大气影响相关距离内时，可以认为它们的大气影响相位相等，即α(xm，ym，Ti)≈α(xn，yn，Ti)，由于线性速度和高程误差是常数，而非线性形变相位和噪声相位都是随机信号，因此(6)式可以演变为线性模型:
|
(7) |
其中，Δv和Δξ 分别是速度和高程误差增量.为了求得这两个参量，建立如下的目标函数方程，通过最优化方法求解:
|
(8) |
式中，γ 为整体相位相干系数，其最大值在[0,1]之间，反映了速度和高程误差增量对差分相位的拟合程度.当对所有的边完成最大化求解后，选择γ≥0.7的边作为可靠的连接，从某一参考点开始通过区域增长的方法对速度和高程误差增量积分，得到各点目标的线性形变速率和高程误差绝对量.
(4) 非线性形变反演
从原始差分干涉图相位减去模型相位(7)式后，得到点目标上的残余相位φres, 它包括大气影响相位φatm, 非线性形变相位φnon-linear 以及噪声相位φnoise.为了得到完整的形变信息，需要对解缠后的残余相位进行时空频谱特征分析，以分离出非线性形变相位.
残余相位的三个分量中，大气影响相位在时间序列上是不相关的，为高频信号，在空间分布是相关的，为低频信号；非线性形变相位在时间序列上是低频信号，在空间上是不相关的，为高频信号；噪声相位则是时间和空间都不相关的随机高频信号.利用这些表现特征，可以将三者分离出来.对位于x处的点目标，首先在时间序列上做频域低通滤波，提取出低频的非线性形变相位，然后利用最小二乘方法及干涉组合关系计算出各时刻的非线性形变量，将其与线性形变叠加可得到全部形变信息.对低通滤波后的残余相位，以1km(大气相关距离)的窗口大小作空间低通滤波消除噪声相位影响，从而提取大气相位.
在具体实现中，关键步骤是时间序列的频域低通滤波.我们发现采用基于Kaiser窗函数的有限长单位冲激响应数字滤波器(Finite Impulse Response, FIR)对残余相位进行卷积滤波可以达到理想的效果[20]:
|
(9) |
式中，h(·)为时间低通滤波核.因为Kaiser窗可选择主瓣宽度和旁瓣衰减，其适应性强，而FIR 滤波器结构简单且稳定，二者组合可以很好地抑制高频信息，提取出非线性形变相位.但是由于不同成像时刻､不同成像地区的大气性质差异较大，因此需要对残余干涉相位进行频谱分析，以选定合适的截止频率､过度带宽和阻带衰减，构造出最佳参数的Kaiser窗函数.
2.2 数据准备本文选择时间跨度从2003年8月17日~2009年2月1日的23景欧空局ENVISATASAR 影像，数据格式为原始信号数据(RAW)，经距离多普勒成像算法生成单视复影像(Single Look Complex, SLC).ASAR 的中心工作频率采用C 波段(频率为5.331GHz)，成像中心入射角为23°，具体成像参数见表 1.研究区范围为太原市，覆盖面积约为40×30=1200km2(图 1).为消除干涉相位中的地形影响，需要研究区的数字高程模型数据(Digital Elevation Model, DEM)模拟地形相位.为此，我们下载了研究区约90m 分辨率的SRTM DEM 数据.为了对地面沉降反演结果定标和验证，我们从太原市水务局获取了研究区2006~2007年11个水准点实测数据(图 1).
|
表 1 本文所用ENVKSAT ASAR数据成像参数 Table 1 ENVISAT ASAR images used in this research |
|
图 1 研究区范围及水准点空间分布背景为23景ASAR影像平均幅度图，红色点为水准点. Fig. 1 Study area and leveling data location Background is the average amplitude of 23 ENVISAT ASAR |
利用2.1节介绍的方法，在小基线(时间基线距不超过3a, 垂直基线距不超过400m)干涉组合原则下，由23景ASAR影像生成93组干涉像对.并以0.3为相干系数阈值，从平均相干图上提取了9940个高相干的点目标.在此基础上，我们首先得到太原市2003年8 月~2009 年2 月平均沉降速率图(图 2).可以看出在该时间段内，太原市地面沉降主要分布在中部和南部地区.北部平原地区地面沉降速率较小，大多不超过-10mm/a, 部分地区地面还有反弹回升.东西两侧太行山､吕梁山山缘洪积扇冲击平原的沉降速率也较小，多数在-25mm/a以内.总的来看，太原市的地面沉降从下元沿汾河往南呈喇叭状展开.区内最大沉降速率位于孙家寨，达-77.28 mm/a.主要沉降中心有4个，从北往南依次为:下元，最大沉降速率为-55.83 mm/a；吴家堡，最大沉降速率为-51.79mm/a；小店镇，最大沉降速率为-64.20mm/a；孙家寨，最大沉降速率为-77.28 mm/a.此外，从孙家寨沿汾河谷地至小店镇，发育了一条沉降速率较大的沉降带，沉降速率为-47.5~-77.28mm/a.
|
图 2 2003〜2009年太原市平均沉降速率 Fig. 2 Average subsidence rate of Taiyuan City from 2003 to 2009images, the red points are leveling data location. |
与2000年前相比，该阶段的太原市地面沉降分布出现明显变化.2000年以前太原市地面沉降漏斗中心主要有西张､万柏林､下元和吴家堡，沉降中心沿汾河呈串状分布，最大沉降中心位于吴家堡一带，南部地区如小店镇则没有发生明显沉降[4].2003~2009年，沉降中心则明显南移，呈喇叭状分布，且集中在小店区一带.在该阶段，过去的西张沉降中心较为稳定，平均沉降速率在-5mm/a以内，部分地区还有反弹，最大反弹速率达5.7mm/a；万柏林沉降中心沉降速率也大为减小，最大速率约为-16.21mm/a.下元和吴家堡虽仍为沉降中心，但沉降速率与1989~2000年的相比都明显减缓，说明太原市近年采取的“关井压采"等措施已取得初步成效.该阶段最显著的地面沉降发育特点是:小店镇至孙家寨一带的沉降速率快速上升，成为最大的沉降中心；而2000年前该地区沉降现象并不明显，没有形成沉降漏斗.表 2对比了太原市各沉降中心在1989~2000年与2003~2009年两时间段沉降速率变化情况.
|
|
表 2 不同时间段太原地面沉降中心沉降速率对比(mm/a) Table 2 Comparison of subsidence rates between different stages for subsidence centers in Taiyuan |
在反演线性形变速率的基础上，继续反演非线性形变，将线性形变量和非线性形变量结果叠加，对离散点空间插值，得到太原市各时间段的累计形变量(图 3).从图中可以看出，随着时间的延长，南部沉降中心逐渐显现出来，最大沉降中心出现在孙家寨附近，累计沉降量达-422.40 mm；而北部地区则出现一定反弹，从20030817 到20090201 最大累计反弹量为74.13 mm.为了定量地显示各时间段累计形变量，我们从11 个水准点中选取7 个点位(P0､P1､P2､P3､P5､P7､P10)，分析这些点2003~2009年的累计沉降序列(图 4).可以看到，太原南部4个沉降中心附近水准点(P3､P5､P7､P10)的累计沉降量较大，如杜家寨(P10)累计沉降量超过-350mm, 小店地下(P7)超过-300 mm；从时间上来看，杜家寨(P10)和小店(P7)两处大致以20050821为节点，沉降速率开始增加；吴家堡(P5)则从20050717 开始，沉降速率减小；下元地下(P3)的沉降速率基本没有变化.对于北部水准点(P1､P2)则表现为地面反弹回升，上兰地下(P0)由于地处基岩区，总体较为稳定，没有明显的形变，这与实际情况相吻合.
|
图 3 2003〜2009年太原市各时间段累计沉降量 Fig. 3 Accumulative deformation of Taiyuan City in different stages from 2003 to 2009 (a) 20030817〜20040104；(b) 20030817〜20050123；(c) 20030817〜20060212；(d) 20030817〜20070128； (e) 20030817〜20080113； (f) 20030817〜20090201. |
|
图 4 7个水准点位置的累计沉降序列(2003〜2009年) Fig. 4 Accumulative deformation on 7 leveling points from 2003 to 2009 |
本实验结果的验证采用太原市水务局提供的2006~2007 年11 个水准点数据(图 1).表 3 显示了相应时间段内DInSAR形变反演结果与水准实测形变值的误差对比.可以看出，除水准点P1和P10外，其余水准点的相对误差都在(5 mm 以内，误差的均方差为2.90 mm, 该结果满足我国对城市地面沉降5mm/km 的偶然中误差和10mm 的全中误差[21]精度要求.这说明小基线DInSAR 技术的反演精度较高，可以应用于地面沉降监测.
|
|
表 3 水准形变值与IHnSAR反演结果对比 (2006〜2007 年)(mm) Table 3 Comparison of deformation from 2006 to 2007 between DInSAR result and levelling measurement (mm) |
本文研究了小基线DInSAR 技术的处理算法，深入分析了线性形变和非线性形变反演方法.该方法的主要特点是:(1)与PS 法反演形变相比，它不需要大量的SAR 影像；(2)通过多主影像进行干涉组合，可以生成更多的干涉图，为高精度的形变反演提供更多的源信息；(3)相对传统DInSAR 技术而言，它可以提供更高时间分辨率的地表形变监测.利用该技术，本文对太原市2003~2009年的地表形变场进行了反演，水准实测数据验证了监测结果.实验结果表明，该阶段太原市存在有4 个明显的沉降中心，即下元､吴家堡､小店､孙家寨，其中最大沉降中心位于孙家寨，平均沉降速率为-77.28mm/a.
综合考虑地面沉降的原因，不外乎三个方面:城市建筑压载､构造活动､人为开采地下水､气､油､煤等资源.太原市地面沉降严重区主要分布在郊区而非城市中心区，这就排除了由于城市建筑压载而产生地面沉降的可能.太原地处内陆新生代断陷盆地，目前该区处于活动期的断裂主要有晋祠大断裂.但相关研究表明[4]，晋祠断裂与地面沉降变化趋势上无明显的相似性或关联，因此可以认为构造活动对该地区地面沉降没有较大影响.太原虽是我国的产煤重镇，但煤矿开采集中在西山和东山两地区1)，在中部和南部并没有煤矿开采区，因此煤矿开采不是太原市地面沉降的主要原因.由于太原地处温带半干旱大陆性气候，降雨量小，人均水资源量少，是我国严重缺水城市之一.加之城市工业的快速发展和人民生活水平的不断提高，长期以来超采地下水导致太原市地下水降落漏斗不断加剧，地面沉降也不断加深，形成西张､万柏林､下元､吴家堡四个沉降中心.
1)太原市国土资源局.《太原市矿产资源规划》(2001~2010年)
2003~2009年，随着“关井压采"措施的实施和小店区的太原经济技术开发区的建立，使得太原市地面沉降中心的空间分布发生了显著变化，呈现南移趋势.在“关井压采"控制区内地面沉降趋势已经趋于减缓，如西张和万柏林地区沉降速率大为减小，甚至出现反弹回升，说明控制沉降措施取得显著成效.而小店区，由于经济技术开发区的建立，大量工厂入驻，经济快速发展，且该区农业灌溉用水需求很大，导致该地区存在大量自备深水井，地下水超采异常严重，地面沉降快速发展，成为太原市新的沉降中心.图 5显示了开发区内沉降速率较大相干点对应的工业厂房分布位置(由于部分厂房建于2005年之后，所以地表相干性不高，难以提取高相干点目标).为了全面控制太原市地面沉降发展，遏制南部沉降区扩大，需要地方政府部门加大“关井压采"措施，从而促进该地区社会经济的可持续发展.
|
图 5 太原市经济技术开发区内高相干点目标对应的厂房分布位置(黄圈表示，背景为GoogleEarth提供的卫星影像) Fig. 5 Distribution of factory buildings (yellow circles) corresponding to high coherent points in Taiyuan Economic and Technological Development Zone (Background is satellite image provided by Google Earth) |
| [1] | Whittaker B N, Reddish D J. Subsidence: Occurrence, Prediction and Control. Amsterdam: Elsevier Science Publishers, 1989 |
| [2] | 郑铣鑫, 武强, 侯艳声, 等. 城市地面沉降研究进展及其发展趋势. 地质评论 , 2002, 48(6): 612–618. Zheng X X, Wu Q, Hou Y S, et al. Advances and trends in research on urban land subsidence. Geological Review (in Chinese) , 2002, 48(6): 612-618. |
| [3] | 孙自永, 马腾, 马军, 等. 太原市地层空间异质性对地面沉降分布的影响. 岩土力学 , 2007, 28(2): 399–405. Sun Z Y, Ma T, Ma J, et al. Effect of strata heterogeneity on spatial pattern of land subsidence in Taiyuan City. Rock and Soil Mechanics (in Chinese) , 2007, 28(2): 399-405. |
| [4] | 闫世龙, 王焰新, 马腾, 等. 内陆新生代断陷盆地区地面沉降机理及模拟—以山西省太原市为例. 武汉: 中国地质大学出版社, 2006 . Yan S L, Wang Y X, Ma T, et al. Mechanism and simulation of land subsidence in the Cenozonic inland faulted basin: a case study of Taiyuan City, Shanxi Province, China (in Chinese). Wuhan: China University of Geosciences Press, 2006 . |
| [5] | 王超, 张红, 刘智. 星载合成孔径雷达干涉测量. 北京: 科学出版社, 2002 . Wang C, Zhang H, Liu Z. Spaceborne synthetic aperture radar (SAR) interferometry (in Chinese). Beijing: Science Press, 2002 . |
| [6] | Gatelli F, Guarnieri A M, Parizzi F, et al. Wavenumber shift in SAR interfereometry. IEEE Transactions on Geosciences and Remote Sensing , 1994, 32(4): 855-864. DOI:10.1109/36.298013 |
| [7] | Zebker H A, Rosen P A, Hensley S. Atmospheric effects in interferometric synthetic aperture radar surface deformation and topographic maps. Journal of Geophysical Research , 1997, 102(B4): 7547-7563. DOI:10.1029/96JB03804 |
| [8] | Ferretti A, Rocca F, Prati C. Permanent scatterers in SAR interferometry. Proceedings of IGARSS 1999, Hamburg, Germany, 28 June~2 July, 1999 |
| [9] | Ferretti A, Prati C, Rocca F. Nonlinear subsidence rate estimation using permeanent scatterers in differential SAR interferometry. IEEE Transactions on Geosciences and Remote Sensing , 2000, 28(5): 2202-2212. |
| [10] | Ferretti A, Prati C, Rocca F. Permanent scatterers in SAR interferometry. IEEE Transactions on Geosciences and Remote Sensing , 2001, 39(1): 8-20. DOI:10.1109/36.898661 |
| [11] | Kim J S, Kim D J, Kim S W, et al. Monitoring of urban land surface subsidence using PSInSAR. Geosciences Journal , 2007, 11(1): 59-73. DOI:10.1007/BF02910381 |
| [12] | Vilardo G, Ventura G, Terranova C, et al. Ground deformation due to tectonic, hydrothermal, gravity, hydrogeological, and anthropic processes in the Campania Region (Southern Italy) from Permanent Scatterers Synthetic Aperture Radar Interferometry. Remote Sensing of Environment , 2009, 113: 197-212. DOI:10.1016/j.rse.2008.09.007 |
| [13] | 张永红, 张继贤, 龚文瑜. 基于SAR干涉点目标分析技术的城市地表形变监测. 测绘学报 , 2009, 38(6): 482–487. Zhang Y H, Zhang J X, Gong W Y. Monitoring urban subsidence based on SAR interferometric point target analysis. Acta Geodaetica et Cartographic Sinica (in Chinese) , 2009, 38(6): 482-487. |
| [14] | Berardino P, Fornaro G, Lanari R, et al. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Transactions on Geosciences and Remote Sensing , 2002, 40(11): 2375-2383. DOI:10.1109/TGRS.2002.803792 |
| [15] | Berardino P, Casu F, Fornaro G, et al. Small baseline DIFSAR techniques for earth surface deformation analysis. FRINGE 2003 workshop, Franscati, Italy, 1~5 December, 2003 |
| [16] | Berardino P, Casu F, Lanari R, et al. On the exploitation of the SBAS algorithm for the analysis of the deformations detected from the ERS and ENVISAT DInSAR data. ERS-ENVISAT DInSAR Symposium, Salzburg, Austria, 6~10 September, 2004 |
| [17] | Casu F, Manzo M, Pepe A, et al. SBAS-DInSAR analysis of very extended areas: First results on a 60000 km2 test site. IEEE Transactions on Geosciences and Remote Sensing Letters , 2008, 5(3): 438-442. DOI:10.1109/LGRS.2008.916199 |
| [18] | Mora O, Mallorqui J, Broquetas A. Linear and nonlinear terrain deformation maps from a reduced set of interferometric SAR images. IEEE Transactions on Geosciences and Remote Sensing , 2003, 41(10): 2243-2253. DOI:10.1109/TGRS.2003.814657 |
| [19] | Wu T, Wang C, Zhang H, et al. Deformation retrieval in large areas based on multi-baseline DInSAR algorithm: a case study in Cangzhou, northern China. International Journal of Remote Sensing , 2008, 29(11-12): 3633-3655. |
| [20] | 吴宏安. 基于最小二乘DInSAR技术反演大区域地面沉降研究. 北京:中国科学院遥感应用研究所, 2009. Wu H A. Retrieval of ground subsidence over larger areas based on advanced least squares DInSAR technique (in Chinese). Beijing: Institute of Remote Sensing Application, Chinese Academy of Sciences, 2009 |
| [21] | 北京市测绘设计研究院. 城市测量规范(CJJ 8-99). 北京: 中国建筑工业出版社, 1999 . Beijng Institute of Surveying and Mapping. Code for Urban Survey (CJJ 8-99) (in Chinese). Beijing: China Building Industry Press, 1999 . |