陕西安康机场始建于2016年，竣工于2020-05，正式通航于2020-09-25。机场位于以膨胀土为主要土体的安康月河盆地^[1]。膨胀土是自然形成的一种多裂隙地质体，具有显著胀缩性、超固结性和多裂隙性等不良性质^[2]。此外，安康机场场址地形复杂，机场建设中填土3 000万方，挖填垂直高差94 m，目前是国内最大体量的膨胀土高填、挖方机场^[3]，大规模填方土的固结压缩性也是地面沉降的常见诱因^[4]。因此，安康机场区域存在填方土固结压缩、膨胀土潜在形变等灾害隐患，时刻威胁着安康机场安全通航，十分有必要对安康机场实时进行监测和稳定性分析，获取机场正式运营前期形变监测及稳定性分析评价则是保证机场安全通航的关键组成部分。

目前针对机场区域地表形变问题的常规监测方法有全球定位系统(global positioning system，GPS)测量、水准测量、分层标定技术等。但上述监测方法人工设备成本较高，只能获取较为离散的点位形变数据，可能遗漏形变信息，较难对研究区域进行整体形变监测。合成孔径雷达干涉测量(interferometric synthetic aperture radar，InSAR)技术作为一种新型对地观测手段，可以克服传统观测手段的不足，通过强度和相位两种观测值，获取高分辨率地表形变信息，实现空间大范围探测和长周期精密监测^[5]。目前，时序InSAR技术已广泛应用于以缓慢形变为主的各类机场进行长时序形变监测与机理分析，如以填海土为沉降主因的中国上海浦东机场^[6]、厦门翔安机场^[7-8]和香港赤腊角机场^[9]，沉降与高填方强相关的康定机场^[10]和九寨黄龙机场^[11]，以及热带泥炭特殊性质导致沉降的马来西亚吉隆坡国际机场^[12]等机场区域。近年来，通过分析InSAR观测结果与土壤环境之间的规律性，证明了时序InSAR技术用于监测与土壤性质相关的地表形变的可靠性^[13]，时序InSAR技术的特点同样被证明契合于监测膨胀土特殊水力学性质导致的周期性缓慢形变问题^[14]。

本文采用一种基于高程改正的InSAR形变监测策略，结合永久散射体干涉测量(persistent scatterer InSAR，PS-InSAR)技术与小基线集干涉测量(small baseline subset InSAR，SBAS-InSAR)技术对膨胀土机场区域进行运营前稳定性分析，并结合膨胀土胀缩特性、填方层厚度以及区域降水量分布对膨胀土形变规律进行分析。

1 研究区概况与数据来源

1.1   研究区域

安康机场位于秦岭南麓月河盆地北缘台地区域，存在大量膨胀土，成分包括粘性土、粉质粘土、碎石类土、砂岩、泥岩与砾岩等。膨胀土特殊水力作用是膨胀土土体失稳的根本原因，主要包括水分和应力变化造成胀缩变形、裂隙发育和强度降低的作用^[15]。这些特性极易对土体上的基础设施造成潜在的地质危害。

综合研究区域水文条件与地质环境，首先，安康机场兼备高强降雨期集中与膨胀土大量分布等膨胀土灾害发生必要条件。强降雨期，膨胀土土体吸水，发生膨胀现象；旱期，膨胀土土体失水，发生收缩现象，雨期与旱期反复胀缩，造成安康机场区域表层土体裂隙发育。其次，在罗家河、俭沟等复杂地貌条件下进行高填方工程会导致两类潜在地质病害，即原有地基固结与构筑体自身压缩导致的高填方地基沉降变形，以及填方体施工后变形和施工进度差异造成的不均匀沉降^[16]。因此，有必要结合降水量分布对安康膨胀土机场地面形变进行持续监测。

1.2   实验数据来源

