采煤引起的地面沉降是矿区常见的一种破坏性灾害，这种渐变、长期的过程会严重影响人类的生存环境和当地经济发展^[1]。星载雷达干涉测量(InSAR)作为一种有效的大区域高空间分辨率地表形变监测技术^[2]，是对逐点测量、以点代面的局部测量手段(如水准测量、GPS等)的补充。在常规InSAR技术基础上发展而来的PSI时序分析技术^[3]，解决了时间和空间失相干及大气波动两大问题^[4-5]，具有对微小形变敏感、mm级监测精度等优势^[6]。

在唐山市地面沉降监测的相关研究中，晏明星^[7]利用2004~2006年ENVISAT影像监测了唐山市地面沉降，结果表明，南湖公园附近沉降趋于稳定，唐山矿沉降在2004年之后有所减缓；张学东等^[8]利用2004-11~2010-07的ENVISAT数据对唐山地区进行形变时序分析，发现2009年以后唐山矿的沉降量明显增大。除此之外，未见2010年以后其他利用InSAR监测该地区地面沉降的相关研究。煤矿区地面沉降具有形变速率较快的特点，前人的研究工作主要利用中等分辨率雷达影像研究面上目标形变。由于InSAR技术最大可监测梯度与雷达天线波长、SAR影像空间分辨率和重访周期有关，以往的中等分辨率数据限制了其可探测的最大形变梯度，也缺乏对特定目标的精细监测。因此，本文使用高分辨率TerraSAR-X影像改善中分影像在应用中的不足，基于PSI技术通过两次回归分析的方法降低地形残余误差，进而估算每个相干目标的形变信息。

1 PSI技术

PSI技术是Ferretti等^[3]提出的PS技术和随之发展起来的其他技术方法^[9-11]的统称。该技术与DInSAR技术的主要不同之处在影像数量和数据模型。PSI技术以差分干涉图上相干特性保持良好的相干目标为研究对象，以避免地物散射特性变化引起相干性的降低。然后，利用干涉图上两个相邻PS点的相位之差进行时序分析^[3]，以减小大气相位和轨道误差的影响。由于两个相邻点的非线性形变速率之差非常小，因此可忽略不计。第k幅差分干涉图上相邻两个PS点x和y的相位差为：

(1)

式中，Δφ_{topo(x, y)}^k为地形残余误差，Δφ_{def(x, y)}^k为相干点在LOS向的位移之差，Δφ_{res(x, y)}^k为残余相位，主要包含非线性形变相位和未去除干净的轨道误差、大气延迟及噪声相位。

地形相位与地物高程之间呈线性相关，斜率大小受垂直基线影响^[12]，为了获取较多的干涉对，垂直基线阈值设置不能过小，因此地形残余相位不可忽略。多数PSI方法的形变估算均为线性模型，本文干涉对的垂直基线在时间域上呈随机变化，因此该模型是对高程误差和形变速率同时估算的二维回归分析：

(2)

式中，λ为雷达信号波长，R为传感器到主影像的斜距，θ为天线入射角，Δφ_mod和Δφ_res分别为k幅干涉图上N个PS点的差分相位矩阵和相对残余相位矩阵，B和T为k幅干涉图的空间垂直基线矩阵和时间基线矩阵，ΔH和ΔV为N个PS点的相对高程误差矩阵和相对形变速率矩阵。其中，Δφ_mod、B和T为已知参数，ΔH和ΔV为待求参数，与水准测量平差类似，利用最小二乘方法进行求解。

为降低地形残余误差，需进行两次回归分析。首次估算得到高程改正值，然后从差分干涉相位减去第一次估算的高程误差，最后再次进行二维回归分析求解形变速率。为了求得形变速率增量，基于整体相位相干系数建立最优化方程，利用已知的少数稳定点上的线性形变速率通过区域增长的方法得到各相干目标上的绝对形变速率^[13]。然后，利用空间低通滤波的方法得到非线性形变相位和大气延迟相位，再利用时域高通滤波的方法去除大气延迟得到非线性形变相位，每个PS点的形变即为线性和非线性部分的总和。PSI技术获取形变速率的处理过程见图 1。

2 研究区概况和数据处理 2.1 研究区概况

