矿山地下煤炭资源开采会造成土地塌陷及采空区积水，进而导致耕地面积减小和农作物产量受损，地表房屋、建筑物及交通运输道路扭曲变形，严重时还会导致山体开裂、崩塌、滑坡、泥石流等重大自然灾害。合成孔径雷达差分干涉测量(differential interferometric synthetic aperture radar，DInSAR)技术是基于面的大范围测量，具有空间分辨率高、覆盖范围广、成本低、快速准确及大尺度连续覆盖等优点^[1-2]，能获取地表雷达视线向(light of sight，LOS)的高精度小量级形变。但矿山形变往往具有沉降速度快、形变梯度大、多量级和多方向等特点，而DInSAR技术仅能监测地表发生的缓慢、小量级形变，因此监测结果不能反映矿山地表真实形变。基于SAR影像的偏移量跟踪(offset-tracking，OT)技术^[3]具有不受云雾遮挡、抗噪声能力强、不受形变梯度制约等特点，能有效监测大范围、大梯度的地表形变，但其监测精度通常在1/10~1/30像元，并不能有效监测下沉盆地边缘附近的小量级形变，在精确识别下沉盆地边缘时会受到影响。因此，有学者联合DInSAR和OT(简称OT DInSAR)技术分别获取了矿山地表LOS向不同量级的形变^[4-8]。但OT DInSAR技术仅能提取矿山的LOS向非连续多量级形变，不能提取矿山连续的三维多量级形变，而传统的三维形变监测方法对SAR影像数量和获取时间有较高要求^[9-12]。基于单轨道InSAR影像提取矿山三维形变能很好地解决这一问题，但该方法目前仅实现了矿山三维非连续大量级形变提取且水平方向位移存在方向偏差的问题^[13-16]。为此，本文联合OT DInSAR技术、多量级形变先验融合模型及概率积分法，提出一种基于单轨道OT DInSAR的矿山时序三维多量级形变重建方法，以获取高精度、完整和真实的矿山地表形变信息。

1 矿山三维多量级形变重建模型 1.1 LOS向多量级形变先验融合模型

基于矿山地理信息和InSAR影像信息的形变先验知识，可构建矿山地表多量级形变融合模型，对DInSAR技术和OT技术获取的大量级和小量级的LOS向形变进行矿山多量级形变融合。基于形变先验知识的融合模型可表示为^[7]：

