随着人类活动的影响，机场及其附近地区出现不均匀形变，对机场安全运营构成了一定的安全隐患。首都国际机场位于北京市顺义区，目前，针对首都机场地表形变的研究基本认为地下水开采是导致地表形变的主要原因^[1]。合成孔径雷达干涉测量（interferometric synthetic aperture Radar，InSAR）技术具有全天候、全天时对形变信息敏感的优势，已被广泛应用于地表形变监测，其中，InSAR时序分析方法主要包括永久散射体技术（persistent scatterer InSAR，PS-InSAR）、小基线集技术（small baseline subsets InSAR，SBAS-InSAR）和相干点目标分析技术（interferometry point target analysis InSAR，IPTA-InSAR）等 ^[2]。基于InSAR技术的北京市地下水开采和地表形变演变现已开展了时空定性定量对比分析的研究^[3]，但未分析地下水和InSAR地表形变观测存在时间滞后对形变预测结果的影响。

随着人工智能技术的发展，深度学习模型也被广泛应用于时序InSAR地表形变监测和预测方面。常用于InSAR地表形变趋势预测的深度学习模型包括：长短期记忆网络（long short-term memory，LSTM）、循环神经网络（recurrent neural network，RNN）和多层感知机（multilayer rerception，MLP）等。其中，LSTM模型为解决传统深度学习模型受初始权值和值域影响较大，且难以准确捕获时序波动的问题而被提出。该模型在学习时间序列数据特征方面具有独特的优势^[4]。目前，利用InSAR数据和深度学习模型实现了地表形变趋势的预测，并不断改进模型结构以提升预测精度。但形变预测研究局限于模型算法本身精度的提升，忽略了地表形变的特征因子对形变趋势预测结果的影响，并且研究模型对时序信息的特征提取不足，因此模型的预测精度和泛化性仍有待提高。

本文采用干涉图叠加（Stacking-InSAR）和IPTA-InSAR技术对首都国际机场进行地表形变监测，获取2018—2022年InSAR形变监测数据，深入分析地表形变与地下水变化的相关性和滞后时间间隔。引入地下水位和降雨量为地表形变特征因子，结合排列组合方法和LSTM神经网络模型，提出C-LSTM地表形变趋势预测模型。将地下水位、降雨量及时序InSAR作为特征因子进行单因子和多因子地表形变趋势预测，获取2022年地表形变趋势预测数据，并将预测数据与InSAR形变数据对比分析，确保实验的准确性。

北京首都国际机场位于平原地区，地势相对较低，夹于温榆河与潮白河之间。主要含水层包括砂、卵和砾石，地下水资源主要来源于河流渗透、降水和水库水下渗。本文的研究区位于顺义区天竺镇（图1），该地区内第四系含水层自上而下分为潜水层和第一至第三承压水层。该地区处于温榆河和潮白河冲洪积扇下方，为北京地区地表形变显著发育的典型区域。研究区地下水资源总体呈超采状态^[5]，地下水位大幅度下降，地表形变范围不断扩大。因此，地表形变成为附近高频率的地质灾害之一^[6]。

本文采用的研究数据包括：Sentinel-1A卫星C波段SAR数据，采用干涉宽幅（IW）成像模式，VV极化方式，如表1所示；欧洲航天局（European Space Agency，ESA）提供的精密轨道星历数据；欧盟地球观测计划提供的30 m分辨率的数字高程模型（digital elevation model，DEM）；北京市水务局提供的地下水数据和降雨量数据（https://swj.beijing.gov.cn/swdt）。

本文研究过程包括3部分：首先，进行InSAR数据处理获取InSAR时序形变信息；其次，分析时序InSAR数据的地表形变过程及与地下水位的相关性；最后，结合时序InSAR、地下水位和降雨量，构建C-LSTM形变预测模型。基于模型预测进行单因子和多因子预测，并完成预测精度的全面评价。本文技术路线如图2所示。

Stacking-InSAR技术由Sandwell于1998年提出，本质上是干涉图的叠加^[7]。原理是将N景影像根据一个配对阈值，利用差分干涉测量技术（differential InSAR，D-InSAR）获取多幅干涉图，然后将干涉对进行相位上的加权平均^[8]，可有效降低大气对形变速率相关相位的影响，同时提升形变信号的信噪比，进而快速获取形变速率，其公式如下：

式中：Δt_j为第j对干涉图间隔的时间；Φ_j为第j对差分干涉相位。

