西安市与咸阳市地质构造复杂，自20世纪60年代以来一直受到地面沉降和地裂缝的影响^[1]。为对西安-咸阳地区的地面沉降进行防治和预警，需要开展地面沉降监测研究。祝意青等^[2]使用精密水准测量研究西安地面沉降演化特征及成因；张勤等^[3]使用GPS技术监测西安市地面沉降与地裂缝。但这2种方法均只能对地面有限的离散点进行观测，覆盖区域较小。InSAR是近年来迅速发展的空间大地测量技术，具有全天时、全天候特点，可以获取大范围、高精度的地表高程与形变信息，已成为大范围地面沉降监测的主要手段^[4]，其衍生技术PS-InSAR、SBAS-InSAR等也已成功应用于很多城市的地面沉降监测^[5-12]。

本文利用时序InSAR技术研究西安-咸阳地区2015-06~2023-01的地面沉降，并结合有关资料分析该地区的沉降特征及时空演变规律。

1 研究区及数据来源 1.1 研究区概况

西安-咸阳地区(图 1)位于我国西北部关中盆地，属于温带大陆性季风气候，春秋短暂，夏季炎热多雨，冬季寒冷干燥。研究区水资源匮乏，相比于我国东部同纬度季风地区地表径流较少，主要的供水源是地下水。但长期大量开采地下水导致区域出现严重的环境地质问题，如地下水位持续下降、地面沉降、地裂缝和断层发育等。为此，政府已采取一系列措施，如限制地下水开采、人工回灌地下水等，以减缓地面沉降的速度。

1.2 数据来源

Sentinel-1A卫星于2014-04发射，重访周期为12 d。本实验采用覆盖西安-咸阳地区的342景升轨Sentinel-1A数据，成像模式为干涉宽幅，极化方式为VV。同时，采用NASA发布的30 m分辨率SRTM DEM去除地形残差相位，使用欧空局提供的精密定轨星历数据提高Sentinel-1A卫星影像的轨道精度。

2 研究方法

PS-InSAR和SBAS-InSAR是2种常用的时序InSAR技术，可以降低InSAR反演形变过程中时空失相干的影响，是监测长期、缓慢形变的有效方法。本文采用这2种方法监测西安-咸阳地区地面沉降，分析2种方法所得结果，其中PS-InSAR使用StaMPS-InSAR(Stanford method for PS-InSAR)技术^[13]。

2.1 SBAS-InSAR

SBAS-InSAR技术^[6]通过设定的时空基线阈值将所有SAR数据组成若干个小集合，集合内基线较小。对每个干涉像对进行差分干涉处理，获取正确的解缠相位，最后采用最小二乘法或奇异值分解法得到整个时间序列地表形变参数的估计。

假设有同一研究区的N+1景SAR影像，则依据设定的时间基线和空间基线阈值可生成M幅干涉图，M需满足：

$ \frac{N+1}{2} \leqslant M \leqslant \frac{N(N+1)}{2} $

(1)

第k(k=1, …, M)幅差分干涉图任意像素(x, r)的干涉相位φ_int可表示为：

$ \varphi_{\mathrm{int}}=\varphi_{\mathrm{def}}+\varphi_{\mathrm{dem}}+\varphi_{\mathrm{atm}}+\varphi_{\mathrm{orb}}+\varphi_{\mathrm{noi}} $

(2)

式中，φ_int为去除平地相位和高程相位后干涉图的差分相位；φ_def为地面沿雷达视线方向的形变相位；φ_dem为高程误差引起的地形残差相位；φ_atm为干涉像对的大气延迟量之差；φ_orb为轨道误差造成的残差相位；φ_noi为噪声相位。

假定φ_noi很小，在去除大气、轨道和地形相位后，式(2)可表示为：

$ \varphi_{\mathrm{int}} \approx \frac{4 \pi}{\lambda}\left(r_B-r_A\right)=\frac{4 \pi}{\lambda} \delta_r $

(3)

式中，λ为雷达波长；r_A、r_B分别为第k幅差分干涉图中A、B两景SAR影像对应雷达视线向的地表形变；δ_r为该时段的地表形变量。

