基于SBAS-InSAR 和PSO-BP 神经网络算法的矿区地表沉降监测及预测

周定义，左小清

(昆明理工大学 国土资源工程学院，云南 昆明 650093)

矿山开采引起的地表沉降是矿区经常发生的破坏性灾害，是我国最主要的地质灾害之一[1]. 由于采矿活动的性质，煤矿采区地下开采活动引起的地表变形通常呈现非线性特征，可能导致地表变形较快[2]. 因此，有效地预测和管理采矿引起的地表变形，是采矿工业和政府部门关注的一个关键问题[3].

传统的矿区沉降监测方法主要采用精密水准测量、全球定位系统（Global Positioning System，GPS）等[4]. 这些监测方法存在变形监测工作量大、费时、费财、测点难以保存等缺陷，同时传统的监测方法无法对未知区域的沉降进行监测[5]. 与传统的监测方法相比，利用合成孔径雷达干涉（Interferometric Synthetic Aperture Radar，InSAR）技术实施对矿区沉降的监测，不仅可以降低成本，而且可以准实时动态监测整个雷达图像覆盖范围内地表沉降位移量[6].

InSAR 技术是近几年发展起来的一种新型大地测量手段[5]. 理论上可以获得非常精确的数字高程模型和毫米量级的地表形变信息[6]. 目前，基于InSAR 技术对矿区地表沉降监测的研究取得了大量成功范例，但基于InSAR 技术对矿区地表沉降预测的研究相对较少. 陈银翠等[7]利用合成孔径雷达差分干涉（Differential InSAR，D-InSAR）和地理信息系统（Geographic Information System，GIS）对矿区地表沉陷预测，并利用三次样条插值与GM（1，1）组合模型预测，研究结果表明预测模型的精度能够满足D-InSAR 定量分析的精度要求. 李勇发等[3]利用永久散射体合成孔径雷达干涉测量（Persistent Scatterer Interferometric Synthetic Aperture Radar，PSInSAR）与遗传神经网络算法对矿区地表沉降预计，研究结果表明GA-BP(Genetic Algorithm Back Propagation )能有效预测矿区地表沉降. 范洪冬等[8]利用时序合成孔径雷达（Synthetic Aperture Radar，SAR）技术对采空区上方的高速公路变形预测，建立了基于最小二乘支持向量机（Least Squares Support Vector Machines，LS-SVM）的时序沉降预计模型，研究结果表明其最大下沉绝对误差3 mm，最大相对误差2.2%，取得了较为可靠的预计成果. 杨俊凯等[9]利用D-InSAR 技术和灰色Verhulst 模型对矿区地表沉降预计，研究结果表明结合灰色Verhulst模型预测的绝对和相对误差分别为3.4～18.8 mm、1.2%～5.7%，均可满足实际工程需求，有一定的实用性. 以上方法能够有效监测和预测矿区地表沉降，但文献[7-9]提出的方法只适用于短期矿区沉降监测，无法对长时间的矿区地表进行监测，其预测方法存在过度依赖于获取的沉降数据对同一沉降点预测，难以对不同沉降区域预测. 李勇发等[3]提出的PS-InSAR 虽然能有效满足长时间地表监测，但预测上过度依赖于沉降数据，以沉降数据作为遗传神经网络的输入层和输出层进行预测，构建单一的预测模型，存在一定的弊端.

矿区地表变形通常呈现非线性特征，在新采空区，地表以快速、非线性沉陷，老采空区则会在相当长一段时间内产生相对缓慢的残余沉降[10]. 神经网络具有强大的非线性映射能力，可用来处理信息十分复杂、各因子相联规则不明确、背景知识不清楚等非线性问题[11]. 在矿区地表形变监测中，DInSAR 技术受到时空失相关和大气延迟的影响，使得其精度和适用性显著降低[12]. 相比PS-InSAR 技术，小基线集合成孔径雷达干涉（Small Baseline Subsets InSAR，SBAS-InSAR）技术能够较好地克服时空失相关的影响，可以显著增加相干影像对的数目，一定程度上解决了大范围时间、空间失相关和大气效应对InSAR 的影响[13].

综上，本文采用SBAS-InSAR 技术和粒子群优化-反 向 传 播（Particle Swarm Optimization -Back Propagation，PSO-BP）神经网络算法对矿区地表沉降进行监测和预测，从多因子角度构建预测模型，利用SBAS-InSAR 技术对研究区数据进行处理获取到矿区沉降区域及年沉降速率图，将年沉降值作为神经网络输出层，选取矿区地表沉降区域相应影响因子，并将其作为输入层数据，构建矿区影响因子和沉降速率之间的函数模型.