IPTA-InSAR技术融合了最经典的2种时序InSAR技术的优势（PS-InSAR和SBAS-InSAR），不仅可以只对相干点目标进行差分干涉，又可依据小基线集原则来选择干涉对，保证数据的时间和空间密度。其原理为，根据影像特征完成目标点的选取，将目标点组合生成点对，提取相位信息，并用回归模型解算相位信息。

形变相位由线性变形和非线性变形2部分组成。2次获取的SAR影像生成的第i幅差分干涉图中任一像素（x，r），其干涉相位表示为^[9]：

式中：λ为中心波长；d（t_B，x，r）和d（t_A，x，r）分别为t_B与t_A相对于t₀的形变量；t_A和t_B时刻分别为从影像和主影像获取时间。

选择一个相干点（x_i，y_i）作为参考点，计算干涉点目标（x，y）的线性变形相位为：

式中：m为与干涉点目标相连的点数；γ_model（x，y；x_i，y_i）为相干函数模型；Δv（x，y；x_i，y_i）为（x，y）相对于参考点（x_i，y_i）的形变增量。

非线性形变相位存在于相位残差中，由非线性形变相位、大气延迟相位和噪声相位组成。最后，将非线性形变相位与线性相位结合，得到完整且准确的相干点目标的形变过程。

LSTM神经网络^[10]是传统RNN的一种流行变体，可以解决传统RNN无法传输长期信息造成梯度爆炸或消失的问题，因此被广泛应用于涉及序列建模的任务中^[11]。其中，地表形变时间序列分析作为LSTM应用的重要方向，目前已经取得了较大的进展^[6]。但为减少人为因素的影响，保证模型预测的客观性，本文在LSTM模型的基础上，使用排列组合的方法选取所有地下水水层和降雨量的特征组合，结合时序形变信息构建C-LSTM模型。基本结构包括特征组合、输入门、遗忘门和输出门（图3）。

C-LSTM模型通过挖掘所有特征组合预测信息，较为全面地分析了地下水、降雨和时序InSAR对于地表形变的影响程度，进而更加准确地预测地表形变趋势。模型所用方程如下：

式中：σ为sigmoid激活函数；f_t为遗忘门；h_t-1为t-1时刻的隐藏状态；x_t为t时刻的输入向量；W_f和b_f分别为遗忘门的权重和偏差；i_t为输入门；o_t为输出门；${\tilde{C}}_{t}$为候选单元状态；C_t为当前时间步的单元状态；h_t为新的输出值。

为验证时序InSAR地表形变监测方法的可靠性，本文对比分析了Stacking-InSAR和IPTA-InSAR结果，分别用disp_stacking和disp_IPTA表示。从2种方法的结果中选取了46 635个时间序列分析点进行相关性统计，得到决定系数R²为0.88。这一高相关性指数表明2种监测技术的结果具有很好的一致性，从而验证了时序InSAR地表形变监测实验二者结果的一致性，如图4所示。

IPTA-InSAR结果如图5所示，研究区最大形变速率达到-139.39 mm/a，整体呈下降的趋势。由于受到地裂缝的影响，导致机场西北部和东南部的地表形变表现出明显的差异，其中西北部沉降更为显著。2部分沉降之间存在清晰可见的分界线，该线自东北向西南贯穿首都国际机场。这一分界线的存在，进一步证实了地裂缝影响地表形变的空间分布趋势。

基于以上研究结果，选择研究区内4个地下水监测站：小天竺、王家场、米各庄-2和沙陀，分析地下水位变化与时序地表形变的关系（图6）。结果表明，除小天竺略有抬升外，其余监测站都呈现下降趋势，详细信息如表2所示。

通过深入分析长时序的累计形变量变化趋势，发现研究区地表形变有明显的周期性规律。具体来看，4个监测站点的累计形变量均在每年的6—7月达到阶段性的极值，说明地表形变趋势与季节性因素存在显著的相关性。

北京地区的地表形变受到多种因素的综合影响，包括地下水的过度开采^[12]、城市建设活动以及复杂的地质构造，因此地表形变的成因呈现出复杂多样性。根据已有的研究^[13]，北京地区的地下水位下降形成的凹陷锥与地表沉降漏斗区在幅度和位置上具有较高的一致性，从而推断长期的地下水高强度开采是造成该地区地表形变的主要原因。因此，本文选取2022年的地表形变数据，结合同期的地下水监测数据，进行相关性分析。通过线性插值法和三次样条插值法，确保时序InSAR和地下水位获得同一时间轴。如图7所示，二者在时间序列上的演变趋势是一致的。当地下水位下降时，地表沉降量也随之下降，并在6—7月达到全年最低点。随后，随着地下水位的逐渐回升，地表也逐渐抬升。

