全球卫星导航系统(GNSS)是一种传统的大地测量技术手段，具有全天候、高精度、实时服务能力强的特点，可直接获取地面监测点的形变信息，但空间分辨率低，难以对整个区域进行监测。近年来新兴的合成孔径雷达干涉测量(interferometric synthetic aperture radar, InSAR)技术具有低成本、大范围、高空间分辨率的特点，可以得到cm级甚至mm级地表形变监测结果。InSAR技术主要有合成孔径雷达差分干涉测量(DInSAR)技术、永久散射体干涉测量(PS-InSAR)技术和短基线集干涉测量(SBAS-InSAR)技术。其中，SBAS-InSAR采用多主影像，时空基线为短基线，可有效解决时空失相干的问题，适用性更强。GNSS与InSAR两种监测方法各有优缺点，将两种方法相结合可以得到数据质量更优且更可靠的地面监测数据。杨国创等^[1]对InSAR和GNSS技术所得的形变量进行分析，结果表明监测的沿海地区地表形变幅度和趋势高度一致，说明两者具备数据融合的条件。周文韬^[2]提出利用方差分量估计方法融合两者数据的三维形变监测并取得良好效果。但该方法缺少InSAR监测量的优化处理，由于InSAR技术监测的地表形变量包含自然或人为破坏、地表气温变化、降雨等气象因素引起的误差，无法直接用于与CORS数据融合处理。

本文对两种技术所得的数据进行优化处理，首先通过时序InSAR处理获取研究区地表形变量，将其分离出数米以下的岩土层形变量；然后以卫星导航定位基准站(CORS)数据为基准，对InSAR监测结果进行整体平差，修复其差分干涉尺度误差，补偿空间中长波误差影响，控制InSAR监测量随时间的累积误差等，从而实现高分辨率、高精度的地面形变监测。

1 研究区和实验数据 1.1 研究区概况

北京地处我国华北平原北部，自20世纪70年代以来，由于城市的开发建设，对地下水需求加大，导致地下水过度开采，同时城市高层建筑不断增多，造成朝阳、通州等地出现严重的地面沉降^[3]。对地面沉降进行大范围、高效的时空监测不仅能为城市建设、经济发展提供可靠的科学依据，也有助于做好地质灾害防范工作。

1.2 实验数据

时序InSAR处理使用31景Sentinel-1A卫星影像数据，起止时间为2018-01-03~2020-11-12。使用POD精密定轨星历数据修正轨道误差，采用空间分辨率为30 m的SRTM1 DEM数据共同参与本次处理。本文所用的CORS站数据为昌平站(CHPN)、牛栏山站(NLSH)、东三旗站(DSQI)、朝阳站(CHAO)、西集站(XIJI) 2018~2020年连续观测数据。

2 数据预处理与融合方法 2.1 时序InSAR处理

短基线集干涉测量(SBAS-InSAR)原理为：假设研究区内影像有N+1幅，选择一幅影像为超级主影像，并且任意一幅影像至少可以与其他影像形成一个干涉对，共生成M幅差分干涉图，M满足以下条件：

$ \frac{N+1}{2} \leqslant M \leqslant N \cdot \frac{N+1}{2} $

(1)

假设在t_A、t_B(t_A < t_B)时刻生成第j幅干涉图，则该干涉图上任意像素点(x, r)相对应的干涉相位可表示为：

$ \begin{gathered} \delta \varphi_{j}(x, r)=\varphi\left(t_{B}, x, r\right)-\varphi\left(t_{A}, x, r\right) \approx \\ \frac{4 \pi}{\lambda}\left[d\left(t_{B}, x, r\right)-d\left(t_{A}, x, r\right)\right]+ \\ \Delta \varphi_{\text {topo }}^{j}(x, r)+\Delta \varphi_{\mathrm{APS}}^{j}\left(t_{B}, t_{A}, x, r\right)+\Delta \varphi_{\text {noise }}^{j}(x, r) \end{gathered} $

(2)

式中，j∈(1, …, M)；λ为雷达信号中心波长；d(t_A, x, r)和d(t_B, x, r)为t_A、t_B时刻的雷达视线向(LOS)累积形变量；Δφ_topo^j(x, r)为差分干涉图中地形残余相位；Δφ_APS^j(t_B, t_A, x, r)为大气延迟相位；Δφ_noise^j(x, r)为噪声相位。式(2)可进一步写成：