1 矿区地表沉降监测及预测相关技术

1.1 SBAS-InSAR 技 术SBAS-InSAR 是 在 差 分InSAR 基础上发展起来的一种新的时间序列分析方法，能够降低相位噪声和误差[14].

假定在时间t1至ts内获取同一地区的S幅SAR 影像，然后根据干涉组合条件，在短基线距的条件下形成N幅干涉条纹图，且满足

1.2 PSO-BP 算法原理PSO-BP 神经网络将粒子群算法（Particle Swarm Optimization，PSO）和反向传播（Back-Propagation， BP）神经网络相结合，弥补BP 神经网络收敛速度慢、容易陷入局部极值等缺点. 基本原理为PSO 算法与BP 神经网络的误差反向传播训练方法结合，以误差为适应度函数，采用PSO 算法对BP 神经网络的初始权重、阈值进行迭代优化，得到适应度大的初始权重和阈值，再将这些参数用于BP 神经网络，最后输出满足精度要求的值[16]. 其学习步骤如下[17]：

步骤7计算全局误差并判断网络误差是否满足要求.

2 结合SBAS-InSAR 技术和PSO-BP 神经网络算法的方法

本文利用SBAS-InSAR 技术和PSO-BP 算法监测和预测矿区地表长时间沉降，方法的重点在于监测和构建预测模型. 而这种方法存在两个关键问题：①如何获取矿区长时间沉降监测值；②如何克服现有模型存在的弊端.

2.1 基于SBAS-InSAR 长时间沉降监测值的获取矿区地表沉降具有成生缓慢、持续时间长、影响范围广等特点. 传统测量无法满足长时间矿区地表监测，目前矿区广泛应用D-InSAR 技术进行短期沉降监测，取得显著效果，但其不适用于长时间监测. 矿区地表沉降预测需要大量数据进行训练、预测， 短期获取的时序沉降数据无法满足预测要求，亟需一种能够获取长时间沉降监测值的方法，SBAS-InSAR 技术既能获取长时间沉降值，又能为后期的预测奠定基础，为此选择SBAS-InSAR 作为矿区长时间地表沉降监测值获取的手段.

2.2 基于PSO-BP 模型的构建现有预测模型过度依赖沉降数据、模型单一等问题. 其原理仅把获取的沉降数据输入到某一种模型之后得出一个预测值，并未真正意义上实现预测. 能否找到一种方法克服现有模型存在的弊端，真正实现预测? 神经网络具有强大的非线性映射能力，可用来处理信息十分复杂、各因子相联规则不明确、背景知识不清楚等非线性问题[11]. 因此，选用PSO-BP 模型对矿区进行长时间地表沉降预测. 主要实现方法为获取矿区地表沉降的影响因子，对所获取的影响因子与监测值做灰色关联分析，得到最佳的影响因子并将其作为该预测网络的输入层，将沉降值作为输出层，通过该预测网络学习，得到最佳参数并保存网络.需要预测下一年的矿区地表沉降仅需输入与沉降相关的因子，应用其保存的网络，可以实现真正意义上的矿区地表长时间预测.

3 实验分析

3.1 研究区概况和数据来源本文选取云南省昆明市东川区验证该模型的可行性. 东川采矿历史悠久，矿产资源丰富且分布广泛，是中国六大产铜基地之一，享有“天南铜都”的美誉，位于东经102°50′43″～103°16′26″，北纬25°58′14″～26°30′13″，最高海拔4 344.1 m，最低海拔696 m，年平均降水量为1 000.5 mm，东川小江流域为世界典型暴雨泥石流区，被称为“泥石流的天然博物馆”，研究区域如图1所示.

图1 研究区域位置Fig. 1 Location of study area

