基于时序InSAR的改进型PIM矿山采空区沉降模型构建与应用

毛嘉骐,李素敏,成 睿,张龙宇,沈显名

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

0 引言

地下矿体开采后易引发地表塌陷等地质灾害,获取沉降区具体位置对减少灾害事故至关重要[1]。传统的获取方式以物探为主,虽然精度较高,但工作范围小,且需投入大量人力物力[2]。合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar,InSAR)具有全天时、全天候和覆盖范围广等优点,为采空区地表沉降监测提供了有力的技术支持[3]。

利用InSAR对采空区地表进行监测时,由于矿山开采引起的地表沉降速度快、形变量大,且SAR卫星存在一定的重复周期,导致沉降区中心附近出现低相干、失相干现象,形变结果缺失,无法获取准确的沉降区模型[4]。为解决上述问题,当前采用的主要方法是概率积分法(PIM)[5],PIM是由刘宝琛基于随机介质理论而提出的,其主要应用于煤矿规则采空区。PIM利用高精度监测点对相关参数进行反演,以解决采空区地表沉降中心相干性较低的问题。为得到准确度更高的PIM参数,目前应用最广泛的参数优化反演方法有最小二乘法[6]、遗传算法[7]、经验设计值法等。然而不同的反演算法均有其固有缺陷,最小二乘法和经验设计值法的抗粗差能力较弱,反演结果波动性较大;遗传算法全局搜索能力不足,易导致局部最优解而使收敛速度过快,其在采空区沉降监测中的精度较差[8]。

本文提出了一种将时序InSAR与PIM相融合的方法,旨在获取准确的采空区地表沉降信息。该方法首先利用时序InSAR技术获取矿区工作面沉降结果,并对相干性进行分析;之后选取InSAR沉降区边界相干性大于0.3的点以及沉降中心线性插值后的少量水准数据作为特征点,通过改进型PSO算法对PIM参数进行反演,以跳出局部最优解,增强其抗粗差能力,得到可靠的PIM预计参数以及沉降区;提取时序InSAR沉降区边界形变值,与由PIM得到的沉降中心进行插值融合,完成工作面地表下沉的信息提取,在采用少量地表监测数据的情况下构建完整的地表沉降模型。

1 研究区概况及前期研究存在的问题

1.1 研究区概况

以云南省玉溪市大红山铁矿地下采区某工作面为研究对象,其地理位置和范围见图1。矿区地属滇中台南坳南端以及红河断裂与绿汁江断裂夹持的三角地带,西侧出露有变质较深、混合岩化较强的太古代哀牢山群。图1c表示该矿区某地下开采工作面的某巷道开采位置,该工作面对应的地表植被覆盖情况一般[9]。工作面矿层平均厚度72.58 m,主要采用无底柱分段崩落法开采,矿体的分布标高为25.72～945.00 m。该工作面走向开采长度为719 m,倾向开采长度为505 m,平均采深500 m,矿层倾角10°,开采厚度5 m。图1c中:黑色线表示工作面上方走向线,走向线自东南向西北方向布设有H1-H5点位;白色线表示倾向线,自西南向东北方向布设有J1、J2、J3点位。测量工作自2021年6月15日开始,至2022年7月18日结束,共7期。

图1 研究区地理位置

1.2 前期研究存在的问题

前期采用时序InSAR方法对地表形变进行了监测,其原理是利用同一地区的多幅时间基线较短的SAR影像形成的干涉对,通过解缠、滤波等方式去除轨道误差、噪音以及地形的残余相位,然后采用奇异值分解(SVD)的方法,将多个基线集联合求解,并对时间域和空间域的滤波进行分析,分离出残余相位中的大气相位和非线性形变误差,得到目标区内覆盖整个观测时间的地表形变信息[10]。本文选取IW模式下的2021年6月10日至2022年7月23日的34景VV极化方式的Sentinel-1A降轨数据,重访周期为12 d,波长5.6 cm,中心入射角39.5°,使用NASA提供的30 m分辨率SRTM DEM去除干涉计算地形相位对结果的影响。

干涉对参数信息见表1。设置最大临界空间基线占比10%,最大时间基线为60 d,按照基线配对的方式两两差分,按照相干系数法设置0.3为平均相干系数及最小相干系数阈值选取高相干点,之后利用自适应滤波的方法去除噪声对干涉相位的影响,得到的沉降结果见图2。

表1 干涉对参数信息

图2 时序InSAR视线向沉降结果