将所有差分干涉图进行自由组合，φ为N景SAR影像上的干涉相位值组成的矩阵，δφ为M幅差分干涉图上相位组成的矩阵，则：

$ \boldsymbol{\varphi}=\left[\begin{array}{lll} \varphi\left(t_1\right) & \cdots & \varphi\left(t_N\right) \end{array}\right]^{\mathrm{T}} $

(4)

$ \delta \boldsymbol{\varphi}=\left[\begin{array}{lll} \delta \varphi\left(t_1\right) & \cdots & \delta \varphi\left(t_M\right) \end{array}\right]^{\mathrm{T}} $

(5)

将上式写成矩阵形式：

$ \boldsymbol{A} \boldsymbol{\varphi}=\delta \boldsymbol{\varphi} $

(6)

式中，A为M×N维矩阵。

小基线集合数L=1时，矩阵A列满秩，采用最小二乘法可求解累积形变相位的估值$\hat{\boldsymbol{\varphi}} $；小基线集合数L>1时，矩阵A秩亏，秩亏数为N-L+1，需对系数矩阵A进行奇异值分解，以求出累积形变相位φ最小范数意义上的最小二乘解。

2.2 StaMPS-InSAR

StaMPS-InSAR^[13]是PS-InSAR的改进，其使用时间相干性对点目标进行识别。首先在N+1景时间序列影像中选取1幅影像作为主影像，其余N景作为辅影像进行差分处理，使用振幅离差选择PS候选点，每个点的干涉相位表达式与式(2)相同。

φ_def、φ_atm、φ_orb在空间一定范围内相关，而φ_dem和φ_noi在空间不相干。根据该特点，对每个干涉图中各像素计算周围相干像素的空间平均值φ_int：

$ \bar{\varphi}_{\mathrm{int}}=\bar{\varphi}_{\mathrm{def}}+\bar{\varphi}_{\mathrm{atm}}+\bar{\varphi}_{\mathrm{orb}}+\bar{\varphi}_n $

(7)

式中，φ_n为φ_dem与φ_noi的均值之和。假定φ_n很小，式(2)与式(7)相减得：

$ \varphi_{\mathrm{int}}-\bar{\varphi}_{\mathrm{int}}=\varphi_{\mathrm{dem}}+\varphi_{\mathrm{noi}}+\bar{\varphi}_n^{\prime} $

(8)

式中，$\bar{\varphi}_n^{\prime}$为相位噪声与原始φ_def、φ_atm和φ_orb与其均值之差的综合值。假定φ_noi和$\bar{\varphi}_n^{\prime}$很小，则式(8)可表示为：

$ \varphi_{\mathrm{int}}-\bar{\varphi}_{\mathrm{int}} \approx \varphi_{\mathrm{dem}}=B_{\perp} K_{\mathrm{dem}} $

(9)

式中，B_⊥为垂直基线；K_dem为比例常数。

地形残差相位只和基线相关，二者呈线性关系，因此可以基于所有可用干涉图使用最小二乘平差估计地形残差相位。$\hat{\varphi}_{\mathrm{dem}}$使用式(10)计算每个点的时间相干性γ，假定$ \bar{\varphi}_n^{\prime}$很小，γ可以用来评价PS候选点的相位噪声水平，以确定最终的PS点：

$ \gamma=\frac{1}{N}\left|\sum\limits_{i=1}^N \exp \left[\sqrt{-1}\left(\varphi_{\mathrm{int}}-\bar{\varphi}_{\mathrm{int}}-\hat{\varphi}_{\mathrm{dem}}\right)\right]\right| $

(10)

对式(2)去除地形残差相位得：

$ \varphi_{\mathrm{int}}=\varphi_{\mathrm{def}}+\varphi_{\mathrm{atm}}+\varphi_{\mathrm{orb}}+\varphi_{\mathrm{noi}} $

(11)

最后对干涉相位进行三维解缠，并采用时空滤波方法去除大气与轨道误差，分离出形变相位。

3 西安市地表形变时空分布特征及成因分析 3.1 西安市地表形变时空分布特征

