合成孔径雷达差分干涉测量(differential interferometric synthetic aperture radar，D-InSAR)技术由于具备高精度、全天时、面测量、低成本(特别是免费数据源雷达卫星的发射)等特点，引起矿山测量学者的广泛关注，已经成为当前矿山变形监测领域的研究热点。当前提出的三维变形监测方法有多视线向D-InSAR^[1-3]、D-InSAR+Offset tracking^[4-5]、D-InSAR+MAI^[6]、融合GPS和D-InSAR法^[7-8]、地表先验模型+D-InSAR^[9-10]等，主要用于解决地震、滑坡、火山喷发、冰川移动等特点的移动与变形监测问题，并不适用于矿山开采沉陷监测。

近年来，针对煤矿开采地表移动与变形特点，有学者提出了开采沉陷预计模型+单视线向D-InSAR的矿山开采沉陷融合监测方法。文献[11]提出了概率积分法预计模型+ D-InSAR方法进行开采沉陷监测，该方法首先利用D-InSAR技术近似获取地表下沉场，反演出与下沉有关的概率积分参数；然后在假定水平移动系数已知的条件下，利用概率积分法模型和已反演的参数体系预计出开采引起的地表移动与变形，进而实现了对开采沉陷的监测。针对文献[11]方法要求水平移动系数为已知的条件，文献[12]对其进行了改进，首先利用概率积分法模型预计任意像元下沉、SN和EW向水平移动，并投影合成D-InSAR视线(line of sight，LOS)向变形；然后以D-InSAR LOS向实测值和预计值之差的和最小为准则，基于模矢法对所有概率积分法参数进行全局寻优；最后利用概率积分法模型和全局最优参数预计出地表移动与变形，即实现了对开采沉陷的监测。然而文献[12]方法在应用时仍然存在以下两个难点：①求参模型为高度非线性函数，如何保证基于模矢法求取的概率积分参数为全局最优解；②要求掌握工作面的面长、推进尺寸、采高、采深、工作面和煤层的走向、倾向等精细地质采矿条件以及概率积分参数的量值范围。文献[13]提出了基于开采沉陷规律+单视线向D-InSAR的融合监测方法，即依据D-InSAR LOS向变形为下沉、SN和EW方向水平移动沿LOS向投影之和的关系，融合水平煤层开采沉陷规律(下沉、倾斜、水平移动之间的关系)，并在假设采动区地表LOS向变形场中，最后一行、最后一列像元水平移动为0的条件下，首次从LOS向变形中提取了下沉、SN和EW方向水平移动。显然文献[13]方法在适用条件、准确性、计算机编程难易程度上具有一定优势。但文献[13]方法由于仅融合了水平煤层开采沉陷规律，对于倾斜煤层开采沉陷监测并不适用。而工程实践中煤系地层多为倾斜地层，这也使得该方法的推广应用存在较大困难。

鉴于文献[13]方法存在不适于倾斜煤层开采深陷监测的问题，本文基于单视线向D-InSAR技术，拟结合倾斜煤层开采沉陷规律，构建适用于水平、倾斜煤层开采沉陷监测的一般模型。

1 监测模型 1.1 监测模型构建

在倾斜煤层开采时，由于煤层倾角α的存在，地表任意点除了发生水平煤层开采引起的指向采空区的水平移动之外，在煤层倾向方向上还将附加因煤层倾斜引起的地表水平移动分量ΔUy′，文献[14]给出了倾斜煤层开采引起的地表任意点S(x, v)沿走向x′和倾向y′的水平移动U、倾斜T和下沉W的关系：

$ \Delta U y^{\prime}=W(x, y) \cot \theta $

(1)

$ U_{x^{\prime}}(x, y)=b \cdot r \cdot T_{x^{\prime}}(x, y) $

(2)

$ U_{y^{\prime}}(x, y)=b \cdot r \cdot T_{y^{\prime}}(x, y)+\Delta U{y^{\prime}} $

(3)

式中，b为水平移动系数；r为主要影响半径，其中r=H/tanβ，H为采深；tanβ为主要影响角正切值；θ为最大下沉角。

如图 1假设煤层倾角与SN方向存在φ夹角，当计算SN、EW向水平移动时，需要将ΔUy′向SN、EW方向分解，通过推导可得到SN、EW方向水平移动U与倾斜T、下沉W的关系模型：

$ \begin{array}{c}{U^{\mathrm{SN}}(x, y)=b \cdot r \cdot T^{\mathrm{SN}}(x, y)+} \\ {W(x, y) \cot \theta \cos \varphi}\end{array} $

(4)

$ \begin{array}{c}{U^{\mathrm{EW}}(x, y)=b \cdot r \cdot T^{\mathrm{EW}}(x, y)+} \\ {W(x, y) \cot \theta \sin \varphi}\end{array} $

