基于InSAR技术和GS-SVR算法的矿区地表开采沉陷预计

马 飞，隋立春，姚顽强，汤伏全

(1.长安大学 地质工程与测绘学院,陕西 西安 710054;2.西安科技大学 测绘科学与技术学院,陕西 西安 710054;3.地理国情监测国家测绘地理信息局 工程技术研究中心,陕西 西安 710054)

煤炭开采引起的矿区地表沉陷灾害已经引起了人们越来越多的关注。然而矿区沉降的常规监测方法大多是重复采集被测地区的数据[1-5]，通过计算得出被监测地区在不同时期的沉降速率。但是西部黄土高原地区，土质疏松，塬梁破碎，沟壑纵横，难以寻找稳定的控制点，使得常规的监测方法耗费人力、物力巨大，效率低下。

近年来，随着卫星遥感领域的不断发展，随着大量的商用SAR数据的获取，合成孔径雷达理论技术也在不断地完善，差分雷达干涉测量技术(D-InSAR)在地表形变监测方面的研究应用得到很大发展。D-InSAR技术所获取的不是离散点的信息，而是大面积区域连续的地表形变信息，与传统测量技术相比，其覆盖范围大、成本低、空间分辨率高，可以全天候工作，已成为获取地表形变信息的新手段[6-10]。至今，刘广、黄宝伟、范洪冬等学者利用该技术在城市地表形变监测、矿区沉陷监测等方面获得了大量成果[13-15]。

支持向量机回归(SVR)是支持向量在函数回归领域的应用，它可以高效处理小样本、非线性、高维数据问题。相比于最小二乘法、灰色模型、神经网络等算法具有更高的运算效率和预测精度。因此，本文通过InSAR技术获得彬长矿区的沉降区域范围，监测随着时间推移矿区地表累计沉降变化量。将沉降结果作为SVR算法的训练学习样本，建立矿区地表沉降预测模型，利用格网搜索法(GS)优化选取模型参数，对矿区地表沉降进行预测估计；结果表明，InSAR技术能够准确获取矿区的沉降范围，监测得到的沉降值与常规水准数据用于矿区沉陷预计可以获取相当的预测精度。

1 InSAR技术获取矿区地表沉降量

InSAR(合成孔径雷达干涉测量技术)是以波的干涉为基础，用卫星两次飞过同一地区(重复轨道方式)上空所获得的两幅微波图像，两幅图像满足相干条件，对其进行相位干涉处理，产生干涉条纹，反映出相位的变化，这种图像叫做干涉图。理想状态，如果地面没有变形或其它影响，通过解缠处理，解算出每一点的相位值，从而计算得出地面点到雷达的斜距以及地面点的高程。

二轨法首先利用一对跨越形变期的SAR图像进行干涉处理(见图1)，获得的干涉相位可以表示为

Δφ=Δφtopo+Δφdef+

Δφflat+Δφatm+Δφnoise.

(1)

上式中Δφ表示干涉图的缠绕相位，Δφtopo表示地形相位，Δφdef表示形变相位，Δφflat表示平地相位，Δφatm表示大气影响，Δφnoise表示噪声相位(包括系统热噪声、时间与空间失相关噪声等)。为了得到准确的形变相位Δφdef，(1)式右边其余各项应当逐项消除。利用二轨差分法结合已有的外部DEM模拟地形信息从而实现地形相位的去除[1]。大气影响Δφatm是最主要的误差源之一，尤其是电离层延迟误差对L波段的SAR数据干涉测量影响不容忽视[3]，本文利用GPS数据进行电离层延迟误差改正；平地相位Δφflat和噪声相位Δφnoise分别通过基线估计和自适应滤波进行去除。各项误差消除之后通过相位解缠得到形变相位Δφdef，进而计算出地表的形变量。

图1 InSAR原理图

由图1中几何关系推导可知，两幅影像的相位差为

(2)

其中λ为雷达波长，ΔRd为地形变化在卫星视线向上的投影值；

由已知DEM反演得到的地形相位为

(3)

得到形变相位

(4)

式(4)给出了差分相位对地形形变的敏感度，ENVISAT/ASAR C波段，波长为5.6 cm，一般情况下，2.8 cm的斜距向变化即可引起一个2π的相位变化，ALOS/PALSAR L波段，波长为23.5 cm，11.7 cm的斜距向变化即可引起一个2π的相位变化。

在二轨法差分中，地形误差对差分相位的影响主要取决于外部DEM的精度。

(5)

由上式可知，地形误差对干涉相位的影响为

(6)

其中B⊥为基线垂直于卫星视线方向的分量，θ为主图像视角。

2 基于GS-SVR的预计模型基本思想

SVR的基本思想是通过一个非线性函数映射将数据转换到高维特征空间，然后对其进行线性回归处理，转化为求解高维特征空间的最优决策函数：

(7)

其中，X→Rn,w∈f，b为阈值。

上述问题通过经验风险和VC维理论分析转化为

