李如仁 孙加瑶

(沈阳建筑大学交通与测绘工程学院,辽宁 沈阳 110168)

尾矿库作为一种具有人造高势能泥石流危险源的特殊工业建筑,一旦发生溃坝或泄露,会对社会安全造成巨大威胁,同时因尾矿中含有重金属物质,也将严重污染环境[1]。据统计,在2001—2015年全国尾矿库共发生99起溃坝和泄漏事故[2]。2008年9月8日山西新塔矿业有限公司尾矿库溃坝直接导致泥石流发生,造成下游居民277人死亡、4人失踪,直接经济损失9 619.21万元[3]。2015年11月5日,巴西淡水河谷公司尾矿库发生溃坝,造成19人死亡,数千人失踪[4]。因此对尾矿库进行长期的形变监测及预警具有十分重要的意义。

目前,针对尾矿库的监测手段主要依赖于传统的水准测量、全球定位系统(Global Positioning System,GPS)、测量机器人等。通过该类技术手段只能实现对监测站点绘制出沉降剖面曲线图并进行相关分析,如需对监测区域进行大范围、全方位的分析则必须大量布设沉降监测点,并通过拟合、插值等技术手段绘制沉降图,这种大范围沉降监测方式不仅存在工作量大、费时、费力等不足,还会由于测点难以保存产生更大的误差[5]。近年来,合成孔径雷达干涉测量(Interferometric Synthetic Apeture Rader,InSAR)技术作为一种新型的对地观测技术,可以根据卫星的轨道信息与观测点的几何关系,测算出观测点的空间信息和形变,理论上可以获取高精度的数字高程模型(Digital Elevation Model,DEM)和毫米级的地表形变信息[6]。为获取连续形变值,提高测量精度,有学者提出了时间序列InSAR技术,很好地解决了时空失相关和大气延迟等问题,例如永久散射体合成孔径雷达干涉测量(Permanent Scatterer Synthetic Aperture Radar,PSInSAR)技术[7],小基线集合成孔径雷达干涉测量(Small Baseline Subsets Interferometric Synthetic Aperture Radar,SBAS-InSAR)技术[8]等,目前这些方法被广泛用于尾矿库形变监测。例如,王丽娟[9]利用DInSAR技术对万年沟尾矿库进行监测,通过实测数据证明了该方法的可靠性;陈娅男等[10]利用SBAS-In-SAR技术获取了卡房尾矿库的形变演化规律,研究发现卡房尾矿库形变受人为工程和雨季影响较大;吴昊等[11]利用DS-InSAR技术对巴西Brumadinho尾矿库进行监测,分析了其演化过程及溃坝原因。虽然InSAR技术在尾矿库形变监测方面的应用研究已有较多案例,但是针对尾矿库形变预测模型的研究涉及较少,在很大程度上影响了对尾矿库灾害的实时预警。

目前,常用的预测方法主要有BP神经网络[12]、支持向量机回归(Support Vector Regression,SVR)[13]、循环神经网络(Recurrent Neural Network,RNN)[14]等。然而上述预测模型存在过拟合、梯度消失、梯度爆炸[15]等问题,近年来通过RNN改进的长短期记忆(Long Short-Term Memory,LSTM)神经网络在各方面时间序列预测中被广泛应用[16],并表现出良好的长期及短期预测能力。本研究在现有成果的基础上融合SBAS-InSAR技术、LSTM神经网络和网格搜索(Grid Search,GS)算法,对鞍山市西果园尾矿库进行形变监测及分析,总结形变规律,将所获沉降结果作为学习训练样本,建立预测动态模型,实现监测与预测的一体化,为尾矿库灾害防治提供参考。

1 形变监测与预测原理

1.1 SBAS-InSAR基本原理

SBAS-InSAR作为一种在InSAR基础上发展的新型时序分析手段,极大地改善了常规InSAR技术的空间、时间失相关等问题,可以有效降低相位噪声和误差[17]。其基本思路是对空间基线与时间基线进行限制,设置阈值,将SAR影像按规定组合生成若干个小基线集合从而获取干涉对,使同一子集内的基线相对较短,集合之间的基线相对较长。在各短基线子集内部利用最小二乘法(Least Squares,LS)计算时间序列形变,各集合之间使用奇异值分解法(Singular Value Decomposition,SVD)联合求解具有高精度的完整时间序列形变值。