$ \boldsymbol{D}_{\mathrm{LOS}}=\left\{\begin{array}{l} \left(\operatorname{ISNAN}\left(\boldsymbol{D}_{\mathrm{LOS}}^{\mathrm{DInSAR}}\right) \cdot \mathrm{FILTER}\left(\boldsymbol{D}_{\mathrm{LOS}}^{\mathrm{OT}}\right)+\right. \\ \ \ \ \ \left.\boldsymbol{D}_{\mathrm{LOS}}^{\operatorname{DInSAR}}\right)+\boldsymbol{D}_{P}^{\mathrm{IDW}}, \bar{c}<0.3 \\ \boldsymbol{D}_{\mathrm{LOS}}^{\mathrm{DInSAR}}+\boldsymbol{D}_{P}^{\mathrm{IDW}}, \bar{c} \geqslant 0.3 \end{array}\right. $

(1)

式中，$\boldsymbol{D}_{\mathrm{LOS}}^{\mathrm{DInSAR}}$为DInSAR技术获取的LOS向形变矩阵，其中包括每个像素的LOS向形变值$d_{\mathrm{LOS}}^{\mathrm{DInSAR}}$和空值NaN；$ \boldsymbol{D}_{\mathrm{LOS}}^{\mathrm{OT}}$为利用OT技术获取的LOS向形变矩阵，其中包括每个像素点的形变值$ d_{\mathrm{LOS}}^{\mathrm{OT}}$；函数ISNAN(·)的作用是将形变矩阵中的空值NaN转换为1，非空值转换为0；$ \bar{c}$为目标研究区的平均相干值；函数FILTER(·)的作用是滤除由SAR影像配准和其他噪声引起的失真值和低精度值，其滤波范围为d_{LOS, M}≤$ d_{\mathrm{LOS}}^{\mathrm{OT}}$≤$ d_{R}^{\mathrm{OT}}$，其中d_{LOS, M}为LOS向的最大形变，$ d_{R}^{\mathrm{OT}}$为OT技术可测量地表形变的最高精度值^[7]；$\boldsymbol{D}_{p}^{\mathrm{IDW}} $为利用反向距离加权法(inverse distance weighted，IDW)对空值区域插值后的矩阵；D_LOS为矿山地表多量级形变融合后的结果。对于平均相干值低于0.3的区域，利用DInSAR技术获取沉降边缘的高精度小量级形变，利用OT技术获取沉降中心的大量级形变，然后对其进行LOS向形变融合处理；对于平均相干值高于0.3的区域，将直接以DInSAR技术获取的形变作为最终LOS向的形变监测结果。

1.2 三维多量级形变重建模型

根据随机介质理论推导出的水平移动与倾斜成正比的关系，对于相同的采深H，比例系数是一个常数，则可得水平位移与沉降梯度的关系为：

$ \left\{\begin{array}{l} U_{e}(x, H)=B(H) \cdot \frac{\partial W_{e}(x, H)}{\partial x}= \\ \ \ \ \ b \cdot \frac{H}{\tan \beta} \cdot P(x, H) \\ U_{e}(y, H)=B(H) \cdot \frac{\partial W_{e}(y, H)}{\partial y}= \\ \ \ \ \ b \cdot \frac{H}{\tan \beta} \cdot P(y, H) \end{array}\right. $

(2)

式中，b为水平移动系数，H为采深, β为主要影响角，这些参数均可通过采矿资料获取。P(x, H)为沿x方向的单元梯度，P(y, H)为沿y方向的单元梯度。

水平方向真实的叠加开采工作面的矿山地表形变网格如图 1(a)所示，图中给出了开采工作面的位置和观测线AA′方向的5个形变点，横坐标对应东西向，纵坐标对应南北向，i和j代表形变图的行和列，用来表示形变点的地理位置，i < m、j < n，ΔE和ΔN表示地理编码后的形变图在东西向和南北向的空间分辨率。图 1(b)给出5个形变点在垂直向、东西向和南北向的移动规律，其中水平移动的正值与x轴正向一致，负值与x轴负向一致，形变图中点(1，1)对应计算过程中用到的真实形变矩阵，这对重建真实水平方向的形变尤为重要。

利用矿山地表监测的水平移动数据、下沉数据、采深H及采空区分布，基于概率积分法模型可求取水平移动系数b和主要影响角正切tanβ。Yang等^[15-16]基于单轨道InSAR数据利用李志伟等^[13]的方法提取了同一矿山的三维形变，但测量结果中水平位移的正向和负向移动出现颠倒的情况，造成严重的方向偏差，大大降低了水平向的测量精度。本文以图 1的几何图为参考，给出修正后的矿山地表形变点(i, j)在水平向的位移，具体公式为：

$ \left\{\begin{array}{l} d_{E}(i, j)=\mu_{E}\left[d_{W}(i, j)-d_{W}(i, j+1)\right] \\ d_{N}(i, j)=\mu_{N}\left[d_{W}(i+1, j)-d_{W}(i, j)\right. \end{array}\right] $

(3)

式中，μ_E=bH/(tanβ·ΔE)，μ_N=bH/(tanβ·ΔN)；d_W(i, j)为垂直向形变，d_E(i, j)为东西向形变，d_N(i, j)为南北向形变。

根据SAR传感器的成像几何原理，LOS向形变可分解为三维方向的投影，具体公式为：

$ \begin{gathered} d_{\mathrm{LOS}}(i, j)=d_{W}(i, j) \cos \theta-d_{N}(i, j)\cdot \\ \cos \left(\alpha-\frac{3 {\rm{ \mathsf{ π} }}}{2}\right) \sin \theta-d_{E}(i, j) \sin \left(\alpha-\frac{3 {\rm{ \mathsf{ π} }}}{2}\right) \sin \theta \end{gathered} $

(4)

式中，α为SAR传感器的轨道方位角，θ为SAR传感器发射电磁波的入射角。将式(3)代入式(4)得到：

$ \begin{gathered} d_{\mathrm{LOS}}(i, j)=\left[L_{1}(i, j) \quad L_{2}(i, j) \quad L_{3}(i, j)\right]\cdot \\ {\left[d_{{W}}(i, j) \quad d_{{W}}(i, j+1) \quad d_{W}(i+1, j)\right]^{\mathrm{T}}} \end{gathered} $

(5)

其中，

$ \left\{\begin{array}{l} L_{1}(i, j)=\cos \theta+\mu_{N}(i, j) \cos \left(\alpha-\frac{3 {\rm{ \mathsf{ π} }}}{2}\right)\cdot \\ \ \ \ \ \sin \theta-\mu_{E}(i, j) \sin \left(\alpha-\frac{3 {\rm{ \mathsf{ π} }}}{2}\right) \sin \theta \\ L_{2}(i, j)=\mu_{E}(i, j) \sin \theta \sin \left(\alpha-\frac{3 {\rm{ \mathsf{ π} }}}{2}\right) \\ L_{3}(i, j)=-\mu_{N}(i, j) \sin \theta \cos \left(\alpha-\frac{3 {\rm{ \mathsf{ π} }}}{2}\right) \end{array}\right. $

(6)

由此可得每个像元的垂直向形变，然后代入式(3)可得东西向和南北向的形变。矿山三维多量级形变重建的步骤(图 2)主要为：

1) 输入覆盖研究区的高分辨率InSAR影像和数字高程模型(digital elevation model，DEM)数据。

2) 对时序InSAR影像对进行DInSAR处理，分析相干图和干涉图，得到形变区域的平均相干值$ \bar{c}$。当$ \bar{c}$≥0.3时，利用式(1)进行多量级形变提取，设定相位解缠阈值为0.3，低于0.3的相位信息将被掩膜为空值NaN，通过IDW对其插值处理，最后获取目标区域的LOS向形变结果D_LOS。

3) 当$ \bar{c}$＜0.3时，设定相位解缠阈值为0.3，首先对该区域进行DInSAR处理，获取高精度的LOS向小量级形变结果$\boldsymbol{D}_{\mathrm{LOS}}^{\mathrm{DInSAR}} $；再利用OT技术获取该区域LOS向的大量级形变结果$\boldsymbol{D}_{\mathrm{LOS}}^{\mathrm{OT}} $，使用d_{LOS, M}≤$d_{\mathrm{LOS}}^{\mathrm{OT}} $≤$d_{R}^{\mathrm{OT}} $对其进行滤波处理，对形变区域的空值点进行IDW处理，最终获取连续的LOS向多量级形变D_LOS。

4) 将矿山地表形变点(i, j)对应的LOS向形变值D_LOS(i, j)、α、θ、b、H及tanβ代入式(3)和式(4)，利用式(5)和式(6)进行回代处理可得d_W(i, j)。

