Print

出版日期: 2016-07-25
点击次数：
下载次数：
DOI: 10.11834/jrs.20165134
2016 | Volumn20 | Number 4





                              上一篇|





下一篇
遥感应用
短基线INSAR相干点探测及应用
expand article info 刘利敏1 , 宫辉力2,3 , 余洁1,2,3 , 李小娟2,3 , 柯樱海2,3
1. 武汉大学遥感信息工程学院, 湖北 武汉 430079
2. 首都师范大学资源环境与旅游学院, 北京 100048
3. 首都师范大学资源环境与地理信息系统北京市重点实验室, 北京 100048
收稿日期: 2015-06-02; 修改日期: 2016-03-09; 优先数字出版日期: 2016-07-25
基金项目: 国家自然科学基金（编号：41130744，41171335，41401493）；国家重点基础研究发展计划（973计划）前期研究专项课题（编号：2012CB723403）
第一作者简介: 刘利敏(1988-),女,博士研究生,从事SAR/INSAR数据处理与分析研究。E-mail:xiaoliulimin@163.com
通讯作者简介: 余洁(1964-),女,教授,从事遥感影像处理与应用及空间数据分析研究。E-mail:yuj2011@whu.edu.cn
中图分类号: TP97
文献标识码: A
文章编号: 1007-4619(2016)04-0643-10

摘要

提出一种结合子孔径相关测度的时序高相干点探测方法，首先对时序 SAR 影像进行谱分解获得子孔径视图，通过时序子孔径相关测度进行强散射点筛选，然后分别根据振幅离差和干涉相位空间相关性，对目标点进行相位稳定性分析，探测出既满足强散射且在时间序列上散射稳定的高相干点。利用改进的点探测方法和短基线INSAR技术，对北京2003-2009年间40景ASAR影像进行相位建模，获取地面沉降时空分布特征，分析典型地物时序沉降过程，研究北京地面沉降与地下水开采关系。结果表明：相比已有方法，本文结合子孔径相关测度的高相干点探测结果更准确可靠；利用本文点探测方法反演的历史累积沉降信息，与水准结果一致，平均速率之差在3.69 mm/a以内，中误差为1.36 mm/a；研究区地面沉降最大速率达92.25 mm/a，空间不均匀分布明显，地面沉降量与地下水开采量呈分段的非线性相关。

关键词

子孔径分解 , 高相干点探测 , 短基线INSAR , 时空分布 , 地下水

Stable pointwise target detection method and small baseline subset INSAR used in beijing subsidence monitoring
expand article info LIU Limin1 , GONG Huili2,3 , YU Jie1,2,3 , LI Xiaojuan2,3 , KE Yinghai2,3
1.School of Remote Sensing and Information Engineering, Wuhan University, Wuhan 430079, China
2.College of Resources Environment and Tourism, Capital Normal University, Beijing 100048, China
3.Beijing Key Laboratory of Resource Environment and Geographic Information System, Capital Normal University, Beijing 100048, China

Abstract

Identification of stable pointwise target is an important procedure in multi temporal multi-temporal INSAR analysis and application in monitoring regional surface deformation. The accurate identification result helps to improve the land subsidence inversion precision. Various methods for pointwise target detection have been proposed during the past two decades from different respective. The methods can be divided into some main categories according to the criteria for coherence point selection, such as amplitude dispersion index DA, signalto-Clutter ratio and phase stability. The DA method performs a time series analysis on amplitude instead of phase, and reflects the stability of series amplitude. The advantage of coherence point selection by phase stability can identify some special objects with a stable phase, which further increases the density of the stable pointwise target points, but ignores the highly scattering reflection characteristics of the coherent point. The existing methods take insufficient account of the overall features of stable pointwise targets. For ensuring stable scattering mechanism and temporal stability of pointwise target, an improved method with subaperture correlation was proposed in this paper. First, the subaperture correlation properties IHP of SAR images were obtained by spectral decomposition. Then the stability of targets is evaluated based on series two-aperture spectral coherence, by which the coherence points with high scattering could be identified and detect as PSC1. The DA threshold is utilized as the second criterion, which means any pixel in PSC1 with amplitude dispersion less than 0.4 can be determined as PSC2Then phase stability analysis was carried out to screen out the true stable points from PSC2 with the Characteristics of high scattering mechanism and temporal stability. The experiments of stable pointwise target detection were performed using 33 high resolution SAR images collected by the TerraSAR X-band radar sensor covering Beijing. The detection results demonstrated that the improved method can detect more accurate and reliable pointwise targets than those identified by traditional methodsTo further confirm the effectiveness of the proposed method, the small subset INSAR technique based on the proposed coherence point detection method was adopted to retrieve the ground deformation by 40 scenes dataset acquired from 2003 to 2009 in Beijing. The vertical surface displacement rates during this period was validated by the leveling observations, with RMSE=1.36 mm/a, indicating two types of subsidence matched very well. The maximum subsidence rate of Beijing in investigated area has reached -92.25 mm/a, with an obvious uneven spatial distribution. Subaperture correlation is sensitivity to the high scattering body and can ensure stable scattering mechanism and temporal stability of pointwise target. Both coherence point detection results and the primary surface deformation proved the effectiveness of the proposed method. The deformation result during 2003-2009 has undergone severe land subsidence with high spatial aggregation characteristic, and the regional subsidence and the groundwater exploitation reveal good corresponding relationship, the more exploitation, the higher deformation rate.