假设在时间t0至时间tn内获取同一研究区的N幅SAR影像,选择一幅影像作为主影像,然后根据所设组合阈值,形成M幅干涉图,且满足下式:

对于ta与tb(ta＜tb)时刻生成的第j幅干涉图像,在去除平地效应和地形相位的影响后,若方位向坐标以x表示,距离向坐标以y表示,则干涉图j中任意一点的相位可表示为

式中,φtb(x,y)、φta(x,y)分别表示坐标为(x,y)的像元在其对应时刻的相位;λ为雷达波长;dtb、dta分别代表在tb和ta时刻该像元在雷达视线方向上的形变值;φtopo、φatm、φnoise分别为干涉图中残余的地形相位、大气延迟相位、噪声相位。

暂不考虑大气延迟相位与噪声相位,地形相位可通过高精度DEM消除。因而式(2)可以简化为

因此,解缠后的所有差分干涉图可以表示为

式中,φ为所有影像数据相位值所组成的向量;δφ为所有差分干涉相位值所组成的向量;A为M×N阶矩阵。当所有数据均属于一个基线集时,M≥N,矩阵A的秩为N,此时可以通过最小二乘法求出φ的估计值(AT·A)-1·AT·δφ。

当数据集分散在几个子集中时,矩阵A的秩小于N,为秩亏矩阵,方程组有无穷多组解。为求得唯一解,SBAS-InSAR技术采用奇异值分解法,将多个基线集联合求解,求出最小范数意义上的最小二乘解,从而获取累计形变量。

1.2 LSTM算法原理

传统的RNN模型存在梯度消失、梯度爆炸等问题,会遗忘距离较远的信息,而LSTM模型作为RNN的一种特例[18],在细胞状态也称作记忆单元中增添遗忘门、输入门和输出门[19],可以有效控制信息的保留、遗忘和传递更新,使其具有长期记忆的能力,大幅提升了预测精度。图1所示为一个LSTM单元结构[20]。

图1 LSTM单元结构Fig.1 Structure of LSTM cell

当上一时刻隐藏层输出的ht-1和当前时刻输入的xt流入记忆单元后,遗忘门可以控制信息的取舍,通过Sigmoid函数得到当前时刻遗忘门的输出值ft:

式中,ft取值为0～1,表示对上一细胞状态输出信息的遗忘程度;Wf为这一步骤的权重系数矩阵;ht-1为前一时刻输出值;xt为当前时刻输入值;bf为偏置项;σ·()为Sigmoid函数;t表示时刻。

输入门接受到当前时刻输入的信息后,通过Sigmoid层和Tanh层共同筛选将新信息储存到细胞状态中。Sigmoid函数决定将要更新的输出值it作为长期记忆。Tanh函数创建一个候选向量作为短期记忆。公式为

式中,Wi、Wc为当前步骤Sigmoid层与Tanh层的权重系数矩阵;bi、bc为Sigmoid层与Tanh层的偏置项。

为了忘记部分之前的信息和记住部分新信息,当前细胞状态Ct可用上一状态Ct-1与ft、it、共同表示:

输出门确定细胞输出的有效信息Ot,而后当前细胞状态Ct通过Tanh函数进行处理生成新的隐含层输出结果:

式中,Wo为当前步骤输出门权重系数矩阵;bo为偏置项。

2 形变监测结果与分析

2.1 研究区概况与试验数据

西果园尾矿库位于鞍钢集团鞍山矿业公司东鞍山烧结厂东南11 km处的西果园村[21]。库区中心地理坐标为东经122°59′52″、北纬40°58′42″,东侧距离鞍隆公路1.9 km,西侧紧邻排土场,库南北长2.3 km,东西长1.6 km。该尾矿库由3座库组成,分别为西沟、南峪沟和井峪沟尾矿库,西沟和南峪沟尾矿库始建于20世纪70年代,于1979年投入使用,原设计标高为260.0 m。1990年在南峪沟东侧新建了井峪沟尾矿库,该尾矿库原设计堆积标高260.0 m。尾矿库采用上游法筑坝,堆积坝平均坡比为1∶5,库中尾矿砂主要成分为石英、磁铁矿、赤铁矿、钛铁矿、石榴子石、白云母、方解石,还有少量稳定和非稳定矿物。在标高225.0 m之前3座库为独立的尾矿库,225.0 m之后三者堆积坝连成一体。后续为满足生产需要加高扩容后,设计最终标高为275.0 m,总库容为7 349.5×104m3,总坝体最高处为158 m,按照相关标准,该尾矿库为二等库,研究区域如图2所示。

