延安新区位于湿陷性黄土沟壑地带，土地结构疏松，极易发生阻碍工程进展的地质灾害^[1]，因此有必要对该地区地表形变时空变化进行分析。目前，有学者对延安新区沉降机理进行研究^[2]，但是对沉降与抬升区域的时空分析还不够深入。本文使用小基线集合成孔径雷达干涉测量SBAS-InSAR技术^[3]，利用覆盖延安新区2016-07~2020-11的43景Sentinel-1A升轨雷达影像数据，获取延安新区的地表形变信息，并运用经验正交函数(empirical orthogonal function，EOF)对SBAS点进行分解，挖掘延安新区地面变化的主要时空演化特征，为该地区的城市建设、灾害防控等提供监测依据。

1 研究区概况及数据来源

延安新区(北区)施工时间为2013-04~2018-08，最大挖方厚度与最大填方厚度分别为118 m和112 m，总施工面积22.3 km²。降雨是延安新区地下水的主要补给方式^[4]。本文使用Sentinel-1A卫星2016-07-03~2020-11-27共43幅升轨影像，数据类型为单视复数，模式为干涉宽幅，极化方式为VV极化。利用每幅影像对应的精密轨道数据校正卫星轨道，使用30 m分辨率的SRTM DEM作为地形数据反演形变速率。

2 研究方法与数据处理 2.1 SBAS-InSAR技术

SBAS-InSAR技术适用于分布式目标的时序分析。相比于传统的D-InSAR技术，SBAS-InSAR技术能够在一定程度上克服空间失相干和时间失相干，在最小形变速率标准的基础上对数个小基线数据集进行组合，得到所需的干涉图。相比于PS-InSAR技术，SBAS-InSAR技术能够同时监测出线性和非线性形变，具有更强的适用性。

2.2 经验正交函数分析

EOF分析方法的分析对象为场的时间序列，该方法能够在先验条件未知的前提下获取研究区域时空数据中的空间分布特征和时间序列。随时间变化的变量场可以被分解为只随时间变化的时间函数以及不随时间变化的空间函数，可以在没有事先人为规定的前提下，由变量场序列本身特征来确定场的结构，并且从分解出来的结构中解读出物理意义^[5-6]。

2.3 数据处理

首先进行SBAS-InSAR数据处理，将2016-10-19的影像设置为超级主影像，设置最大空间基线为临界基线的37%，最大时间基线为175 d，生成的时空基线如图 1所示。对所有干涉像对进行干涉处理前首先需要进行多视处理，根据以往经验，设置多视视数为4 ∶1。然后使用Goldstein方法对干涉结果进行滤波处理，使用Delaunay最小费用流法(minimum cost flow，MCF)进行3D相位解缠。在选择地面控制点时，应该避开形变区域，使用奇异值分解法(singular value decomposition，SVD)反演获得研究区域可靠的形变序列。最后，将SBAS-InSAR处理后得到的栅格结果转为矢量格式，方便接下来进行EOF分析。

将获得的矢量结果构建成i个观测点、j个观测时间的矩阵X：

$ \begin{array}{l} \mathit{\boldsymbol{X}} = \left[ {\begin{array}{*{20}{c}} {{x_{11}}}&{{x_{12}}}& \cdots &{{x_{ij}}}\\ {{x_{21}}}&{{x_{22}}}& \cdots &{{x_{2j}}}\\ \cdots & \cdots & \cdots & \cdots \\ {{x_{i1}}}&{{x_{i2}}}& \cdots &{{x_{ij}}} \end{array}} \right]\\ i = 1, 2, \cdots , m, j = 1, 2, \cdots , n \end{array} $

(1)

方程求解前，需要先对矩阵内的数值进行正则化，将原始数据转化成无量纲的形式。计算得到矩阵X的协方差矩阵A为：

$ \mathit{\boldsymbol{A}} = \frac{1}{N}\mathit{\boldsymbol{X}}{\mathit{\boldsymbol{X}}^{\rm{T}}} $

(2)

式中，N为SAR影像的观测期数。根据实对称矩阵分解定理得：

$ \mathit{\boldsymbol{A}} = \mathit{\boldsymbol{V \boldsymbol{\varLambda} }}{\mathit{\boldsymbol{V}}^{\rm{T}}} $

(3)

式中，对角矩阵Λ(λ₁, λ₂, …, λ_i)为矩阵X的特征值，V为特征向量，每个非零特征根对应的特征向量值为EOF模态。时间系数T可表示为：

$ \mathit{\boldsymbol{T}} = {\mathit{\boldsymbol{V}}^{\rm{T}}}\mathit{\boldsymbol{X}} $

(4)

特征根越大表示模态越重要，方差贡献率也越高。每个模态的贡献率P可表示为：

$ P = {\lambda _i}/\sum\limits_{i = 1}^m {{\lambda _j}} $

(5)

一般情况下，通过分析前几个占比较高的EOF模态得分和对应时间系数，就可以解释原始数据的主要时空变化规律。

3 结果与讨论 3.1 SBAS-InSAR形变结果

从Sentinel-1A雷达数据集中可以得到延安新区视线向年平均地面形变速率。已知雷达波入射角，就可以把视线向的形变转换为垂直向的形变(图 2)，其形变结果和空间分布与前人所得结果基本一致^[7]。由图可见，研究区域整体形变速率为-56~32 mm/a，绝大部分区域形变比较稳定，形变量为-5~5 mm/a。西北至东南方向有一条带状沉降区域，东南方向沉降最严重，沉降速率达-56~-45 mm/a。东北部可以看出有明显沉降，沉降速率为-45~-20 mm/a。该区域同时存在明显的地表抬升，抬升速率为15~32 mm/a，但由于该区域失相干严重，SBAS点并没有将其完全覆盖。