哨兵1号A星(Sentinel-1A)提供覆盖安康地区、重访周期为12 d的卫星数据，满足对安康机场区域进行时序InSAR监测的数据条件。实验数据如图 1所示，通过对图1(b)、1(c)中不同时期高分系列光学影像目视解译发现，直至2020-04，机场区域仍有较为明显的施工痕迹；2020-05后，机场地形基本固定，跑道施工已经完成，满足InSAR监测的必要条件，而后续合成孔径雷达(synthetic aperture radar，SAR)影像中干涉对的相干性也证明了这一点。因此，本文选择2020-05-2021-01共17景Level-1 Sentinel-1A升轨影像单视复数(single look complex，SLC)数据进行实验，研究区域及数据覆盖范围如图 1(a)所示。

图  1  实验数据

Figure  1.  Study Data

下载: 全尺寸图片 幻灯片

使用Sentinel-1A卫星SAR数据成像模式为干涉宽幅(interferometric wide swath，IW)，极化方式为VV与VH双极化，入射角为39.5°，方位向×距离向分辨率为2.33 m×13.99 m，影像数据如下：轨道号100，轨道配置ASCENDING，获取日期分别为2020年5月24日，6月5、17、29日，7月11、23日，8月4、16、28日，9月21日，10月3、15、27日，11月8、20日，12月2、26日。

在InSAR数据处理中，用SRTMGL1^[17](shuttle radar topography mission global 1 arc second)扣除地形相位。SRTM计划提供的分辨率为1″的SRTMGL1产品相较于3″分辨率的数字高程模型(digital elevation model，DEM)数据细节信息更为丰富，更符合真实地貌；相较于其他分辨率为1″的DEM产品，如ASTER GDEM(advanced spaceborne thermal emission and reflection radiometer global digital ele‑vation model)颗粒感更弱，丢失信息更少。因此，在DEM偏差估计与改正中，采用NGA自2015年7月提供的分辨率为1″(约30 m) 的SRTMGL1作为外部DEM产品。在有关DEM的相关研究中，SRTMGL1被验证在平地区域垂直精度小于9 m^[18]，在多山丘陵区域，垂直精度为4.31~14.09 m^[19]。

2 数据处理与方法

2.1   InSAR数据处理策略

针对安康膨胀土机场施工历史中挖填方程度较大，而外部DEM数据采用的SRTMGL1时效性较差的问题，在InSAR数据处理中需要重点考虑DEM误差的改正。对此，本文利用GAMMA与StaMPS软件，同时采用基于高质量干涉对并顾及DEM误差的SBAS-InSAR与PS-InSAR进行形变监测。整体数据处理流程如图 2所示，具体步骤为：(1)利用精密轨道数据以及SRTMGL1作为外部DEM产品对覆盖研究区域的Sentinel-1A数据进行预处理，包括影像裁剪、多普勒效应的消除以及影像配准；(2)选取短时间基线、长空间基线的干涉对子集来计算DEM误差，将得到的DEM误差迭代更新到原始DEM数据中^[20]；(3)选取短时间基线、短空间基线的干涉对子集，使用SBAS-InSAR技术计算研究区域的形变时间序列；(4)使用PS-InSAR技术计算研究区域累计形变与形变时间序列。

图  2  数据处理流程

Figure  2.  Flowchart of Data Processing

下载: 全尺寸图片 幻灯片

2.2   数据预处理

IW成像模式下的Sentinel-1A影像数据以循序扫描地形观测(terrain observation by progressive scans，TOPS)成像实现250 km幅宽的扫描成像，但在TOPS模式下，不同burst之间存在较大多普勒中心频率变化，配准误差将会引起干涉相位的偏移^[21]，这导致不同burst之间容易出现明显的相位跳变现象，因此需要利用GAMMA软件进行数据预处理。由于安康机场位于两景Sentinel-1A数据接缝处，因此先将对应时间两景影像拼接并进行burst选取。所有影像需添加精密轨道参数，并精确配准到某一参考影像雷达坐标下。为实现SLC影像间的精确配准，利用主副影像参数文件和DEM计算初始查找表，将辅影像重采样到参考影像的几何框架下。结合距离向和方位向配准偏移多项式对影像进行重采样。在此基础上，有少量数据需进一步基于频谱多样性进行重叠区配准，消除差分干涉图中burst间的相位不连续造成的误差。

