黄河流域是我国重要的生态屏障和经济地带，在推动我国社会经济发展与生态建设方面都具有非常重要的作用。黄河三角洲作为我国最大的三角洲，蕴含着丰富的卤水、石油、天然气等自然资源。随着人们对地下卤水的抽取及石油的不断开采，黄河三角洲存在多处不同程度的地表沉降区，严重影响当地居民的生产生活。因此，准确监测黄河三角洲地区地表沉降对于减少该地区的自然灾害和经济损失具有重要意义^[1]。

目前已有研究人员利用水准测量和全球定位系统对黄河三角洲地区地表沉降进行监测，这些测量方法虽有较高的精度，但存在费时费力、工作效率较低、地面测量点位易被破坏、测量点密度低、测站点不一定与沉降中心重合、难以进行大范围监测等诸多不足。2000年以来，永久散射体干涉测量(persistent scatterers InSAR, PS-InSAR)、短基线集干涉测量(small baseline subset, SBAS)等时序InSAR技术已得到较为广泛的应用。有学者利用PS-InSAR或SBAS方法对黄河三角洲区域进行地表沉降监测^[2]，虽然研究人员将InSAR监测结果与CORS站结果进行了精度验证^[2-3]，但该类时序InSAR方法多依赖于具有强散射特性的地物目标，而在裸地和多植被覆盖区域，其选点密度不高，难以反映沉降区域详细的沉降特征，降低了监测结果的准确性。此外，现有研究监测时段较早，无法反映黄河三角洲沉降的现状。

为解决常规时序InSAR的局限性，Ferretti等^[4]提出SqueeSAR方法，将相干性高的DS点与PS点联合进行形变解算，自此揭开了分布式目标InSAR(distributed scatterers InSAR，DS-InSAR)方法研究的序幕。在此基础上DS-InSAR不断发展，2015年Cao等^[5]详细论述各种相位优化方法的原理，并将相干矩阵采用特征值分解方法实现相位优化过程；2016年蒋弥等^[6]提出建立在快速同质点选取下的干涉处理框架，利用置信区间估计来判断像元的同质点，得到美国加州Lost Hills油田2002~2004年间的地表形变监测数据，表明该方法在运算效率和增加空间观测密度上具有明显优势；2019年Zhao等^[7]利用非线性相位优化，并基于StaMPS算法将PS点与优化后DS点联合处理，对美国图森地区进行地表形变监测，有效提高了相干点密度，但计算效率方面优势不大。目前DS-InSAR方法在山体滑坡、煤矿开采区^[8]、铁路沿线、水利设施^[9]等领域均得到成功应用，这些成功应用案例表明，DS-InSAR能够提高地表形变监测点的密度及低相干性区域形变监测的精度。

考虑黄河三角洲大范围区域存在大量非城市区域的状况，本文利用DS-InSAR方法获取2019-12~2020-12期间黄河三角洲区域的地表形变监测数据，并与PS-InSAR结果进行对比分析和交叉验证，证明结果的可靠性；最后，根据该地区人类活动特征及沉降中心分布情况分析黄河三角洲地表沉降的影响因素。

1 研究区域及数据概况 1.1 区域概况

黄河三角洲位于我国山东省东营市，北倚渤海湾，东邻莱州湾。现代黄河三角洲是1855年黄河改道携泥沙汇入渤海后，泥沙淤积形成的扇形冲积平原，面积约为5 450 km²，是我国最大的三角洲。该地区地势平坦，海拔高度多在6 m以下。现代黄河三角洲植被覆盖范围广，湿地覆盖面积约4 500 km²，且具有丰富的石油、天然气等自然资源。在特殊的地质构造下，现代黄河三角洲存在约60%面积的软土层，且分布不均匀。软土分布区含水量高、承载力低，因此易发生地表沉降。