第1 类数据用于获取研究区域的沉降值，其数据是从欧洲航天局(European Space Agency, ESA)下载的20 景C 波段Sentinel-1A 升轨影像，时间跨度为2018 年4 月25 日至2019 年4 月20 日，极化方式为VV，成像模式为IW. 实验数据参数如表1所示. 为了提高影像轨道精度，引入了POD 精密定轨星历数据，使用美国航空航天局提供的SRTM数据，精度为30 m 的DEM 数据，用于去除地形相位影响，拼接后的DEM 如图2 所示. 第2 类数据用于构建PSO-BP 模型的影响因子，包括研究区域的：（Ⅰ）高程；（Ⅱ）坡度；（Ⅲ）土地利用类型；（Ⅳ）年降雨量数据；（Ⅴ）土壤侵蚀数据；（Ⅵ）土壤质地（粉砂土和粘土）数据；（Ⅶ）地层岩性数据；（Ⅷ）地下水数据；（Ⅸ）NDVI 数据. （Ⅰ）、（Ⅱ）来源于美国航空航天局提供的DEM，根据DEM 制作坡度图；（Ⅲ）～（Ⅵ）来源于中国科学院资源环境数据中心；（Ⅶ）来源于全球岩性图数据库；（Ⅷ）来源于美国航空航天局；（Ⅸ）依据Landsat 8 数据制作.

图2 研究区域DEMFig. 2 DEM of study area

表1 Sentinenl-1A 数据参数Tab. 1 The data parameters of Sentinenl-1A

3.2 矿区地表 沉降监测实验选取2018 年4 月25 日 至2019 年4 月20 日 的20 景Sentinel-1A 斜距单视复数（Single Look Complex，SLC）影像，利用SARscape 软件进行处理，选取日期为20180928的影像作为超级主影像，通过设置临界基线和时间基线，生成91 对像对，设置多视数为1∶4，这样设置可以较好地抑制斑点噪声，采用Minimum Cost Flow 解缠方法和Goldstein 滤波方法做干涉工作流，最终生成干涉图，调整删除不理想的数据，在研究区生成较为理想的部分干涉图（图3）. 从图3 可以看出，针对山区地区，运用SBAS-InSAR 处理后得到的干涉图相干性较好.

图3 研究区生成的部分干涉效果图Fig. 3 Partial interference image generated in the study area

经过轨道精炼和重去平，首先估算和去除残余的恒定相位和解缠后还存在的相位坡道；然后进行第1 次反演、第2 次反演；最后对序列信息进行地理编码获得研究区域2018 年4 月25 日至2019 年4 月20 日雷达视线方向（Los）的形变速率值，如图4 所示. 从图4 可以看出，矿区主要沉降区为A、B、C，其中区域A 位于老采空区茂麓厂、大朵、尖子洞、九龙冶铜遗址和新采空区落雪老背冲矿段、石将军矿洞、沙坝中路铁矿、双龙铅锌矿等；B 区域位于会泽大村子磷矿；C 区域位于东川小水井铁铜矿区；研究区域最大沉降速率为-46.965 4 mm/a，最大抬升速率为179.428 mm/a，最大抬升速率位于小江断裂带中北部，该区域断裂带宽5～20 km，呈现弱剪切强挤压活动特征[18]，以挤压穹起隆升变形为主[19]，致使断裂带抬升. 同时，该区域共有泥石流沟113 条，泥石流频发，大量泥沙堆积导致监测抬升速率过大，如图5 所示，SBAS-InSAR 监测结果与之相吻合.

图4 研究区年沉降速率图Fig. 4 Annual settlement rate of the study area

图5 抬升区泥石流堆积物Fig. 5 Debris flow deposits in uplift area

3.3 矿区地表沉降预测实验本实验选取A 区域作为沉降预测的研究区域，为了解9 个影响因子是否适合作为矿区地表沉降预测的因子，运用SPSSAU 软件对9 个因子分别与沉降值做灰色关联分析，得到灰色关联值，其排序为高程（0.808）＞土壤质地（0.803）＞地下水数据（0.803）＞年降雨量数据（0.802）＞土壤侵蚀数据（0.797）＞土地利用类型数据（0.796）＞地 层 岩 性 数 据（0.792）＞NDVI 数 据（0.776）＞坡度（0.769）. 可以看出选取的9 个因子与矿区沉降值关联程度高，但还不能选取出最终的因子. 运用因子分析方法，通过因子分析确定上述9个因子之间存在的共性，分别对9 个因子进行不同组合，得到最优组合（高程、土地利用类型、年降雨量、地下水数据、土壤质地、土壤侵蚀），其KMO检验统计量为0.631. 当KMO＜0.6 时，说明不适合进行因子分析，因子之间共性极差；当KMO＞0.6 时，说明可以进行因子分析. 由于进行因子分析的目的是为了从与沉降值高关联度的9 因子中选出更为合适的因子，所以KMO 统计量不要求最优，只是起到一种筛选因子的作用.

