李鸣庚 张书毕 高延东 李世金 张艳锁 贾义琨 孔令辰

(1.中国矿业大学环境与测绘学院,江苏 徐州 221116;2.自然资源部国土环境与灾害监测重点实验室,江苏 徐州 221116)

露天开采是现阶段我国资源获取的主要手段之一[1],然而随着露天矿开采深度的增加,工作面滑坡、排土场边坡失稳、矿区地表沉陷等重大灾害对露天矿安全生产造成的威胁也随之增加[2-3]。辽宁省鞍山市是我国东北老工业基地重要的铁矿资源城市,其周边的西鞍山、东鞍山、大孤山、齐大山、眼前山、鞍千矿等十几个露天矿区不可避免地面临由于长期大规模露天开采引发的地表形变问题[4-5]。因此,对鞍山市周边露天矿区进行地表形变监测对于城市健康发展具有重要意义[6]。

传统地表形变监测方法如水准测量、全站仪、GPS (Global Position System)测量等,虽然精度较高,但也存在耗时长、成本高、且难以实现长时序、大范围的地表形变监测等不足[7]。时序InSAR技术是在合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar,InSAR)技术的基础上发展而来,如小基线集(Small base lines subset InSAR,SBAS-InSAR)技术与永久散射体 (Permanent Scatters InSAR,PS-InSAR)技术,已在大范围形变监测中得到了广泛应用,具有高效率、低成本、长时序、大范围等优势[8]。时序In-SAR技术虽然在城市地区具有较高的监测精度与点位密度,但在裸地、植被覆盖区与非城区则存在因相干性较低引起的监测点密度不足的问题[9]。在此基础上,学者们提出了分布式目标InSAR (Distributed Scatterers InSAR,DS-InSAR)技术,可有效提高低相干区域的监测点密度[10],在长时序、大范围的地表形变监测方面具有明显优势。

近年来,不少学者将DS-InSAR技术应用于大范围地表形变监测中,取得了较为理想的结果。王波等[11]利用DS-InSAR技术获取了西部矿区时序形变结果,并将DS-InSAR监测结果与水准数据进行对比,相关性系数为0.97,证明该技术具有较高的可靠性。李柱等[12]通过DS-InSAR技术对乌达煤田火区进行了长时序地表形变监测分析,通过与TCP-InSAR监测数据的对比分析,反映出DS-InSAR技术可靠性较好。吴昊等[13]通过DS-InSAR技术对巴西Brumadinho矿区尾矿库进行了形变监测,结果表明:DSInSAR获取的尾矿库监测点密度为PS-InSAR的2.2倍,两者获取的形变速率具有较强的相关性。王新玲等[14]通过DS-InSAR技术对郓城矿区进行了地表形变监测,结果表明:DS-InSAR技术可显著提升裸地和稀疏低矮植被区的监测点数量,获取的形变场更为完整。赵立峰等[15]针对常规InSAR技术在矿区受失相干影响严重、监测点数量低的问题,通过DS-InSAR技术提取了张双楼煤矿长时序地表形变结果,并与传统SBAS-InSAR技术进行对比分析,验证了DS-InSAR监测精度。杨嘉威等[16]基于DS-InSAR技术对矿区铁路线进行了形变监测,通过相位优化显著提升了高相干像元数量与解算质量。

上述成果分析表明:虽然DS-InSAR技术在大范围地表形变监测方面已有较多应用,但研究区以煤矿或露天矿尾矿库、排土场居多,对于该技术在露天矿形变监测方面的应用研究相对较少。为此,本研究提出了一种基于FaSHPS同质像元提取与复相干矩阵特征值分解相位优化的DS-InSAR时序地表形变监测技术,获取了2020年1月—2021年8月辽宁省鞍山地区齐大山矿、大孤山矿、鞍千矿3个矿区的时序地表形变特征,结合常规时序InSAR技术监测结果对比验证了DS-InSAR技术的可靠性,丰富了DSInSAR技术在露天矿形变监测方面的研究,提高了DS-InSAR技术露天矿形变监测的精度,对城市周边大范围露天矿形变分析具有一定的参考意义。

1 研究区概况与数据来源

1.1 研究区概况