由图2可知:研究区出现了明显的失相干现象,导致形变结果缺失;边缘部分的形变结果显示,2022年1月24日采区沉降面初步形成,说明此时该工作面的开采活动对上覆岩层稳定性产生了影响,但由于该工作面开采时间较短,沉降量较小,截至2022年7月23日,利用时序InSAR监测到的最大累积沉降量为48 mm。随着开采沿东南向西北方向移动,采区地表逐渐形成了较大的沉降区。

由研究区相干系数图(见图3)可知,在靠近采空区地表中心附近的少部分区域相干系数低于0.3,相干性较差[11],除失相干区域外,其余大部分区域相干性良好。结合InSAR监测结果对实地进行了踏勘,发现地表出现了明显的局部沉降现象(见图4a),且周边存在较大的贯穿裂缝(见图4b、图4c),说明该区域地表产生了较大形变。

图3 研究区相干系数图

图4 采空区地表现场踏勘

2 研究方法

2.1 PIM模型

PIM是一种非连续的随机介质模型,由于岩层在移动过程中会破坏介质的连续性,介质单元发生变化,导致沉积环境也发生变化,从而产生不同的新岩层,在地质构造作用下,新岩层中形成节理、断层等弱面,弱面即为非连续的随机介质[11-12]。矿层开采后岩层与地表移动是一种随机的过程,形成的地表塌陷区见图5。

图5 沉降区在某一点处下沉示意图

由PIM给出任意点T(x,y)的下沉预计模型:

(1)

(2)

(3)

式中,Wi为最大沉降值,Wi(x)、Wi(y)分别为某一点在主断面上的走向、倾向沉降量,u为PIM开采单元,H为工作面采深,m为矿层厚度,q为下沉系数,α为矿层倾角,θ为开采影响传播角,tanβ为主要影响角的正切值,D1、D3分别为工作面的走向、倾向长,s1、s2分别为走向上的左右边界拐点偏移距,s3、s4分别为上山下山方向的拐点偏移距[11]。

2.2 改进型PSO算法

PSO(Particle Swarm Optimization)算法[13]是对鸟类集体觅食行为的模仿,不同个体之间通过集体协作以及信息共享,使种群的移动方向保持一致,而个体间保持一定距离来寻找最优位置,是一种基于粒子群的优化算法。经典的粒子群算法随着种群进化会使其多样性呈递减趋势,从而过早地快速收敛到局部最优解。本文主要以较大的惯性权重和学习因子为开局,有利于增强初始粒子群的全局迭代寻优能力,最后使用较小的惯性权重和学习因子结束,以使粒子群跳出局部最优解。

在对种群进行初始化处理的基础上,选择合适的参数,迭代计算每个粒子的适应值,将个体极值和种群极值代入式(4)、式(5),不断更新每个粒子的位置和移动方向,最终得到所需的寻优解[13]。

Vi=ω×Vi+c1×rand()×(pbesti-Xi)
+c2×rand()×(gbesti-Xi),

(4)

Xi=Xi+Vi,

(5)

式中,Vi为第i个粒子的当前速度,Xi为第i个粒子的当前位置,c1、c2为学习因子,rand()是(0,1)之间的随机数,ω为整体惯性权重。

2.3 工作流程

a.对覆盖矿区的SAR影像进行时序InSAR处理,获取矿区的时序形变信息。收集矿区与SAR影像同时段内的走向和倾向线上的水准数据;

b.在工作面坐标系下,选取时序InSAR沉降区边缘附近高相干点和水准数据作为特征点,利用改进型PSO算法对PIM参数进行寻优反演;

c.根据改进型PSO算法反演出的PIM参数,建立矿区的PIM沉降模型,并将结果统一至WGS-84坐标系,利用克里金插值法,将InSAR沉降区边缘形变值与由PIM获得的中心形变值结合,获取连续且完整的沉降结果。具体的工作流程见图6。

图6 工作流程

3 基于融合方法的矿区沉降结果提取及分析

3.1 基于时序InSAR的特征点选取

为构建概率积分沉降模型,需先将InSAR的LOS向形变投影至垂直方向,之后提取InSAR沉降结果中大于等于0.5的点位进行反演分析[14]。选取InSAR沉降盆地边缘的30个相干性较高的点位,并对该工作面上方的走向和倾向线的7个水准点进行线性插值,使InSAR沉降结果与水准数据在时间上统一,共选取37个特征点,特征点位置见图7(白色圆点表示InSAR高相干性点,黑色三角点表示走向和倾向线水准监测点,黄色圆点表示边界GPS点,黄色三角点表示其余水准点)。

图7 特征点分布

3.2 基于改进型PSO算法的PIM参数反演