2.3   DEM误差改正

在InSAR数据处理中，引入了1″(约30 m)分辨率的SRTMGL1作为外部DEM剔除地形相位。考虑研究区域进行了高挖、填方工程，以及外部DEM产品相对于Sentinel-1A影像数据时效性较低，这会导致差分干涉图中出现与垂直基线相关的DEM误差相位^[22]，此外地形坡度、地表植被覆盖等因素也会影响SRTMGL1的精度质量^[23]，因此本文对DEM偏差值进行了估计并改正。

时序InSAR技术从覆盖同一地区多幅SAR影像中提取地形相位，再根据高程与相位间的关系解算出DEM误差。基于文献[24]提出的SBAS算法对DEM误差进行估计并校正。假设地表形变为完全线性，则图上像素点的线性形变及其DEM误差$∆z$的求解函数模型如下：

式中，$\mathit{A}$为系数矩阵，反映了待求参数向量$\mathit{X}$与观测值$\mathrm{\delta }\mathit{\phi }$间的关系；$\mathit{X}$由平均速率$\overline{v}$、平均加速度$\overline{a}$、平均加速度的增量$∆\overline{a}$以及待求DEM误差$∆z$共同组成。$\mathit{A}$的计算式为：

其中$,{\mathit{B}}_{M\times N}{\mathit{M}}_{N\times \mathrm{3}}$为$\mathit{X}$的形变部分对应的系数矩阵，反映形变速率与差分干涉图解缠相位之间的关系；${\mathit{C}}_{M\times \mathrm{1}}$为$\mathit{X}$的DEM误差$∆z$对应的系数矩阵，反映了DEM误差与差分干涉图解缠相位之间的关系，表达式为：

式中，$\lambda$表示雷达波长；$B_{\perp \mathrm{1}}$表示垂直基线长度；$r$表示轨道半径；$\theta$表示卫星入射角。将式(3)与式(1)联立，基于奇异值分解，即可求取最小二乘准则下的DEM误差值。式(3)中反映的DEM误差与干涉对垂直基线之间的关系说明，干涉对垂直基线越长，差分干涉解缠相位对DEM误差越敏感。此外，大气伪影和相位噪声的影响在时间上是不相关的，因此具有长垂直基线子集的多干涉图的加权平均可以大大减少最终DEM上的相位伪影。综上，为更加准确地解算DEM误差，选取兼顾高质量与长垂直基线的干涉对子集。

2.4   时序形变监测

相比于差分InSAR方法，为了更好地降低空间失相干影响，SBAS-InSAR^[24]依据短基线原则选取影像对。求解形变时，利用奇异值分解方法(singular value decomposition，SVD)将由于基线限制而分开的独立数据集联系起来，解出基于最小二乘准则的地表形变解^[25]。SBAS-InSAR基线图如图 3(a)所示。

图  3  SBAS-InSAR和PS-InSAR基线图

Figure  3.  SBAS-InSAR and PS-InSAR Baseline Maps

下载: 全尺寸图片 幻灯片

针对小区域、高填方、膨胀土机场的特殊性质，通过迭代更新外部DEM改正DEM误差后，采用SBAS-InSAR技术进行地表形变监测，主要处理步骤如下：(1)生成干涉图，本文选取空间基线阈值为±100 m，时间基线阈值为36 d，共生成56对干涉对；(2)计算相干性，相干性是衡量InSAR干涉图质量的一个重要指标，通过对相干性的计算，不依赖外部测量数据^[26]，结合地质踏勘资料以及DEM误差分布，得到较为可靠的安康机场InSAR形变监测结果；(3)选取地面相干性稳定且良好的解缠点，采用基于最小费用流的相位解缠；(4)去除趋势向与大气误差；(5)通过SVD方法解算形变时序相位的加权最小二乘解。