5) 将d_W(i, j)代入式(3)可得d_N(i, j)和d_E(i, j)，最终得到基于单轨OT DInSAR重建的矿山三维多量级形变监测结果。

2 研究区分析 2.1 区域概况

大柳塔煤矿地处内蒙古自治区与陕西省的交界处，地理位置为39°13′53″~39°21′32″N、110°12′23″~110°22′54″E。井田位于黄土高原北侧和毛乌素沙漠东南侧，地势中部高、东西低，南低北高，平均海拔为1200m，矿山内大部分地貌属风沙堆积地貌，植被稀少。

研究区为52303工作面，该区域煤层储量大、埋藏深度浅且煤炭质量优良，工作面及观测线地理位置如图 3所示，图中包含52303工作面和52304工作面的地理位置及沿工作面走向(2′2″)和倾向(1′1″)布设的地表移动观测站位置，研究区煤层倾角小于2°，观测站下方煤层倾角近似于水平煤层。52303工作面开采方向为自东向西沿煤层走向回采，工作面为刀把式，采煤方法为走向长壁式采煤法，全部垮落法顶板管理。设计采高为6.5 m，52303工作面面宽为301 m，走向推进长度为4 443.3 m，开采始于2012-12-06，止于2014-03，水平移动系数b=0.32，平均采深H=230 m，主要影响角正切tanβ=1.24。