图2 研究区及尾矿库范围Fig.2 Scope of the study area and tailings pond

本次试验采用60景Sentinel-1A雷达卫星升轨数据,入射波为C波段,微波波长为5.6 cm,入射角为39°,重访周期为12 d,时间跨度为2018年7月18日—2020年6月25日,极化方式为VV同极化,同时结合SRTM精度为30 m的DEM以及对应时期的精密轨道文件进行地形相位和轨道误差的去除。

2.2 时序InSAR监测结果

利用SARscape软件,采用SBAS-InSAR技术对所获取的影像进行处理,选取日期为2019年4月8日的影像作为超级主影像,生成了318个干涉对;利用精密轨道数据进行去平地效应,利用DEM数据消除地形相位,采用最小费用流法(Minimum Cost Flow)进行相位解缠;在进行轨道精炼与重去平后,进行形变速率反演与地理编码,以获取雷达视线方向的年平均形变速率值。

2018年7月—2020年6月研究区年平均形变速率分布如图3所示。由图3可知:尾矿库堆积坝形变量较小,初期坝及坝底较为平稳,主要沉降区为尾矿库西侧排土场与尾矿库顶部干滩区域,最大沉降速率为95 mm/a,年平均形变速率为-95～30 mm/a。图中,A、B点为研究区内2个GPS观测站,G1、G2、G3点为后期建立预测模型时选取的点位,这3个点分别分布在3座尾矿库顶部位置,且十分靠近最上层堆积坝,一旦发生剧烈形变,其监测数据最具参考意义。

图3 年平均形变速率Fig.3 Annual average deformation rate

由于本研究使用的影像较多,虽然已经获取了研究区整体形变速率,但为了更加直观地展示出整个监测时间段内研究区的形变演化过程,绘制了等间隔的12幅时间序列形变图,间隔时间为60 d,所有图像均显示以2018年7月18日为参考时间的时序累计形变值,具体如图4所示。由图4可知:随着时间推移,尾矿库区域的坝体顶部沉降量不断增大,沉降范围开始从中心向外扩散,最终覆盖整座尾矿库,3座尾矿库沉降趋势及变化规律大致相同;在2 a内,研究区最大的地面累计沉降量达到-181 mm,最大累计抬升量为60 mm。

图4 时序累计形变值Fig.4 Time series cumulative deformation values

为了检验SBAS-InSAR技术的监测精度,本研究在A、B两观测站点中获取了与SAR影像时间相互对应的10期GPS观测数据,将其与SBAS-InSAR监测结果进行对比验证。但由于SBAS-InSAR技术得到的形变为雷达视线方向的形变,属于真实形变在入射角上的投影,可以根据入射角的余弦值对视线方向形变进行解算,得到垂直方向的形变值。GPS观测值和SBAS-InSAR垂直方向沉降值结果对比如图5所示。由图5可知:A点最大绝对误差、均方根误差分别为2.8 mm、1.6 mm;B点最大绝对误差、均方根误差分别为2.3 mm、1.5 mm。两种监测方法所呈现的形变趋势总体上较为一致,拟合度较高,说明SBAS-InSAR技术在该区域监测精度较理想,适用于尾矿库沉降监测。

图5 SBAS-InSAR与GPS形变监测结果对比Fig.5 Comparison of the deformation monitoring results of SBAS-InSAR and GPS

2.3 时序InSAR形变分析

为了清晰地观察西果园尾矿库沉降趋势,绘制出G1、G2、G33个点包含60个沉降记录的沉降曲线,可以看出库内处于下沉状态,且形变量存在周期性波动,如图6所示。

图6 西果园尾矿库沉降曲线Fig.6 Subsidence curves of Xiguoyuan Tailings Pond

在图6所示的沉降曲线中可以发现其剧烈波动时间段恰好与降雨时间吻合。本研究将沉降曲线进行差分处理,得到形变的周期项,其波动原因主要受外部因素影响[22],所以结合研究时段内的月降水量进行分析,结果如图7所示。