利用PS-InSAR^[26]技术获取地表形变信息的数据处理流程可分为永久散射体(persistent scatterer，PS)点选取、相位解缠、时空滤波和形变量计算4个部分。基于StaMPS选取长时间间隔内保持良好稳定性的PS点，所选点受外界噪声影响较小，所以具有较高相干性。这样既可以避免不稳定地面点参与数据解算，又削弱了时间、空间失相干和大气延迟对形变监测结果的影响。为生成具有较短时空基线的干涉对，选取2020-09-21的影像为主影像，时空基线图图 3(b)所示。PS-InSAR技术适用于研究长时间缓慢的地表形变变化规律，同时其监测精度更高，可用于与SBAS-InSAR获取的研究区域形变时间序列进行精度检核。

3 安康膨胀土机场算例分析

3.1   DEM误差分析

将预处理后17幅影像数据都进行差分干涉，经过干涉对质量挑选后，共生成56对干涉对。由§2.3中所述DEM误差与垂直基线的关系，为了估计更准确的DEM误差相位，以100 m为垂直基线阈值，将干涉对分为两个子集：(1)长垂直基线(> 100 m)和短时间基线的DEM子集；(2)短垂直基线(< 100 m)和短时间基线的形变子集，子集干涉对组合见图 3(a)。得到了13对干涉对计算DEM误差和42对干涉对计算地表形变。

利用SBAS-InSAR技术，将用于DEM误差的干涉对子集，经过滤波、解缠去除干扰相位等处理步骤后，基于最小二乘准则，得到DEM误差并将其对应至机场挖填方分布，如图 4所示。

图  4  DEM误差分布图

Figure  4.  DEM Error Distribution

下载: 全尺寸图片 幻灯片

结合图 4中阴影区域所示施工挖填方资料分析，SBAS-InSAR方法获取的DEM误差分布与施工期挖填方资料基本符合，这从一定程度上说明了本文InSAR数据处理结果的可靠性。图 4中安康机场由西向东依次分布两个较大的挖方区(DEM误差最大值达到了52 m)和3个较明显的填方区，填方造成DEM误差最大值达到61 m。这类高挖填方区域存在的潜在形变问题需要在InSAR结果中重点关注。

综上，虽然使用SBAS-InSAR技术估计的DEM误差精度有一定局限，但其反映的地表变化的空间分布情况较为准确，可以作为一种辅助判别资料选取潜在形变区域，从而进行重点分析。

3.2   时序InSAR形变分析

使用小基线集InSAR技术进行DEM误差改正后，迭代更新DEM产品，分别采用SBAS-InSAR与PS-InSAR进行处理，计算地表形变时序分布情况，其累计形变情况如图 5所示。

图  5  时序InSAR累计形变情况

Figure  5.  Time-Series InSAR Cumulative Deformation

下载: 全尺寸图片 幻灯片

图 5(a)中结果表明，安康膨胀土机场竣工后7个月累计形变量小于26 mm，机场内部最大形变量为－15.9 mm，未出现较大形变，基本稳定。监测区域内较大形变均出现在机场周围膨胀土填方边坡区，分别为图 5(a)中P₁、P₂、P₃、P₄。其中最大形变量出现在机场北侧膨胀土边坡(图 5(d))上P₂点，为－21.7 mm。计算SBAS-InSAR结果标准差作为精度评定标准之一。标准差分布于0~－4 mm，最大值为－3.9 mm，说明SBAS-InSAR获取研究区域累计形变结果较为可靠。由图 5(c)可以看出，累计形变分布与SBAS-InSAR高度吻合，最大形变点行列号一致，均出现在机场北侧膨胀土边坡影响区，PS-InSAR解算出该点最大形变量为－19.4 mm，分别提取PS-InSAR与SBAS-InSAR累计形变结果中机场跑道两侧膨胀土边坡上形变最大两个点P₁、P₂的形变时间序列(见图 6)。

图  6  膨胀土边坡区特征点位时序形变图

Figure  6.  Time-Series Deformation of Point Position in Terminal Area

下载: 全尺寸图片 幻灯片

