[0001] 本发明涉及地表形变监测领域，具体涉及一种基于InSAR和双向门控循环单元的地表形变预测方法。

[0002] 地面沉降是一种累进性且发育过程缓慢的地质灾害，一旦形成难以恢复。随着地面沉降的持续发展，极易造成城市路面和路基形变，甚至引起建筑物受损、城市基础设施破坏、地下工程等多种危害，严重威胁城市居民的生命财产安全，制约社会经济的可持续发展，地面沉降作为地质灾害类型之一，严重影响基础设施的结构健康，可见，探索行之有效的地面沉降预测手段对于基础设施损害的早期预警和及时采取补救措施具有重要意义。

[0003] 根据是否考虑环境因素，目前的SD预测研究可分为两类，第一类是仅根据GPS、水准测量数据或多时相干涉合成孔径雷达（Multi‑Temporal Interferometric Synthetic Aperture Radar ，MT‑InSAR）地表形变（SD）监测数据构建SD预测模型，例如，已有研究学者基于永久散射体干涉合成孔径雷达（Persistent scatterer interferometric synthetic aperture radar，PS‑InSAR）SD监测数据构建了长短期记忆（Long short‑term memory，LSTM）SD预测模型，结果表明，LSTM 的短期预测性能更好；已有研究学者基于小基线子集干涉合成孔径雷达（Small baseline subset interferometric synthetic aperture radar，SBAS‑InSAR）技术获得的SD监测结果，提出了一个CNN‑PhLSTM模型，结果表明，CNN‑PhLSTM预测模型在多项评价指标上优于支持向量回归模型（Support vector regression，SVR）和CNN‑LSTM模型，其长时序预测结果更符合InSAR原始SD趋势；已有研究学者利用独立分量分析（Independent component analysis，ICA）从原始InSAR序列中分离出独立的形变信号，并构建了LSTM预测模型，其在两个形变区域的预测结果都比LSTM更准确；然而，上述研究并未考虑自毁监测数据与环境因素（EF）之间的长时序依赖关系，使得预测结果缺乏可信度和有效性，此外，未考虑环境因子变化而构建的SD预测模型在预测SD时假设EFs（归一化差异植被指数（Normalized difference vegetation index，NDVI）、地表温度（Land surface temperature，LST）和降水量（Precipitation）等）和沉降区地质条件不变，导致预测精度不高；第二类SD预测模型是在第一类预测模型之上，考虑了形变区域中EFs与SD之间的响应关系；例如，已有研究学者基于transformer中的自注意（Self‑attention，SA）机制，综合考虑了EFs与SD之间的关系，较为准确地模拟了非季节性信号和季节性信号，实现了盐湖周边SD的预测；已有研究学者利用地理和时间加权回归方法，结合SD监测数据得到了EFs的相对权重，同时，构建了基于LSTM和注意力机制（Attention mechanism，AM）的SD预测模型，该模型可考虑时空非平稳关系，并获得了比第一类预测模型更好的预测结果；但是，上述SD预测模型没有考虑InSAR SD和EFs的多尺度特征，使得他们的SD预测结果缺乏可信度和有效性。

[0016] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

[0017] 本发明实施例的一种基于InSAR和双向门控循环单元的地表形变预测方法流程示意图如图1所示，包括：利用InSAR方法获取地表形变数据和环境因子特征数据并进行时序差分，根据得
到地表形变时序序列和环境因子特征时序序列建立原始数据集；
在本发明一个具体实施例中，为实现有效可靠的SD（地表形变）预测，本发明以青藏高原冻土区四个典型区域为例，首先通过时间序列插值、按月填充等方法构建数据集，如表1所示：
表1 输入模型的数据信息