使用SBAS-InSAR进行数据处理，获得西安市2015-06~2023-01地表形变速率分布(图 2)。由图可知，西安市地表形变表现出明显的空间差异性，大多数形变区域受到地裂缝和断层的影响，通常沿东北-西南方向延伸。为更好地研究西安市地表形变的时空分布特征，分别绘制每年形变速率分布，结果如图 3所示。可以看出，西安市共有5个主要的形变区域，分别为鱼化寨(P₁)、城墙南部(P₂)、余家庄(P₃)、电子城(P₄)以及凤栖原(P₅)-等驾坡(P₆)区域。其中，鱼化寨区域在2015~2017年为西安市沉降最剧烈的区域，沉降速率最大达到-150 mm/a。2018年年底，鱼化寨区域由沉降转为快速抬升，至2019年年中，其抬升速率逐渐趋于平缓，之后未发生较大形变。城墙南部区域为西安最大的形变区域，主要受f₅、f₆和f₇三条地裂缝的影响，呈东北-西南走向，形变最大区域位于石油大学生活区，最大平均抬升速率达18 mm/a，整个区域在监测期间一直为抬升趋势。余家庄区域位于鱼化寨南部，形变速率较缓，平均抬升速率为5~10 mm/a。电子城区域形变过程与鱼化寨区域类似，在2021年之前一直处于沉降阶段，2021年年初开始由沉降转为快速抬升，最大速率达100 mm/a。凤栖原-等驾坡区域为东北-西南走向，其中有2个沉降剧烈的地区，分别为凤栖原和等驾坡，二者一直在沉降，但沉降范围不断缩小。

3.2 成因分析

地下水开采与地面沉降关系密切。含水层系统是一个饱和非均质、互层渗透和低渗透的水文地质单元，由水力连接的含水层和弱透水层组成。地下水抽水量的增加会导致地下水水位下降，从而使上方物质的部分压力从填充孔隙空间的加压流体转移到含水层系统的颗粒骨架，致使孔隙压力下降，最终导致地表下沉。大量研究表明，西安市地面沉降主要是由地下水抽取造成，地面沉降中心与承压水开采区大致吻合，历史上西安市地下水开采阶段与地面沉降阶段基本一致^[14]。西安市于1990年实施黑河引水计划，此后于2007年通过《西安市地下水资源管理条例》，2018年停止开采地下水，并开始人工回灌地下水。对比西安市沉降时空分布特征可以发现，2018-11西安市实施鱼化寨地区地面沉降应急治理地下水回灌专项工程，随后鱼化寨地区地表由沉降转为抬升，说明地下水开采与地面沉降关系密切。

此外，地裂缝和断层也与地面沉降有一定的关系。西安地区地裂缝是由地下水开采引起的地面沉降所致。西安市地表沉降带、地裂缝和长安-临潼断层(CAF)均大致呈东北-西南走向，这表明西安市地面沉降的空间分布在一定程度上受地裂缝和断层控制。绘制沿剖面线AA′的累积形变量，同时将主要地裂缝和断层叠加在剖面上，结果如图 4所示。可以看出，地裂缝和断层两侧的沉降梯度多数都发生变化，这表明西安地裂缝和地面沉降之间存在明显的相关性。

综上可知，地下水的开采与回灌是目前引起西安市地表形变的主要因素，同时地裂缝和断层也会影响西安地表形变的走向与时空分布。

4 咸阳市地表形变时空分布特征及成因分析 4.1 咸阳市地表形变时空分布特征

图 5为使用SBAS-InSAR获得的咸阳市2015-06~2023-01地表形变情况。由图可知，与西安市相比，咸阳市地表形变总体上并不剧烈，大部分地区的形变速率在-5~5 mm/a之间，仅有局部地区沉降明显。选择6个沉降明显的区域绘制累积沉降序列(图 6)，分别为咸阳机场南部(Q₁)、艺术中心(Q₂)、周陵(Q₃)、西咸广场(Q₄)、两寺渡(Q₅)和马村(Q₆)地区。由图可知，咸阳机场南部、艺术中心和马村区域累积沉降量均呈匀速下降状态，累积沉降量分别为-135 mm、-76 mm和-56 mm。周陵区域在2017年之前未发生明显形变，至2017-01累积沉降量仅为-1 mm，之后沉降速率明显加快，至监测末期累积沉降量达到-133 mm，仅次于咸阳机场南部区域。西咸广场区域在大部分时间沉降速率均为-9 mm/a，仅在2021年出现沉降加速的情况，累积沉降量为-82 mm。两寺渡区域在监测早期处于稳定阶段，至2018-01累积沉降量仅为-2 mm，之后沉降加速，沉降速率在2021年达到最大，为-15 mm/a，整个监测期间累积沉降量为-79 mm。