本文研究区覆盖河口区、利津县及垦利县的大部分区域，包含胜利油田(中国第二大油田)、孤东油田、东营港口、东营机场及大量盐田区。由于黄河三角洲地区特殊的地质构造特点及大量开采油田、地下卤水等人类活动导致地层压力发生变化，该地区存在多个沉降区，且沉降量级大，严重影响了当地人民的生产生活。

1.2 实验数据

实验选取覆盖研究区的26景Sentinel-1A SAR影像，时间跨度为2019-12-19~2020-12-13，成像模式为IW模式，极化方式为VV极化，距离向和方位向的像元尺寸分别为2.33 m和13.93 m。为去除地形相位，本文采用90 m分辨率的SRTM数字高程模型(DEM)作为外部DEM数据。实验以2020-04-29影像作为公共主影像，分别与其他影像生成25个干涉对。

2 DS-InSAR方法原理 2.1 FaSHPS同质点识别

FaSHPS(fast SHP selection)算法^[6]的核心思想是将假设检验问题转化为置信区间估计，通过影像像元数据在时间序列上与参考点的相似性来判断是否为同质像元点。本文利用FaSHPS方法逐个计算出像元周围的同质点数量，以20个同质点作为阈值，将同质点数量高于该阈值的参考像元作为DS预选点。

由中心极限定理可知，时间序列上SAR影像数量越多，平均振幅A(S) 将越接近于高斯分布，此时，A(S) 的区间估计可表示为^[6]：

$ \begin{aligned} P\left\{\mu(S)-z_{1-\alpha / 2} \cdot \sqrt{{\rm{var}}({\boldsymbol{A}}(S)) / N}<\overline{{\boldsymbol{A}}}(S)<\right. \\ \left.\mu(S)+z_{1-\alpha / 2} \cdot \sqrt{{\rm{var}}({\boldsymbol{A}}(S)) / N}\right\}=1-\alpha \end{aligned} $

(1)

式中，P{·} 表示概率，z_1-α/2表示在标准正态分布中1-α/2分位点，μ(S) 为像元S的期望，var(A(S)) 表示像元S在时间序列上振幅的真实方差，N表示SAR影像的个数。

在同质区域可认为SAR影像的振幅服从瑞利分布^[10]，则变异系数CV可由下式计算：

$ {\mathrm{CV}}=\sqrt{{\rm{var}}(\cdot)} / E(\cdot)=\sqrt{4 / \pi-1} \approx 0.52 $

(2)

式中，E(·) 表示期望。假定像元的散射特性在观测时间范围内基本稳定，则式(1)可转换为：

$ \begin{aligned} P\left\{\mu(S)-z_{1-\alpha / 2} \cdot \frac{0.52 \cdot \mu(S)}{\sqrt{N}}<\overline{{\boldsymbol{A}}}(S)\right. \\ \left.<\mu(S)+z_{1-\alpha / 2} \cdot \frac{0.52 \cdot \mu(S)}{\sqrt{N}}\right\}=1-\alpha \end{aligned} $

(3)

将A(S) 作为像元S在时间维度上振幅均值的真值，根据式(3)可求出振幅均值的置信区间。计算在时间序列上相邻像元的平均振幅，若其处于参考像元平均振幅的置信区间内，则认为二者是同质点。

2.2 特征值分解相位优化

分布式目标像元内地物应具有相同的散射特性，由于受到噪声及时空失相干的影响，时间序列SAR影像干涉图每个像元的信号矢量中包含多种类型的散射信号，使得分布式目标像元的相位稳定性较差。通过对像元的相位优化可达到对干涉图像元相位降噪的目的。像元相干矩阵对应的特征值代表不同的后向散射特性，将相干矩阵进行特征值分解^[5]即可得到最大特征值及其对应的特征向量，认为此特征向量对应的像元散射特性是稳定的，并将其作为优化后的相位。

假设在匀质区域Ω内部存在N_P个后向散射特性相近的像元，其相干矩阵T可表示为：