[0018] 表1中SD数据于2022年4月24日的图像存在缺失，本发明首先利用InSAR观测数据（如图2（b）所示）提取了原始InSAR SD图像和EFs（环境因子）图像（如图2（a）所示），本发明实施例中，环境因子包括归一化差异植被指数（Normalized difference vegetation index，NDVI）、地表温度（Land surface temperature，LST）以及降水量（Precipitation）），分别对InSAR SD、LST和NDVI进行时序差分，并对降水量进行了月度填充，上述操作使EFs和InSAR SD在时间步长上严格对齐，从而使MHCA（多头交叉注意）机制能够在每个时间步长上获得InSAR SD特征与EFs之间的依赖关系；月度填充是将月平均降水量填充到当月SD时间步长的操作；本发明实施例中时序差分公式如下：
；
式中， 和 分别是前一个时间点和相应的LST， 和 分别是后一
个时间点和相应的LST， 是需要插值的时间点。

[0019] 由于选取的四个典型区域位于冻土带，具有冻胀融沉的特点；因此，本发明实施例为确保时序重叠分割结果具有一定的周期性和足够的模型训练样本，采用快速傅立叶变换（fast Fourier transform，FFT）计算时间步长，选择较小的频率作为时间步长，并确保所选频率对应的相位角绝对值最小。此外，时间步长的经验值一般在15左右。

[0020] 构建基于双向门控循环单元的多输出地表形变预测网络；所述多输出地表形变预测网络包括全局‑局部特征提取模块、多因素交叉注意残差模块和嵌入自注意机制的局部残差模块；如图3所示，本发明构建了一个双向门控循环（Bi‑directional gated 
recurrent，BiGRU）单元多输出SD预测网络（GLER‑BiGRUnet），该网络包含创新的全局‑局部特征提取（Global‑local feature extraction，GLFE）模块、多因素交叉注意残差（Multi‑factor cross‑attention residual，MCAR）模块和嵌入自注意机制的局部残差（Local residual module embedded in self‑attention mechanism，RCSA）模块，此外，维度提升（Dimension raising，DR）模块利用时间分布dropout（TD‑Dropout）和时间分布密度（TD‑Dense）的串联结构来抑制模型过拟合，恢复RCSA模块输出结果的维度，GLER‑BiGRUnet的主要模块详情如下：
GLFE模块：在神经网络训练过程中，选择合理的特征提取方法对于提高模型特征
学习效率和增强泛化能力至关重要，在SD时序预测任务中，现有研究一般直接将原始特征维度通过 Dense层进行降维和全局特征提取，然而，原始特征维度与RNN层输入维度之间的差异过大，可能会导致局部形变特征的遗漏，这会导致模型拟合不足，影响预测精度；因此，本发明结合密集连接层（Dense）和一维卷积层（Conv1D）构建了一个GLFE模块，如图4所示，该模块旨在捕捉时间序InSAR SD和EFs的全局和局部特征；在Dense和Conv1D层中，InSAR SD和LST使用的激活函数为TanH，而降水和NDVI则使用ReLU，此外，当将Conv1D叠加为两层时，模型有过拟合的趋势，这是因为两个Conv1D层提取了离群特征或噪声，导致模型在验证集上的泛化能力下降。

[0021] 模块作用于形变特征维度，模块公式如下式：；
其中， 是时序InSAR数据和EFs的特征向量， 是权重矩阵， 是偏置
向量。

[0022] ；其中， 表示输入密集连接层原始数据集的特征向量， 为权重矩阵，
为偏置向量， 表示一维卷积层输出结果的相应元素， 表示一维卷积核对应的位置，为一维卷积核的大小， 为卷积核在位置 的权重， 为
的滑动窗口；本发明中， 的值为3。

[0023] MCAR模块：模型表达式为：；
上式中， ， 和 分别表示LST，NDVI和
precipitation经过GLFE模块后的输出。