唐山矿位于河北省唐山市西南部，井田范围由市中心经西南郊延至丰南市，整个井田长14.55 km，宽3.5 km，自1878年开凿1号竖井至今已有130多年的开采历史。唐山矿属于典型的“三下”(建筑物下、铁路下、水体下)开采，地面建筑物众多，地下管线复杂，原京山铁路(货运线)横贯井田中央，铁路压煤线路长达10 km。

2.2 PSI数据处理

TerraSAR-X卫星载有X波段雷达天线，重访周期为11 d^[14]。本文选取13幅3 m分辨率的TerraSAR-X影像，接收时间为2012-12~2014-02，影像覆盖唐山市区及周边煤矿区。利用资源三号立体像对提取精细DEM作为辅助高程数据，其空间分辨率为15 m，高程精度为10 m。为了提高干涉图的相干性和时间采样密度，选择空间基线小于220 m，时间间隔小于400 d的干涉像对组成47幅干涉图序列，图 2为干涉图组合模式。

PSI技术的核心为选取相干性保持良好的地面目标。本文采用强度稳定性和频谱相干性相结合的方法，共提取出92万多个相干目标，约1 200个/km²。PS密度相比于中等分辨率影像大幅增加^[8]，主要原因是影像空间分辨率较高，天线发射信号的波长较短^[15]。然后，利用Delaunay三角网建立相邻相干目标的关系以进行空间相位解缠。

对式(2)进行二维回归分析，然后从差分干涉相位减掉回归分析获得的高程改正值，以降低地形残余误差的影响。将所有PS点的高程误差进行统计，得到地形残余误差去除前后的直方图(图 3)。图中红色为地形残余误差改正前的直方图，高程误差平均值为7.24 m；绿色为改正后的直方图，高程误差靠近0 m的频率大幅增高，平均值降低至1.51 m，说明该方法可有效改善地形残余误差。

3 唐山矿地面沉降时序分析

图 4为利用PSI二维回归分析方法得到的唐山矿2012-12~2014-02的地面沉降速率图。从图 4(a)可以看出，唐山矿地面沉降主要分布在东北和西部地区，沉降中心与前人研究的结果^[8]一致。研究区内的主要沉降中心为中国唐山地震博物馆(最大沉降速率为58.24 mm/a)和北新庄子村(最大沉降速率为34.96 mm/a)。图 4(a)中红色矩形框为研究区较稳定的区域，该区域提取的PS点约3万个，经统计得出其形变速率均值为1.41 mm，标准差为3.33 mm，验证了本文使用的PSI技术mm级的测量精度。

实地调查发现，该地区已出现因地面不均匀沉降而引发的建筑和设施损坏的现象。如图 4(b)所示，左上图和左下图拍摄于唐山地震博物馆，左上图在纪念墙附近因不均匀沉降导致地面拱起，左下图为不同方向的应力作用导致玻璃开裂；右上图为地面塌陷导致地面积水，右下图为建筑出现了墙体裂缝。

此外，通过提取PS点的沉降时间序列，可进一步分析沉降特征。图 5为东北沉降区A点的沉降时间序列。该点在2013-06~09有回弹现象，其他时间段近似线性下沉，累积沉降量达65 mm。杜兴明等^[16]认为，目前唐山市地基土的稳定性主要取决于采空区顶部的土体而不是岩体，土体中的土-水比例变化经各土层压缩固结表现为地面沉降或回弹。从河北省网站统计数据得知，2013-06~09唐山市总降水量达486 mm，占全年总降水量的87%，本文得到的地面回弹现象可能与降水量较大有关。

4 结语

本文利用2012-12~2014-02获取的13景高分辨率TerraSAR-X影像，采用PSI时序分析技术对唐山矿进行地面沉降监测，得到了该地区的地面沉降速率图和沉降时间序列图。研究结果表明：

1) 高分辨率影像提取的PS点密度与中等分辨率影像相比大大提高。

2) 本文通过两次回归分析的方法降低了地形残余误差，改正前后的直方图显示，高程误差平均值从7.24 m降低到1.51 m，说明该方法效果显著。

3) 两个沉降中心分别位于唐山矿的东北和西部地区，最大沉降速率达58.24 mm/a。实地调查发现，该区域地面和建筑均有不同程度的破坏，主要受地面不均匀沉降影响。

4) 沉降中心A点的沉降时间序列图表明，2013-06~09地面有所回弹，可能与该时间段降水量较大有关。