(5)

工作面沿煤层走向布置，o′-x′y′、o-x′y′、o-xy分别表示工作面坐标系、地面坐标系、测量坐标系，o′-x′y′与o-x′y′坐标轴平行，且x′、y′分别与工作面走向和倾向一致。P为过S(x, y)点且与倾向平行的竖直面。

考虑到地表下沉W与倾斜T的关系，

$ T^{\mathrm{SN}}(x, y)=\frac{\mathrm{d} W(x, y)}{\mathrm{d} x} $

(6)

$ T^{\mathrm{EW}}(x, y)=\frac{\mathrm{d} W(x, y)}{\mathrm{d} y} $

(7)

综合上述式(4)~(7)，可得到:

$ \begin{array}{c}{U^{\mathrm{SN}}(x, y)=b \cdot r \cdot \frac{\mathrm{d} W(x, y)}{\mathrm{d} x}+} \\ {W(x, y) \cot \theta \cos \varphi}\end{array} $

(8)

$ \begin{array}{c}{U^{\mathrm{EW}}(x, y)=b \cdot r \cdot \frac{\mathrm{d} W(x, y)}{\mathrm{d} y}+} \\ {W(x, y) \cot \theta \sin \varphi}\end{array} $

(9)

当前SAR卫星多数属于极轨飞行、侧视成像，D-InSAR变形场相邻4个像元排列示意图如图 2所示，假设任意像元的下沉值为W_ij，则式(6)~(9)写成差分形式可表达为:

$ T_{i j}^{\mathrm{SN}}=\frac{W_{i-1, j}-W_{i, j}}{\Delta x} $

(10)

$ T_{i j}^{\mathrm{EW}}=\frac{W_{i, j-1}-W_{i, j}}{\Delta y} $

(11)

$ U_{i j}^{\mathrm{SN}}=b \cdot r \cdot\left(\frac{W_{i-1, j}-W_{i, j}}{\Delta x}\right)+W_{i j} \cot \theta \cos \varphi $

(12)

$ U_{i j}^{\mathrm{EW}}=b \cdot r \cdot\left(\frac{W_{i, j-1}-W_{i, j}}{\Delta y}\right)+W_{i j} \cot \theta \sin \varphi $

(13)

根据开采沉陷移动与变形的定义，同理可得到SN、EW方向水平变形ε、曲率变形K的表达式为:

$ \varepsilon^{\mathrm{SN}}(x, y)=\frac{\mathrm{d} U^{\mathrm{SN}}(x, y)}{\mathrm{d} x}=\frac{U_{i-1, j}^{\mathrm{SN}}-U_{i, j}^{\mathrm{SN}}}{\Delta x} $

(14)

$ \varepsilon^{\mathrm{EW}}(x, y)=\frac{\mathrm{d} U^{\mathrm{EW}}(x, y)}{\mathrm{d} y}=\frac{U_{i, j-1}^{\mathrm{SN}}-U_{i, j}^{\mathrm{SN}}}{\Delta y} $

(15)

$ K^{\mathrm{SN}}(x, y)=\frac{\mathrm{d} T^{\mathrm{SN}}(x, y)}{\mathrm{d} x}=\frac{T_{i-1, j}^{\mathrm{SN}}-T_{i, j}^{\mathrm{SN}}}{\Delta x} $

(16)

$ K^{\mathrm{EW}}(x, y)=\frac{\mathrm{d} T^{\mathrm{EW}}(x, y)}{\mathrm{d} y}=\frac{T_{i . j-1}^{\mathrm{EW}}-T_{i, j}^{\mathrm{Ew}}}{\Delta y} $

(17)

对于D-InSAR监测的任意像元LOS向变形L_ij与下沉W_ij、SN向水平移动U_ij^SN、EW向水平移动U_ij^EW存在如下投影关系：

$ \begin{array}{c}{L_{i j}=-U_{i j}^{\mathrm{SN}} \sin \varphi_{i j} \cos \left(\alpha_{i j}-\frac{3}{2} \pi\right)-} \\ {U_{i j}^{\mathrm{EW}} \sin \varphi_{i j} \sin \left(\alpha_{i j}-\frac{3}{2} \pi\right)+W_{i j} \cos \varphi_{i j}}\end{array} $

(18)

式中，φ_ij为卫星入射角；α_ij为卫星飞行方向方位角。

联合式(10)~(13)以及式(18)可得到：