4.2 成因分析

咸阳市地面沉降的发展过程与地下水开采关系密切。咸阳市作为陕西省重要的工农业基地，生产生活用水量巨大，早期基本只能依靠地下水供给，地下水过度开采引发地面沉降。咸阳市地面沉降的发展过程可以分为发生、发展和稳定3个阶段^[15]。20世纪80年代，由于地下水开采导致地下水位下降，咸阳市开始出现地面沉降。1991~2006年随着社会发展，地下水抽水量逐年增加，地下水水位快速下降，地面沉降速率明显增加。2006年开始，咸阳市采取部分措施，如关停自备井、实施“引石过渭供水工程”等，降低地下水开采量。这些措施使地下水水位逐渐趋于稳定，并有所回升，地面沉降也基本稳定，不再出现区域性地面沉降。

另外，城市建设也是咸阳市地面沉降的重要影响因素。自2011年咸阳市人民政府编制《咸阳市城市总体规划(2011~2030)》后，咸阳市开启新一轮的城市化进程。咸阳市先后对机场、铁路、地铁、工厂、小区等基础设施进行大规模的改造与建设，城镇化率由2015年的49%增长到2021年的57.25%。这些工程的施工对地面沉降造成一定影响，如咸阳国际机场三期扩建工程导致咸阳国际机场南部沉降，城际铁路车辆段和地铁14号线的修建引发艺术中心沉降等。这是因为土体上方的荷载逐渐增大，土体所受应力状态发生变化，致使土层发生压缩变形，从而导致地面沉降。这种方式引发的地面沉降较开采地下水更为实时、短暂，并且只对局部地区有影响。

综上，咸阳市地面沉降受地下水开采、城市建设等多方面因素的共同影响。

5 交叉验证

由于本次实验缺乏西安-咸阳同期水准和GNSS数据，故采用StaMPS-InSAR技术获取西安-咸阳地区地表形变，并与SBAS-InSAR的结果对比，以进行交叉验证。图 7和8分别为使用StaMPS-InSAR技术获取的西安市和咸阳市地面形变速率图，通过绘制西安-咸阳各特征区域在监测期间的累积形变量来对比分析2种技术获取的形变结果(图 9、10)。

从图 9和10可以看出，StaMPS-InSAR所得的形变曲线波动较大，而SBAS-InSAR所得的形变曲线更为光滑，趋势也更为稳定。2种方法在多数区域的形变曲线趋势和形变结果基本一致，仅在个别点上存在差异，如西安地区P₁点，2种方法获得的累积形变量相差近60 mm，其原因可能为模型参数以及环境差异所致。

计算在相同经纬度下StaMPS-InSAR与SBAS-InSAR的形变速率差值，结果如图 11所示。由图可知，2种方法的形变速率差值99%在-5~5 mm/a之间，95%在-3~3 mm/a之间，说明2种方法具有相似的形变速率。总体而言，2种方法获得的结果具有较强的一致性，形变区域大致重合，实验结果具有较高的可靠性。

6 结语

本文利用SBAS-InSAR与StaMPS-InSAR技术对2015~2022年西安-咸阳地区地表形变进行监测，得到以下结论：

1) 西安地区地表形变呈现明显的空间差异性，北部较为稳定，中部和南部发生大规模形变；在监测期间，中部形变由沉降变为抬升，而南部一直在沉降。西安地区地表形变的主要原因为地下水开采与回灌，同时地裂缝和断层也会影响西安市地面形变的走向与时空分布。

2) 咸阳市大部分地区处于稳定状态，形变速率在-5~5 mm/a之间，仅局部区域沉降较为明显。咸阳市地面沉降是地下水过量开采、大规模城市建设等多方面因素影响的结果。

3) SBAS-InSAR与StaMPS-InSAR获得的地表形变结果具有较高的一致性，形变趋势与形变范围基本相同，证明了本文结论的可靠性。