将选择出来的6 个影响因子与A 区域沉降值构建BP 模型，如图6 所示. 其中，输入层节点数为6，输出节点数为1，隐含层节点数为9. 随机选取740 个点作为PSO-BP 算法的训练样本，对其中的25 个点进行预测. 在学习过程中，为了验证训练样本数量对预测结果的影响，将选出来的740 个数据点分成3 组，其中第1 组训练样本数为275，第2 组训练样本数为475，第3 组训练样本数为715. 利用Matlab 软件构建PSO-BP 算法，各层网络参数如表2所示. 该模型训练参数如表3 所示.

表2 PSO-BP 模型各层网络函数Tab. 2 The each layer network function of PSO-BP model

表3 PSO-BP 模型训练参数Tab. 3 The training parameters of PSO-BP model

图6 本文构建的BP 神经网络模型Fig. 6 BP neural network model constructed in this paper

将3 组样本分别输入PSO-BP 模型里，预测结果与SBAS-InSAR 监测结果对比如图7 所示. 从图7 可以看出，预测曲线和实际曲线大致一致，存在个别“异常”数据（紫色区域）误差过大. “异常”数据并不影响其整体的预测[20]，求各组残差，如表4所示. 为了有效比较各组残差，剔除异常数据，其中，第1 组、第2 组和第3 组的残差区间分别为0.27～6.12 、0.33～4.01 、0.07～3.31 mm，随着训练样本数据的增加，预测效果越好.

表4 各组预测数据残差Tab. 4 Residual data of each group mm

图7 各组数据基于PSO-BP 模型的预测效果图Fig. 7 Prediction effect chart of each group of data based on PSO-BP model

为验证PSO-BP 算法预测性能，根据各组数据产生的性能指数（Performance）验证PSO-BP 算法预测性能，如图8 所示. 各组的收敛最佳迭代次数分别为40、3、2 次，均方误差分别为0.020、0.017、0.016 mm. 随着训练样本增大，收敛迭代次数加快，均方误差减小；反之，收敛迭代次数慢，均方误差大. 从图8 可以看出，第2 组和第3 组迭代快，误差小. 为验证PSO-BP 算法预测精度的可靠性，采用误差绝对值均值（MAE）和平均绝对百分比误差（MAPEA）评价指标，各组MAE 分别为2.46、2.11、1.94 mm；各组MAPEA 分别为25.79%、22.44%、20.94%. 通过对比可以发现，随着训练样本数据的增加，MAE 和MAPEA 逐渐减小，预测精度逐渐增高. 至此，验证了利用PSO-BP 算法对矿区地表沉降预测性能和预测效果较好.

图8 各组数据产生的性能指数Fig. 8 Performance index of each group of data generation

3.4 模型验证分析

3.4.1 模型对比分析实验 目前，对矿区地表沉降监测主要通过D-InSAR 技术获取[7-9]. 为验证SBAS-InSAR 技术在矿区地表沉降监测中的优势，选 取2018 年4 月25 日 和2019 年4 月20 日 的2景Sentinel-1A 斜距单视复数（Single Look Complex，SLC）影像，以20180425 作为主影像，20190420 作为辅影像，采用二轨法进行D-InSAR 处理. 设置视数为4∶1，这样设置能有效抑制斑点噪声，采用Minimum Cost Flow 解缠方法和Goldstein 滤波方法做干涉工作流，最终生成干涉图（图9）. 从图9可以看出，采用D-InSAR 技术处理获得的干涉图效果差，相干性低，仅在城区有一定的干涉效果，在矿区、山区干涉效果极差. 选择控制点（Ground Control Points，GCP）进行轨道精炼和重去平，将得到的相位转形变以及地理编码，最终得到雷达视线方向（Line Of Sight，LOS）方向的形变速率值，如图10 所示. 基于D-InSAR 技术得到研究区年沉降速率范围与SBAS-InSAR 沉降范围大体一致，但在沉降速率上有很大差异，从图10 可以看到研究区最大年沉降速率为-206 mm/a，最大年抬升速率为331 mm/a，这与SBAS-InSAR 结果差距巨大. 这主要同影像的时空分辨率、相干性等因素有关. 在时间和空间上欠采样的雷达影像是无法满足实际的矿区形变监测需要的，只能得到开采沉陷的范围和趋势，对地表破坏的程度分析则显得无能为力[21].D-InSAR 技术易受大气相位延迟和失相关的影响，矿区地表具有散射性不稳定的特点，容易造成失相干，用传统D-InSAR 技术进行长时间监测比较困难[3]. 由此验证了SBAS-InSAR 技术在长时间矿区沉降监测中的优势.