$ {\boldsymbol{T}}=E\left[{\boldsymbol{y}} {\boldsymbol{y}}^{{\mathrm{H}}}\right] \approx \frac{1}{N_P} \sum\limits_{{\boldsymbol{y}}=\Omega} {\boldsymbol{y}} {\boldsymbol{y}}^{{\mathrm{H}}} $

(4)

式中，y =[y₁, y₂, …, y_N] 为分布式目标的同质点在N景SAR影像上的复数观测量经归一化处理后的复数向量，(·)^H表示矩阵的共轭转置。上式得到的相干矩阵为半正定Hermitian矩阵，经特征值分解可得：

$ {\boldsymbol{T}}={\boldsymbol{U}} {\boldsymbol{\varLambda}} {\boldsymbol{U}}^{-1}={\boldsymbol{U}} {\boldsymbol{\varLambda}} {\boldsymbol{U}}^{{\mathrm{H}}}=\sum\limits_{i=1}^N \lambda_i \cdot {\boldsymbol{u}}_i \cdot {\boldsymbol{u}}_i^{{\mathrm{H}}} $

(5)

式中，Λ为非负实数特征值λ_i的对角矩阵，U为不同特征值对应的正交特征向量。λ_i越大，对应的地物散射机制就越占优势，因此，可将最大特征值λ₁对应的特征向量μ₁作为主散射体对应的相位。由此可计算主散射体信号T_signal，其余的作为去相干噪声信号T_noise。将去相干噪声信号的相位分量去除，保留主散射体信号的相位分量，即可达到优化相位的目的：

$ \begin{gathered} {\boldsymbol{T}}={\boldsymbol{T}}_{\text {signal }}+{\boldsymbol{T}}_{\text {noise }}= \\ \lambda_1 \cdot {\boldsymbol{\mu}}_1 \cdot {\boldsymbol{\mu}}_1^{{\mathrm{H}}}+\sum\limits_{i=2}^N \lambda_i \cdot {\boldsymbol{\mu}}_i \cdot {\boldsymbol{\mu}}_i^{{\mathrm{H}}} \end{gathered} $

(6)

相位优化完成后，还需对相位的优化质量进行评估。将优化前后的干涉相位进行差值拟合计算得到拟合度^[4]：

$ \gamma_{{\mathrm{DS}}}=\frac{2}{N(N-1)} {\rm{Re}} \sum\limits_{n=1}^N \sum\limits_{m=n+1}^N \exp ^{j\left(\varphi_{m n}-\left(\varphi_m-\varphi_n\right)\right)} $

(7)

式中，γ_DS为拟合度，还可看作为分布式目标的时序相干性；φ_mn表示两幅SAR影像优化之前的干涉相位；φ_m和φ_n分别表示两幅SAR影像优化之后的相位。本文以0.4作为阈值，将相位优化质量高于阈值的像元选定为最终的DS点。

2.3 处理步骤

本文利用FaSHPS算法获得DS候选点，经特征值分解相位优化及优化质量评估后确定最终DS点，将得到的DS点联合振幅离差阈值初选得到的PS点，构建用于相位分析的Delaunay三角网。由于DEM误差的不确定性，利用StaMPS软件首先对视角误差的空间不相干部分进行相位校正，然后在时间和空间两个维度进行干涉相位的3D解缠，利用时间上的高通滤波器及空间上的低通滤波器对具有时间和空间相关的相位进行滤波，以估计剩余的空间相关干扰项，如大气相位和轨道相位，在减去上述分量之后最终可提取黄河三角洲区域的地表时序形变。具体的处理流程如图 1所示。

3 InSAR监测结果分析 3.1 地表形变监测结果