3.2 EOF分析地面沉降时空分布特征

对2016-07~2020-11共43景Sentinel-1A影像得到的SBAS点进行EOF分析，结果如表 1所示。由表可见，第1模态的方差贡献率高达85.24%，第2模态为12.52%，第3模态为0.71%。前2个模态的累积方差贡献率为97.76%，超过总体的95%，因此前2个模态就可以解释延安新区的地表形变时空演化特征。

3.2.1 第1模态特征

第1模态的方差贡献率为85.24%，特征向量值有正有负，正值和负值覆盖区域与年平均形变速率图大致相同。作为占比最高的一个模态，第1模态反映的是延安新区(北区)地表形变的整体空间分布特征。使用反距离权重法(inverse distance weighted，IDW)对SBAS点和第1模态的得分进行插值，结果如图 3所示。由图可见，两者的空间分布高度相似，相关性为90.43%，因此能准确反映出延安新区的沉降和抬升区域。

由图 3(a)可见，地表形变区域与填挖方区域有着高度一致的对应关系。为进一步证实这一结论，选择覆盖延安新区(北区)的2018-08-22~2020-11-27共22景Sentinel-1A数据进行SBAS-InSAR分析，得到的年平均形变速率如图 4所示。选择2个横穿填挖方区域的剖面CC′和DD′，对比剖面上各点的形变速率和2012年施工前的高程数据(图 5)。由图可见，CC′剖面点号0~30、56~93、120~130，DD′剖面点号48~73、88~103为填方区域，对应的形变速率均为负值。根据以往的研究可知，造成地表沉降的主要原因有3个：1)填方体在填充后自身发生压缩下沉；2)原地基在填方后被施加大量载荷，破坏其内在应力结构，发生压缩下沉；3)原地基和填方体的土质均属于湿陷性黄土^[8]，容易受到土体湿度增加的影响从而发生压缩下沉。研究资料显示，延安新区土体的含水率与填方深度成正比^[9]，CC′与DD′剖面各填方区域中填方厚度随高程的降低而增加，沉降速率的绝对值相比其他剖面点较大。CC′剖面点号36~50、100~114、132~142，DD′剖面点号25~40、114~123为挖方区域，对应的形变速率为正值。在挖山造地的过程中，地表荷载量的减少打破原有的应力平衡，为达到新的平衡，挖方区域的土体会出现一定程度的回弹^[10]。由图 5可见，挖方区域中高程值越大的点号挖方越深，卸荷量越大，抬升速率相对较高。在挖方厚度较小的区域土体回弹量小，且后期施工建设必然会对该区域重新加载，所以没有出现明显的抬升。

由图 6可见，第1模态的时间系数均为负值。2016-07~2017-01时间系数绝对值快速上升时段正是延安新区快速建设时期。2017-01~2020-11时间系数绝对值维持在0.9~1的范围内，绝对值在2018-10达到最大，说明此时空间分布最为典型，而且这个时间点正是延安新区(北区)建设完成的节点，随后整体的地表形变逐渐趋于稳定状态。通过分析延安新区抬升与沉降的形成机理和发生变化的时间节点可知，挖山造地正是影响延安新区(北区)地表形变的最主要原因。

3.2.2 第2模态特征

第2模态的方差贡献率为12.52%，特征向量值有正有负。由图 7(b)可见，第2模态得分为正的区域分布在东北方向，得分为负的区域分布在西南方向，二者形成较为明显的分界线。B区域正值较大，这与B区域较晚的开发建设时间有一定的关系。据资料显示，特征点P₁、P₂、P₃所在区域的工程任务完成时间为2014年初到2016年底，而特征点P₄所在的区域直到2018-08才彻底开发完成。

第2模态的时间系数和特征点累积形变量如图 8所示。由图可见，第2模态的时间系数从2016-07开始一直处于增长状态，对应选取的P₁、P₂、P₃、P₄特征点一直处于不断下沉的状态，直到2018-08后时间系数变为正值，此时所选特征点的累积形变量也由正值变为负值。随后时间系数曲线逐渐趋于稳定，特征点持续下沉的趋势也得到缓解，说明此时第2模态的分布特征基本成型。第2模态的分布特征与填挖方区域和时间对应较好，延安新区(北区)西南方向从2012年开始施工，2016年初整体沉降速率已经趋于稳定；2016-02~2018-08东北方向开始施工，包括该区域内的挖方填方以及地表建筑建设，上述工程使得该区域在这一时段内的形变速率加快，项目整体完工后，地表形变速率不断下降直至达到稳定状态。

4 结语

1) 根据获得的地表形变信息可知，发生沉降的区域大多集中在A、B区域，最大沉降速率为-56 mm/a。发生抬升的区域覆盖范围较广，最大抬升速率为32 mm/a。工程前期整体形变速率较快，工程结束时形变速率逐渐趋于稳定。

2) 通过EOF分解SBAS点得到空间分布特征和时间系数，第1模态的空间分布与地表形变的整体空间分布相似，相关性高达90.43%。对比施工前的高程与施工结束后的形变速率，证实工程填挖方与形变速率有关。第1模态时间系数趋于稳定的时间点与工程结束的时间点相吻合，证明延安新区(北区)的工程建设是影响其地表形变的主要原因；第2模态在空间分布上具有明显的分界线，分界线的西南方向开发时间为2012~2016年，后期地表形变速率已经趋于稳定，分界线东北方向从2016年开始施工，直到2018年才建设完成，这一点在第2模态时间系数的变化中也有所体现。