图9 研究区生成的干涉效果图Fig. 9 Interference image generated in the study area

图10 D-InSAR 处理获得年沉降速率图Fig. 10 The annual subsidence rate chart by D-InSAR

为验证本文PSO-BP 模型在矿区地表沉降预测中的优势，选取LS-SVM[8]和灰色Verhulst 模型[9]与本文PSO-BP 模型在矿区地表沉降预测中预测精度进行对比实验. 控制沉降数据来源、样本量等对预测精度的影响，选取SBAS-InSAR 获取的沉降值作为沉降数据，随机选取740 个样本数据，其中735 个数据作为训练样本，对剩余的5 个同名点数据（S1，S2，S3，S4，S5）进行预测，预测结果如表5 所示. 由表5 可以看出LS-SVM 模型预测值、灰色Verhulst 模型预测值和本文PSO-BP 模型预测值与监测值的绝对误差范围为1.148～6.933、1.093～4.922、0.043～4.104 mm. 为有效评估3 种不同预测模型在矿区地表预测中的预测精度，通过获取3 种预测模型预测值与监测值的平均绝对误差，结果表明，本文PSO-BP 预测模型平均绝对误差优于灰色Verhulst 预测模型和LS-SVM 预测模型，其平均绝对误差分别为2.169、2.937、4.131 mm.至此验证了本文PSO-BP 预测模型在矿区地表沉降预测中的有效性和优势.

表5 不同预测模型下的精度对比Tab. 5 Accuracy comparison under different prediction models mm

3.4.2 PSO 和SBAS-InSAR 对BP 网络性能的影响分析 SBAS-InSAR 沉降值作为PSO-BP 模型和BP 模型输出层的训练数据和预测对比数据，对两种模型预测起到决定性的作用. 为有效分析PSO 对BP 网络性能的影响，随机选取740 个点作为样本，其中训练样本为715，对其中的25 个样本点用标准BP 神经网络进行预测. 这样设置的目的是为了避免参数设置不同致使预测结果产生较大差异. 将BP 网络参数与PSO-BP 网络参数设置一致，BP 预测结果与PSO-BP 预测结果（图7 训练样本751）进行对比，如图11 所示. 通过计算绝对误差，评估两种预测模型与监测值的误差，如图12 所示. 从图12 中可以看出，BP 模型的绝对误差的点线图高于PSO-BP 模型. BP 预测模型和PSO-BP 预测模型的平均误差绝对值分别为2.418、1.942 mm.根据BP 和PSO-BP 网络模型显示其收敛最佳迭代次数分别为8、2 次，均方误差分别为0.017、0.016 mm，MAE 分别为2.42 、1.94 mm，MAPEA 分别为26.15%、20.94%. 通过平均误差绝对值、收敛最佳迭代次数、均方误差、MAE 和MAPEA 验证了本文PSO 能够克服传统BP 神经网络导致收敛速度慢、学习速率慢等弊端，有效提高网络预测性能及预测精度，至此验证了PSO-BP 模型的合理性和有效性.

图11 两种模型的预测值和监测值比较Fig. 11 Comparison of predicted and monitored values of two models

图12 两种预测模型的绝对误差点线图Fig. 12 The absolute error point lines of the two prediction models

4 结论

与传统测量相比，SBAS-InSAR 可有效运用于矿区长时间地表沉降监测，能更直观的反应矿区地表沉降区域范围及沉降值. 利用SBAS-InSAR 监测沉降值与影响因子构建PSO-BP 矿区地表预测模型，有效解决现有预测模型过度依赖沉降数据且模型单一等问题. 同时，随着训练样本的增加，PSOBP 预测值与SBAS-InSAR 沉降值残差逐渐减少，算法收敛迭代加快，均方误差降低，虽然存在个别“异常”数据，但并不影像整体的预测. 实验证明，该模型能有效地监测和预测矿区地表沉降，可作为矿区地表沉降监测和预测的有效手段. 矿区地表沉降与矿区开采深度、工作面等有关，受条件限制，无法获取这些因子，如果加入这些因子，可有效提高该模型的预测精度.