2019, Vol. 39
2. 湖南科技大学资源环境与安全工程学院，湖南省湘潭市桃园路2号，411201;
3. 武汉大学测绘学院，武汉市珞喻路129号，430079
鄂尔多斯市的红庆河矿区位于人口密度低、人工建筑物少的沙漠、半沙漠地区，若采用PS-InSAR技术[1]监测该地区地面沉降，将难以获得足够密度的点。本文拟利用DS-InSAR技术[2]对覆盖该地区的33景Sentinel-1A数据进行处理，以获取更详细的形变信息，并与PS-InSAR结果进行对比分析，证明DS-InSAR技术在矿区形变监测领域的优势[3-5]。
1 DS-InSAR基本原理DS提取技术通过比较干涉相位时间维样本的相似度，选取统计学上具有相同散射特性的像素进行同质滤波，并在复数域对其相位进行优化估计，可有效提高干涉图相位的质量，有利于相位解缠和大气信号分离，提高每个目标的时间序列位移估计精度，且相较于其他时序处理方法，提高了干涉相位信息的利用率[6]。与PS不同，DS是分辨率单元中所有较小散射体回波信号的相干累积。
DS-InSAR技术主要包括空间维的同质滤波技术和时间维的缠绕相位平差技术[7-8]。Ferretti等[2]提出利用Kolmogorov-Smirnov(KS)检验基于振幅数据选取统计同质像元，假设给定2个向量d(P1)和d(P2)，如果在特定的显著水平下，来自同一概率分布函数的像素点P1和P2零假设不被拒绝，则2个像素在统计上具有一致性。实际数据处理过程中会以每幅影像的每一个像素点为中心，以合适的估计窗口利用统计检验对其统计一致性进行判别。在N幅SAR图像中，选出1景主影像，将剩余(N－1)幅辅影像在主图像上进行重采样，则：
| $ \boldsymbol{d}(P)=\left[d_{1}(P), d_{2}(P), \cdots, d_{N}(P)\right]^{\mathrm{T}} $ | (1) |
式中，d为复数向量，P为任意图像像素，di(P)为像素点P在干涉图像堆栈中第i幅图像的复反射率值。复数向量d的振幅值| d |序列可以用累积分布函数的无偏估计SN(| d |)表示：
| $ {S_N}(X) = \left\{ {\begin{array}{*{20}{l}} {0, |\mathit{\boldsymbol{d}}| < {d_1}}\\ {\frac{k}{N}, {x_k} < |\mathit{\boldsymbol{d}}| < {d_1}}\\ {1, |\mathit{\boldsymbol{d}}| \ge {d_N}} \end{array}} \right. $ | (2) |
式中，xi为振幅排序序列的第i个元素。为判断干涉图像上P1、P2点是否具有统计一致性，计算其累积分布函数SNP1和SNP2绝对差的最大值DN：
| $ {D_N} = \sqrt {N/2} {\sup _{x \in R}}\left| {S_N^{{P_1}}(|\mathit{\boldsymbol{d}}|) - S_N^{{P_2}}(|\mathit{\boldsymbol{d}}|)} \right| $ | (3) |
式中，DN的分布可以用KS分布近似替代，DN的累积分布函数为：
| $ P\left( {{D_N} \le t} \right) = H(t) = 1 - 2\sum\limits_{n = 1}^\infty {{{( - 1)}^{n - 1}}} {{\rm{e}}^{ - 2{n^2}{t^2}}} $ | (4) |
从式(4)可以看出，其结果不受向量具体分布函数的影响。c为阈值，当DN≤c时，认为2个向量属于同一分布，定义c的固定显著水平α为：
| $ \alpha=1-H(t) $ | (5) |
KS检验得到同质目标后，需要根据同质目标的协方差矩阵进行参数估计。对于SAR图像中的像素点，其协方差矩阵A的数学表达式为:
| $ \mathit{\boldsymbol{A}} = E\left[ {\mathit{\boldsymbol{y}}{\mathit{\boldsymbol{y}}^{\rm{H}}}} \right] \approx \frac{1}{{{N_P}}}{\sum\limits_{y \in \mathit{\Omega} } {\mathit{\boldsymbol{yy}}} ^{\rm{H}}} $ | (6) |
式中，E为期望，
| $ \mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} 1&{{{\tilde \gamma }_{1, 2}}{{\rm{e}}^{{\rm{j}}{\phi _{1, 2}}}}}& \cdots &{{{\tilde \gamma }_{1, 2}}{{\rm{e}}^{{\rm{j}}{\phi _{1, N}}}}}\\ {{{\tilde \gamma }_{1, 2}}{{\rm{e}}^{{\rm{j}}{\phi _{2, 1}}}}}&1& \cdots &{{{\tilde \gamma }_{1, 2}}{{\rm{e}}^{{\rm{j}}{\phi _{2, N}}}}}\\ \vdots & \vdots & \vdots & \vdots \\ {{{\tilde \gamma }_{1, 2}}{{\rm{e}}^{{\rm{j}}{\phi _{N, 1}}}}}&{{{\tilde \gamma }_{1, 2}}{{\rm{e}}^{{\rm{j}}{\phi _{N, 2}}}}}& \cdots &1 \end{array}} \right] = |\mathit{\boldsymbol{T}}|\phi $ | (7) |
对相干矩阵特征进行分解，可获取不同散射机制的最佳相位[6]。一般认为，前n个特征矩阵为有效相位，剩余的(N－n)个为相位噪声，则：
| $ \mathit{\boldsymbol{T}} = \sum\limits_{i = 1}^N {{\mathit{\boldsymbol{\lambda }}_i}} {T_i} + \sum\limits_{i = n + 1}^N {{\mathit{\boldsymbol{\lambda }}_i}} {T_i} = {\mathit{\boldsymbol{T}}_{{\rm{signal }}}} + {\mathit{\boldsymbol{T}}_{{\rm{noise }}}} $ | (8) |
主导散射机制T1作为形变信号，特征向量相位λ1用于估计DS的地面变形，理论上最大特征值所占比例越高，相位估计结果越可靠。一般认为，最大特征值对应的特征向量为其相位最优解。DS提取完成后，再联合PS对展开的干涉相位进行2D线性回归分析，校正外部DEM带入的误差，然后利用时空特性的不同去除残余的大气相位和非线性形变，最后迭代求出最终的形变速率[8]。
2 实验数据及结果分析 2.1 研究区介绍红庆河煤矿所在地区地质结构稳定性好，但顶板是砂质泥岩，稳定性较差。煤矿共有10个可采煤层，煤层平均厚度为6.9 m，倾角为1°~3°，埋深700 m左右，采用综合机械化一次采全高采煤法[9-10]。矿区分为南翼、西翼、北翼3个盘区，2016-11南翼开始开采，2017-03北翼开始开采。到2018-04的开采进度如图 1所示，其中，白色矩形区域为规划巷道，斜线填充区域为已经掘进的巷道，箭头为开采起点及开采方向，左下角为南翼，累计推进3 210 m，右上角为北翼，累计推进2 489 m。
|
图 1 红庆河煤矿采掘工程平面图 Fig. 1 Hongqinghe coal mine mining project plan |
Sentinel-1A卫星是C波段合成孔径雷达，卫星轨道为太阳同步近极地轨道，轨道高度为693 km, 倾角为98.15°，卫星的重返周期为12 d，理论上卫星沉降监测的最大速度为42 cm/a[11]。本文选用覆盖红庆河煤矿的33景Sentinel-1A数据进行干涉处理。由于Sentinel-1A卫星进行干涉处理时方位向配准精度需要达到0.001像元，为满足配准要求，本文在精密轨道几何配准基础上进行增强谱分集配准，并以多主影像的方式生成干涉图[12-13]，影像时间范围为2017-03-12~2018-04-24，干涉对空间基线范围小于100 m，时间基线小于50 d，如图 2所示。时序处理时选取覆盖研究区的burst进行DS-InSAR处理，原始差分干涉图及DS-InSAR处理后的差分干涉图见图 3。从图中可以看出，经过同质滤波后的干涉图(图 3(b))质量明显好于原始干涉图(图 3(a))。
|
图 2 时空基线组合 Fig. 2 Space-time baseline combination |
|
图 3 差分干涉图 Fig. 3 The differential interferometry |
图 4和5中椭圆区域是红庆河煤矿。图 4为采用PS-InSAR技术获取的沉降结果，图中PS分布大致均匀但十分稀疏，密度为0.6个/km2，不能满足Colesanti等[14]提出的最小有效密度(3~4个/km2)，难以获取实验区沉降的具体信息。图 5为采用DS-InSAR获取的结果，图中矿区沉降范围大致呈椭圆形，沉降区边界明显，DS的密度为12.8个/km2，远大于PS的密度。红庆河煤矿从2016-11开始开采，虽然开采年限较短，但煤层顶板岩性多为砂质泥岩，结构较软，顶板稳定性差，在开采过程中易垮塌，因此沉降区中心影像干涉图出现严重失相干。DS-InSAR结果显示，沉降区呈月牙形，这是由于沉降中心无有效结果。另外，由图可知，红庆河煤矿累积沉降面积约为32.3 km2，最大沉降速率为41 mm/a，但实际塌陷中心沉降速度应远大于该值。
|
图 4 PS-InSAR技术获取的红庆河煤矿沉降结果 Fig. 4 Settlement results of Hongqinghe coal mine obtained by PS-InSAR |
|
图 5 DS-InSAR获取的红庆河煤矿年沉降结果 Fig. 5 Settlement results of Hongqinghe coal mine obtained by DS-InSAR |
图 6为研究区局部放大图，其中PS-InSAR结果叠加了已开采区域的采掘工程平面图。可以看出，沉降区呈NE-SW走向，沉降范围与煤矿采掘工程平面图的空间分布相吻合。一般情况下，水平煤矿开采过程中，采矿区上方会逐渐形成比采空区略大的椭圆形地表移动盆地，盆地中间沉降值最大，四周沉降逐渐减小[15]，由于煤矿南翼开采长度比北翼多721 m，故南翼生成的沉降面积更大。从结果可以看出，未来沉降趋势大致符合采空区地表移动规律。红庆河煤矿附近人工建筑稀少，且该地区地质结构稳定，故认为煤炭资源开采是引发红庆河矿区沉降的主要原因。
|
图 6 时序InSAR局部图 Fig. 6 The local figure of time series SAR interferometry |
本文对覆盖红庆河煤矿的Sentinel-1A数据进行DS-InSAR处理，并与PS-InSAR结果进行对比分析，证明DS-InSAR技术更适合煤矿区形变监测，获取的煤矿累积沉降面积约为32.3 km2。DS-InSAR技术克服了时间去相关和大气干扰对干涉结果的影响，实现了大面积长时间序列的监测，为煤矿开采对周边影响的研究提供了一定的参考，同时也为地面沉降监测布点提供了有效依据。另外，Sentinel-1A卫星具有重返周期短的优势，因而可以获得相干性较好的干涉图，能为矿区沉降监测提供较好的数据源。
致谢: 感谢天地科技股份有限公司开采设计事业部煤炭科学研究总院开采研究分院的冯友良博士提供红庆河煤矿的数据资料。
| [1] |
Ferretti A, Prati C, Rocca F. Nonlinear Subsidence Rate Estimation Using Permanent Scatterers in Differential SAR Interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2000, 38(5): 2 202-2 212 DOI:10.1109/36.868878
(  0) |
| [2] |
Ferretti A, Fumagalli A, Novali F, et al. A New Algorithm for Processing Interferometric Data-Stacks: SqueeSAR[J]. IEEE Transactions on Geoscience and Remote Sensing, 2011, 49(9): 3 460-3 470 DOI:10.1109/TGRS.2011.2124465
( 0) |
| [3] |
Graniczny M, Colombo D, Kowalski Z, et al. New Results on Ground Deformation in the Upper Silesian Coal Basin(Southern Poland) Obtained during the DORIS Project(EU-FP 7)[J]. Pure and Applied Geophysics, 2015, 172(11): 3 029-3 042 DOI:10.1007/s00024-014-0908-6
( 0) |
| [4] |
Du Z Y, Ge L L, Ng A H M, et al. An Innovative Distributed Scatterer Based Time-Series InSAR Method over Underground Mining Region[C]. IEEE International Geoscience and Remote Sensing Symposium, Fort Worth, 2017
( 0) |
| [5] |
Esfahany S S. Exploitation of Distributed Scatterers in Synthetic Aperture Radar Interferometry[D]. Delft: Delft University of Technology, 2017
( 0) |
| [6] |
汪慧.星载合成孔径雷达相位萃取算法及应用研究[D].重庆: 重庆大学, 2017 (Wang Hui. Study on Phase Reconstruction Algorithm of Spaceborne Synthetic Aperture Radar and Its Application[D]. Chongqing: Chongqing University, 2017)
( 0) |
| [7] |
Cao N, Lee H, Jung H C. A Phase-Decomposition-Based PS-InSAR Processing Method[J]. IEEE Transactions on Geoscience and Remote Sensing, 2016, 54(2): 1 074-1 090 DOI:10.1109/TGRS.2015.2473818
( 0) |
| [8] |
Guarnieri A M, Tebaldini S. On the Exploitation of Target Statistics for SAR Interferometry Applications[J]. IEEE Transactions on Geoscience and Remote Sensing, 2008, 46(11): 3 436-3 443 DOI:10.1109/TGRS.2008.2001756
( 0) |
| [9] |
冯友良.大断面煤巷开挖卸荷帮部破坏机制与控制技术研究[D].北京: 中国矿业大学, 2017 (Feng Youliang. Study on Rib Breakage Mechanism and Control Technology of Large Cross-Section Coal Roadway for Excavation Unloading[D]. Beijing: China University of Mining and Technology, 2017)
( 0) |
| [10] |
杨晓杰, 毛文彬, 张星宇, 等. 红庆河煤矿厚煤层巷道切顶卸压关键参数研究[J]. 中国煤炭, 2018, 44(10): 76-80 (Yang Xiaojie, Mao Wenbin, Zhang Xingyu, et al. Research on Key Parameters of Roof Cutting and Pressure Releasing of Roadway in Thick Coal Seam of Hongqinghe Mine[J]. China Coal, 2018, 44(10): 76-80 DOI:10.3969/j.issn.1006-530X.2018.10.014)
( 0) |
| [11] |
吴岳.TOPS模式Sentinel-1数据地面沉降监测研究[D].北京: 中国矿业大学, 2017 (Wu Yue. Study on Ground Subsidence by Sentinel-1 TOPS Mode Data[D]. Beijing: China University of Mining and Technology, 2017)
( 0) |
| [12] |
吴文豪, 周志伟, 李陶, 等. 精密轨道支持下的哨兵卫星TOPS模式干涉处理[J]. 测绘学报, 2017, 46(9): 1 156-1 164 (Wu Wenhao, Zhou Zhiwei, Li Tao, et al. A Study of Sentinel-1 TOPS Mode Co-Registration[J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(9): 1 156-1 164)
( 0) |
| [13] |
吴文豪.哨兵雷达卫星TOPS模式干涉处理研究[D].武汉: 武汉大学, 2016 (Wu Wenhao. TOPS Interferometry with Sentinel-1[D]. Wuhan: Wuhan University, 2016) http://cdmd.cnki.com.cn/Article/CDMD-10486-1016124261.htm
( 0) |
| [14] |
Colesanti C, Ferretti A, Prati C, et al. Monitoring Landslides and Tectonic Motions with the Permanent Scatterers Technique[J]. Engineering Geology, 2003, 68(1-2): 3-14 DOI:10.1016/S0013-7952(02)00195-3
( 0) |
| [15] |
盛耀彬.基于时序SAR影像的地下资源开采导致的地表形变监测方法与应用[D].北京: 中国矿业大学, 2011 (Sheng Yaobin. Method for Land Deformation Monitoring Caused by Underground Resources Mining Activities Based on Time-series SAR Images and Its Application[D]. Beijing: China University of Mining and Technology, 2011) http://cdmd.cnki.com.cn/article/cdmd-10290-1011281027.htm
( 0) |
2. Taoyuan Road, Xiangtan 411201, China 2 School of Resources Environment and Safety Engineering, Hunan University of Science and Technology, 2 Taoyuan Road, Xiangtan 411201, China;
3. School of Geodesy and Geomatics, Wuhan University, 129 Luoyu Road, Wuhan 430079, China