鞍山市位于辽宁省中部,是我国重要的钢铁生产基地,其周边主要矿区分布如图1所示。由图1可知:东鞍山、齐大山、大孤山、眼前山等大型露天矿采场从东北至东南向南环绕鞍山市城区,最远距市中心12 km,最近距市中心区仅2 km,形成一个宽约6 km、长为40 km的矿区生产活动带[17]。鞍千矿业铁矿位于胡家庙村和肖家堡村之间,据统计原有19个小型铁矿采场和2家小型硅石矿采场,占地面积41.18 km2[18]。东鞍山、大孤山、齐大山各矿区占地面积均在1 000 km2以上,且距离市区较近,东鞍山铁矿距离市中心仅7 km,距离最近的建成居民区仅2 km。由此可见,露天矿的安全监测对于鞍山市健康发展具有极其重要的意义。

1.2 数据源

为对鞍山市周边矿区进行大范围、长时序的地表形变监测,本研究选取欧空局2014年发射的Sentinel-1A卫星、单视复数(Single Look Complex,SLC)、升轨SAR影像,时间跨度为2020年1月9日—2021年8月31日,共50景,覆盖范围为41.013°～41.203°N,123.000°～123.224°E,如图1所示。数字高程模型(Digital Elevation Model,DEM)采用SRTM (Shuttle Radar Topography Mission) DEM,分辨率为30 m。为削弱大气误差,本研究引入了合成孔径雷达通用性大气改正在线服务(Generic Atmospheric Correction Online Service for InSAR,GACOS)的天顶对流层延迟(Zenith Troposphere Delay,ZTD )产品[18]。另外,试验还使用了Sentinel-1A POD精密定轨星历(Precise Orbit Ephemerides,POD)数据,定位精度优于5 cm。

图1 研究区地理位置Fig.1 Location of the study area

2 试验方法

为实现对鞍山地区齐大山矿、大孤山矿、鞍千矿3个矿区大范围、高精度的时序地表形变监测,本研究采用了一种基于FaSHPS同质像元提取与复相干矩阵特征值分解相位优化的DS-InSAR技术。该技术与传统时序InSAR技术的区别在于增加了同质像元提取与相位优化2个步骤,能够有效解决传统时序InSAR技术在低相干区域监测点数量不足且空间分布不均匀的问题,在露天矿区长时序、大范围地表形变监测方面具有明显优势。

2.1 FaSHPS同质像元提取

本研究采用快速同质像元识别算法FaSHPS(Fast

Statistically

homogeneouspixels

selection,FaSHPS)进行同质像元提取。该算法基于SAR影像任一像元的振幅均值服从高斯分布的假设,通过采用置信区间估计实现快速高效的同质像元识别,有效解决传统的基于假设检验的同质像元识别算法效率不高的问题,且该算法提取的同质像元间异质性较低,在露天矿等高相干点分布不均且数量不足的区域具有明显优势[19]。设有N景已配准的SAR影像,对于影像上任一像元I,其时间维度上的振幅A(I)为

式中,振幅A(I)的均值(I)服从高斯分布,可由式(2)获取(I)的区间估计:

式中,P{·}为概率;E(·)为期望;Var(·)为方差;Z1-α/2表示标准正态分布中置信度为α时所对应的分位点。

单视复数SAR影像上的任一像元在时间维度上服从瑞利分布,则变异系数Cv为一常数,计算公式为

假设在均质地区像元I后向散射稳定,则式(2)可转化为

通过确定置信区间α,若待估计像元I的振幅平均值(I)在目标像元置信区间内,则为同质像元。

2.2 复相干矩阵分解法相位优化

为提高获取的同质像元密度及可靠性,本研究采用复相干矩阵分解法进行相位优化。由于分布式目标在SAR影像上具有相似的后向散射性,因此复相干矩阵可通过特征值分解实现相位优化。该方法在运算过程中无需进行迭代,相位优化效率较高,可获取高质量的优化相位[20]。分布式目标的相干性矩阵可表示为

式中,向量y=[y1,y2,…,yN]为分布式目标的同质像元在N景影像上的复数观测量经过时间维振幅归一化后的复数观测值;NSHPs为分布式目标同质点数量;Ω为分布式目标同质点集;yH表示向量y的Hermitian共轭转置。

由于分布式目标的相干性矩阵为半正定Hermitian矩阵,因此可根据特征值分解原理将式(5)转化为