图7 尾矿库形变量波动与降水的关系Fig.7 Relationship between deformation fluctuation and precipitation of tailings pond

由图6及图7可知:形变量的波动与降水量呈正相关,在2018年及2019年8月前后降水量最大时,周期项波动最为剧烈,波动最大值达到12 mm,在2020年6月雨季开始时,曲线开始出现抬升现象,而在10月至6月非雨季时间段内,其周期项在零值上下小幅波动。这主要是因为排土场及库区内主要物质为废土及尾矿砂,大多含有黏土颗粒[23],在雨季期间降水量大时,地表吸水膨胀,并且西果园尾矿库为山谷型尾矿库,大量的降水会降低坝体的黏结系数,使得坝体稳定性降低,短期内出现抬升现象,随后由于质量增加,沉降速率较降水前也会变大,所以沉降曲线会出现先抬升、再加速下降的现象,在降水减少时,库内水分由于蒸腾作用逐渐流失,沉降速率及其波动曲线也趋于平稳。由此可见,降水对库区沉降影响较大,在后续的预测中应考虑这一波动现象。

3 GS-LSTM模型构建及预测结果分析

在构建LSTM网络模型时,众多超参数的选取对模型的预测精度会产生较大影响,为了提高预测精度,本研究采用网格搜索算法搜索全局最优解[24]。网格搜索算法是穷举搜索的一种,在所有候选的超参数选择中,通过循环遍历,将超参数互相组合,尝试每一种可能性,以得到全局最优的超参数组合。

本研究将前50期时序沉降数据作为训练样本,后10期数据作为测试样本进行沉降预测。由于数据集不够庞大,数据复杂度相对较低,为避免模型过拟合,搭建了单层LSTM神经网络,在使用GS算法调参后,最终设置批次大小batch_size=2,窗口大小window_size=50,神经元个数num_units=128,学习率为0.001,迭代次数为4 500次,将均方根误差作为损失函数,训练过程中采用Adam优化算法更新网络权值。为验证该模型预测结果的准确性,将其与传统的SVR模型和灰色GM(1,1)模型结果互相比对,逐一分析,结果见表1至表4。

表1 G1点各预测模型结果对比Table 1 Comparison of prediction results of each model at G1 point mm

表2 G2点各预测模型结果对比Table 2 Comparison of prediction results of each model at G2 point mm

表3 G3点各预测模型结果对比Table 3 Comparison of prediction results of each model at G3 point mm

表4 各模型精度对比Table 4 Comparison of accuracy of each modelmm

由表1至表4可知:3种预测方法与沉降值均较为一致,为更直观地分析及评价各模型的精度,分别计算各预测模型的平均绝对误差(MAE)和均方根误差(RMSE),其中GS-LSTM模型最大MAE及RMSE为2.51 mm、2.90 mm,SVR模型分别为3.58 mm、4.48 mm,GM(1,1)模型分别为4.01 mm、4.95 mm。在两种评价指标上GS-LSTM模型预测结果误差均为最小,优于传统的SVR模型及GM(1,1)模型。

3种模型预测结果与实测值对比如图8所示。由图8可知:SVR模型和GM(1,1)模型可以较好地预测沉降趋势,但是尾矿库的沉降特点及波动性无法得到有效预测,在G3点表现得尤为突出,沉降波动剧烈时其误差也随之增大,而GS-LSTM模型在预测时即使波动剧烈,仍能保持较高的精度,预测结果十分理想,为后续的预测模型持续优化与应用提供了思路。

图8 不同监测点处3种模型预测精度对比Fig.8 Comparison of the prediction precision of three models at different monitoring points

4 结 论

(1)为有效预防和减少尾矿库事故的发生,以鞍山西果园尾矿库为例,利用SBAS-InSAR技术对库区进行沉降监测,提取连续沉降信息,在此基础上使用GS-LSTM模型对库区沉降值进行预测,实现了尾矿库沉降监测与预测的一体化。研究表明:本研究所提方法的监测值与实测值的误差较小,预测值与传统模型相比不仅精度更高,并且能记忆尾矿库的沉降波动特点,监测精度和预测精度均满足工程需要,可为尾矿库沉降规律研究及沉降预测提供依据。

(2)后续工作中考虑使用层析SAR技术结合更高空间和时间分辨率雷达数据获取库区三维立体形变信息,进一步提升尾矿库监测与预测模型的精度。