从图 6中可以看出，SBAS-InSAR与PS-InSAR获取同一像元点的形变时间序列分布基本吻合，将数据解算方法作为影响因子，对选取两个点时间序列分别进行方差分析，结果见表 1。通过一元方差分析可以看出，在99%的置信区间下，SBAS-InSAR与PS-InSAR监测结果不存在显著性差异，结合表 1与图 6说明了监测方法的可靠性与有效性。

表  1  方差分析表

Table  1.  Variance Analysis

统计项	P1点					P2点
	离差平方和	自由度	均方离差	F值		离差平方和	自由度	均方离差	F值
组间	12.644	1	12.644	0.68		7.314	1	7.314	0.2
组内	596.589	32				1 159.334	32	36.229

下载: 导出CSV 

| 显示表格

根据文献[27]，下面主要对机场膨胀土填方边坡、机场飞行区以及航站区稳定影响区进行稳定性分析。

1) 膨胀土边坡稳定影响区稳定性分析。填方膨胀土边坡稳定影响区分布于机场飞行区两侧，均属于高填方区域。施工过程中进行了混凝土浇筑以及建设防滑挡墙等防护性措施，但由于其地基与填方土均以膨胀土为主要成分，该区域可能出现与膨胀土水力作用相关的形变。图 6所选取的形变特征点均位于安康膨胀土机场膨胀土填方厚度最大的区域。两个点位形变趋势相似，最大累计形变达到－21.7 mm，但仍然处于民用机场岩土工程设计规范规定范围内。结合现场踏勘与丰水期形变趋势分析，特征点位于机场外侧的膨胀土边坡，顶部存在未进行混凝土浇筑防雨防水处理的裸露土壤区域。雨季出现的持续抬升形变趋势，判断为大量降雨后，水分下渗导致膨胀土发生胀缩性形变，最终造成膨胀土边坡稳定影响区发生抬升形变。

2) 飞行区与航站区稳定性分析。从图 6可以看出，安康机场飞行区基本稳定，累计形变量较小，区域形变分布在－15.9~3 mm范围内。根据§3.1中估算的DEM误差分布情况与挖填方资料，对安康机场飞行区的跑道区和机(站)坪区的挖填方施工区域，分别提取特征点进行时序形变分析，如图 7所示。

图  7  跑道及航站区点位时序形变图

Figure  7.  Deformation Diagrams of Point Position of Runway and Terminal Area

下载: 全尺寸图片 幻灯片

图 7(a)为提取的机场跑道挖填方厚度最大区域特征点形变时间序列，可以看出，跑道填方区在丰水期出现小范围内的地面回弹，最大形变值为3 mm，挖方区形变时间序列基本未出现明显形变，形变分布在－3 mm以内，均属于合理工后沉降范围内，跑道形变趋势平稳，整体稳定，足以满足通航等功能性需求。图 7(b)为机场机(站)坪内挖填方最大区域特征点形变时间序列，可以看出，机场机(站)坪区填方区与跑道填方区形变趋势相似，竣工后发生较小形变后于丰水期出现回弹现象，最大形变量达到3.6 mm，累计形变量较小，为2 mm，基本处于稳定状态；从机(站)坪内挖方区点位时序形变信息可以看出，该点发生持续形变，累计形变量小于－3 mm。机(站)坪区整体形变较小，趋于稳定。

根据文献[27]对飞行区规定使用年限内的形变标准(跑道：工后沉降2~3 cm；滑行道和机坪：工后沉降3~4 cm)，安康膨胀土机场飞行区累计形变远小于使用年限内的形变标准的该时间跨度的允许平均形变值，且考虑到工后沉降随时间逐渐减小，机场飞行区稳定性将进一步增强，可知安康膨胀土机场运营前夕飞行区较为稳定。

3.3   膨胀土填方区形变规律分析

一般情况下，高填方建筑体发生形变的主要原因有填方土下地基的固结沉降、填方部分自身压缩导致的工后沉降以及施工不连续导致的不均匀沉降。此外，在膨胀土区域需要考虑干湿循环对压实膨胀土微观结构的影响造成的地面形变。综合考虑上述潜在影响因素，分别结合填方层深度及降水量对形变进行分析。安康膨胀土机场最大填方厚度超过了56 m，填土超过3 000万方，结合§3.1中DEM误差估计结果表明造成DEM误差达到了64 m。为分析安康膨胀土机场填方体与地表形变的关系，结合DEM变化，提取研究区域监测点的累计形变及其所处位置的填方体深度进行分析(见图 8)。