式中,λi为特征值;μi为与λi对应的特征向量;μ1为优化后的相位;为μ1的Hermitian共轭转置;signal为主散射体信号;noise为噪声信号。

完成同质像元提取与相位优化后,可利用干涉相位在经过相位优化前后的拟合度进行相位优化质量评价。通过设置相干性阈值筛选分布式目标中具有高相干性与高信噪比的部分,将筛选出的高相干点融合PS点进行时序InSAR处理,实现对研究区高精度形变监测。

本研究方法的数据处理流程如图2所示。为保证运算效率,本研究对SLC采用2×1多视处理。

图2 DS-InSAR数据处理流程Fig.2 Flow of DS-InSAR data processing

3 监测结果与分析

3.1 监测结果分析

利用DS-InSAR技术,基于2020年1月—2021年8月共50景Sentinel-1A影像,通过如图3所示的时空基线组合,获取鞍山市周边矿区年形变速率变化特征。在缺少水准数据的情况下,采用交叉验证[21]方法将DS-InSAR技术与常用的StaMPS-InSAR技术获取的年沉降速率值进行对比分析,以验证DS-In-SAR技术的可靠性。本研究以StaMPS监测点为基准,通过经纬度搜索最邻近DS点为同名点,共获取了16 890对同名点,占StaMPS总监测点数量的73.53%。通过对同名点的形变速率值进行相关性分析,获取二者形变速率相关性与形变速率差值直方图。

图3 时空基线分布Fig.3 Layout of spatial-temporal baseline

鞍山市周边矿区形变速率分布如图4所示。由图4可知:DS-InSAR与StaMPS监测结果具有较好的空间一致性。在研究区内, StaMPS技术共选出22 971个监测点;DS-InSAR技术共选出28 804个监测点,为StaMPS技术的1.253倍。DS-InSAR技术在尾矿库、丘陵等低相干区域的点位数量明显提升,且DS点在具有较高点位密度的同时,在空间分布上也较为均匀,在一定程度上提高了形变区的监测精度,提供的形变信息更加准确、完整。

图4 鞍山市周边矿区卫星视线向年平均沉降速率Fig.4 Annual average subsidence rates of LOS direction in Anshan City

StaMPS-InSAR与DS-InSAR监测结果交叉验证获取的形变速率相关图及差值分布如图5所示。

由图5(a)可知:StaMPS技术与DS-InSAR技术获取的同名点地表形变速率具有较高的相关性,通过Pearson函数计算出的两者相关系数R为0.855 6。由图5(b)可知:2种技术获取的同名点形变速率差值较小,标准差为13.91 mm/a。在缺乏实测水准数据的情况下,通过DS-InSAR技术与StaMPS技术监测结果的交叉验证,可有效反映出DS-InSAR技术的可靠性。

3.2 鞍山市周边矿区大范围地表形变分析

为系统研究鞍山市周边矿区地表形变特征,本研究基于DS-InSAR技术,对鞍山市周边矿区进行了长时序的地表形变监测,获取了各时段相对于2020年1月9日的时序累积形变量,如图6所示。从时间序列上分析,鞍山市周边矿区自2020年1月9日起,出现明显的、连续的地表形变现象。2020年1月9日—2020年9月17日,各形变区的形变范围呈缓慢扩大趋势,形成了A～E5个主要形变区。2020年11月4日—2021年8月31日,各形变区的扩大趋势逐渐减缓,形变区内的累积形变量呈增加趋势。从空间分布上分析,在监测时段内,鞍山市周边矿区地表形变现象较为严重,出现了A～E5个主要形变区。其中,形变区A、B位于齐大山矿范围内,形变区C、D主要位于鞍千矿范围内,形变区E则位于大孤山尾矿库以西、露天矿坑以东的区域。齐大山矿区形变范围与形变量均较高,形变区A与形变区B呈自西北到东南的连续分布趋势,与齐大山矿床的NW—SE走向一致。鞍千矿范围内存在多个形变中心,主要形变集中于形变区C,形变区D存在多个形变中心,形变量均较小。形变区E位于大孤山范围内,形变现象较为明显,且大孤山露天矿坑西侧也存在少量的形变现象。

图6 鞍山市周边矿区时序形变累积结果Fig.6 Time series deformation accumulation results of surrounding mining areas in Anshan City