Key words

Subaperture decompositions , high coherence point detection , small subset INSAR , temporal-spatial distribution , ground water

1 引 言

随着城市进程加剧，地下水超量开采，北京市地面沉降呈快速增长趋势，累积沉降大于300 mm的地区达到1300 km2，严重威胁城市安全，因此有必要对北京开展地面沉降监测研究。

基于时序INSAR技术地面沉降反演，国内外学者开展了一系列研究。Ferretti等人(2001)提出的PSInSAR技术通过单一主影像进行干涉计算，在高相干点构成的稀疏图上，根据形变及误差相位时空分布特性进行信号分离，反演地表形变。Berardino等(2002)提出小基线集(SBAS)方法，利用多主影像进行干涉建模，有效减弱时空基线引起失相关，获得更多干涉组合。该方法利用SVD方法将多个小基线集联合求解，有效解决不同SAR数据集空间基线过长造成的时间不连续问题，提高监测时间采样率。后续研究者在PSInSAR及SBAS算法框架下提出一系列时序INSAR分析方法(Strozzi等，2006；Hooper等，2008；Perissin和Wang，2012)，如IPTA、STAMPS、QPS等。不同于传统利用最小二乘或周期图求解线性形变速率的方法，STAMPS在空间侧视角误差校正的基础上通过时空3D相位解缠反演真实相位，基于相位组分空间相关性进行时序建模分析，逐步校正和分离误差组分，不受先验形变模型限制，具有较好的可靠性(Hung等，2011)，因此此本研究采用3D相位解缠的短基线INSAR方法对北京地面沉降进行形变解算。

以上基于时序分析INSAR方法，均对高相干点进行相位建模分析，如何准确高效地识别出真正高相干点，直接影响形变解算和误差组分估计。Ferretti等人提出了相干系数阈值，振幅离差指数DA阈值和相位离差阈值法(Ferretti等，2001；Kampes M. 2006)，后续研究者提出振幅离差相位离差，振幅和振幅离差等双阈值选择方法(罗小军，2007)。但这些方法对高相干点整体特性考虑不足，如相干系数法考虑相干点强散射特性而忽略稳定性(聂运菊等，2013)，且导致影像分辨率降低，易产生错误目标点；DA阈值和相位离差阈值法只考虑目标点的时序稳定性，忽略高相干点应是强散射点的前提。Hooper等人(2007)从相位稳定性分析出发，提出的基于时间相干的相干点探测方法，利用相位空间相关性代替空间上窗口估计，提高相干性估计精度，但对所选点的可靠性和稳定性缺乏分析。

一些学者根据地物目标散射特性，利用谱相关分析法根据目标在不同波谱范围的相关性进行单幅影像相干强散射点提取，已被有效应用于单幅SAR影像目标识别及小数据集的相干目标提取中(Lu等，2008；黄昭权，2010)。通过谱分解获得子孔径影像，根据子孔径影像相关性提取在时序上均表现强散射的相干点。为充分考虑高相干点强散射及相位噪声水平低的特性，本文提出了改进的相干点探测方法，首先基于谱分析进行时序强散射点初选，然后利用DA阈值和干涉相位空间相关性分析评价初始候选点相位稳定性，精选出可靠的相干点。

在验证改进的高相干点探测方法有效的前提下，利用该高相干点探测方法及短基线INSAR技术，对40景ASAR影像反演解算北京2003—2009年间区域地表形变，分析地面沉降的时空分布特征及典型地物时序形变过程，研究地面沉降与地下水的作用关系。

2 改进的高相干点探测方法

2.1 利用子孔径相关性的强散射相干点初选