(8)

式中，C为惩罚因子，ζ为误差，ε为损失函数参数。

损失函数参数ε用于控制支持向量的个数和泛化能力，取值越小，精度越低，则支持向量越少，为了达到最优的拟合效果一般取值为(0.000 1～0.01)；惩罚因子C用于控制模型的复杂度，一般取值为(1～1 000)。为了选择合适的预计模型参数，前人提出了多种参数优化选取方法，本文采用网格优化算法(GS)对模型参数进行寻优。首先确定C和ε的初始值，然后基于网格法全局搜索，获取最优的C和ε值，确定预测方程。

根据于广明[23]、于学义[24]等人的研究成果，通过大量实测资料分析矿区地表沉陷，发现沉陷非线性机理导致地表点的下沉过程在时间上沉陷呈不规则性，地表点下沉量在相等时间内大小不等呈非线性。本文通过InSAR技术获取一系列与时间相关的沉降值。这些沉降值表现为在时间上非线性关系：{xi}={x1,x2,…,xn}，取前n-p个数据作为GS-SVR算法的训练数据，构造如式(7)所示的预计函数，通过对式(8)进行计算分析，求解得到预计函数；对p组数据进行预测，将其与实测值比较分析，从而确定预测模型的精度。

3 实验研究

3.1 实验准备

彬长矿区位于陕西省关中西北部长武和彬县境内，是国家规划的黄陇基地的主力矿区之一。本区地处渭北黄土高原塬梁沟壑区，地势从黄土塬梁向中间泾河谷地倾斜。矿区水土流失较为严重，土壤主要是黑垆土、黄绵土、红土、淤土、潮土等，植被类型以阔叶落叶灌丛和草本植被为主。矿区东西长46 km，南北宽36.5 km，规划面积978 km2，地质储量为67.29亿t，整个矿区生产能力达5 000万t。其中，大佛寺煤矿2006年建成，大佛寺40301工作面为首采区，开采煤层厚度平均为11 m。

本文以彬长矿区工作面为例，选取5景PALSAR数据，组成干涉对进行方法验证，为减少时间去相关的影响，选取的干涉对的时间间隔尽量小。干涉对组成情况如表1所示。

表1 ALOS PALSAR影像对基本参数

3.2 矿区沉降预测

本文数据处理采用二轨法差分干涉的方法，包括影像的预处理、主辅影像配准、生成干涉图、干涉图滤波、去平地效应、相位解缠、基线参数计算、去地形相位、生成差分干涉图、地理编码等环节，最终得到相干性分析图、相位干涉图、差分干涉相位图、形变图，具体处理流程见图2。按照表1中的干涉对可获取4组沉降图，以第一组为参考将后续得到的沉降图依次叠加即可得到在监测时间段内的累积沉降量，具体如图3所示。

由图3可知，随着时间的推移，矿区的累积沉降量越来越大，逐渐形成为一个沉陷盆地， 本文选取沉降中心的一个沉降点(O)以及沿走向线(T1、T2)和倾向线(Q1、Q2)方向各两个点为研究对象，对第3节提出的方法进行验证。取前4组数据作为测试样本，对第5组沉降值进行预测，选取的沉降信息如表2所示。按照式(7)和式(8)建立预测函数，预测结果如表3所示：表3中MAE表示绝对误差，MRE表示相对误差。

由表3可知，预测结果与GPS技术测得的结果较为一致，最大绝对误差为3 mm，最大相对误差为5.9%，其预测精度满足工程实践的需求，将InSAR技术与SVR算法相结合可以用于矿区沉降预测的实际应用。分析表3可知在沉降量大的区域其预测精度相对较高，所以在绝对误差相等时，O点处的相对误差最小。

图2 二轨法差分干涉处理流程图

图3 时间序列沉降图

获取时间时间间隔/dO/mmT1/mmT2/mmQ1/mmQ2/mm2007年7月0000002007年8月46-15-10-12-12-102007年10月46-26-17-20-18-142007年11月46-57-37-39-37-322008年1月46-81-52-51-48-46

4 结 论

将InSAR用于矿区开采沉陷可以获得矿区的整个开采沉降影响范围和发展趋势，且其监测技术获得的沉降值与常规水准测量方法的精度相当，可以为沉陷预计模型提供良好的预测数据。

本文提出用InSAR技术与GS-SVR算法相结合预测矿区沉降值，在样本数据量较少的情况下也可以准确的预测矿区的沉降值。本实验采用5组数据进行实验，结果表明：本文提出的预测模型获得的沉降值与GPS技术监测得到的沉降值最大误差为3 mm，最大相对误差为5.9%，其预测精度符合工程的应用需求。

表3 预测结果与实测结果对比

分析表3中的绝对误差和相对误差可知，二者在绝对误差相同的情况下，沉降量较大值的相对误差反而小。而表3中的S0点为实验区的沉降中心点。所以该方法用于沉降中心的沉降预测比非沉降中心可获得更高的预测精度。