分别采用PS-InSAR和DS-InSAR两种方法提取黄河三角洲2019-12-19~2020-12-13时序地表沉降信息(图 2)。PS-InSAR方法共选取425 729个监测点，最大沉降速率为-227 mm/a；而DS-InSAR方法共选取2 365 545个监测点，约是PS-InSAR方法选点个数的5.56倍，最大沉降速率为-238 mm/a。可见，DS-InSAR方法显著增加了地面监测点的空间分布密度。图 2表明，在该研究时间段内黄河三角洲区域存在严重的地表沉降现象，两种方法均监测到4处明显的大量级沉降区域，且形变结果接近，具有较好的一致性。本文监测得到的沉降区域与文献[11]得到的2015~2018年黄河三角洲的沉降区域相吻合，且通过实地调研发现，这些沉降区域存在大量的盐田池和采油机，表明本实验可准确识别沉降区域。由于DS-InSAR方法选点数量增多，在沉降区域其监测效果要优于PS-InSAR方法，沉降边界更加明显。

3.2 监测结果分析

由于黄河三角洲区域的地物覆盖类型多以农田、盐田与养殖池为主，导致PS-InSAR方法选点数量较少，使得地表形变的分布、特征、范围不够显著；而融合分布式目标的DS-InSAR方法有效提高了非城市区域的地面监测点密度。本文监测到的4个沉降区域主要分布在盐田和油田附近，其中A、B、C区域均位于东营市河口区，且分布较为集中；D区域位于东营市垦利县南部。

为分析PS-InSAR与DS-InSAR所获取的地表沉降的差异性，对4个沉降区域分别选取了剖面线(图 3)，并绘制了对应的地表沉降散点图(图 4)。可以看出，DS-InSAR有效提高了测量点空间分布密度，能够更加直观地反映沉降区域地表沉降变化趋势；同时也明显看出，两者具有较高的一致性，验证了DS-InSAR结果的可靠性。

A沉降区中心位于东营市河口区广河村国星盐场^[12]，其最大年沉降速率为-134 mm/a。由于存在过度开采地下卤水的情况，导致盐场附近出现较大范围的沉降漏斗，图 4(a)显示沿A′-A″线距A′点约5~12 km的范围内，地表沉降速率达到-75 mm/a。

B沉降区中心位于盐田与油田混合区域，最大年沉降速率达-238 mm/a。图 4(b)为B′-B″线沉降速率散点图，该区域年沉降速率较大，沉降现象最为严重。该区域被黄河故道、东营港疏港高速及孤北水库包围，自2015年就已出现沉降现象，由于盐田分布较为集中且开采力度大，如今地表沉降速率明显加快。

C沉降区中心位于桩西油田、孤东油田及盐田油田综合开采区，其最大年沉降速率达-155 mm/a。由图 4(c)、(d)可知，在C₁-C′₁剖面方向沉降现象较为严重，沉降速率可达-150 mm/a；沿C₂-C′₂线由西北向东南方向地表存在持续缓慢沉降现象，其沉降速率约-60 mm/a。由于该区域存在大面积的盐田区域，并且零星分布着许多游梁式抽油机，因此在地下卤水开采及油气开采的双重影响下，该地区发生了一定程度的地表沉降，并且其影响范围接近100 km²。

D沉降区中心位于垦利县东营机场附近的永丰盐场和东营景洪盐化有限公司，最大年沉降速率达-122 mm/a。由图 4(e)可看出，剖面线上的3个沉降漏斗对应于永丰盐场周围3个集中晒盐区域，其中两个区域地表沉降速率均已超过-80 mm/a。另外，D′-D″剖线端点D′附近位于东营机场，其年沉降速率达到-30 mm/a。

3.3 DS-InSAR可靠性分析

为进一步验证DS-InSAR方法的可靠性，本文通过对比分析PS-InSAR与DS-InSAR在一定范围内同名监测点的沉降速率进行交叉验证。将PS-InSAR监测结果作为基准，以每个PS点为中心，15 m为半径，找到与PS点同名的DS点。一共得到144 772个PS-DS同名点，利用皮尔逊相关性分析算法得到同名点的沉降速率相关性，并进一步计算误差分布(图 5)。