图  8  累计形变与填方体厚度关系

Figure  8.  Relationship Between Cumulative Deformation and Filling Thickness

下载: 全尺寸图片 幻灯片

从图 8中可以看出，随着填方体厚度的增加，沉降观测点的数量及形变程度明显增加，填方土的分布及其厚度明显影响着地表形变的分布和大小。形变最明显的区域均分布在填方层深度 > 30 m的区域，该区域的形变监测点占全部形变沉降监测点的70%，尤其只有在填方层深 > 20 m的区域，才出现累计形变 > 20 mm的监测点，相反，填方层深 < 15 m的区域，沉降点较少，累计形变值较小。由于填方土不同深度的压缩变形以及对下方原状土的压力随填方层深增加而增大^[28]，因此，安康机场填方土的压缩变形是地表形变的诱因之一。

安康地区处于亚热带季风气候，降水量资料表明，机场所处月河流域降雨较多且集中分布在6月下旬至8月下旬。对研究区域时序形变的分析中发现，整个机场区域填方区形变趋势较为相似，具有一定规律性，即形变速率均在雨季逐渐减缓，发生回弹。为分析膨胀土填方区形变规律与降水量分布之间的关系，在监测区域中出现较大形变的膨胀土填方边坡上(图 5(a)中P₁~P₄)提取特征点时序形变分布与降水量分布(见图 9)。

图  9  填方区点位时序形变图

Figure  9.  Time-Series Deformation Diagram of Fill Area Points

下载: 全尺寸图片 幻灯片

从图 9中可以看出，安康机场膨胀土填方区形变时间序列分布普遍在雨季出现减缓，甚至在7月降雨最丰沛时发生回弹。结合膨胀土填方区特征点时序形变分布与降水量分析，推测为安康6月-8月，降雨量增加，安康机场膨胀土吸水膨胀的水力学性质抵偿了机场地表由于工后形变、填方土固结压缩等因素引起的微小形变。这种现象与膨胀土和土壤湿度相关的胀缩性形变特征相吻合，因此分析其特殊形变分布是由于填方土工后沉降和膨胀土胀缩性的联合作用导致。

4 结语

本文结合估计高程改正的SBAS-InSAR与PS-InSAR技术，基于Sentinel-1A与高分系列影像数据，通过适宜的数据处理策略进行DEM误差校正并获取区域时序形变分布，结合区域时序形变分布与降水量数据进行分析。结果表明：

1) 利用SBAS-InSAR技术估算的DEM误差结果较为可靠，可以在一定程度上作为筛查潜在形变区域的辅助性资料。

2) 在SBAS-InSAR数据处理过程中，采用构建不同用途的干涉对子集的干涉对组合策略，可以在进行高挖填方区域形变监测时校正DEM误差造成的相位噪声。

3) SBAS-InSAR监测结果中，监测点密度明显多于PS-InSAR监测结果，两种方法监测精度相当，SBAS-InSAR监测结果能够更完整地显示机场形变信息的优势。

4) 安康膨胀土机场地表形变是由于填方土以及膨胀土湿胀干缩特性的联合作用导致。

5) 安康膨胀土机场跑道以及航站区地面形变量较小，整体稳定满足设计要求，可以满足机场正式运营需要。

总体上，基于高程改正的时序InSAR监测技术对安康机场区域的监测结果表明，机场区域形变规律与膨胀土胀缩特性高度吻合。机场虽然暂时未发生显著的地表形变，但其与降水量相关的时序形变变化趋势说明，机场开始运营后，对安康膨胀土机场进行长期形变监测是十分必要的，相关研究需要进一步结合更长时序的多源遥感数据与实时降水资料，进行更长时间的形变监测与机理分析。