2.2 InSAR实验数据和DEM数据

选用TerraSAR-X(TSX)雷达卫星的单视斜距复影像数据产品，雷达波入射角为42.43°，轨道方位角为189.53°，像元尺寸为0.91×0.86。表 1为TSX影像参数，为去除地形相位的影响，DEM采用SRTM的30 m分辨率数据产品。

3 实验结果及分析 3.1 OT DInSAR处理

利用OT DInSAR技术获取52303工作面开采造成的LOS向不同量级形变。在DInSAR数据处理中，距离向和方位向采用3×5的多视因子来消除斑点噪声，网格尺寸约为5 m。为提高SAR影像对的相干性，滤波方法选取Goldstein滤波，由于X波段的SAR影像对去相关效应高度敏感，故滤波窗口设为32×32。解缠方法采用最小费用流法，解缠阈值设为0.3，对相干性大于0.3的点进行解缠处理，反之则舍弃。在OT数据处理中，为获取与DInSAR处理结果相同的像素尺寸，多视因子同样采用3×5；综合考虑计算的性能和效率，互相关窗口大小取128×128，过采样因子取4。

52303工作面的上覆岩层为中坚硬岩层，煤层平均厚度为6.85 m，根据矿区地理经验模型计算垂直向最大沉降为-5.75 m，则估计LOS向最大形变d_{LOS, M}约为-4.24 m^[7]。为了滤除OT监测结果中的低精度测量值，$d_{R}^{\mathrm{OT}} $选为-0.25 m, 则滤波范围为-4.24 m≤$d_{\mathrm{LOS}}^{\mathrm{OT}} $≤-0.25 m。图 4为DInSAR和OT的观测结果。可以看出，沿工作面走向存在时间序列相干性较低的黑色区域及沉降中心混叠且不连续的干涉条纹，这是由于矿区地表沉降漏斗中心发生大梯度形变，进而造成严重的失相干。从DInSAR的监测结果可看出，该技术可有效监测矿区沉降漏斗边缘的小量级形变，而沉降漏斗中心出现较多的空值区域；从OT监测结果可看出，该技术可有效监测沉降漏斗中心的大量级形变，对于沉降漏斗边缘的小量级形变监测精度较低。此外，图中还包含一些由SAR影像配准和其他噪声引起的失真值和低精度值，需要进行滤除。

图 5为52303工作面开采造成的矿山地表LOS向多量级形变，其中图 5(a)~5(f)为DInSAR测量的小量级形变与OT测量的大量级形变的融合结果。可以看出，多量级形变不连续且存在较多空值点，这些点多是由DInSAR测量结果中的失相干点和OT测量结果中低精度测量点构成。图 5(g)~5(l)为利用IDW估计后的LOS向连续多量级形变测量结果。可以看出，沉降的方向与该工作面由东向西的开采方向一致，沉降范围与开采周期直接相关。但由于LOS向的测量结果只能反映雷达视线向的形变方向和大小，不能精确获取三维方向的多量级形变，基于LOS向的多量级形变还需重建矿山三维多量级形变。

3.2 矿山三维多量级形变重建结果分析

图 6为52303工作面在2012-12-13~2013-06-07开采造成的时序三维形变，其中2012-12-13~2013-01-04垂直向的最大沉降值为2.31 m，东西向及南北向的最大形变值分别为3.3 m和2.27 m；2013-01-04~02-06垂直向的最大沉降值为3.73 m，东西向及南北向的最大形变值分别为2.24 m和6.42 m；2013-02-06~03-22垂直向的最大沉降值为3.79 m，东西向及南北向的最大形变值分别为1.86 m和5.74 m；2013-03-22~04-02垂直向的最大沉降值为2.75 m，东西向及南北向的最大形变值分别为1.65 m和4.23 m；2013-04-02~05-05垂直向的最大沉降值为4.74 m，东西向及南北向的最大形变值分别为2.55 m和5.59 m；2013-05-05~06-07垂直向的最大沉降值为2.31 m，东西向及南北向的最大形变值分别为2.33 m和2.95 m。针对该工作面的6组观测周期内，最大下沉值为4.74 m，最大下沉速率为0.144 m/d。