$ \begin{array}{c}{L_{i j}=-\left[b \cdot r \cdot\left(\frac{W_{i-1, j}-W_{i, j}}{\Delta x}\right)+W_{i, j} \cot \theta \cdot\right.} \\ {\cos \varphi]_{-} \sin \varphi_{i j} \cos \left(\alpha_{i j}-\frac{3}{2} \pi\right)-\left[b \cdot r \cdot\left(\frac{W_{i, j-1}-W_{i, j}}{\Delta y}\right)+\right.} \\ {W_{i, j} \cot \theta \sin \varphi ] \sin \varphi_{i j} \sin \left(\alpha_{i j}-\frac{3}{2} \pi\right)+W_{i, j} \cos \varphi_{i j}}\end{array} $

(19)

由式(19)推导得到：

$ W_{i, j}=\frac{L_{i, j}+b \cdot r \cdot W_{i-1, j} \cdot P_{1} / \Delta x+b \cdot r \cdot W_{i, j-1} \cdot P_{2} / \Delta y}{b \cdot r \cdot P_{1} / \Delta x-\cot \theta \cdot \cos \varphi \cdot P_{1}+b \cdot r \cdot P_{2} / \Delta y-\cot \theta \cdot \sin \varphi \cdot P_{2}+\cos \varphi_{i j}} $

(20)

式中，

根据开采沉陷规律，开采沉陷移动盆地下沉范围往往大于水平移动范围，据此可假设D-InSAR反演的变形场中第一行、第一列像元水平移动量为零，此时式(18)退化成式(21)，即：

$ W_{i j}=L_{i j} / \cos \varphi_{i j} $

(21)

基于上述假设，利用式(21)可求取采动区地表变形场第一行、第一列像元的下沉值，然后在此基础上基于递推关系式(20)，可求取剩余(m－1)×(n－1)像元下沉值(m为影像行数，n为影像列数)，最后利用任意像元下沉值，根据式(10)~(17)求取任意像元沿SN和EW向水平移动值T^SN、T^EW、U^SN、U^EW、ε^SN、ε^EW、K^SN、K^EW。

当最大下沉角取θ=90°，上述倾斜煤层开采地表沉陷监测模型退化成水平煤层开采地表沉陷监测模型。

1.2 监测模型应用流程

基于单视线向D-InSAR技术的倾斜煤层开采地表沉陷监测模型应用流程为：首先利用D-InSAR技术获取采动区地表LOS向变形场；然后根据研究区域地质采矿条件获取采深H、水平移动系数b、主要影响角正切tanβ、最大下沉角θ，以及影像的EW、SN方向分辨率Δx和Δy，并假设影像第一行、第一列水平移动为零；最后利用本文模型计算采动区地表的下沉、SN和EW方向移动与变形。监测模型应用流程如图 3所示。

2 模拟监测实验 2.1 模拟地质采矿条件概况

以淮南矿区地层条件为背景，模拟工作面上覆岩层岩性为中硬，煤层倾角为10°，釆高M=2.5 m，釆深H=300 m，工作面沿倾向和走向水平投影尺寸为D₁×D₃=300 m×600 m，工作面倾向方位角为180°，工作面沿走向布置，顶板管理方法为长壁垮落法。工作面尺寸满足D₁/H=1.0 < 1.2，D₃/H=2>1.2，工作面采动程度为倾向非充分采动，走向超充分采动，总体为非充分采动。工程实践表明淮南矿区沉陷盆地形态可用概率积分法模型进行预计，假设模拟工作面的开采沉陷预计参数为：下沉系数q=0.8，水平移动系数b=0.3，主要影响角正切tanβ=2.0，最大下沉角θ=85°，拐点偏移距S₁=S₂=S₃=S₄=0 m。

2.2 模拟实验及结果分析

以TerraSAR卫星为例进行卫星参数模拟。假设卫星回访目标区几何成像时雷达卫星的入射角为Φ=44.494°，卫星视线向的方位角α=350.596°，分辨率为6 m。模拟工作面开采引起的地表LOS向变形方法如下：首先将采动区地表沿SN、EW方向划分间隔为6 m的方格网，模拟空间分辨率为6 m的影像；然后利用概率积分法模型及选定的参数值预计所有格网点(模拟像元)的下沉值W_{i, j}、SN方向水平移动U_ij^SN和EW方向水平移动U_ij^EW，依据式(18)可计算目标像元沿LOS方向的模拟移动与变形，模拟D-InSAR监测的采动区地表LOS向变形场如图 4(a)所示。

基于上述模拟的采动区地表D-InSAR LOS向变形场，利用本文建立的基于单视线向D-InSAR技术的倾斜煤层开采地表沉陷监测模型，提取了目标区下沉场W、SN和EW方向水平移动场U，提取的采动区地表移动与变形场如图 4(b)-4(d)所示。