$ \boldsymbol{B} \boldsymbol{v}=\delta \boldsymbol{\varphi} $

(3)

式中，B为系数矩阵，v为变化率，δ表示估计值，φ为相对差分相位。B在满秩时可利用最小二乘法求解，秩亏时需用奇异值分解法求得最小范数解，进而得到研究区形变速率^[4]。由式(4)可将LOS向形变量时间序列转换为大地高方向形变量时间序列：

$ D_{U}=D_{\mathrm{LOS}} / \cos \theta $

(4)

式中，D_U为大地高方向形变量，D_LOS为雷达视线向形变量，θ为卫星雷达方向入射角。

本次实验采用SARscape软件，选择SBAS-InSAR处理方法，时间基线阈值为120 d，共生成87对干涉对，时空基线分布如图 1所示。

2.2 InSAR监测量优化方法

由于InSAR监测量中存在野值、粗差和突变等非动力学形变信号，需要将其进行探测与分离。以InSAR监测量采样历元时刻为单元，由高斯低通滤波器构造低通监测量参考面，再以3倍标准差为阈值进行粗差探测。高斯低通滤波器^[5]定义如下：

$ H(u, v)=\mathrm{e}^{-D(u, v)^{2} / 2 \sigma^{2}} $

(5)

式中，e为自然常数；σ为每单元InSAR监测量的标准偏差；(u, v)为图像点坐标；H(u, v)为滤波后的结果；D(u, v)表示从点(u, v)到频率平面原点的距离，即$D(u, v)=\sqrt{u^{2}+v^{2}}$。

由于在每个采样历元中所分离的粗差点不同，为保证监测量在整个时序上的时空分布，依据动力学地面垂直形变量与动力源作用点距离或距离平方近似成反比的空间变化性质，以InSAR监测量采样历元时刻为单元，按高斯滤波算法计算粗差点的监测量，对粗差点进行修复，从而保证监测点的时空分布不变。通过对InSAR监测量进行优化处理，以抑制或削弱非地质动力学作用的极浅地表局部变化影响，从而得到地表数米以下的深岩土层垂直形变^[6]。

2.3 融合CORS网与时序InSAR协同监测

采用GAMIT/GLOBK软件进行CORS站连续观测数据的基线解算，获得单日解文件，计算CORS站大地高周变化时间序列。对时间序列进行粗差探测，并分离线性项和常数项，重构非线性项。采用连续Fourier和切比雪夫组合基函数，由不规则采样时序构造低通滤波曲线参考基准，计算采样值与参考值之差，并对残差时序进行统计，将大于指定倍数残差标准差的采样记录作为粗差并分离。剔除粗差后，利用连续切比雪夫基函数分离CORS站大地高周变化时间序列的常数项和线性项^[7]，并对非线性项进行低频重构。将重构后的非线性项与线性项、常数项相加，得到新的CORS站垂直方向周变化时间序列。

去除InSAR监测量中极浅地表局部变化影响，得到地表数米以下的深岩土层垂直形变量，与CORS站数据进行融合。选取一个采样历元作为CORS数据与InSAR数据融合的参考历元，将两部分数据的参考历元进行统一。由CORS站周边InSAR监测量内插得到CORS站处的InSAR监测量，对CORS站大地高变化时序按照时间内插得到InSAR采样历元时刻的CORS站大地高变化。以CORS站为基准，根据每期InSAR散点间监测量的相对变化量，采用间接平差方法对全部InSAR监测量进行整体平差，从而实现CORS网与InSAR监测时序的高度统一与高精度传递。平差模型^[8]为：

$ \underset{n \times 1}{\boldsymbol{V}}=\underset{n \times {u} }{\boldsymbol{B}} \underset{u \times 1}{\hat{\boldsymbol{x}}}-\underset{n \times 1}{\boldsymbol{l}} $

(6)

式中，V为每期InSAR监测量间相对变化量的改正数向量，B为系数矩阵，$\hat{\mathit{\boldsymbol{x}}}$为每期InSAR监测量的平差值，l为误差方程常数项。根据最小二乘原理可求得$\hat{\mathit{\boldsymbol{x}}}$，即平差后的每期InSAR监测量。

通过融合CORS网与时序InSAR数据，可以将时序InSAR高空间分辨率与CORS站数据高精度的优点结合起来，构建具有高分辨率、高精度的地面垂直方向形变量时间序列。