理想相干点在雷达回波信号获取时间内应具有恒定雷达响应，即在任意频谱范围内雷达响应相同，可通过目标在谱分解获取的子孔径影像间的相关性进行强散射点判断，如子孔径相关系数、子孔径共轭内积(Souyris等，2003)及子孔径信息熵(Schneider等，2006)等。共轭内积Internal Hermitian Product(IHP)通过计算两子孔径影像未归一化的相关系数，如式(1)，来衡量子孔径影像在不同频谱上雷达响应的相似性，可有效区分相干强散射目标和分布目标。

$ r = \left| {{s_1} \cdot s_2^*} \right| $ (1)

为尽量保持原影像的数据分辨能力，本文选择将SAR影像进行两视分解，分为前后2个子孔径视图，并利用IHP计算子孔径相关性测度。

本文基于子孔径分解的时序SAR影像强散射点探测，主要包括以下步骤：

(1) 对时序SAR影像进行频谱分析，获取两个子孔径视图

首先对方位向进行FFT频率变换，在距离向上对方位向频谱进行平均，根据得到的频谱估计权重，并对原始频谱进行去加权校正；然后对校正后的方位向频谱进行分割，并对子频谱进行Hamming窗口等的滤波，抑制影像的旁瓣影响；最后将分割后的两个子孔径频谱分别进行逆FFT变换，转换到空间域，完成子孔径分解。

(2) 计算时序SAR影像各自子孔径IHP

计算各SAR影像的两子孔径影像的共轭内积，获得各影像子孔径相关性测度图；为消除不同时相SAR影像受成像星地距离、入射角及天线方向图等参数差异导致的灰度基准不一致问题，需对时序子孔径相关测度进行归化校正，矫正同到同一尺度基准。设有N景影像，首先计算各影像子孔径测度图的均值；并以N景子孔径测度图均值的均值，作为影像时序均值，论文以各影像子孔径测度图均值与时序均值之比，作为该影像的共轭内积校正系数。

(3) 设定相关测度阈值，选择时序相干点

根据矫正后时序子孔径相关测度图[r1，r2，…，rn]，通过设定阈值识别在时序上均保持强散射的相干点。计算各点时序平均相关测度$\bar r$和最小相关测度${r_{\min }}$；然后根据各目标点平均测度$\bar r$均值μ和标准差δ，用于目标检测的阈值的经验值为${T_{{\text{sub}}}} = \sqrt {\mu \delta } $；之后根据时序最小相关测度${r_{\min }}$与阈值${T_{{\text{sub}}}}$的关系进行强散射候选点选取，获得初始候选点集PSC1。

2.2 振幅离差及相位空间相关性的高相干点精选

从高相干目标在时序上保持稳定的特性出发，进行高相干点的精选。对于高信噪比下，振幅离差指数与影像干涉相位噪声水平成正比，可利用振幅离差指数进行时序相干点识别。利用影像上各像元时序振幅标准差及均值的比值，即${{D}_{A}}={{\delta }_{A}}/{{\mu }_{A}}$，计算各像元的振幅离差，通过设定${T_{{\text{DA}}}}$进行相干候选点的探测。研究表明，当阈值取0.4时，可保证约90%具有稳定相位的点被选出(卢丽君等，2008)，本文选0.4作为${T_{{\text{DA}}}}$阈值。将同时满足${T_{{\text{sub}}}}$和${T_{{\text{DA}}}}$阈值的目标点作为候选点PSC2，通过相位空间相关性进行候选点精选。

第x点在第i幅差分干涉图中缠绕相位为：

$ \begin{array}{l} {\varphi _{x，i}} = W \displaystyle \left \{ {{\phi _{{\rm{def}}，x，i}} + {\phi _{{\rm{aps}}，x，i}} + {\phi _{{\rm{ort}}，x，i}}} \right. + \\ \quad \;\; \ \quad\quad\quad\quad\quad\\ \left. { {\phi _{{\rm{top}}，x，i}} + {\phi _{n，x，i}}} \right\} \end{array} $ (2)

式中，${\phi _{{\rm{def}}，x，i}}，{\phi _{{\rm{aps}}，x，i}}，{\phi _{{\rm{ort}}，x，i}}，{\phi _{{\rm{top}}，x，i}}，{\phi _{n，x，i}}$分别表示x在第i幅干涉图两次成像时的形变相位，大气延迟误差、轨道误差、参考DEM误差及时空去相干和系统热噪声误差。${\phi _{{\rm{def}}，x，i}}，{\phi _{{\rm{aps}}，x，i}}，{\phi _{{\rm{ort}}，x，i}}$在一定空间范围内相关，主要分布在低频部分，而${\phi _{{\rm{top}}，x，i}}，{\phi _{n，x，i}}$空间非相关。

在一定空间范围内，根据差分干涉相位各组分空间分布特性差异，利用格网化及频域的自适应带通滤波，粗略估计并滤除差分相位的低频部分$\varphi _{x，i}^c$，则残余相位为：