为了验证模型的正确性，沿SN、EW方向分别设置了两条移动与变形检验剖面线A-B(点间隔6 m，累计148个点)和C-D(间隔6 m，累计101个点)。沿着剖面线A-B和C-D分别提取了下沉和水平移动，并计算了下沉和水平移动计算相对误差。利用本文模型，沿剖面线A-B、C-D的实际下沉、提取下沉以及下沉提取相对误差绝对值如图 5(a)-5(b)所示；在SN方向水平移动场中沿剖面线C-D的实际水平移动、提取水平移动、水平移动提取相对误差绝对值如图 5(c)所示；在EW方向水平移动场中沿剖面线A-B的实际水平移动、提取水平移动、水平移动提取相对误差绝对值如图 5(d)所示。

在去除移动与变形零边界影响后，从图 5可以看出：①沿剖面A-B提取的下沉相对误差绝对值小于9.53%，平均为1.58%，沿剖面C-D提取的下沉相对误差绝对值小于6.01%，平均为0.92%，沿剖面A-B、C-D提取的下沉相对误差绝对值平均为1.31%；②沿剖面C-D提取的SN方向水平移动相对误差绝对值小于9.78%，平均为4.61%；沿剖面A-B提取的EW方向水平移动相对误差绝对值小于9.22%，平均为3.11%，沿剖面A-B、C-D提取的水平移动相对误差绝对值平均为3.71%；③移动盆地边界附近(下沉0~50 mm区域)，下沉提取相对误差稍大，水平移动零边界附近(水平移动-50~50 mm区域)，水平移动提取误差稍大，这主要因为零边界附近移动与变形量小，较小的提取误差可导致较大的相对误差所致，但是下沉、水平移动提取相对误差极值不超过10%。从实验结果可以看出，基于单视线向D-InSAR的倾斜煤层开采地表移动与变形提取模型可以从D-InSAR监测的LOS向变形场中较为精准地提取下沉和SN、EW方向水平移动。

3 工程应用实验 3.1 实验概况

实验区位于中国山东省兖州矿区南屯煤矿九采区上方地表。实验数据为TerraSAR-X卫星获取的单视复数雷达影像，分辨率为3 m，影像大小为16 224×10 710像素，实验数据覆盖了整个南屯煤矿。主辅影像成像时间分别为2012-01-27和2012-02-07，时间基线为11 d，雷达卫星入射角为θ=44.494°，卫星视线向方位角α=350.596°。在D-InSAR监测期间(2012-01-27—2012-02-07)，南屯煤矿正在回采3上煤层的9310工作面西段，工作面北侧9312工作面已回采结束形成老采空区。9310工作面面长197 m，实验期间工作面推进长度约为26 m，工作面平均采深约588 m，煤层倾角为8°，煤层倾向方位角约为133°，工作面沿走向布置，采用综放采煤方法。

3.2 下沉与水平移动提取实验

通过对主辅影像进行InSAR差分干涉测量，获取了9310工作面上方地表在2012-01-27—2012-02-07时段LOS向移动与变形场如图 6(a)所示。

根据南屯煤矿开采沉陷预计参数经验值，取水平移动系数b=0.30, 主要影响角正切tanβ=2.0，最大下沉角θ=86°，利用本文建立的基于单视线向D-InSAR的倾斜煤层开采地表移动与变形提取模型，从LOS向移动与变形场中提取了监测时段的下沉、SN方向的水平移动、EW方向的水平移动，其中下沉等值线、SN方向的水平移动等值线、EW方向的水平移动等值线如图 6(b)-6(d)所示。

从图 6(b)-6(d)可以看出，最大下沉约74 mm，SN、EW方向最大水平移动约为33 mm；下沉和水平移动分布形态略偏向9312老采空区上方，这主要由于9310工作开采时诱发了相邻的9312采空区边界上方中的悬臂梁、砌体梁结构失稳“活化”导致地表二次移动与变形所致。利用本文模型监测的地表移动与变形特征符合南屯煤矿开采沉陷规律，进一步验证了基于单视线向D-InSAR技术的倾斜煤层开采地表移动与变形提取方法的科学性。

4 结语

本文针对现有的基于单视线向D-InSAR技术的监测方法不适用倾斜煤层开采沉陷监测问题，依据D-InSAR监测的采动区地表LOS向变形与下沉、SN和EW方向水平移动的投影关系，融合倾斜煤层开采地表沉陷规律，通过理论建模，构建了基于单视线向D-InSAR技术的倾斜煤层开采地表沉陷监测模型。模拟实验表明，本文模型下沉监测相对误差绝对值小于9.53%，平均为1.31%，SN、EW方向水平移动监测相对误差绝对值小于9.78%，平均为3.71%，满足开采沉陷监测精度要求。将本文监测模型用于中国山东省兖州矿区南屯煤矿9310工作面开采地表沉陷监测，获取了2012-01-27—2012-02-07时段采动区地表下沉、SN和EW方向水平移动，验证了模型的科学性。