为更详细地分析鞍山市周边矿区地表形变,结合图4的DS-InSAR监测结果对各形变区进行了定量分析。综合图4、图6可知:A区最高形变速率为-260.43 mm/a,B区为-258.49 mm/a;A区的形变中心存在因形变梯度较大导致失相干严重而出现DS点数量不足的现象;C区最高沉降速率为-229.05 mm/a,D区为-175.92 mm/a;E区最高沉降速率为-312.18 mm/a,高于上述各区。

3.3 鞍山市周边矿区时序地表形变分析

为分析鞍山市周边矿区时序形变特征,在如图6所示的A～E5个主要形变区内,根据DS-InSAR监测结果选取了Q1～Q8共8个监测点,位置分布如图7所示。提取上述8个监测点对应的时序形变结果,绘制时序形变曲线。同时选取Q4、Q6、Q7、Q84个DS监测点对应的同名StaMPS监测点Q4′、Q6′、Q7′、Q8′绘制时序形变曲线进行对比分析,以验证DS-InSAR技术的可靠性。最终获取的时序形变曲线如图8所示。由图8(c)、图8(d)可知:DS监测点与对应的同名StaMPS监测点在时序形变趋势上具有较高的一致性,但整体上形变量存在部分差异。考虑到DS监测点与StaMPS监测点间存在一定距离与误差的影响,认为形变量的差异处于合理范围内,可以使用选取的DS监测点进行时序形变分析。

图7 监测点位置Fig.7 Location of the monitoring points

图8 各监测点的时序形变曲线及StaMPS对比结果Fig.8 Time series deformation curves and stamps comparison results of each monitoring point

图8(a)为监测点Q1～Q3时序形变结果。监测点Q1形变速率最高,为-260.43 mm/a,累积沉降量为397.48 mm,且Q1点周围出现因相干性较低而引起的DS点缺失现象,说明A区存在更为严重的沉降趋势。Q2点沉降速率为-254.66 mm/a,累积沉降量为-407.35 mm。Q3点沉降速率为-258.49 mm/a,累积沉降量为-426.99 mm/a。3个监测点均位于齐大山矿采区内,自2020年1月起呈连续沉降趋势。Q1点自2020年1月开始呈现较高的沉降速率,在2020年6月—2020年8月沉降速率出现减缓,2020年8月后沉降速率再次提升至较高水平。Q2、Q3点位于齐大山矿东南部,地处同一采区,整体形变趋势上具有一定的相似性。Q2点的形变速率与形变量均低于Q3点,且Q3点在2020年全年均保持较高的沉降速率,进入2021年则明显放缓。推测这种沉降速率的变化可能与监测时段内齐大山铁矿的扩帮生产[22]规划有关。

图8(b)中,Q4点位于C区,为鞍千矿主要形变区。Q5、Q6位于D区。Q4点形变速率为-229.53 mm/a,累积形变量为-315.86 mm。Q5、Q6点的形变速率分别为-175.92 mm/a、-166.76 mm/a。累积形变量分别为-240.12 mm、-256.03 mm。分析图8(b)可知:2020年1—8月,监测点Q4～Q6的形变速率均保持较低水平,在2020年9月后形变速率出现明显提升,且呈连续的沉降趋势。Q4点的沉降速率与沉降量均高于Q5与Q6点。Q5、Q6点在形变速率的变化趋势与累积形变量方面具有一定的一致性。鞍千矿区内小型采区较多,因此沉降中心数量相对较多、空间分布上较为分散。C区推测可能为监测时段内鞍千矿的主要采区,沉降面积、沉降速率与沉降量均高于D区。

图8(c)为监测点Q7、Q8时序形变结果。Q7、Q8点主要处于大孤山矿区内,Q7点位于大孤山露天矿坑以东、尾矿库以西,形变速率为-312.18 mm/a,累积沉降量为-493.43 mm。Q8点位于大孤山露天矿坑西北侧,形变速率为-110.61 mm/a,累积形变量为-169.88 mm。Q7点自2020年1月即呈连续沉降趋势。除了在2020年10月、2020年12月、2021年5月等少数时段沉降速率有所降低外,其余时段内均保持较高的沉降速率。Q8点沉降速率与沉降量均较低,推测是由自然原因引发的地表沉降。2020年,大孤山铁矿正处于露天开采转为地下开采的过渡期[23],原本作为主要开采区的大孤山露天矿坑内部沉降速率相对较小,而矿坑西部、尾矿库东部区域发生明显的地表沉降,这种变化可能是由于大孤山矿由露天开采转为地下开采引起的。