3 实验结果与分析 3.1 InSAR结果分析

将北京地区SBAS-InSAR处理结果转化为垂直方向形变量时间序列，绘制年平均形变速率，结果如图 2所示。由图可知，沉降比较明显的地区为朝阳区与通州区交界处、大兴区南部、海淀区北部，沉降速率超过-60 mm/a，最大沉降速率约为-110 mm/a；抬升比较明显的地区为昌平区中西部、门头沟区和海淀区交界处以及房山区中北部，这些地区地表年平均形变速率约为10~20 mm/a。该结果与文献[9]中结论大致相同。

将昌平站(CHPN)、牛栏山站(NLSH)、东三旗站(DSQI)、朝阳站(CHAO)、西集站(XIJI)5个CORS站在大地高方向的年平均变化率与SBAS-InSAR结果在大地高方向的年平均变化率进行比较，发现位于昌平区中心位置的CHPN站与位于顺义区北部的NLSH站略微抬升，而位于昌平区东南部的DSQI站、朝阳区中心位置的CHAO站和通州区东部的XIJI站发生沉降，且CHAO站沉降最为明显。CORS站数据与SBAS-InSAR结果相吻合，表明数据具有可靠性。

3.2 InSAR优化方法结果

对SBAS-InSAR结果进行优化处理后，选取采样历元进行对比。表 1(单位mm)为2018-07-02监测量进行优化处理前后的数据统计，由表可知，处理后的标准差更小，平均值几乎不变，并剔除了极值，表明InSAR数据去除极浅地表局部变化影响后获得的数据质量更可靠。图 3为该历元下监测量优化前后结果对比，由图可知，处理后的误差点得到明显修复，进一步说明了优化方法的有效性。

3.3 CORS站数据处理结果

表 2(单位mm/a)为2018~2020年5个CORS站在大地高方向的年平均变化率，其中CHPN和NLSH站表现为略微抬升，DSQI、CHAO和XIJI站呈现沉降趋势，并且DSQI站沉降最为显著，3 a间累积沉降量为87 mm。图 4为北京地区5个CORS站大地高周变化原始数据与其线性变化和非线性项重构结果。

3.4 融合CORS网与时序InSAR监测结果分析

以CORS数据为基准，对InSAR监测量进行间接平差，计算SAR影像覆盖范围内的散点在2018~2020年间年平均形变速率，结果如图 5所示。

由图 5可知，朝阳区与通州区交界处地面沉降最严重，最大沉降速率达到-90 mm/a；昌平区中西部、海淀区西部、门头沟区东部、石景山区、顺义区北部地面抬升比较明显，年平均形变速率约为10~20 mm/a；其他地区沉降变化不明显，沉降速率约为-10~10 mm/a。

对比SBAS-InSAR结果可知，地面沉降与抬升情况一致。选取2018-03-04、2018-11-11和2020-09-13三个采样历元，利用距离反比方法获取参考历元为2019-04-21 00:00的SBAS-InSAR监测量以及数据融合后的监测量在昌平站等5个CORS站处的形变量，与CORS站大地高周变化时序相比并分别作差，结果见表 3(单位mm)。由表可知，融合CORS网与时序InSAR的监测量更接近CORS站的监测量。

4 结语

本文基于SBAS-InSAR处理获取的地表形变时间序列和CORS站大地高周变化时间序列，经过优化处理后，以CORS网为基准对InSAR数据进行间接平差，得到融合后的地面形变时间序列。对北京地区2018~2020年地面沉降进行分析，得到以下结论：

1) 通过去除InSAR监测时间序列中非动力学信号形变，可以得到更加可靠且稳定的地面形变数据。

2) 融合CORS网与时序InSAR进行监测整体可行，可将InSAR技术的高空间分辨率特点与CORS网的高精度特点结合起来，实现高分辨率、高精度的地面形变监测。融合CORS网与时序InSAR在监测方面具有很大优势。

3) 融合后的地面沉降监测结果显示，朝阳区、通州区、大兴区中南部地面沉降最为明显，最大沉降速率达到-90 mm/a；昌平区中西部、海淀区西部、门头沟区东部、石景山区、顺义区北部存在较为明显的抬升，年平均形变速率约为10~20 mm/a；其他地区地面形变相对较小，年平均形变速率约为-10~10 mm/a。