$ \begin{array}{l} W\left\{ {{\varphi _{x，i}} - \varphi _{x，i}^c} \right\} = W\left\{ {\phi _{{\rm{def}}，x，i}^{uc} + \phi _{{\rm{aps}}，x，i}^{uc}} \right. + \phi _{{\rm{ort}}，x，i}^{uc} + \\ \quad \;\; \ \quad\quad\quad\quad\quad\\\quad\quad\quad\quad\ \left. { \phi _{{\rm{top}}，x，i}^{uc} + \phi _{n，x，i}^{uc}} \right\} \end{array} $ (3)

式中，$\phi _{{\rm{def}}，x，i}^{uc} + \phi _{{\rm{aps}}，x，i}^{uc} + \phi _{{\rm{ort}}，x，i}^{uc}$表示空间相关相位残余高频成分之和，值很小；空间非相关$\phi _{{\rm{top}}，x，i}^{uc}$与基线成正相关：$\phi _{{\rm{top}}，x，i}^{uc} = 4\pi {{\rm B}_{ \bot x，i}}\theta _x^{nc}$，可通过空间搜索估计求解。去除空间相关相位$\varphi _{_{x，i}}^c$及非相关相位$\varphi _{_{{\rm{top}}，x，i}}^{uc}$后，最终残余相位约为噪声相位。

噪声水平越低，被选为相干点的可能性越大，定义像元相位残差变化时间相干模型γ对噪声水平进行评价：

$ {\gamma _x} = \frac{1}{N}\left| {\sum\limits_{i = 1}^N {\exp \left({j\left({{\varphi _{x，i}} - \varphi _{x，i}^c - \varphi _{{\text{top}}，x，i}^{uc}} \right)} \right)} } \right| $ (4)

利用时间相干系数γ对候选点PSC2点集进行精选，最终探测出在满足强散射体特性且在时间序列上保持相位稳定的高相干点。

3 短基线INSAR技术形变信息解算

基于短基线的INSAR方法，首先根据N景影像间的时空基线距及多普勒频率差构建的相干性测度，进行短基线组合选取，形成M幅差分干涉图。针对本文方法探测得到高相干点，利用空间非相关相位 对初始缠绕相位进行矫正，保证相邻相干点满足Nyquist采样定理，满足可相位解缠先验条件。通过高相干点构网，利用时空3维解缠，恢复网络上各相干点在时序真实相位，并根据误差组分空间分布特性通过时空高低通滤波逐步分离大气延迟和轨道误差AOE及空间相关侧视角误差SCLA，获得各点真实形变相位。对于由ta和tb时刻影像构建的干涉对i，其累计形变如式(5)。利用M个干涉对，构建如式(5)的，包含N-1个待求参数的M个观测方程组，通过最小二乘求解各时间段的形变速率及序列。

$ \begin{aligned} {\varphi _{{\rm{def}}，x，i}} & =\varphi \left\{ {tb} \right\} - \varphi \left\{ {ta} \right\} = \\ & \;\;\; \frac{{4\pi }}{\lambda }\sum\limits_{j = ta + 1}^{tb} {{v_{k，k - 1}}\left({{t_k} - {t_{k - 1}}} \right)} \end{aligned} $ (5)

4 实验结果及分析

4.1 改进的高相干点探测方法有效性验证

为验证本文高相干点探测方法的有效性，同时为了便于高相干点细节展示，采用北京地区33幅分辨率为3 m的TSX影像进行方法验证，研究区大小为2000×1226像素，约6×3.6 km2。

对时序TSX影像进行谱分解，获取对应子视图，以子视图的共轭内积IHP作为强散射点选取指标，为保证识别结果整体可靠的前提下获得更多的候选点，选择${T_{{\text{sub}}}} = 60.8$作为子孔径阈值，共选出250848个点，占总像素数的10.23%，如图 1(a)，相干点主要分布在建筑物及道路上，且以建筑物为主，在农田和河流上几乎未有分布，相干散射点和分布散射目标得到较好的区分，缺点是相干点空间分布不够均匀，而且该方法只从相干点的强散射特性出发，而未考虑相干点稳定性，可能导致部分在时序上不稳定点被选出；以DA=0.4为阈值，利用振幅离差指数进行相干点识别，识别出163774个点，约占总像素数6.68%，如图 1(b)。相比图 1(a)，DA阈值法结果，虽然只识别到6.68%的点，但是依旧有大量河流及植被上的点被误选，说明该方法对水体及植被等弱散射体的区分能力不如子孔径相干测度法，但该方法识别的相干点空间分布相对均匀。在振幅离差初选的基础上，利用干涉相位的空间相关性分析对初选相干点进行精选，结果如图 1(c)，该方法与振幅离差阈值探测结果在空间分布上类似，但可有效筛选掉农田与河流中的离散目标点，但由于它也是从时序相位稳定出发进行相干点探测，依旧会导致部分在时序上表现稳定的非强散射点被误识。本文首先根据子孔径相关测度筛选强散射体点，再结合振幅离差指数和空间相位稳定性分析对强散射体点的时序稳定性进行分析，精选出10507个既满足强散射体特性又在时序上散射稳定目标，相干点在研究区上均匀分布，结果如图 1(d)。