PIM预计的参数分别为:下沉系数q,主要影响角正切值tanβ,开采影响传播角θ,拐点偏移距s1、s2、s3、s4。设M=[q,tanβ,s,θ],W为M的搜索空间,若地面点的测量沉降量为wi,PIM计算沉降量为wi′,以二者的误差平方和最小为准则,将PIM参数的反演计算过程表示为式(6)的约束优化问题,即在一个确定的空间W中找到一个向量M0,使得目标函数的值达到最小[15],即

(6)

式中,fp为目标函数,m为引入的观测站点数量。

在算法设计过程中,需反复变换参数进行调试分析,以确定改进的PSO算法参数。设置初始的种群规模为300,最大迭代次数为2 000,学习因子c2=c3=0.5,惯性系数ω=0.9,粒子飞行速度v=1.5,寻优的空间维度为4。利用InSAR计算出沉降区边缘高相干点的形变值,并用位于走向和倾向线上的水准数据对PIM模型参数进行反演。参数反演结果见表2。

表2 改进PSO算法反演结果准确性验证

地下采区经验设计值根据同类型的大红山铜矿开采经验、矿山实际开采方式和矿山地质岩性等资料获得,相比于经验设计值法,利用改进型PSO算法反演得到的开采影响传播角相对误差最小,主要影响角正切值相对误差最大,反演得到的参数值与相关经验设计值存在一定误差,还需进一步分析。

3.3 融合结果分析

由于矿区在开采后地表形变规律易受地质、地形等因素影响,为表示该工作面开采后的地表形变情况,将SBAS-InSAR与PIM进行融合,融合后的下沉曲面见图7。由图7可知,拟合效果较好,能够反映出地表的沉降情况,沉降规律符合PIM模型[16]。将得到的沉降结果转至WGS-84坐标系,采用克里金插值法获得具有连续性的结合PIM的矿区沉降结果(见图8)。

从图7、图8中可以看出,在研究时段内,探测到该工作面上方存在漏斗状沉降区,最大沉降值为175 mm,周边沉降值递减,最终在沉降区边缘逐渐收敛,其收敛速度在走向线上沿西北方向较慢、东南方向较快,在倾向线上均以较快的速度收敛。其原因是矿区沉降受地形、断层及其他外部因素的影响,而在PIM中仅考虑开采引起的沉降,故监测到的地面沉降集中在开采工作面的上方[3,11,17]。

图8 基于改进型PSO算法的预计下沉曲面

图9 插值后连续沉降区

3.4 融合法监测结果精度分析

为了对融合法得到的监测精度进行验证,将其与该工作面上走向及倾向线上8个水准点的监测数据进行对比分析,结果见图10、图11。由图10、图11可知,采用融合方法获取的形变量与水准监测得到的形变量较接近,而采用经验设计值法得到的形变量与水准监测得到的形变量存在明显差异。为进一步判断模型的准确性,利用图8中未参与模型反演的黄色点监测值进行精度分析,结果见图12。图12中,K1-K3表示边缘部分GPS监测点,H6和G5表示走向和倾向线附近水准点,通过对比发现,构建出的沉降模型符合现场监测值的变化趋势。

图10 走向线沉降量对比图

图11 倾向线沉降量对比图

图12 现场监测点沉降量对比图

在相同的条件下,对不同方法在走向和倾向线上的形变结果误差进行统计,结果见表3。从表3中可以看出,融合方法的平均绝对误差为43 mm,最大绝对误差为70 mm,最大沉降量绝对误差为52 mm,均方根误差为46 mm,联合方法的拟合效果明显优于经验设计值法,说明融合时序InSAR与改进型PSO算法可以用于开采沉陷预计参数的求取,并能获得精度较高的沉降区模型。误差产生的主要原因是PIM为一个理想化的模型[18],未考虑断层、排土场和降雨等因素对地表沉降的影响。

表3 实测值与融合值、经验设计值的总体误差对比

4 结论

本文针对时序InSAR在矿山采空区地表沉降监测中的不足,提出了一种融合时序InSAR与改进型PIM的矿区地表沉降监测方法。以大红山矿区某采空区为研究对象,利用该方法获取在研究时段内的沉降区模型,并与时序InSAR在同一时段内的水准数据以及由经验设计值法得到的沉降结果进行对比,得到以下主要结论:

a.采用时序InSAR技术与改进型PIM融合的方法,弥补了InSAR在沉降中心附近由于失相干而导致的形变结果缺失的问题,且仅结合少量水准数据便获取了矿山地下开采引起的地表沉降区完整模型。

b.融合改进型PSO算法与时序InSAR技术反演出的PIM参数效果较好,在走向线和倾向线上误差结果均优于由经验设计值法得到的沉降结果,且各项误差均不超过100 mm,说明该方法满足实际工程需求,为时序InSAR技术在矿山开采沉陷监测中的应用提供了新思路。