小天竺和王家场地下水位在第156天达到最低点，分别为-21.44 m和-25.47 m，累计形变量在第184天达最低点，分别为-2.70 mm和-7.70 mm，两者相差28 d。米各庄-2地下水位最低点出现在第156天，为-44.22 m，累计形变量最低点则在第177天，为-6.61 mm，两者相差21天。沙陀地下水位最低点在第163天，为-17.43 m，累计形变量最低点在第184天，为-9.20 mm，相差21天。实验表明，地下水位的最低点与累计沉降量的最低点之间存在一定的时间滞后。

小天竺形变趋势拟合程度最高，其地下水对该监测站的影响程度显著大于其他3个监测站。值得注意的是，随着地下水位上升，该站点处的地表也呈现为轻微抬升趋势。与此同时，其他3个监测站也显示出与地下水位变化的趋势一致性，表明地下水位变化显著影响着研究区的地表形变过程。

计算地下水位和地表形变量的相关系数r、线性拟合函数斜率、R²和均方根误差（root mean square error，RMSE）（表3），并绘制散点图（图8）。小天竺的相关系数高达0.76，R²为0.54，RMSE为0.18，这表明地下水位和累计形变量具有一定的相关性，拟合程度中等偏上。相比之下，王家场、米各庄-2和沙陀的相关系数大约为0.5，R²在0.2左右，这表明这些监测点地下水位和地表形变的相关性较弱，难以用线性回归模型解释它们之间的相关性。

地表形变是一个复杂的现象，受到多因素影响，包括地下水位的变化、地质构造特性和可压缩黏性土的厚度等。研究区域土壤深度大约为65 m，主要由粉土、黏土和砂土构成。其中，黏性土层的渗透性较弱，因此在水位变化时，水分的释放速度较慢，从而表现出滞后性^[14]。

通过相关性分析可知，地下水位与地表累计形变量随时间变化虽总体趋势一致，但两者之间存在一定的滞后量。特别在降雨量集中的6—8月，滞后现象较为明显，表明地下水位的变化对地表形变的影响并非即时发生，而是存在一定的延迟现象。

为分析这种延迟现象，本文针对季节性变化显著的6—8月期间，对地下水位和地表累计形变量进行了相关性分析并计算两者间的滞后量，结果如图9所示。4个站点的相关性峰值分别为0.86，0.55，0.84和0.74，通过峰值推算得到地表形变相对于地下水位变化滞后分别约为6 d，5 d，12 d和14 d。以上结果表明，地表形变与地下水位变化之间不仅存在相关关系，而且这种关系在不同监测站的时间滞后量不同。

本实验采用了LSTM模型，对InSAR数据进行拟合，并对结果进行精度验证和可视化。LSTM模型输入层具备接收多维特征信息的能力，而隐藏层由3层LSTM层组成，并设置了2个Dropout层防止模型过拟合。输出层则为1个全连接层，专门用于输出一维的InSAR数据预测结果。该实验的超参数设置如表4所示。模型的评估指标综合了5个关键参数，包括RMSE、平均绝对误差（mean absolute error，MAE）、平均绝对百分误差（mean absolute percentage error，MAPE）、R²和相关系数r，以全面衡量模型的预测精度。实验选取了6个位于研究区内的站点进行模型预测。6个站点分别为顺义沟北村（D1）、顺义柳衍村（D2）、佟辛庄（D3）、马圈村（D4）、高丽营（D5）和张镇赵各庄（D6），如图10所示。

1）单因子特征预测

利用LSTM模型进行单因子特征地表形变趋势预测，InSAR形变数据作为单一输入特征，对6个站点进行实验，详细预测结果如表5所示。所有站点的r均超过0.97，R²均高于0.88，这表明模型具有较高拟合度。特别是站点D6的R²高达0.973 9，表明预测结果准确程度较高。

然而，单因子特征对于某些站点预测结果的精度评价结果较差。D4的R²仅为0.887 4，这表明在使用单一特征输入的情况下，模型整体预测能力还有待提升。此外，D4，D1及D3等的RMSE，MAE和MAPE结果也较差。

综上所述，仅使用时序InSAR作为单一输入特征在一定程度上能够对某些站点形变趋势进行准确预测，但模型的拟合精度和预测能力仍存在明显的提升空间。这说明单一特征并不足以全面捕捉地表形变所有的复杂影响因素，需要考虑更多相关特征，以提高预测准确性和可靠性。

2）多因子特征预测

