2021, Vol. 41
表1(Table 1)
2. 内蒙古煤炭地质勘查(集团)测绘院有限公司，呼和浩特市展东路25号，010052
矿山地下煤炭资源开采会造成土地塌陷及采空区积水，进而导致耕地面积减小和农作物产量受损，地表房屋、建筑物及交通运输道路扭曲变形，严重时还会导致山体开裂、崩塌、滑坡、泥石流等重大自然灾害。合成孔径雷达差分干涉测量(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) |
式中，
根据随机介质理论推导出的水平移动与倾斜成正比的关系，对于相同的采深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)对应计算过程中用到的真实形变矩阵，这对重建真实水平方向的形变尤为重要。
|
图 1 矿山三维地表形变曲线 Fig. 1 Curves of three-dimensional multi-level mining deformation |
利用矿山地表监测的水平移动数据、下沉数据、采深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)；dW(i, j)为垂直向形变，dE(i, j)为东西向形变，dN(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)主要为：
|
图 2 矿山三维多量级形变解算流程 Fig. 2 Calculating flowchart of three-dimensional multi-level mining deformation |
1) 输入覆盖研究区的高分辨率InSAR影像和数字高程模型(digital elevation model，DEM)数据。
2) 对时序InSAR影像对进行DInSAR处理，分析相干图和干涉图，得到形变区域的平均相干值
3) 当
4) 将矿山地表形变点(i, j)对应的LOS向形变值DLOS(i, j)、α、θ、b、H及tanβ代入式(3)和式(4)，利用式(5)和式(6)进行回代处理可得dW(i, j)。
5) 将dW(i, j)代入式(3)可得dN(i, j)和dE(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。
|
图 3 工作面位置 Fig. 3 Location of the workface |
选用TerraSAR-X(TSX)雷达卫星的单视斜距复影像数据产品，雷达波入射角为42.43°，轨道方位角为189.53°，像元尺寸为0.91×0.86。表 1为TSX影像参数，为去除地形相位的影响，DEM采用SRTM的30 m分辨率数据产品。
|
表 1 TSX影像参数 Tab. 1 The parameters of the TSX images |
利用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向最大形变dLOS, M约为-4.24 m[7]。为了滤除OT监测结果中的低精度测量值，
|
图 4 2012-12-30~2013-06-07大柳塔矿区DInSAR和OT观测结果 Fig. 4 Results of DInSAR and OT in Daliuta mining area from 2012-12-13 to 2013-06-07 |
图 5为52303工作面开采造成的矿山地表LOS向多量级形变，其中图 5(a)~5(f)为DInSAR测量的小量级形变与OT测量的大量级形变的融合结果。可以看出，多量级形变不连续且存在较多空值点，这些点多是由DInSAR测量结果中的失相干点和OT测量结果中低精度测量点构成。图 5(g)~5(l)为利用IDW估计后的LOS向连续多量级形变测量结果。可以看出，沉降的方向与该工作面由东向西的开采方向一致，沉降范围与开采周期直接相关。但由于LOS向的测量结果只能反映雷达视线向的形变方向和大小，不能精确获取三维方向的多量级形变，基于LOS向的多量级形变还需重建矿山三维多量级形变。
|
(a) 和(g) : 2012-12-13 ~2013-01-04; (b) 和(h) : 2013-01-04~2013-02-06; (c) 和(i) : 2013-02-06 ~2013-03-22; (d) 和(j) : 2013-03-22 ~2013-04-02; (e) 和(k) : 2013-04-02~2013-05-05; (f) 和(l) : 2013-05-05 ~2013-06-07 图 5 2012-12-30~2013-06-07大柳塔矿区LOS向多量级形变融合结果 Fig. 5 Multi-level deformation fusion results in Daliuta mining area from 2012-12-13 to 2013-06-07 |
图 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。
|
(a) ~ (c) : 2012-12-13 ~2013-01-04; (d) ~ (f) : 2013-01-04 ~02-06; (g) ~ (i) : 2013-02-06 ~03-22; (j) ~ (l) : 2013-03-22 ~04-02; (m) ~ (o) : 2013-04-02~05-05; (p) ~ (r) : 2013-05-05 ~06-07; 第1列为垂直向形变, 第2列为东西向形变, 第3列为南北向形变 图 6 2012-12-13~2013-06-07矿山时序三维多量级形变监测结果 Fig. 6 Three-dimensional multi-level mining deformation from 2012-12-13 to 2013-06-07 |
分析每组影像对的三维多量级形变最大值后发现，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的时间间隔。
|
图 7 2012-12~2013-06累积形变监测结果 Fig. 7 The cumulative deformation from December2012to June 2013 |
用本文方法计算水准观测线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，满足矿山三维多量级形变监测的精度要求。
|
图 8 本文方法监测结果与水准测量结果对比 Fig. 8 Comparison of the deformation obtained by proposed method and leveling |
针对矿山地表三维多量级形变获取难的问题，本文提出一种矿山三维多量级形变重建方法。利用该方法对研究区2012-12-13~2013-06-07期间7景X波段的TerraSAR-X数据进行时序处理，获取大柳塔矿山52303工作面的时序三维多量级形变，与水准测量结果对比可知，垂直向、东西向及南北向的RMSE分别小于0.1880m、0.2243m和0.2074m，监测结果与水准测量数据形变趋势一致。本文提出的方法可满足监测矿山沉降漏斗边缘地带cm级到沉降中心m级的多量级形变监测需求，将为矿山三维多量级形变测量提供一种大范围、高分辨率、高精度、三维方向、低成本的监测技术，同时为InSAR技术在矿山灾害监测应用的推广提供技术保障。
| [1] |
Massonnet D, Feigl K L. Radar Interferometry and Its Application to Changes in the Earth's Surface[J]. Reviews of Geophysics, 1998, 36(4): 441-500 DOI:10.1029/97RG03139
(  1) |
| [2] |
孙赫, 陈巍然, 牛玉芬, 等. 基于InSAR的矿区沉降中心动态变化监测与分析[J]. 大地测量与地球动力学, 2020, 40(3): 276-280 (Sun He, Chen Weiran, Niu Yufen, et al. Monitoring and Analysis of Dynamic Change of Mining Subsidence Center Based on InSAR[J]. Journal of Geodesy and Geodynamics, 2020, 40(3): 276-280)
( 1) |
| [3] |
Michel R, Avouac J P, Taboury J. Measuring Ground Displacements from SAR Amplitude Images: Application to the Landers Earthquake[J]. Geophysical Research Letters, 1999, 26(7): 875-878 DOI:10.1029/1999GL900138
( 1) |
| [4] |
Chen B Q, Deng K Z, Fan H D, et al. Combining SAR Interferometric Phase and Intensity Information for Monitoring of Large Gradient Deformation in Coal Mining Area[J]. European Journal of Remote Sensing, 2015, 48(1): 701-717 DOI:10.5721/EuJRS20154839
( 1) |
| [5] |
Fan H D, Gao X X, Yang J K, et al. Monitoring Mining Subsidence Using a Combination of Phase-Stacking and Offset-Tracking Methods[J]. Remote Sensing, 2015, 7(7): 9166-9183 DOI:10.3390/rs70709166
( 0) |
| [6] |
Wang Z W, Yu S W, Tao Q X, et al. A Method of Monitoring Three-Dimensional Ground Displacement in Mining Areas by Integrating Multiple InSAR Methods[J]. International Journal of Remote Sensing, 2018, 39(4): 1199-1219 DOI:10.1080/01431161.2017.1399473
( 0) |
| [7] |
Zhang L L, Cai X X, Wang Y, et al. Long-Term Ground Multi-Level Deformation Fusion and Analysis Based on a Combination of Deformation Prior Fusion Model and OTD-InSAR for Longwall Mining Activity[J]. Measurement, 2020, 161: 107911 DOI:10.1016/j.measurement.2020.107911
( 3) |
| [8] |
Ou D P, Tan K, Du Q, et al. Decision Fusion of D-InSAR and Pixel Offset Tracking for Coal Mining Deformation Monitoring[J]. Remote Sensing, 2018, 10(7): 1055 DOI:10.3390/rs10071055
( 1) |
| [9] |
Gray L. Using Multiple RADARSAT InSAR Pairs to Estimate a Full Three-Dimensional Solution for Glacial Ice Movement[J]. Geophysical Research Letters, 2011, 38(5): 132-140
( 1) |
| [10] |
Hu J, Li Z W, Ding X L, et al. Resolving Three-Dimensional Surface Displacements from InSAR Measurements: A Review[J]. Earth-Science Reviews, 2014, 133(2): 1-17
( 0) |
| [11] |
祝传广, 邓喀中, 张继贤, 等. 基于多源SAR影像矿区三维形变场的监测[J]. 煤炭学报, 2014, 39(4): 673-678 (Zhu Chuanguang, Deng Kazhong, Zhang Jixian, et al. Three-Dimensional Deformation Field Detection Based on Multi-Source SAR Imagery in Mining Area[J]. Journal of China Coal Society, 2014, 39(4): 673-678)
( 0) |
| [12] |
范洪冬, 高晓雄, 邓喀中, 等. 基于多轨道SAR的老采空区地表三维形变监测[J]. 采矿与安全工程学报, 2017, 34(6): 1156-1161 (Fan Hongdong, Gao Xiaoxiong, Deng Kazhong, et al. Detection of the Three-Dimensional Surface Deformation above the Abandoned Goaf Based on SAR[J]. Journal of Mining and Safety Engineering, 2017, 34(6): 1156-1161)
( 1) |
| [13] |
Li Z W, Yang Z F, Zhu J J, et al. Retrieving Three-Dimensional Displacement Fields of Mining Areas from a Single InSAR Pair[J]. Journal of Geodesy, 2015, 89(1): 17-32 DOI:10.1007/s00190-014-0757-1
( 2) |
| [14] |
Yang Z F, Li Z W, Zhu J J, et al. Deriving Time-Series Three-Dimensional Displacements of Mining Areas from a Single-Geometry InSAR Dataset[J]. Journal of Geodesy, 2018, 92(5): 529-544 DOI:10.1007/s00190-017-1079-x
( 0) |
| [15] |
Yang Z F, Li Z W, Zhu J J, et al. Retrieving 3-D Large Displacements of Mining Areas from a Single Amplitude Pair of SAR Using Offset Tracking[J]. Remote Sensing, 2017, 9(4): 338
( 1) |
| [16] |
Chen B Q, Li Z H, Yu C, et al. Three-Dimensional Time-Varying Large Surface Displacements in Coal Exploiting Areas Revealed through Integration of SAR Pixel Offset Measurements and Mining Subsidence Model[J]. Remote Sensingof Environment, 2020, 240: 111663 DOI:10.1016/j.rse.2020.111663
( 2) |
2. Surveying and Mapping Institute of Inner Mongolia Coal Geological Exploration (Group) Co Ltd, 25 Zhandong Road, Hohhot 010052, China