分析每组影像对的三维多量级形变最大值后发现，2013-01-04~05-05南北向的形变最大值不符合开采沉陷规律，出现了失真。造成这种现象的原因是由于地表塌陷引起了地裂缝，在形变重建算法中计算水平向移动形变时利用了相邻两点的垂直向形变梯度(式(3))，而严重的地裂缝会产生过大的非线性垂直向形变梯度，进而造成水平向的形变失真。这样的失真点多发生于沉降漏斗底部和塌陷式沉降的边缘，但失真点的数量是有限的，通过对形变图的分析可以人为地识别出这些失真点，不影响对水平向地表形变的分析。

图 7(a)为52303工作面在2012-12-13~2013-06-07垂直向的累积沉降量，随着时间的推移，沉降区域的面积逐渐增大，最大沉降量为5.1 m，最大沉降速度为0.066 m/d。综合开采条件下，由于工作面快速推进，使得煤层上覆岩层各层次下沉速度均加快，相对悬空的时间减少，移动变形集中，引起严重的地表多量级形变，对当地环境造成严重破坏。图 7(b)为52303工作面垂直向形变的三维显示结果，能够直观反映该工作面地表每个点的累积沉降情况，便于对研究区地表损坏程度进行评估与治理。为验证52303工作面三维多量级形变监测结果的精度，与同期地面在观测线2′2″和1′1″的水准监测数据进行对比。由于实际情况中很难获取同步的InSAR与水准数据，其中高程水准数据获取时间为2012-08-19~2013-03-19，水平位移水准数据获取时间为2012-03-20~2013-03-28，而52303工作面开采时间为2012-12-06，本文假设2012-12-06之前52303工作面地表不发生地表位移，因此水准数据与InSAR影像获取时间分别在垂直向和水平向有3 d和6 d的时间间隔。

用本文方法计算水准观测线1′1″垂直向、东西向及南北向形变，并与水准测量结果进行对比，结果如图 8(a)~8(c)所示。通过对比可知，该观测线垂直向、东西向及南北向的多量级形变与水准测量结果的形变趋势较为一致，RMSE分别为0.1880m、0.2243m和0.2074m。沉降漏斗中心有部分观测点与水准测量结果相比整体偏大，主要原因有：1)由于水准测量数据和SAR影像获取时间存在6 d的时间间隔，该监测期间是矿山的主要沉降阶段，矿山形变速率快、形变梯度大，导致与水准测量结果偏差大。2)沉降中心的形变是由OT技术获取的，该技术虽可获取大量级形变，但精度比DInSAR技术低，故会导致大量级形变点存在误差。用本文方法计算水准观测线2′2″垂直向、东西向及南北向形变，并与水准测量结果进行对比，结果如图 8(d)~8(f)所示。可以看出，观测线2′2″上垂直向、东西向及南北向的多量级形变与水准测量结果的形变趋势较为一致，RMSE分别为0.0709m、0.1346m和0.0816m。水准观测线2′2″针对52303工作面相邻的52304工作面布设，在观测周期内只能监测到沉降边缘的小量级形变，在实际情况中，由于52304工作面在2013-03开采完毕，前期开采期间势必会对相邻的52303工作面上方的地表造成影响，会有小量级甚至中量级形变发生，因此本文方法重建的垂直向、东西向和南北向的多量级形变监测结果的RMSE要小于0.1880m、0.2243m和0.2074m，满足矿山三维多量级形变监测的精度要求。

4 结语

针对矿山地表三维多量级形变获取难的问题，本文提出一种矿山三维多量级形变重建方法。利用该方法对研究区2012-12-13~2013-06-07期间7景X波段的TerraSAR-X数据进行时序处理，获取大柳塔矿山52303工作面的时序三维多量级形变，与水准测量结果对比可知，垂直向、东西向及南北向的RMSE分别小于0.1880m、0.2243m和0.2074m，监测结果与水准测量数据形变趋势一致。本文提出的方法可满足监测矿山沉降漏斗边缘地带cm级到沉降中心m级的多量级形变监测需求，将为矿山三维多量级形变测量提供一种大范围、高分辨率、高精度、三维方向、低成本的监测技术，同时为InSAR技术在矿山灾害监测应用的推广提供技术保障。