由图 5(a)可知，本文采用的DS-InSAR与PS-InSAR方法监测得到的地表沉降速率具有良好的一致性，其同名点沉降速率的皮尔逊相关系数为0.727。图 5(b)表明，沉降速率误差分布趋近于高斯分布，且均值接近0，其中误差值小于20 mm/a的同名点个数为128 149，占所选同名点个数的88.5%。因此，DS-InSAR在提高地面目标点密度的同时，具有很好的可靠性。

4 黄河三角洲地表形变影响因素分析

黄河三角洲地区地下卤水丰富，且含盐量高于海水，研究区域内的河口区、垦利县、东营区已建成多处盐场并投入生产，这些工厂开采地下卤水用于工业制溴、制盐和制造相关化工产品等，但长期不合理开采地下卤水将会导致地下卤水水位下降、局部海水入侵和地面裂缝，严重时则会导致地表沉降、发生地质灾害。众多企业大量开采地下卤水不仅致使其资源减少，而且形成浅层地下卤水沉降漏斗。因此，过度开采地下卤水用于制盐和化工产业是造成黄河三角洲地区地面沉降的主要原因之一。

本文监测到的A、B、C、D沉降区中心均为盐场及盐田开采区，其中C沉降区的中心为盐田和油田综合开采区。图 6为叠加了盐田区和居民区大致范围的典型沉降区域的卫星影像图。可以看出，本实验监测得到的沉降区大都分布在盐田区域，说明黄河三角洲区域开采地下卤水制盐是造成地面沉降的主要因素之一。A沉降区域(国星盐场)附近有少量的居民区，未发现该区域地表形变对居民区造成明显影响；最大沉降区域B紧邻东营疏港高速公路和仙河镇居民区，长期过度开采地下卤水必将对高速公路和居民区造成损害；C沉降区域为盐田和油田的综合开采区，从五号桩以南沿兴港路两侧分布有大量的盐田，且同时存在许多采油机，在二者的共同影响下造成大面积地表沉降，且附近建构筑物出现裂缝现象；D沉降区域紧邻红光新村和东营机场，目前其沉降速率缓慢，但永丰盐场的卤水开采制盐活动对东营机场附近的地基稳定性造成一定程度的影响。

东营市胜利油田作为中国第二大油田，其石油开采活动从上世纪开始一直在进行，长期的深层石油开采活动会增加储油层压力，使得储油层压实收缩，从而导致缓慢的地表沉降，而人为进行地下水回注能够起到减缓地表沉降的作用，甚至会发生小幅度地面抬升现象^[13]。此外，人工回注过程中，往往需要从油田周围地区抽取地下水，这又导致油田周围地区进一步发生地面沉降。沉降中心C的孤岛-孤东油田便是胜利油田的重要采油区之一，由图 6(a)可看出，油田中心区域的沉降量级较大，地面沉降影响范围广，周围地区也发生了缓慢的地表沉降现象。由此可见，黄河三角洲地区的地表沉降现象与油气开采活动也有一定关系。

5 结语

本文采用DS-InSAR方法对2019-12~2020-12期间26景Sentinel-1A影像进行数据处理及结果分析，得到黄河三角洲区域的地表沉降速率，结论如下：

1) 与传统PS-InSAR方法相比较，融合分布式目标的DS-InSAR方法有效提高了地面监测点的密度，能够详细描述沉降区域特征，实现了地表形变的精细化监测，并且该方法具有很好的可靠性。

2) 在监测时间段内，监测到黄河三角洲存在多处不同程度的地表沉降现象，最大沉降速率达-238 mm/a。同时利用DS-InSAR方法探测到黄河三角洲地区较为准确的时序地表形变趋势，可为地表形变规律及地表沉降诱发因素的研究与分析提供参考。

3) 导致黄河三角洲地表形变的因素主要有过度开采地下卤水、超负荷开采油气资源等。由于缺乏实测数据及水文地质资料等，本文对黄河三角洲地表形变影响因素的分析稍有欠缺，之后可根据实地调研情况结合相关资料对地表形变机理作进一步分析与完善。