[0024] 如图3中（c）所示，在MCAR模块的Line 1中，由GLFE模块处理的SD和EFs数据被一起输入MHCA，以获取不同注意头SD和EFs之间的关系，经GLFE模块处理的EFs数据被输入BiGRU，在Line 2中进行建模，本发明实施例中MHCA的结构示意图如图5所示，MHCA首先通过八个注意头的权重矩阵（ ， 和 ）将SD和EFs映射到多个不同的表示子空间，得到相应得查询（Quary， ）、键（Key， ）和值（Value， ）向量，其计算表达式为：；
接下来，对每个注意力头使用基于点积交叉注意力的高度优化矩阵乘法，计算
和 之间的注意力得分矩阵（ ），注意力得分矩阵（ ）还按 向量维数
的平方根缩放，以确保梯度的稳定性，计算表达式如下：
；
随后， 通过 函数转换为交叉注意权重，并通过 向量的加权
求和得到每个注意头的输出计算公式如下：
；
然后，将所有注意力头的输出进行连接，并使用权重矩阵 映射连接结果，确保
连接结果与标清具有相同的维度，最后得到注意力权重特征矩阵（ ）：
；
本发明实施例中MHCA的输出进一步通过TD‑Dense层和TD‑Dropout层，以防止模型过度拟合，如图6所示，同时，减少了输入BiGRU层的参数数量，以提高模型的泛化能力；随后，将这些结果输入3层BiGRU，以进一步探索InSAR SD序列正向和反向的时间序列依赖关系，并更好地提取InSAR SD特征，最后，Line 2和Line 1的BiGRU输出分别通过一个TD‑Dropout层，以抑制模型的过拟合。

[0025] 需要说明的是，本发明实施例中图3的（c）部分，Line 1中BiGRU的输入是MHCA模块输出的注意力权重特征矩阵（Attention weight feature atrix，AWFM），Line 2是GLFE模块中EFs的输出，前向‑后向BiGRU能进一步提高模型对不同时间步之间注意力权重的解释和表达能力，提取出更丰富、更准确的特征表征。

[0026] 假设输入BiGRU的序列长度为时间步长T，当前时间步的输入和输出分别为 、，前向和后向的隐藏状态分别为 、  ，则BiGRU的计算公式如下：前向： ；
后向： ；
BiGRU的输出： ；
其中， 表示前向GRU， 表示后向GRU， 表示更新门， 表示重置门， 表
示点乘， ， ， 和 ， ， 表示当前时间步的权重矩阵和偏置向量， ， ，
表示上一步和下一步的权重矩阵。

[0027] RCSA模块：如图3中（d）所示，该模块中，GLFE模块的输出与BiGRU模块的输出通过残差融合一起输入RCSA模块，这样可以防止梯度消失，保留长距离和多尺度SD特征，提高模型拟合能力；随后，将融合结果输入SA（自我注意）机制进行进一步的特征提取和关联学习，其结构示意图如图7所示，以提高SD预测的准确性和鲁棒性，最后，通过恒等映射残差连接将SA机制的输出和残差添加结果串联起来，以增强模型的特征共享能力、学习效率和参数共享，并减少过拟合。

[0028] 残差融合和SA机制的公式如下：；
；
；
；
上式中， ， 和 代表 的权重矩阵， 代表向量 的维度。

[0029] 损失函数的正确选择对于SD预测模型的训练和性能至关重要。回归任务中常用的两个损失函数是平均绝对误差（MAE）和平均平方误差（MSE），然而，MAE不利于梯度下降，对优化过程中的异常值不敏感，而MSE对异常值敏感，但可能导致优化不稳定，本发明的 SBAS‑InSAR SD结果虽然经过了一系列的时间序列差分和归一化操作，但无法完全避免异常值对模型训练和拟合的影响，此外，在模型训练过程中期望在完全学习到形变特征的同时，还能获得较高的拟合精度，因此，本发明选择对数双曲余弦函数（Logarithmic hyperbolic cosine function，LogCosh）作为基础损失函数；该函数使模型在训练过程中更加稳定（在闭合形式解中平滑且二次微分），并对异常值具有一定程度的鲁棒性；考虑到本发明是一个多输出模型，本发明实施例提出了多目标平均 LogCosh（Multi‑objective average LogCosh，MOA‑LogCosh）和多目标平均决定系数（Multi‑objective average Coefficient of Determination，MOA‑R2）来训练模型，其数学公式如下：；
；
其中， 代表样本数量， 和 代表分别代表第 次观测值的真实值和预测
值， 代表观测值 的平均值。

[0030] 本发明采用均方根误差（Root mean square error，RMSE）、平均绝对误差（Mean absolute error，MAE）、R2和对称平均绝对百分比误差（Symmetric mean absolute percentage error，SMAPE）来评价所构建的GLER‑BiGRUnet模型的性能，其计算方式均可参考现有计算方法，本发明对此不做限制。

[0031] 将所述原始数据集输入所述多输出地表形变预测网络中，获取地表形变预测结果本发明构建的多输出GLER‑BiGRUnet能够考虑EFs的变化并实现长时序预测，基于训练好的SD预测模型，本发明实施例设计了一种考虑EFs的循环预测方法，如图8所示，具体步骤包括：1）利用GLER‑BiGRUnet同时对SD、LST、NDVI和降水进行一步预测（图8（a）），并将预测结果与输入序列进行连接（图8（b））；2）将合并的SD和EFs输入GLER‑BiGRUnet，进行一步预测（图8（b）），并将预测结果与输入序列合并；3）迭代应用上述循环预测模式，并对最终SD进行反时间序列差分和反归一化处理（图8（d））；由于Sentinel‑1A卫星的重访周期为12天，而本发明的长时间序列预测是从测试集开始的，因此本发明的2年SD预测需要72次迭代，即图8（d）中的 P=72。

[0032] 在本发明另一个具体实施例中，为了验证本发明所构建的GLER‑BiGRUnet SD预测模型的预测性能，本实施例拟从模型训练过程，全局‑局部预测结果和不同模型评价指标对比三个方面来进行说明。

[0033] 如图9所示，给出了本发明提出模型的Loss和R2的变化曲线示意图，本发明提出的GLER‑BiGRUnet模型在四个典型区域的Loss和R2都正常收敛到0、1，此外，当epoch位于[4, 40]时，A区域的模型过拟合现象严重，随后有所缓解，当epoch>40时，R2逐渐收敛到1；在此过程中，学习率的调整对抑制过拟合起到了关键作用，说明表3中超参数的设置是合理的；
此外，由于C区域的SD特征数（50596）是A区域（19578）、B区域（24036）和C区域（24726）的2至
3倍，因此Loss和R2最终在迭代55次之后收敛。

[0034] 图10显示了2022年6月23日四个典型区域的SD结果与GLER‑BiGRUnet模型预测结果的对比分析，从图10的第一列和第二列来看，原始SD与预测SD在A、B、C和D地区的空间分布是一致的，在第三列中，四个典型区域的原始SD和预测SD的残差有50%位于[‑1.20，0.81]、[‑1.13，1.66]、[‑1.41，0.89]和[‑1.45，1.71]的范围内，残差的平均值优于±0.5 mm；此外，第三列中A、B、D区域的残差结果更符合正态分布，而C区域的残差结果则集中在[‑
2.00，2.00]区间，上述结果表明，本发明提出的GLER‑BiGRUnet模型能够准确预测不同区域的空间分布和局部剖面。

[0035] 为了验证GLER‑BiGRUnet模型的预测精度，本实施例进一步选择了Conv1D‑MHCA‑BiGRU、MHCA‑BiGRU、MHCA‑GRU、MHSA‑BiGRU、MHSA‑GRU、SA‑BiGRU、SA‑GRU、BiGRU和GRU等模型进行评价指标的比较，比较结果如表2所示：表2 GLER‑BiGRUnet模型与其他预测模型在不同评价指标上的比较

[0036] 可以看到，在仅考虑InSAR SD的模型中，SA‑GRU在A、B和D区域的所有评价指标得分均低于BiGRU，而在C区域的表现优于BiGRU；相反，在SD特征较多的区域，模型复杂度的增加有助于模型更好地学习SD特征，此外，MHSA‑GRU和SA‑GRU的A区和C区也存在类似情况；在考虑了SD和EF的多尺度特征的模型中，本发明提出的GLER‑BiGRUnet模型在所有评价指标上的得分都较高，同时，并没有出现模型复杂度越高，评价指标反而越低的情况；此外，在A区域，GLER‑BiGRUnet模型的RMSE分数比MHCA‑GRU模型（只考虑EFs）高68%；C区域，GLER‑BiGRUnet的RMSE分数比未考虑多尺度SD特征的Conv1D‑MHCA‑BiGRU高24%，此外，Conv1D‑MHCA‑BiGRU、MHCA‑BiGRU和MHCA‑GRU的两两比较表明，考虑局部或双向SD特征可提高预测精度，综上可以看出，本发明所提出的GLER‑BiGRUnet在不同形变区域的预测精度最高。

[0037] 在本发明另一个具体实施例中，获取盐湖至五里段2019年至2022年的88幅Sentinel‑1A图像，这些图像均采用下降轨道、干涉宽（IW）模式、VV极化，入射角和方位角分别为34.69 º和90 º。盐湖至五里段的InSAR形变数据是利用基于可视化图像环境 SARscape 5.6.2平台的SBAS‑InSAR技术获得的；此外，本实施例进一步将视线（LOS）方向的形变监测结果 转换为垂直方向的形变监测结果 并将其用于后续工作，转换表达式为：；
其中， 是入射角度。

[0038] 在对SD监测结果和谷歌地球图像进行分析后，在盐湖至五里段选择了四个典型区域（A、B、C和D）进行SD预测；这些典型区域的特点是地面沉降严重、热岩溶湖面积较大，且被QTR穿越。基于本发明构建的GLER‑BiGRUnet模型，对四个典型区域的时间序列InSAR地表形变进行了模拟和训练，最终的超参数详情如表3所示：表3 典型区域地表形变预测模型的超参数训练

[0039] 本发明利用提出的GLER‑BiGRUnet模型预测了QTR沿线四个典型区域未来两年的SD变化趋势，并将预测结果与其他模型进行了比较，比较结果如图11所示，图中（a）、（c）、（e）和（g）分别是 A、B、C 和 D 地区原始SD平均值和预测SD平均值的变化，（b）、（d）、（f）和（h）分别是四个区域在测试集上原始SD平均值和预测SD平均值的变化；由图11（b）和图11（f），在测试集上，MHSA‑BiGRU，GRU等模型的预测结果与InSAR SD差异较大，预测趋势不稳定；由图11（b）、图11（d）、图11（f）以及图11（h）可以看出，MHCA‑BiGRU等模型在InSAR SD趋势方向发生变化的某些位置明显缺乏特征捕捉能力，相反，所提出的GLER‑BiGRUnet模型考虑了InSAR SD的局部细节和多尺度特征，能够更好地拟合InSAR SD。

[0040] 未考虑EFs变化的模型对未来两年的长时间序列预测结果不稳定，而未考虑多尺度特征的模型预测结果则明显更加稳定，但预测结果并不能准确反映未来SD的变化趋势，特别是在A区域，不同预测模型的预测结果存在较大差异；然而，如图11（a）所示，与其他预测模型相比，所提出的GLER‑BiGRUnet能够学习到所示2021年1月11日、2021年2月28日和2021年4月29日具有陡增或陡减SD趋势的A区域的多尺度形变特征，这使得GLER‑BiGRUnet模型的预测结果与InSAR SD变化趋势的吻合度最高；此外，在B、C和D区域，InSAR SD趋势平稳，不同预测模型的预测结果较为集中，虽然Conv1D‑MHCA‑BiGRU、MHCA‑BiGRU和MHCA‑GRU模型的预测结果与InSAR的形变趋势基本一致，但仍无法预测某些细节位置的SD，相反，所提出的GLER‑BiGRUnet模型在不同区域的预测结果最为稳定，最能代表典型区域未来SD的发展趋势。

[0041] 以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。