为进一步提升预测模型的准确性，本文又对6个监测站点进行了多因子特征的预测实验，详细预测结果如表6所示。选取的数据包括时序InSAR、潜水层水位、第一至第三承压水层水位及降雨量。每个站点以InSAR作为固定输入数据，将其余5组特征数据进行排列组合，构建C-LSTM形变趋势预测模型。采用组合的方法旨在通过整合不同类型的影响因子，以期获得更为精确和全面的预测结果。

在站点D1的实验中，以InSAR为固定输入值，结合潜水层水位及第一和第二承压水层水位作为多因子特征组合构建模型，预测结果R²达到0.984 3为最佳，相较于单因子预测R²提高2.45%。对于站点D2和D3的实验，仅将InSAR和潜水层水位作为特征组合构建的模型预测效果最佳，R²分别达到0.981 4和0.981 7，相较于单因子R²分别提高1.52%和4.16%。在站点D4和D6实验中，以InSAR和降雨量组成的特征组合取得了最优预测效果，R²分别为0.959 1和0.988 0，比单因子提升了8.08%和1.45%。对于站点D5的实验，采用InSAR、潜水层水位和降水量作为特征组合预测结果，R²达到0.988 2为最佳，比单因子提高5.08%。

多因子特征模型的预测准确度相较于单因子特征模型有显著的提升。实验分析得到，站点D1，D2，D3和D5的InSAR形变值和潜水层水位间存在明显的相关性。其中，D1站点的InSAR形变数据不仅与潜水层水位相关，还与第一和第二承压水层水位存在较高的相关性。进一步说明，以上4个站点（D1，D2，D3和D5）的地表形变在很大程度上受到浅部含水层水位变化的影响。而站点D4和D6的InSAR形变值与累计降雨量之间的相关性更为显著，这表明这2个站点的地表形变与降雨量的关系更为紧密，降雨量的增减可能对地表形变有更直接的影响。

图11左列展示了单因子特征模型的实验结果图，图11右列呈现了多因子特征模型的实验结果图，具体为模型预测结果（橙色曲线）与实际观测结果（蓝色曲线）的对比情况。从图中可以明显看出，无论是单因子特征模型还是多因子特征组合模型，其预测结果都与实际形变具有较高的吻合度。D1至D6预测结果的R²分别达到了0.984 3，0.981 4，0.981 7，0.959 1，0.988 2和0.988 0。即在C-LSTM模型中引入地下水层特征后，预测结果与实际观测数据（橙色曲线与蓝色曲线）的吻合程度得到了一定的提升，其提升范围为1.45%至8.08%。这表明，地下水层特征的加入能够有效增强模型的预测精度。

本文采用2018—2022年InSAR时序形变信息，探讨了地下水和地表形变之间的相关性和滞后性，并通过单因子和多因子特征的C-LSTM模型对首都机场地表形变趋势进行了分析及预测。经过实验分析，得出以下结论：

1）本文通过IPTA-InSAR技术得到平均形变速率最大值-139.39 mm/a。除小天竺观测到小幅度抬升外，在王家场、米各庄-2及沙陀监测得到的累计形变量分别为-14.38 mm，-23.46 mm及-27.12 mm；平均形变速率分别为-1.59 mm/a、-3.66 mm/a及-3.69 mm/a。地表形变趋势显示出以年为单位的周期性变化，特别是在每年的6—7月，累计沉降量达到极小值，这表明在这一时期地表沉降最为显著。

2）地表形变与地下水位的相关性分析显示，在小天竺，两者的相关系数达到0.76，R²为0.54。而在王家场、米各庄-2和沙陀，相关系数均在0.5左右，R²在0.2左右，说明其相关性较差。尽管2个变量的整体趋势是一致的，但它们之间存在明显的滞后性。进一步相关性标定，得到4个监测站分别存在6 d，5 d，12 d和14 d的滞后量。研究表明，地表沉降与地下水位变化之间不仅存在相关性，而且这种关系伴随着滞后效应，且滞后时间在不同站点存在差异。

3）本研究通过构建C-LSTM模型，将地下水位、降雨量和时序InSAR形变结果作为特征因子，以提高地表形变预测的准确性。分别对单因子和多因子模型进行精度评价，发现在最优特征组合下，多因子预测R²相比于单因子预测分别提升了2.45%，1.52%，4.16%，8.08%，5.08%和1.45%。研究表明，本文提出的多因子预测在提升地表形变预测精度方面具有显著效果，也说明地下水层特征信息的加入能够有效增强地表形变预测模型的精度。