图 1 不同方法相干点探测结果分布
Fig. 1 Coherence point detection results of different methods

选取图 1中典型目标点A、B、C，统计分析其时序子孔径相关测度及振幅变化。A点为带顶棚的建筑物；B点为河流堤岸的点；C点主要分布在植被区，散射特性较弱。

A点时序最小子孔径相关测度为77.53，大于子孔径阈值，表明在时间序列上具有强散射特性，但是如图 2示，该点在时序上振幅震荡明显，振幅离差达0.557，在子孔径相干测度选择结果中，该点被误选。由B及C点时序变化图发现，B及C振幅在时序上虽然震荡，但总体变化平稳，振幅离差均小于0.4，满足时序相位受噪声影响水平低、散射稳定的需求，但通过时序子孔径相关测度发现，B及C的最小子孔径相干测度均小于阈值60.8，不满足强散射特性，并非真实相干点，但在DA阈值及空间相位稳定性分析的相干点识别方法未能进行有效筛除。而由于本文方法分别通过强散射特性及目标相位稳定的特点进行探测，保证目标在满足强散射特性前提下时序相位变化平稳，在基于本文方法的结果中，A，B，C等只满足一种散射特点的目标，被逐步排除，得到较好的相干点识别结果，且探测的相干点点密度486个/km2，满足密度要求，说明改进的相干点点探测方法更为合理可靠。


图 2 典型点ABC振幅及子孔径相关测度时序变化
Fig. 2 Amplitude and IHP time sequences of typical points

