-
合成孔径雷达差分干涉测量(differential interferometry synthetic aperture radar，D-InSAR)由于具有全天时、高精度、似面测量等特点, 已成为当前开采沉陷监测的研究热点，但如何基于D-InSAR视线(line of sight, LOS)方向变形求取开采沉陷全部预计参数(概率积分参数)是困扰该项技术应用于矿山变形监测的难点。利用D-InSAR监测开采沉陷示意图见图 1，即当LOS方向与采动地表点A的空间移动矢量垂直时，D-InSAR LOS方向变形为零，完全无法监测A点变形；当LOS方向与采动地表点C的空间移动矢量斜交时，D-InSAR LOS方向变形仅能部分反映B点的空间移动；仅当LOS方向与采动地表点B的空间移动矢量平行时，D-InSAR LOS方向变形可以准确地反映C点的空间移动。因此，基于单视线向D-InSAR变形无法准确获取采动地表的三维移动变形，进而无法求取开采沉陷预计的概率积分参数。
Figure 1. Existing Problems in D-InSAR Monitoring of Mining Subsidence
解决上述难题，可以从两个途径去突破：①建立适用于开采沉陷特点的InSAR三维变形监测方法，并进行全面求参；②仍基于单视线向D-InSAR技术，寻找开采沉陷中下沉、水平移动与LOS方向移动变形的关系模型，基于非线性理论[1-2]实施沉陷预计参数反演。对于途径①, 通常的解决方案主要有多视线向D-InSAR[3-5]、D-InSAR+Offset Tracking[6-7]、D-InSAR+MAI(multiple-aperture interferometry)[8]、融合GPS和D-InSAR法[9-10]、地表先验模型+D-InSAR[11-12]，但上述方法主要用于解决地震、滑坡、火山喷发、冰川移动等特点的移动变形问题，且上述对象的移动机理不同于矿山开采沉陷，并不能妥善解决矿山开采地表的三维变形监测问题。文献[13]结合常规开采沉陷模型，利用水平移动与倾斜的关系，提出了基于单视线向D-InSAR的三维变形提取方法，但该方法仅适用于水平煤层的开采沉陷监测，对倾斜与极倾斜煤层、厚松散层等特殊地质采矿条件的开采沉陷监测并不适用。因而, 目前基于途径①进行概率积分求参仍需要进一步研究。对于途径②，研究表明，基于单视线向D-InSAR技术虽然无法直接监测开采沉陷三维移动变形，但LOS方向移动变形的本质为目标点的下沉、水平移动分量沿LOS方向投影之和。理论上可根据这种关系模型，结合目标区域地质采矿条件，融合概率积分法得到沉陷预计模型，利用非线性规划理论反演全部沉陷预计参数。
鉴于此，本文针对途径②重点开展研究，结合采动变形规律和雷达几何成像特点，首先构建融合单视线向D-InSAR和遗传算法(genetic algorithm，GA)的开采沉陷预计全参数反演模型，然后通过模拟实验开展可行性研究，最后将模型应用到兖州矿区南屯煤矿9310工作面沉陷预计参数反演。
HTML
-
根据概率积分法预计基本原理[14](见图 2)，在倾斜煤层中开采某一单元i引起的地表任意点S(x, y)的静态下沉函数模型为：
Figure 2. Principle of Subsidence Prediction of Arbitrary Point in Mining Face
式中，
ri为单元i处的主要影响半径，ri=Hi/tanβ，Hi为单元i处的采深；tanβ为主要影响角正切；θ为最大下沉角；(xi, yi)为单元i中心点的平面坐标；(x, y)为地表任意点S的平面坐标。
根据概率积分法原理，单元i开采引起地表S处沿x、y方向的水平移动函数模型为：
式(3)和式(4)中, b为水平移动系数；
其他参数含义同前。
假设将工作面按照图 2(b)沿走向、倾向剖分成足够小的正方形单元集合，剖分尺寸为Ai=0.1H×0.1H时(H为采深), 可满足工程精度需要[14]。若将任意正方形Ai视为开采单元，并假设该地质采矿条件下最大下沉值为W0，则基于叠加原理可获得工作面开采地表任意点的下沉W和沿x、y方向的水平移动U：
根据式(8)和式(9)可计算工作面开采地表点任意方向φ的水平移动：
-
综合D-InSAR LOS方向变形与三维变形关系、开采沉陷预计模型以及GA算法原理，提出了融合单视线向D-InSAR和GA的概率积分法预计全参数反演模型构建方法，其具体步骤如下：
1) 建立GA适应度函数。根据工作面地质采矿条件选取一组概率积分参数初值，利用开采沉陷预计模型(即式(7)至式(10))，计算目标区域任意像元j的下沉Wj、南北方向水平移动UjSN、东西方向水平移动UjEW。根据D-InSAR的LOS投影原理，可得到像元j的LOS方向预计移动变形r′jLOS为：
式中，ϕj为入射角; αj为卫星飞行方向方位角。
令目标像元j实测的LOS方向移动变形值为rjLOS，则LOS方向预计残差为：
假设选取D-InSAR变形场中n个目标像元参与概率积分参数求参，则根据式(12)可构造遗传算法适应度函数F：
式中，C为使适应度函数F始终大于零的一个系数。
2) 编码和种群生成。根据地质采矿条件和概率积分参数的经验关系确定参数(左拐点偏移距S1、右拐点偏移距S2、上拐点偏移距S3、下拐点偏移距S4、下沉系数q、主要影响角正切tanβ、最大下沉角θ、水平移动系数b)的范围，基于二进制编码规则随机生成各参数对应的二进制编码，建立初始种群。
3) 解码计算适应度函数值。将种群解码还原为各参数，基于解码参数和预计移动残差vj，计算适应度函数值。
4) 由适应度函数值计算个体被选择的概率。计算该个体的适应度与所有个体适应度总和之比，该值为该个体相对于种群的适应度，即被选到的概率。
5) 进行(轮盘赌)选择、(一点)交叉、变异操作，产生新一代种群。
6) 迭代计算。重复步骤2)~6)，直至达到适应度要求，即参数达到足够高的精度。遗传算法的迭代终止条件为满足种群迭代次数或者适应度达到某个要求。本次参数设定为：最大迭代次数maxG=100次，种群初始个体数p_num=100，交叉率p=0.95，变异率为0.05。
7) 最终解码得到最优概率积分参数。
基于上述思想，提出了如图 3所示的概率积分法全参数反演模型技术路线图，并编制了求参程序。
Figure 3. Technology Roadmap of Full Parameter Inversion Model Based on Probability Integral of Fusion of D-InSAR and GA
1.1. 开采沉陷下沉和水平移动预计模型
1.2. 融合单视线向D-InSAR和GA的开采沉陷预计求参模型构建
-
以淮南矿区煤系地层为背景，模拟工作面上覆岩层岩性为中硬，煤层倾角为6°，釆高M=3.0 m，釆深H=400 m，工作面开采尺寸沿倾向和走向水平投影尺寸为D1×D3=200 m×800 m，工作面倾向方位角为180°，顶板管理方法为长壁垮落法。工作面尺寸满足条件D1/H=0.5 < 1.2，D3/H=4 > 1.2，工作面采动程度属于倾向非充分采动，走向超充分采动，总体为非充分采动。工程实践表明，淮南矿区沉陷盆地形态符合概率积分法模型，假设模拟工作面的开采沉陷预计参数为：下沉系数q=0.8，水平移动系数b=0.3，主要影响角正切tanβ=2.0，最大下沉角θ=87°，拐点偏移距S1=S2=S3=S4=10 m。
-
以TerraSAR为例进行雷达LOS方向变形模拟，假设TerraSAR卫星回访目标区进行几何成像时，雷达卫星视线入射角为ϕ=44.494°，视线方向方位角α=350.596°。模拟D-InSAR监测的LOS方向采动变形时，首先在工作面上方地表沿南北和东西方向划分间隔为3 m的方格网，用格网点模拟分辨率为3 m(方位向和距离分辨率为3 m)的TerraSAR影像。根据概率积分法模型和选定的模型参数预计任意格网点的下沉值Wj、南北方向水平移动UjSN和东西方向水平移动UjEW，再根据D-InSAR的LOS三维投影原理(式(11))计算目标点沿LOS方向的累计移动量，并以此模拟D-InSAR的LOS方向实测地表累计移动量，模拟效果如图 4所示。
Figure 4. Simulation Diagram of Moving Deformation Along LOS Derection in Simulated Working Face
-
基于§2.2模拟的D-InSAR LOS方向变形值，利用本文建立的融合D-InSAR和GA的概率积分法预计全参数反演模型进行了概率积分参数反演实验，求参拟合效果如图 5所示，求参结果如表 1所示。
Figure 5. Parameters Fitting Effect Diagram of Simulated Experiment
名称 q tanβ b θ/(°) S1/m S2/m S3/m S4/m 参数真值 0.800 2.00 0.30 87 10 10 10 10 反演参数 0.782 2.03 0.28 87 8 12 9 8 绝对误差 0.018 0.03 0.02 0 -2 2 -1 -2 相对误差/% 2.3 1.5 6.7 0 -20 20 -10 -20 Table 1. Simulation Experimental Results of Parameters Fitting
从表 1可以看出，利用融合D-InSAR和GA方法反演的概率积分参数q、tanβ、b、θ的相对误差不超过6.7%，拐点偏移距S的相对误差不超过20%(参数敏感度低，对求参整体效果影响小)；求参拟合下沉误差为-5.90~6.10 mm，拟合中误差约为±2.20 mm。模拟实验表明，融合D-InSAR和GA方法能够有效准确地反演出全部概率积分预计参数。
2.1. 地质采矿条件概况
2.2. InSAR LOS方向变形模拟
2.3. 实验结果分析
-
实验区位于兖州煤田南屯煤矿九采区，其上方有邹济公路和矿用铁路穿过。雷达数据覆盖时间段为2011-12-25—2012-03-11，在此期间，南屯煤矿正在回3上煤层9310和9308工作面，9310北侧和9308南侧工作面均回采结束形成老采空区, 且地表沉陷基本稳定。9310和9308工作面煤层近水平(倾角平均约为8°)，工作面沿走向布置，采用综放采煤方法。9310工作面采厚4.90 m，工作面采深约588 m，工作面长197 m，实验期间推进尺寸约220 m；9308工作面采厚4.76 m，工作面采深约630 m，工作面长140 m，实验期间推进尺寸约415 m。
实验数据为TerraSAR-X获取的单视复数雷达影像，分辨率为3 m，影像大小为16 224× 10 710像素，实验数据覆盖了整个南屯煤矿。考虑到矿区变形梯度大易导致空间失相干特点，实验选取了2011-12-25—2012-03-11期间的时间基线为11 d和22 d的7景SAR影像进行实验。TerraSAR-X成像的主要几何参数如表 2所示。
干涉对 主影像名称 辅影像名称 间隔/d 垂直基线长度/m 主影像入射角/(°) 主影像方位向/(°) 1 20111225 20120105 11 -25.10 44.497 4 350.596 2 20120105 20120116 11 105.60 44.500 0 350.596 3 20120116 20120127 11 6.43 44.489 3 350.596 4 20120127 20120207 11 90.50 44.494 0 350.596 5 20120207 20120218 11 -141.72 44.483 1 350.596 6 20120218 20120311 22 -19.63 44.496 1 350.596 Table 2. Image Geometric Parameters of TerraSAR-X in the Experimental Area
-
对表 2中的6对主从影像分别进行差分干涉测量，获取了各干涉对沿LOS方向的地表移动变形量。由于各干涉对主影像LOS方位向相同, 入射角近似相等(可按相同考虑，误差可忽略)，利用二轨差分技术分别获得了各个时段LOS方向的移动变形, 如图 6所示；通过将各干涉对LOS方向移动变形量进行累加，获得了2011-12-25—2012-03-11时段内工作面开采引起的LOS方向地表累积移动变形场，如图 7所示。
Figure 6. Accumulated Amount of Surface Movement Deformation of LOS Direction in Each Period
Figure 7. Accumulated Amount of Surface Movement Deformation of LOS Direction During 2011-12-25—2012-03-11 Period
9310工作面于2012-02-28停采，截止到2012-03-11时，地表移动盆地基本稳定，可近似反演截止到2012-03-11时的移动盆地动态概率积分参数。由于9310工作面形状为矩形且煤层为近水平，开采沉陷盆地具有对称性，可沿走向西侧和倾向南侧两个方向分别取半条求参线C、D(见图 7)，并假设2011-12-25之前的工作面采动未对C、D线变形产生叠加影响。基于D-InSAR监测的LOS变形场，利用本文提出的融合D-InSAR和GA方法对9310工作面概率积分预计参数进行了反演。
结果表明，截止到2012-03-11，9310工作面开采沉陷的动态概率积分预计参数为：q=0.172，b=0.13，tanβ=2.08，θ=88°，S1=-10 m，S2=9 m，S3=-85 m，S4=30 m，拟合误差约为-35.00~45.00 mm，拟合中误差为±15.16 mm，拟合效果图如图 8所示。
Figure 8. Parameters Fitting Effect Diagram of Practical Application
从求参结果来看，9310上侧拐点偏移距达到了S3=-85 m(拐点偏移距离出现大幅度外偏)，这主要是因为9312工作面开采时在下山方向上存在较大的悬臂梁或砌体梁，当9310工作面开采时引起了悬臂梁或砌体梁失稳“活化”所致。其他概率积分参数基本正常。
3.1. 地质采矿条件及雷达数据概况
3.2. 概率积分参数反演实验及求参结果分析
-
针对利用常规D-InSAR结果无法求取开采沉陷全部预计参数(特别是水平移动参数b)的问题，本文开展了融合单视线向D-InSAR和GA的概率积分预计全参数反演方法研究。在煤矿采动影响区，D-InSAR LOS方向移动实质为下沉、南北和东西方向水平移动沿LOS方向投影之和，基于GA原理，首次构建了融合D-InSAR和GA的概率积分法预计全参数反演模型，并编制了求参软件。模拟实验结果表明，融合D-InSAR和GA方法能够有效准确地反演出全部概率积分预计参数。利用融合D-InSAR和GA的概率积分法预计全参数反演模型，在一定的假设条件下，对2012年3月11日9310工作面开采沉陷的动态概率积分参数进行了求取，并认为S3=-85 m主要为9310工作面开采导致邻近9312采空区边界的悬臂梁或砌体梁失稳“活化”所致。
DownLoad:
