AM-LSTM网络的北京平原东部地面沉降模拟

曹鑫宇，朱琳，4，7，宫辉力，4，7，郭琳，4，7，尉毓姣，郭涛，陈蓓蓓，4，王海刚，李蕙君

1.首都师范大学 资源环境与旅游学院,北京 100048;2.首都师范大学 水资源安全北京实验室,北京 100048;3.首都师范大学 城市环境过程与数字模拟国家重点实验室培育基地,北京 100048;4.首都师范大学 地面沉降机理与防控教育部重点实验室,北京 100048;5.四川省农业科学院 遥感与数字农业研究所,成都 610066;6.中国地质环境监测院,北京 100081;7.自然资源部京津冀平原地下水与地面沉降野外科学观测研究站,北京 100081

1 引 言

地面沉降是在自然和人为因素作用下，土体压缩而导致地面标高降低的环境地质现象，是一种不可补偿的永久性环境和资源损失（郑铣鑫等，2002）。目前世界上已有60多个国家和地区发生地面沉降，中国有50 多个城市正遭受地面沉降灾害影响（贾三满等，2019）。地面沉降成因复杂，治理难度大，严重制约了区域可持续性发展，成为现代城市的重要安全隐患。

准确模拟预测不同条件下的地面沉降是地面沉降风险预警的基础。地面沉降模拟预测模型主要分为两类：基于物理机制的模型和基于数据驱动的模型。基于物理机制的模型主要依据太沙基有效应力原理和Biot 固结理论建立数值模型（Di Francesco，2013；白林 等，2015；Rogers和Chung，2016；骆祖江等，2018），其优点在于具有明确的物理含义，可解释性好。然而，此类模型需要详细的地质结构信息、清晰的水文地质条件（初边值条件、源汇项等）以及水力学和土力学参数，从而限制了地面沉降模拟预测的准确性（吴吉春和陆乐，2011）。基于数据驱动模型的地面沉降模拟预测方法，得益于统计学原理的支撑和人工智能的发展。基于统计学发展的沉降模拟预测模型主要包括：线性回归模型（范珊珊等，2013；姜媛等，2015）和灰色模型（林跃忠，2000；郭小萌等，2013；吕传振和安动动，2021）。其中，灰色模型能够在少量样本数据条件下进行短期预测，但不适用于中长期预测（闫宏亮和马得花，2020）。基于人工智能技术构建地面沉降预测模型的方法主要包括：BP（Back Propagation）神经网络（周复旦等，2011；Zhu等，2013；成枢等，2015；Wang等，2019）、径向基神经网络（王忠忠和钱为民，2006；谭洋 等，2016）、小波神经网络（蔡东健等，2010；赵凤阳，2016）、循环神经网络RNN（Recurrent Neural Network）（李洛宾等，2020；岳振华 等，2020；刘青豪 等，2021）。RNN 模型近年在地面沉降预测中得到应用。李洛宾等（2016）利用两种RNN 模型（长短时记忆网络LSTM（Long Short-Term Memory）、门控循环单元GRU（Gated Recurrent Unit））和传统BP神经网络模型对隧道沉降进行预测分析，结果显示基于RNN 的沉降预测结果优于传统BP神经网络模型。岳振华等（2020）利用RNN 模型对聚类后的沉降数据分别建立预测模型，结果发现对沉降趋势相同的永久散射体点PS（Permanent Scatterer）时间序列建模结果精度较高。刘青豪等（2021）采用LSTM模型对基于PS观测的地面沉降结果进行预测，并和基于支持向量回归、多元线性回归、二次指数平滑、自回归移动平均的地面沉降预测方法进行比较，发现基于LSTM 模型预测精度最高，结果的平均绝对误差降低27%。上述RNN 模型大多仅从观测点处的累计沉降时序信息中挖掘沉降发展规律，未考虑其他影响因素对地面沉降的影响。Li 等（2020）利用LSTM 模型并结合多影响因素对北京平原沉降进行模拟，结果发现地面沉降严重区域模拟精度低于沉降轻微区域。但是，该研究未考虑不同因素对地面沉降的差异性贡献。为了提高RNN 模型对地面沉降的模拟预测能力，基于历史沉降数据并考虑多影响因素的不同作用程度来构建地面沉降模型是非常有必要的。

LSTM 模型属于循环神经网络的一种，能够记住长时间序列的历史信息，适用于处理非线性变化的地面沉降时序数据。但LSTM 模型对多个输入因子平等对待，不能对输入因子进行有目的的学习（谷丽琼等，2020）。注意力机制AM（Attention Mechanism）源于对人类视觉的研究，人眼会重点关注对图像分类结果影响最大的相关区域，最早在视觉图像领域中提出（Tsotsos 等，1995）。Bahdanau 等（2014）将注意力机制应用于机器翻译任务，之后便成为研究热点，被广泛应用于各种自然语言处理任务中。机器学习领域衡量当前输入因子重要性的方法是权重，而AM 通过关注输入序列本身并发现输入序列之间的内在联系，将注意力集中在多个输入信息中对当前预测结果更关键的信息上，为输入序列分配权重系数，将各序列加权求和，从而准确识别关键特征，提高模型效率和准确性。传统的LSTM 网络中，神经元将输入因子与先前的输出相结合，不能很好地识别关键输入因子对预测结果的影响。谷歌团队（Mnih 等，2014）首次将AM 及LSTM 模型结合起来应用在图像分类领域，结合后的模型能够关注到图像重点区域，有效地减少了图像噪声对分类结果的影响。目前，结合注意力机制的长短期记忆网络AM-LSTM（Attention Mechanism Long Short-Term Memory）已被成功应用到交通流量、水文预报及空气污染等领域的模拟预测研究（Liang 等，2018；Liu 等，2020；Ding 等，2020），但尚未应用在地面沉降领域。

自20 世纪90年代以来，地下水开采量以每年25 亿m3左右持续增加，导致地下水水位平均埋深从1978年的6.4 m 增加至2010年的24.9 m（杨勇等，2013）。长期超采地下水导致水位大幅下降，引发了严重的地面沉降。截至2013年，北京平原区累计沉降量超过50 mm 的区域面积已超4300 km2，约占平原区总面积的68%，最大累计沉降量达到1495 mm（雷坤超等，2016）。合成孔径雷达干涉InSAR（Interferometric Synthetic Aperture Radar）测量技术能够获取毫米级精度的地表形变结果，Ferretti等（2000）提出的永久散射体合成孔径雷达干涉测量PS-InSAR（Persistent Scatterer Interferometric Synthetic Aperture Radar）技术，能够克服传统InSAR技术的时空失相关及大气延迟等方面的不足，有效提高形变监测的精度。很多学者基于PS-InSAR技术对北京地面沉降演化、成因机制、数值模型及风险评估方面进行了深入研究（Zhu 等，2015，2020；Guo等，2019；Lyu等，2020；Li等，2021）。北京平原沉降空间差异性明显，已形成南北两大沉降中心，以朝阳和通州东北部最为严重（Zuo等，2019）。地下水超采是地面沉降的主要因素，第二承压含水层对地面沉降的贡献最大（Zhou 等，2019；Chen等，2020）。

本文以北京平原东部为研究区，基于PSInSAR 技术获取研究区2010年11月—2016年8月地面沉降信息；在此基础上选择不同沉降发育地区的典型位置，结合不同层位的地下水水位时序数据，构建基于AM-LSTM 的地面沉降模拟模型，并与传统LSTM 模型进行对比分析；同时借助注意力机制探究不同含水层水位对沉降模拟结果的影响。研究结果能够为区域地面沉降防控工作提供必要的技术支撑。

2 研究数据及方法

2.1 研究区概况

研究区位于北京平原东部（39°43′N—40°09′N、116°22′E—116°54′E），地处潮白河与永定河冲洪积扇中下部，面积约5301 km2，包括东西城区、朝阳区、通州区、顺义南部及大兴区北部（图1）。区内属大陆性季风气候，1978年—2019年平均气温11.3 ℃，平均降水量约535.2 mm，年内降水分配不均，夏季炎热多雨、冬季寒冷干燥。第四系松散沉积物分布广泛，岩性多为砂砾石、砂、黏土层相互交错出现，并以黏土为主（周毅等，2016）。含水层系统主要包括：潜水含水层和第一、第二和第三承压含水层。其中，潜水含水层深度小于50 m，岩性主要为中细砂、粉砂、砂质黏土，该层地下水接受大气降水、农田灌溉入渗和河流入渗等补给；第一和第二承压含水层深度分别为50—100 m、100—180 m，由多层砂砾石、砂、黏土组成，其中，第一承压含水层主要用于农业开采，第二承压含水层主要用于工业、生活开采，为北京平原地面沉降的主要贡献层（雷坤超等，2016）；第三承压含水层深度大于180 m，该层地下水主要用于生活，少量用于工业，开采程度相对较低，因此本文不考虑第三承压含水层对地面沉降的影响。研究区朝阳—通州一带是北京平原区最大的沉降漏斗，2003年—2015年间沉降速率均大于110 mm/a，最大累计沉降量已经超过1500 mm（Zhou等，2018）。

2.2 数据源与数据处理

本文采用的数据包括：2010年11月22日—2016年8月10日共54景RADARSAT-2影像；2010年10月—2016年12月6 个地下水水位长期观测井监测的潜水、第一承压含水层及第二承压含水层的月监测水位数据（表1）；6 个2012年水准监测点数据。空间位置分布见图1。

图1 研究区分布图Fig.1 Distribution of study area

表1 6个地下水水位监测点监测深度信息Table 1 Monitoring depth information of six groundwater level monitoring points

本文采用PS-InSAR 技术（Ferretti 等，2000）处理SAR 影像数据，通过选取受时空去相关影响较小的PS 点，对这些点进行时序相位分析，从而获取可靠的地表形变信息，测量精度可以达到毫米级。PS-InSAR 技术的主要流程包括：（1）主辅影像配准；（2）生成差分干涉图序列及选取PS点；（3）大气延迟相位估计及去除；（4）PS 点形变速率估计。由于PS-InSAR 测量结果为雷达视线向LOS（Line of Sight）的形变结果，是地表真实三维形变在LOS 向上的投影。考虑到北京平原区水平方向的形变相对于垂直方向的形变小（Ng 等，2012），本文忽略水平位移的影响，依据式（1）将PS-InSAR 技术处理得到的LOS 向形变信息投影至垂向，获取研究区垂向地表形变结果。

式中，Dv为垂向形变速率（mm/a）；Dlos为LOS 向形变速率（mm/a）；θ为雷达入射角。

考虑到水位监测点位置处PS 点较少，本文提取水位监测点位置处300 m 缓冲区内PS 点监测的沉降序列信息（图2），并将80%的PS 点数据作为AM-LSTM 模型的训练样本。各点缓冲区内PS点数量不同：A 点、C 点及D 点周边区域多为农田、果园等种植区，缓冲区内的PS点数量较少，在［15，40］范围内；而B 点、E 点和F 点位于村庄或城镇区域，缓冲区内的PS点相对较多，在［120，185］范围内。B 点和C 点缓冲区内PS 点监测的地面沉降较为严重，最大累计沉降量分别达到768 mm 和870 mm；E点和F点沉降速率较为缓慢，最大累计沉降量分别为112 mm和199 mm。

图2 各水位监测点缓冲区内PS点沉降时间序列Fig.2 Subsidencce time series of PS points in the buffer zone of each water level monitoring point

由于地下水水位数据是逐月观测获得（图2），为保证InSAR 数据与地下水水位数据时间一致，本文对获取的PS 点时序沉降信息进行筛选，针对同一月份多景影像取日期靠后一景。同时为减少数据时间间隔对模型模拟带来的影响，对2010年11月与2010年12月两期形变信息不做考虑。最后选取2011年6月—2016年8月共43 个月的沉降序列数据，利用对应月份的不同层位地下水水位作为模型输入因子。考虑到地下水位和地面沉降的不同量纲，且相应的值域范围相差大，本文采用标准差标准化方法（谷丽琼等，2020）对水位数据和沉降数据进行预处理，将原始数据进行线性变换使其值域映射在［-1，1］范围内。完成模型训练后，对模型的输出值进行反标准化处理，将输出值还原成实际模拟值。

2.3 研究方法

2.3.1 LSTM

LSTM 网络（Hochreiter 和Schmidhuber，1997）是RNN 的改进模型，通过引入门的概念来控制长期状态，增强了对数据权值的控制，从而有效避免传统RNN 中的长期依赖问题，在语音识别及图像处理等方面取得了广泛的应用（Vinyals 等，2015；Soltau 等，2016）。LSTM 包括输入层、隐藏层和输出层（图3）。

图3 LSTM单元结构Fig.3 LSTM cell structure

图3 中ft、it和ot是3 个非线性门函数，称为遗忘门、输入门和输出门。不同于以往的RNN 直接记忆全部历史信息，LSTM 通过3 个门对历史信息进行有选择地更新和删除。计算过程如下（Greff等，2017）：

式中，Wf、Wi、Wc和Wo分别为输入向量xt的权重矩阵；bf、bi、bc和bo为t时刻的偏置向量；C͂t是当前时刻的神经元候选状态，包括t-1 个时刻的隐藏层状态ht-1和神经元状态Ct-1；Ct是当前的神经元状态；ht是当前时刻的隐藏层状态；σ和tanh 是非线性激活函数。

2.3.2 AM-LSTM

不同层位地下水水位数据在时间序列下对沉降的影响程度不同，AM 可以对输入至LSTM 模型中的不同层位地下水水位赋予不同的权重，从而增强对地面沉降影响更大的地下水位序列的关注程度，提高沉降模型的模拟精度。本文构建的AM-LSTM 地面沉降模拟模型（图4），由两部分组成：AM模型和LSTM 模型。AM模型的输入序列是不同层位地下水水位序列，注意力层提取不同层位地下水水位时间序列的注意力权重信息，并对不同层位地下水水位时序数据进行加权平均；然后，使用AM 模型的输出作为LSTM 模型的输入来模拟形变结果，最后通过全连接层输出结果。

图4 AM-LSTM网络框架图Fig.4 AM-LSTM network frame diagram

其中，注意力机制计算过程如下（Lin等，2020）：

式中，向量Vl、Wl和Ul是模型的可学习参数；Bl是偏置向量；表示第k个输入特征序列在时间t处的注意力权重。softmax 函数对标准化以确保各特征序列在t时刻下的注意力权重之和为1。

（2）根据注意力分布对原始输入序列进行加权平均。在时间t获得注意力模型的输出，即加权输入特征序列et如下：

式（8）—（10）对输入地下水水位时间序列赋予相应时间t下的注意力权重值，将得到含有注意力权重的新序列et输入至LSTM 网络中进行训练并模拟。

3 结果与讨论

3.1 PS-InSAR结果验证及空间分布

3.1.1 水准验证

利用收集的6 个2012年水准点数据对PSInSAR 监测结果进行验证，将水准点200 m 缓冲区内的PS 点同期平均形变值与水准测量值进行比较（图5）。验证结果表明，二者形变差值的最大值为5.15 mm/a，最小值为0.47 mm/a，平均误差为0.67 mm/a，R2为0.98，水准监测结果与PS-InSAR结果呈现较好的一致性。此外，本文与以往文献对该地区的监测结果进行对比。InSAR监测结果均表明（Zhou等，2018，2019；Lyu等，2020），北京平原2011年—2015年最大沉降速率大于140 mm/a，沉降主要分布在昌平沉降漏斗、来广营沉降漏斗及朝阳—通州沉降漏斗，与本文PS-InSAR 监测结果相似且空间分布具有一致性。

3.1.2 形变速率空间分布

由研究区形变速率分布图（图6）可以看出，2010年11月—2016年8月研究区地面沉降速率空间差异性明显，地面沉降主要分布在朝阳、通州、海淀、昌平及顺义等部分地区，并且以朝阳东部和通州西北部地区地面沉降最为严重，是研究区内的主要沉降漏斗。该区域地下水开采导致地下水水位下降，最大降幅达到20 m，累计可压缩层厚度超过200 m，使得地面沉降发育严重。研究区内最大沉降速率约153 mm/a，累计沉降量达到1063 mm，位于朝阳三间房附近。区内抬升速率总体较小，最大抬升速率约4 mm/a，位于通州永乐店镇。年沉降速率超过50 mm/a 的PS 点数量占总数的12.4%，面积为1162.9 km2；其中，超过110 mm/a的PS点数量占总数的2.9%，面积为110 km2，主要分布在朝阳东北部（金盏乡附近）和东南部（三间房附近）以及通州西北部（宋庄镇—永顺—梨园镇一带）。东西城区、大兴区北部、顺义区东部及通州区南部沉降速率较小，地面沉降发育缓慢。

图6 研究区形变速率结果Fig.6 Results of subsidence rate in the study area

本文选取的6个水位监测点（典型点）位处于不同发育级别的沉降区域，可用于分析不同沉降区的模型模拟结果。C点处地面沉降最为严重，沉降速率大于110 mm/a；B点处的沉降速率处于80—110 mm/a；A点和D点处的沉降速率为50—80 mm/a；E 点和F 点处的地面沉降相对缓慢，沉降速率为0—30 mm/a。

3.2 AM-LSTM 模型训练与分析

3.2.1 AM-LSTM 模型参数与训练过程

将6 个典型点位缓冲区内的PS 点数据按比例拆分，随机选取80%作为训练集用于AM-LSTM 模型的特征学习，20%作为测试集用于评估模型模拟能力。在训练神经网络时，考虑到每次训练所选取样本数—训练批次（Batch Size）的大小影响模型的训练时间和准确性，数值过大或者过小都会导致模型的结果误差增加，本次设计多个实验以选取最优Batch Size（表2）。经过反复实验，根据模型模拟结果的评定指标RMSE，最终确定模型训练中在各点处的Batch Size分别设置为8、16、8、8、16 和32。利用Adam 优化算法（Adam optimization algorithm）（Kingma 和Ba，2014）对模型参数进行自动调优，在模型迭代（epochs）100 次、学习率0.05、LSTM 神经元数量16 时，得到最优结果。本文分别对6 个典型点300 m 缓冲区的地面沉降观测数据进行建模。首先，AM 作用于每个时间节点，获取不同层位的地下水水位在不同时间节点上的Attention权值，该数值表示水位对当前时刻沉降模拟值影响的大小。然后，通过权值变化来捕捉时间序列建模过程中，不同层位地下水水位对沉降的影响程度，利用AM 模块得到的含有注意力权重信息的地下水水位序列输入至LSTM 神经单元，并对沉降进行模拟。

表2 各典型点不同训练批次大小下AM-LSTM 模型RMSE对比Table 2 RMSE comparison of AM-LSTM model under different training batch sizes at each typical point

为了减少AM-LSTM 训练过程中的学习误差，将前一次迭代的误差反馈到网络中，对权重进行优化，得到的模型训练损失曲线如图7所示。A-D点loss 曲线快速趋于稳定；而E-F 点由于沉降缓慢，形变特征不明显导致收敛较慢。为了能够更快的收敛且防止模型出现过拟合，采用Early-Stopping 技术（Prechelt，1998）对AM-LSTM 模型进行训练，训练集每完成一次epoch 迭代后都用损失函数（loss）来衡量测试集中的误差，并将此时的网络参数保存。

3.2.2 基于AM-LSTM 的沉降模拟结果分析

根据最优参数设置采用未加入注意力机制的LSTM模型与AM-LSTM模型分别对PS点的形变进行模拟，将模拟的形变结果分别与真实PS 点形变结果对比（图8），并统计不同模拟结果误差（表3）。可以看到，利用AM-LSTM 模型模拟的形变值与真实值拟合度更高，能够较为准确地捕捉到各时间节点的形变信息。F 点处模型拟合曲线和PS 点真实形变曲线与其他点处拟合结果对比，具有明显的差异性。原因是由于该点处地面沉降发育缓慢（0—30 mm/a），与其他典型点相比月形变量小，使得毫米级的误差在形变曲线中体现的较为明显。而处于相同速率区间的E 点处的RMSE 低于F 点处，其模拟曲线拟合度优于F点处。虽然无法在每个时间节点上对F点准确模拟，但仍能够模拟得到相似的沉降变化，且拟合程度较好于传统LSTM模型。

图8 AM-LSTM与LSTM对比结果Fig.8 Comparison results of AM-LSTM and LSTM

由各点处AM-LSTM 模型与LSTM 模型的误差（表3） 可以看出，AM-LSTM 模型在各点处的RMSE 误差值均低于LSTM 模型，模拟精度最高可提升22%。数据驱动模型的驱动力是数据集，训练样本的数量越多，有助于模型在训练过程中对时序特征的学习，进而提高模型模拟精度。由表3可知：B、E、F 点处的PS 点数量远高于A、C、D点，模型模拟结果的RMSE 整体低于A、C、D 点，原因可能是B、E、F 点周边的沉降观测值多，能够输入模型的训练样本相应增多，所以模型对沉降时序特征学习能力相对增强，对形变的模拟精度较高。位于沉降速率严重的B、C 点，尽管B 点处沉降速率较大，但训练样本多，模型模拟结果RMSE 低于C 点处。A、C、D 点处PS 点数量整体较少，其中，C点所在区域沉降发育最严重，累计沉降量中位数（780 mm）高于A 点（351 mm）和D 点（398 mm），其模拟结果RMSE（2.31 mm）高于A 点（1.75 mm）和D 点（1.49 mm）。E、F 点处沉降发育缓慢（0—30 mm/a），PS 点数量整体较多，虽然F 点处PS 点数量高于E 点处，但其累计沉降量中位数（158 mm）高于E点处（80 mm），模型模拟结果RMSE（1.23 mm）高于E点（0.35 mm），说明沉降量级在一定程度上影响模型模拟精度。

表3 各典型点LSTM 与AM-LSTM月形变模拟误差对比Table 3 Error comparison of monthly deformation simulation between LSTM and AM-LSTM for each typical point

3.3 不同层位地下水水位对沉降的贡献率

为了研究典型点处不同层位地下水水位变化（潜水、第一承压含水层、第二承压含水层）对沉降模拟结果的影响程度，本文将AM-LSTM 模型训练过程中各时间节点处各层位地下水水位特征的注意力权重值进行加权平均，该数值表示模型训练过程中各层地下水水位特征对沉降模拟结果的影响程度（图9）。发现A—F 点在地面沉降模拟过程中，第二承压含水层较其他层位地下水水位所占注意力权重最高，分别为0.39、0.54、0.59、0.61、0.62 和0.43。说明第二承压含水层对沉降模拟结果影响最大，是现阶段地面沉降的主要贡献层位。第二承压含水层是该地区农村生活用水和工业用水的主要来源，地下水水位多年来呈快速下降趋势。该发现与前人研究结论一致（雷坤超等，2016；Chen等，2020）。

图9 各层地下水位注意力权重Fig.9 Weight of attention of groundwater level at different levels

4 结 论

本文基于PS-InSAR 技术监测北京平原区东部2010年11月—2016年8月地表形变结果，结合同期不同层位的地下水水位时序数据，采用AMLSTM 模型对地面沉降时间序列进行模拟。该模型利用注意力机制对各层地下水水位赋予不同的权重，并将含有注意力权重的地下水水位序列输入至LSTM模型并模拟月沉降变化量。主要结论如下：

（1）研究区2010年11月—2016年8月地面沉降空间差异性明显，区内主要沉降漏斗位于朝阳东部和通州西北部附近，最大沉降速率约153 mm/a，累计沉降量达到1063 mm，位于朝阳区三间房乡附近。

（2）与传统LSTM 模型相比，AM-LSTM 模型模拟结果与真实沉降序列拟合程度更高，模拟精度最高提升22%，说明本文构建的模型能够对时序的地面沉降做出更准确的模拟。

（3）各监测点区域不同层位地下水水位的注意力权重表明，第二承压含水层地下水水位特征所占注意力权重最高，贡献率最大，对地面沉降模拟结果的影响最大。

本文构建的AM-LSTM 模型输入因子仅为3 个层位的地下水水位，然而城市化带来的动静载荷也会影响地面沉降的时序演化。因此，后续研究将增加模型的输入因子。此外，本文仅对地下水水位监测点300 m 缓冲区内的PS 点进行模拟，没有对未知时间段进行预测，后续将采用经典的地下水数值模型模拟未来不同条件下的水流场，之后再利用构建的AM-LSTM 模型对未来时刻的沉降进行预测。