3.4 典型露天矿区形变特征分析

本研究选取A～E5个形变区中监测点位密度较高、分布较为均匀的C区作为典型露天矿区进行形变特征分析。为定量分析C区的形变特征,如图9所示绘制了两条剖面线x、y,并按照图示方向在两条剖面线上各选取了13个监测点,编号分别为X1,X2,…,X13,Y1,Y2,…,Y13绘制了形变区剖面时序下沉曲线,如图10所示。

图9 剖面线与监测点分布示意Fig.9 Schematic of profile lines and monitoring points

图10 形变区剖面线时序下沉曲线Fig.10 Time series subsidence curves of profile lines in deformation area

由图10(a)可知:2020年1月9日—2020年7月7日,C区的主要沉降集中在位于形变区西侧的X3点附近,其余监测点整体呈下沉趋势,但未出现明显的沉降中心。2020年7月7日后,X3发生持续沉降的同时,位于形变区中部的监测点X6、X7点开始出现明显的下沉趋势,2020年12月10日前,X6点沉降量略高于X7点。2021年1月15日后,沉降中心开始向沉降区中部转移,X6、X7点附近地区的沉降量较X3点更为明显,取代X3点成为新的沉降中心。

由图10(b)可知:2020年1月9日—2020年7月7日,位于C区南部的监测点Y10沉降现象明显。8月12日后,位于沉降区中部的Y6、Y7点开始发生明显的沉降现象,在2020年9月29日后沉降量超过Y10点。Y10点在2020年12月10日后形变速率降低,Y6、Y7点则持续保持较高的沉降速率,成为新的沉降中心。综合剖面线x、y的沉降中心变化情况,可以推测2020年1月9日—7月7日,C区的主要开采区可能位于露天矿西南侧。7月7日后,露天矿中心地带的开采强度逐渐增大,而西南侧的开采强度开始减小。进入2021年后,矿区西南侧的开采强度明显减小,主要开采区转至矿区中部,成为新的沉降中心。

4 结 论

为实现对露天矿区的高精度、大范围、长时序的地表形变监测,以鞍山市周边矿区为例,研究了一种适用于露天矿时序形变监测的优化DS-InSAR技术。通过该技术对2020年1月9日—2021年8月31日共50景Sentinel-1A SAR影像进行处理,获取了研究区的时序地表形变特征,并结合监测结果与相关资料进行了地表形变分析。得出如下结论:

(1)与传统的StaMPS-InSAR技术相比,本研究采用的优化DS-InSAR技术在通过FaSHPs算法提取的研究区内部分布式目标的基础上,采用复相干矩阵特征值分解法实现相位优化,再通过时序InSAR处理获取露天矿区时序形变特征。该方法有效提高了获取的同质像元密度及可靠性,在一定程度上解决了传统时序InSAR技术在露天矿区因失相干严重而产生的监测点数量不足、空间分布不均匀的问题,实现了对露天矿区大范围、长时序、高精度的地表形变监测。

(2)研究时段内,通过DS-InSAR技术监测到大孤山矿、齐大山矿与鞍千矿共存在5处沉降量级与形变范围不同的形变区,最大形变速率为-312.18 mm/a。结合相关文献资料分析表明,本研究采用的DS-InSAR技术监测到的形变区位置与范围较为准确,获取的时序地表形变监测结果具有较高的精度和可靠性。该方法可为露天矿区地表形变监测及其形变机理研究提供参考。

(3)本研究仅通过DS-InSAR技术对鞍山市周边矿区进行地表形变监测,结合传统的StaMPS-InSAR监测结果交叉验证其可靠性,并通过相关资料推测各矿区开采活动,缺少具体实测数据资料与矿区开采资料进行辅助验证。后续研究中可结合各露天矿具体开采量、开采时段、开采计划及实测水准数据等资料与DS-InSAR时序监测结果进行对比,进一步分析DS-InSAR技术在露天矿地表形变监测中的适用性与可靠性。