为进一步验证本文方法的有效性，分别基于(振幅离差指数+相位空间分析点识别结果和本文方法相干点识别结果对研究区进行地面沉降形变解算，结果如图 3，与水准对比结果见表 1。


图 3 不同方法识别PS点形变结果解算
Fig. 3 Deformation results caculating by different coherence points

表 1 不同方法形变解算结果对比
Table 1 Deformation results by different Method

下载CSV 
水准点 水准/(mm/a) 传统INSAR /(mm/a) 差异 本文INSAR /(mm/a) 差异
BM1 -21.93 -18.74 3.19 -23.29 1.36
BM2 -18.52 -16.67 1.85 -19.13 0.61

由图 3发现，虽然两组沉降结果形变趋势与沉降中心较为一致、沉降速率较为吻合，但两组结果在部分沉降细节上却存在差异。传统相干点识别方法得到形变图 3(b)，由于未要求相干点满足强散射特性，识别点数目较多，但根据并不完全可靠的相干点解算的形变结果，也存在部分形变细节缺失，如3(b)中椭圆以及方框所示的小型沉降区，在图 3(a)中并未得到有效监测。为评价两种方法形变结果，以研究区范围内2个水准点做可靠性评定，见表 1，结果表明本文相干点识别方法，在满足点密度要求的前提下，获得形变结果更为准确，进一步表明本文点识别方法可靠性。

4.2 改进相干点探测的短基线INSAR地面沉降监测

为对本文点探测方法进一步验证，并分析北京市平原区地表形变状况，利用改进的点探测方法及短基线INSAR技术，采用2003—2009年间40景Envisat ASAR数据集对北京市地面沉降进行形变信息提取及时空演化关系分析。

研究区以(116.60°E，40.06°N)为中心点，覆盖范围达109×104 km2。分别利用SRTM3 DEM数据及ODR精密轨道数据对干涉相位进行地形轨道误差去除。实验为获取该地区可靠的地表形变监测结果，根据时空基线及多普勒频率形成157个连通，且无孤立干涉对的差分干涉图集。采用本文高相干点探测方法，最终识别出394325个高相干点，平均密度为30.38个/km2。根据大气、形变以及噪声在时空域的分布特性，对3维空间解缠后的相位，利用最小二乘及带通滤波逐步去除大气延迟相位、空间相关侧视角误差及轨道误差相位等，获取北京市平原区2003—2009年间的时序形变序列及平均形变速率。图 4及图 5分别为研究区地面沉降年均速率及时序形变速率标准差图。


图 4 研究区年平均沉降速
Fig. 4 Mean deformation rate of study area


图 5 时序速率标准差
Fig. 5 Deformation standard deviation of study area

由图 5发现各相干点速率标准差整体较小，最大为2.75 mm/a，且大部分相干点速率标准差处于较小水平，表明时序形变速率整体较稳定，估算的形变速率值可靠。

4.2.1 沉降结果与水准结果精度验证

为评价基于本文点探测方法和短基线INSAR技术的地表形变反演精度，利用6组2005年—2009年间的水准监测数据与其最邻近高相干点平均垂向沉降速率进行对比，GCP水准点分布见图 6，表 2 为INSAR形变监测结果与水准结果对比。


图 6 GCP点分布位置示意
Fig. 6 Distribution of the GCP points

表 2 时序INSAR技术与水准获取的平均沉降速率对比
Table 2 Deformation results by INSAR and Leveling

下载CSV 
/(mm/a)
水准点 水准 INSAR 差值
GCP1 -32.6 -29.13 3.47
GCP2 -23.6 -23.35 0.25
GCP3 -2 -3.5 1.5
GCP4 -19.8 -17.29 2.51
GCP5 -58.7 -61.4 2.7
GCP6 -13.9 -13.14 0.76

由图 6及表 2发现两种监测方法结果表现一致，沉降速率最大差值为3.69 mm/a，差异均值及中误差分别为2.09 mm/a和1.36 mm/a；对比表明本文通过短基线INSAR技术获取的北京平原区地面沉降结果精度较高，进一步证明本文改进的相干点识别方法的有效性。

4.2.2 地面沉降空间分布及演化

由图 4知，北京市2003—2009年间地面沉降漏斗呈片状分布，最大沉降速率达92.25 mm/a。研究区沉降速率大于30.71 mm/a的相干点约占10%(如图 7)，约70%的相干点沉降速率小于10 mm/a，相对稳定，快速沉降的相干点虽数量不多，但空间分布较广，主要分布在朝阳区、顺义区及通州部分地区，主城区沉降相对不明显。研究区地面沉降差异较大，区域不均匀沉降明显，其中昌平沙河-八仙庄、顺义平各庄、朝阳来广营和通州北部地区沉降较为严重。整体而言，平原区西南部沉降较少，东北部地面沉降严重，且以较大的速率持续沉降，空间聚集特性显著。结合图 5发现，速率标准差相对较大区域，主要分布在地面沉降漏斗中心，沉降过程迅速，速率梯度变化较大，非线性形变特性明显。地面沉降漏斗所在朝阳区，有多个大型人工设施及多条交通线贯穿(如图 8)，如机场、京承高速、首都机场高速、京哈高速、京沈高速以及京津塘高速的部分路段等，交通线运营产生的动载荷会加速地面沉降的发生，同时在不均匀沉降的大背景下，地面沉降会对贯穿其中的交通线安全造成危害。该漏斗区城市近地表空间开发利用加大，在地下水长期过量开采、复杂立体交通网络及高密度建筑群引起的静载荷等因素综合作用，导致该区域地面沉降过程复杂，这也是该区域时序速率标准差相对较大的原因。


图 7 研究区不同形变速率相干点比例
Fig. 7 Rate ratio of coherent point with different deformation


图 8 地面沉降漏斗区与交通线叠加图
Fig. 8 Overlay chart of subsidence funnel and traffic lines

为分析各主要沉降区随时间变化规律，提取研究区3个典型目标进行时序分析，时序变化如图 9。位于主城区的a点，在时序上表现较为稳定；处于沉降漏斗中心的b点，在2003年—2009年间，沉降速率达89.65 mm/a，累计沉降量达574 mm，沉降过程急剧；平均沉降速率为19.52 mm/a的c点，在时间序列上不均匀性明显。从3个典型目标在时序上的沉降过程发现，在2003年12月至2007年8月间，地面沉降不均匀特性显著，之后沉降速率逐渐趋于平稳。


图 9 典型目标a、b、c时序变化过程
Fig. 9 Deformation time sequences of typical points

4.2.3 地面沉降与地下水关系分析

结合2003年—2009年间北京地下水流场数据与本文获取地面沉降结果，进行地面沉降与地下水相关性分析。图 10(a)—(f)分别为2004年—2009年地下水水位以及历年内地下水开采量，图 10(g)为地面沉降结果，图 10(h)为文中使用剖面线P1—P2及黑椭圆典型区在地下水开采图上的放大展示。由图 10发现，地面沉降与历年地下水水位漏斗空间展布具有较好的一致性，地下水流场变化与地面沉降具有较好的响应，计算2004年—2008年历年地下水水位以及6年间地下水水位变化量与地面沉降年均速率的相关系数，结果如表 3。


图 10 不同年份地下水水位及地面沉降漏斗分布
Fig. 10 Overlay chart of subsidence funnel and groundwater contours

表 3 历年地下水水位及地下水开采与地面沉降相关性
Table 3 Correlation of deformation and groundwater exploitation

下载CSV 
2004年水位与
沉降相关性		 2005年水位与
沉降相关性		 2006年水位与
沉降相关性		 2007年水位与
沉降相关性		 2008年水位与
沉降相关性		 地下水累计开采与
沉降相关性
系数 0.30322 0.73374 0.39932 0.46113 0.47388 0.50817

在2003年—2009年6年内，除2004年地下水位与地面沉降相关系数为0.3032外，其余各年水位及2004年—2009年内地下水累计开采与地面沉降的相关系数均约等于或大于0.4，表明地下水开采量与地面沉降整体相关。

以2008年水位为例，沿图 10(g)中P1—P2所示剖线，提取地下水水位及地面沉降量，结果图 11，地面沉降与地下水位变化趋势一致，沉降发生在地下水水位埋深较深的地区。对比图 10(f)，2004年—2009年6年间地下水开采量与图 10(g)地面沉降的空间分布发现，二者空间分布具有一定的一致性，但并不完全吻合；沿P1—P2，分别提取地下水年均沉降速率、地面沉降累计形变与地面水开采的剖面，分别如图 12(a)和(b)，发现地面沉降与地下水开采总体变化趋势一致，地下水开采越大，沉降速率越大，累计沉降越大，对于红框所示特殊位置，地下水水位降落量为正，即2009年地下水水位较2004年地下水有所抬升，但该段剖线在历时6年内却以10—40 mm/a速率持续下沉，其中该段剖线位置为图 10(h)中黑色虚方框所示。这表明地下水开采与地面沉降并非简单的线性相关，而呈现出分段式的复杂非线性变化，这也解释历年地下水水位与地面沉降结果虽表现相关，但是线性相关系数却不够大的原因。同时选取图 10(g)中黑色椭圆所示的朝阳南部沉降漏斗区，结合历年地下水水位图，发现该区域2005年该区域的地下水水位较低，从2006年开始，由于地下水开采中心逐步向东北郊区方向聚集，该区域地下水水位稳定，图 10(f)显示该区在6年内水位变化接近于0，但图 10(g)显示该区地面沉降明显，且4.2.2节中典型目标B正好分布在该区域内，发现该区在2003年—2009年内呈持续快速下沉状态。以上结果表明北京地面沉降与地下水开采具有较好的相关性，地面沉降漏斗与地下水水位及开采量空间展布一致，但并非完全吻合，地下水开采与地面沉降表现为分段的非线性关系，部分地下水水位变化较小区域，也可能发生较大沉降，这是由于地下水超采引发的地面沉降具有滞后性，此外地面沉降除受地下水开采影响，也与水文地质条件、地质构造、可压缩层厚度等存在相关性。


图 11 2008年地下水水位与地面沉降分布
Fig. 11 Spatial variation of deformation and groundwater in 2008


图 12 年均沉降速率及累计沉降与地下水开采量关系
Fig. 12 Time series variation of deformation and groundwater

5 结 论

针对已有相干点识别方法的局限性，本文提出了结合子孔径和相位稳定性分析的子孔径探测方法，利用子孔径共轭内积敏感于强散射点的特性进行时序相干点初选，之后利用振幅离差及空间相位相关性评价候选点受噪声影响水平，最终探测出满足强散射特性且时序散射稳定受噪声影响较小的高相干点，试验结果表明本文方法可获得较普通振幅离差法及相位稳定性分析更好的相干目标探测效果。利用本文相干点探测方法及短基线INSAR技术反演2003—2009年间北京地面沉降时空演化信息，发现北京区域地面沉降空间分布不均匀性明显，在东北部聚集呈片状连通分布，西部城区较为稳定，在时间演化过程中，2003年12月至2007年2月之间地面沉降整体不均匀性较2007年2月至2009年12月更为明显。分析同期历年地下水水位以及地下水降落水位与地面沉降的相关性，发现地面沉降分布受地下水开采影响显著，且地面沉降与地下水开采量呈分段式非线性相关，总体表现为地下水水位开采越大，地面沉降速率越快，但也存在特殊区域如朝阳南部沉降漏斗区，地下水累计开采量减少，但是地面沉降持续明显的情况。此外沉降漏斗区受地下水开采、地表复杂动静载荷的综合作用，形变梯度较大，沉降过程复杂，时序不均匀性较其他区域更为明显。论文定性分析了地下水与地面沉降的时空演化关系，但地下水对地面沉降的具体贡献值和定量作用特点仍有待开展后续研究。

参考文献(References)

  • Berardino P, Fornaro G, Lanari R, Sansosti E.2002.A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms. IEEE Transactions on Geoscience and Remote Sensing, 40 (11) : 2375–2383 . [DOI:10.1109/TGRS.2202.803792]
  • Ferretti A, Prati C, Rocca F.2001.Permanent scatterers in SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing, 39 (1) : 8–20 . [DOI:10.1109/36-898661]
  • Ferretti A, Fumagalli A, Novali F.2011.A new algorithm for processing interferometric data-stacks:SqueeSAR. IEEE Transactions on Geoscience and Remote Sensing, 49 (9) : 3460–3470 . [DOI:10.1109/TGRS.2011.2124465]
  • Huang Z Q. 2010. Monitoring ground deformation with Small Dataset differential InSAR. Zhejiang University, 23-35
  • Hooper A, Segall P, Zebker H.2007.Persistent scatterer interferometric synthetic aperture radar for crustal deformation analysis, with application to Volcán Alcedo, Galápagos. Journal of Geophysical Research, 112 (B07407) : [DOI:10.1029/2006JB004763]
  • Hooper A.2008.A multi-temporal InSAR method incorporating both persistent scatterer and small baseline approaches. Geophysical Research Letters, 35 (16) : L16302
  • Hung W C, Hwang C Chen Y A.2011.Surface Deformation from Persistent Scatterers SAR Interometry and Fusion with Leveling Data:a Case Study over the Chouhui River Alluvial Fan, Tai Wan. Remote Sensing of Environment, 115 : 957–967 . [DOI:10.1016/j.res.2010.11.007]
  • Kampes M. Radar Interferometry:Persistent Scatterer Technique. Germany, German Aerospace Center (DLR). 2006 .
  • Laurent F, Andreas R, Eric P, Wolfgang-Martin B. Scene Characterization Using Subaperture Polarimetric SAR Data. IEEE Transactions ON Geoscience and Remote Sensing. 2004 : 2264 -2276.
  • Luo X J. 2007. Theory of differential SAR interferometry based on permanent scaterer and applications to Shanghai's surface subsidence setection. Southwest Jiaotong University, 52-62
  • Lu L J, Liao M S, Wang T, Tian X.2008.A Multi-Step detection method for extraction of stable pointwise target in long temporal SAR image series. Journal of remote sensing, 12 (4) : 561–567 .
  • ( 卢丽君, 廖明声, 王腾, 田馨. 2008. 遥感学报. 一种在长时间序列SAR影像上提取稳定目标点的多级探测法, 12 (4) : 561–567. )
  • Lu L J, Liao M S, Wang C C, Tian X, Li D R.2008.A new method of identification of stable pointwise target in small SAR dataset. Proceedings of the 2008 Dragon Symposium, ESA SP-655, April 21-25, Beijing, P. R. China Nie Y J, Liu G X, Shi J F, Yu B, Cheng P G. 2013. Improved method for persistent scatterer detection and its application. Jornal of remote sensing, 17 (3) : 632–639 .
  • Nie Y J, Liu G X, Shi J F, Yu B and Cheng P G. 2013. Improved method for persistent scatterer detection and its application. Jornal of remote sensing, 17(3):632-639
  • ( 聂运菊, 刘国祥, 石金峰, 于冰, 程朋根. 2013. 改进的PS探测方法及其应用. 遥感学报, 17 (3) : 632–639. )
  • Perissin D, Teng W.2012.Repeat-pass SAR interferometry with partially coherent targets. IEEE Transactions on Geoscience and Remote Sensing, 50 (1) : 271–280 . [DOI:10.1109/TGRS.2011.2160644]
  • Strozzi T, Wegmuller U, Keusen H R.2006.Analysis of the terrain displacement along a funicular by SAR interferometry. IEEE Geoscience and Remote Sensing Letters, 3 (1) : 15–18 . [DOI:10.1109/LGRS.2005.855072]
  • Souyris J C, Henry C, Adragna F. On the Use of Complex SAR image spectral analysis for target detection:Assessment of polarimetry. IEEE Transactions on Geosience and Remote Sensing. 2003 : 2725 -2734.
  • Schneider R Z, Papathanassiou K P, Hajnsek I. Polarimetric and Interferometric Characterization of Coherrent Scatterers in Urban Areas. IEEE Transactions on Geosience and Remote Sensing. 2006 : 971 -984.