合成孔径雷达干涉测量(InSAR)技术因其监测范围大、时空分辨率高和监测周期长等优势^[1]，在地面高程获取和地表形变监测领域得到广泛应用，但由于受对流层延迟效应的严重影响，其在高精度形变监测领域的应用有限^[2-3]。时序InSAR技术采用时空滤波算法，能在很大程度上削弱对流层延迟误差，但仅能消除大气湍流混合延迟效应，无法削弱大气垂直分层延迟效应^[4]。因此，提高时序InSAR技术的大气延迟校正效果成为国内外研究的重点^[5-6]。

目前我国利用地面气象数据校正InSAR对流层延迟的研究较少。基于此，本文对浙江省宁波市鄞州区2018~2020年58景哨兵影像作干涉处理，得到相位解缠后的干涉图，利用传感器获取的地面气象数据和美国国家环境预报中心(NCEP)提供的气象再分析资料，结合大气校正模型进行对流层延迟误差校正，并利用相位标准差STD评估2种数据的校正效果，以分析不同数据对大气延迟校正模型的适用性。

1 研究区概况和实验数据

鄞州区东起北仑港、宁波保税区，西部与余姚接壤，南濒奉化地区，北临长江三角洲，属于宁波市重要经济枢纽。研究区地理范围为121.13°~121.90°E、29.61°~29.95°N，占地面积约为1 381 km²。

本文选取欧空局提供的2018-01~2020-12共58景极化方式为VV的Sentinel-1A影像，条带为宽幅模式，覆盖范围约为250 km×160 km，空间分辨率为20 m×5 m。为更好地进行干涉测量，外部DEM数据选取美国国家航空航天局提供的空间分辨率为30 m的数字高程模型，去除地形相位的影响；选取欧空局提供的精密轨道数据对轨道信息进行修正，去除平地相位的影响。

选用2种数据集削弱对流层延迟：1)传感器每分钟记录的气压、温度等地面气象数据。测站均匀分布在监测区域，相邻两测站间的距离最大为38.47 km、最小为5.31 km；2)NCEP提供的中国区域大气再分析资料，包含气压、温度、水汽混合比等气象参数，水平分辨率为0.5°×0.5°，垂直方向的气压层共分为37层，时间分辨率为6 h，数据集时间为UTC 00:00、06:00、12:00和18:00。

2 数据处理 2.1 SAR数据干涉处理

本文采用SBAS-InSAR技术监测鄞州区地表形变，数据处理流程如图 1所示。以覆盖研究区的58景SAR影像为数据源，对影像进行裁剪。综合考虑多普勒质心频谱差等因素后，选取2019-05-29最优SAR影像为主影像，其余57景SAR数据为副影像，将副影像分别与主影像进行配准。设置时间基线阈值为50 d，空间基线阈值为临界基线的45%，经SAR影像干涉处理后共生成149对干涉对。为更好地进行干涉测量，分别引入精密轨道数据和外部DEM数据以去除平地相位和地形相位的影响。选用Goldstein自适应滤波去除干涉处理时存在的噪声相位。为更好地处理相干性高的区域，采用Delaunay MCF方法进行相位解缠。

由于鄞州区临近海域，对流层空气湿度较大，与湿度相关的大气延迟较为严重，传统的时空滤波算法只能去除大气湍流混合延迟，对于垂直方向的大气延迟并不适用。因此，必须引入外部气象数据并建立大气校正模型来求解天顶方向的对流层总延迟ZTD。由于微波信号的传播路径为雷达视线方向，故需将雷达视线向的大气延迟投影到天顶方向，对解缠后的干涉图引入改正量天顶总延迟ZTD，以削弱大气延迟误差的影响，提高解缠图的质量，从而更好地反演地面年平均形变速率和地表累积形变量。

2.2 大气延迟相位去除

雷达信号在传播过程中会受大气折射的影响，产生路径延迟和路径弯曲，使得InSAR监测过程中形成大气延迟误差。利用每幅干涉图对应的2期ZTD改正量进行计算，得到对应的大气相位，将解缠后的干涉图与之差分相减，以达到减弱大气延迟误差的目的。天顶总延迟ZTD可表示为：

$ \mathrm{ZTD}=\mathrm{ZHD}+\mathrm{ZWD} $

(1)

式中，ZHD为天顶静力学延迟，ZWD为天顶湿延迟。

天顶静力学延迟ZHD可通过Saastamoinen模型求解得到^[7]。Saastamoinen模型是目前常用的对流层延迟改正模型之一，预测精度可达mm级，模型公式如下：

$ \begin{gathered} \text { ZHD }= \\ \frac{0.0022768 \times e_s}{1-0.00266 \times \cos (2 \varphi)-2.8 \times 10^{-7} \times h} \end{gathered} $

(2)

式中，e_s为气象站点处的气压，φ为站点位置的纬度，h为站点大地高。

本文利用气象传感器获取卫星过境时站点记录的气象数据，可直接用于计算。NCEP提供的气象再分析资料属于格网数据，要获取站点附近的大气参数，首先需对周围格网进行搜索，依次查找距离站点最近的4个格网点，然后利用反距离加权(inverse distance weight，IDW)插值获取站点附近的压强P、温度T等大气参数，计算公式如下：

$ \bar{P}=\frac{\sum\limits_{i=1}^4 \frac{P_i}{d_i^2}}{\sum\limits_{i=1}^4\left(\frac{1}{d_i^2}\right)}, i \in Z $

(3)

$ \bar{T}=\frac{\sum\limits_{i=1}^4 \frac{T_i}{d_i^2}}{\sum\limits_{i=1}^4\left(\frac{1}{d_i^2}\right)}, i \in Z $

(4)

式中，d_i为气象站到最近4个格网点的距离，P_i为格网点处的压强，T_i为格网点处的温度。NCEP再分析数据集时间分辨率为6 h，每个气象站点每日只能记录4个压强数据，为获取卫星过境时站点附近的气象数据，采用三次样条插值法(cubic spline interpolation，CSI)获取特定时刻的压强。同理，站点附近的温度也可通过上述方法求得。

NCEP记录的温度为地球表面数据，而记录的气压为平均海平面数据，式(2)中气象站周围的实际气压值e_s可通过下式求得：

$ \begin{gathered} e_{\mathrm{s}}=\operatorname{Pres}_{\mathrm{MSL}} \times \\ \left(1-\frac{0.006\;5 \times h}{\overline{T}-0.006\;5 \times h+273.15}\right)^{5.257} \end{gathered} $

(5)

式中，Pres_MSL为平均海平面气压，h为气象站点大地高，T为气象站点处的温度。

天顶湿延迟ZWD与大气加权平均温度T_m、大气可降水量PWV密切相关，通过一个转换系数可求得天顶方向的大气水汽量，即

$ \mathrm{ZWD}=\mathrm{PWV} \times \frac{R_v \times\left(k_2-k_1 \times \frac{m_v}{m_d}+\frac{k_3}{T_m}\right)}{10^5} $

(6)

式中，R_v为水汽气体常数，取值461 J·(kg·K)^－1；k₁、k₂、k₃均为大气折射常数，分别取值77.6 K·(hPa)^－1、71.98 K·(hPa)^－1、3.754×10⁵ K²·(hPa)^－1；m_v和m_d分别为水蒸气和干空气的摩尔质量，取值18.015 2 g·(mol)^－1、28.964 4 g·(mol)^－1；T_m为大气加权平均温度。

大气加权平均温度T_m是求解ZWD的关键参数。目前使用最广泛的是Bevis提出的实用模型T_m=70.2+0.72T_s，即T_m与地表温度T_s呈线性关系^[8]。但T_m模型受气候与地形等因素的限制，多适用于欧美区域，在中国区域的应用精度较低。本文利用分布在中国地区的84个探空站资料和对应的气象数据，建立研究区的T_m回归经验模型。同时，为更好地表示加权平均温度，在模型中引入水汽压^[9]，计算公式如下：

$ T_m=\alpha_0+\alpha_1 \times T+\alpha_2 \times e $

(7)

式中，T为气象站点处的绝对温度，e为水汽压。利用84个探空站数据估算多项式系数α₀、α₁、α₂，分别取值92.61、0.634、0.279 7。

水汽压e需通过饱和水汽压公式并结合露点温度间接求得，即

$ e=6.11 \times 10\left(\frac{7.5 \times T_d}{237.7+T_d}\right) $

(8)

式中，T_d为露点温度。

结合美国国家环境预报中心提供的气象再分析资料可知，PWV可通过对地面上空各气压层中的可降水量进行数值积分求得，计算公式为^[10]：

$ \begin{gathered} \mathrm{PWV}=\frac{1}{g} \times \int_{P_t}^{P_s} r \mathrm{~d} p=\frac{1}{g} \times \\ \sum\limits_{i=1}^{i=m}\left(\left(\frac{r_i+r_{i+1}}{2}\right) \times\left(P_i-P_{i+1}\right)\right) \times 0.1 \end{gathered} $

(9)

式中，g为研究区的重力加速度，取值9.793 6 m/s²；P_s、P_t分别为对流层顶气压和地面气压；r为水汽混合比；m为气象再分析资料记录的大气层数；P_i和r_i分别为第i层大气压力和水汽混合比。其中，地面到对流层顶部的气压和水汽混合比可从NCEP数据中获取。

3 实验结果与分析

利用2种数据集对处理后的149幅干涉图进行对流层延迟校正，篇幅有限，本文以2019-04-30~2019-03-25干涉对为例进行分析，校正效果如图 2所示。由图 2(a)可见，原始干涉图的相位标准差为1.68，校正前的大气延迟噪声较多；图 2(b)和2(c)的相位标准差分别为1.37和1.33，大气延迟相位均有所降低，说明2种数据集均能有效改善对流层延迟。但利用NCEP气象再分析资料削弱大气噪声的效果差于实测气象数据，对实测气象数据进行大气校正后噪声相位基本消失，效果最佳。

由图 3可知，经大气校正后的大部分干涉相位标准差都有所降低，其中采用NCEP气象再分析资料进行大气校正后相位标准偏差降低的干涉图共计100幅，占67.1%；采用地面气象信息进行大气校正后相位标准偏差降低的干涉图共计106幅，占71.1%。由此说明，地面气象信息具有更好的校正效果。

3.1 精度评估

地面气象站能够对靠近地面的大气层气象要素进行实时观测，具有较高的时空分辨率。为更好地利用气象再分析资料计算天顶总延迟的精度，本文以实测气象数据为真值，以再分析资料提供的数据为观测值，计算各站点处再分析资料求解ZTD的平均偏差bias和均方根误差RMSE，结果如表 1(单位mm)所示。偏差bias和均方根误差RMSE的计算公式为：

$ \text { bias }=\frac{1}{N} \sum\limits_{i=1}^N\left(X_m^{M_i}-X_{m_i}^{R_i}\right) $

(10)

$ \mathrm{RMSE}=\sqrt{\frac{1}{N} \sum\limits_{i=1}^N\left(X_{m_i}^{M_i}-X_{m_i}^{R_i}\right)^2} $

(11)

式中，N为样本数，X_m^R_i为以实测气象数据为代表的真值，X_m^M_i为NCEP再分析资料观测值。

从表 1可以看出，部分站点求解ZTD的bias和RMSE相对较大，最大偏差和均方根误差分别为15.98 mm、16.05 mm。以最大误差站点58569为例，对该气象站点插值处的气象数据进行分析。如图 4(a)和4(b)所示，2种气象数据的一致性相对较好，证明利用反距离加权插值算法和三次样条插值算法获取的气象数据具有可靠性。由图 4(c)可知，部分数据存在较大偏差。由于NCEP数据自身存在精度误差，因此在解算对流层延迟时会出现误差叠加，导致校正效果降低。

为进一步验证InSAR监测结果的准确性，将雷达视线向的地表形变投影至垂直方向，并引入2018~2020年外部水准数据进行对比验证。选取分布在研究区的7个水准点数据，分别与未加大气校正、利用NCEP气象再分析资料校正和利用地面气象信息校正后的监测数据进行对比，结果如表 2(单位mm)所示。从表可以看出，未进行大气校正的监测结果精度较低，与水准测量的监测结果差值较大，RMSE为5.62 mm；利用NCEP气象再分析资料校正后的监测结果精度有所提高，RMSE为3.86 mm；而经地面气象信息校正后的形变值与水准测量所得的监测数据最为接近，RMSE为2.78 mm。综上所述，经过地面气象信息和NCEP气象再分析资料大气校正后的监测精度都有所改善，但利用地面气象信息进行大气校正的效果最好、可靠性最高。

利用SBAS-InSAR技术对实测气象数据大气校正后的地区进行监测分析，获取鄞州区2018~2020年垂直方向的地表形变结果，如图 5所示(负值代表沉降，正值代表抬升)。由图可见，大部分地区的形变都处于一个相对稳定的状态，即未发生较大形变。但仍存在几个明显的沉降漏斗区，本文选取5个沉降较为严重的区域展开分析。

3.2 沉降区域分析

A区域位于云龙镇南部鄞城大道附近，该道路于2019-12正式通车，年平均形变速率约为-6.7~-55.8 mm/a，与甬新河交界处的累积沉降量最大达到165 mm。引起沉降的原因主要为监测期间该段道路处于修建状态，大量工程材料的堆积和浇筑，施工期间对地面压实不够及修建完成后巨大车流量对道路的长期荷载等，而甬新河附近土质结构疏松，多为沙性土质且不易结成土块，所以在河流附近的沉降量最大。

B区域位于首南街道西北部江山万里小区附近，年平均形变速率约为-3.7~-48.2 mm/a，最大累积沉降量达142.3 mm。该区域形成沉降漏斗的主要原因是在形变监测过程中，江山万里小区六期和七期高层建筑物均处于施工阶段，地基的开挖、地下水的大量抽取和重型载货汽车的过往等均会对地面土质产生扰动，且施工区域临近奉化江，土质相对较软，导致地面的沉降现象更加显著。

C区域位于东钱湖镇西南部钱湖景苑小区附近，年平均形变速率约为-18.8~-54.1 mm/a，最大累积沉降量达147.1 mm。实地踏勘发现，该小区南侧宁波轨道交通正在建设施工，判断该地区沉降主要原因为所在地基受轨道施工影响。

D区域位于邱隘镇北部宁波艺术实验学院教育集团明湖校区附近，年平均形变速率约为-2.6~-60.8 mm/a，最大累积沉降量达179.8 mm。该区域属于新开发区域，整个监测期间处于施工阶段，大量工程建设活动对地表有较大影响，引发地面出现沉降漏斗。

E区域位于五香镇南部四都装村，附近有地铁1号线，年平均形变速率约为-5.8~-56.7 mm/a，最大累积沉降量达164.7 mm。沉降原因主要有2点：1)该区域每天都有列车经过，列车循环碾压轨道的动荷载效应累积引发沉降；2)监测期间正值宝涵路修建，压实度不达标，加上车辆长期荷载对地面带来的影响引起地面沉降。该区域四面环河，土质松软，也加大了沉降风险。

4 结语

1) 利用NCEP气象再分析资料校正后相位标准差降低的干涉对占67.1%，而经过地面气象信息校正后相位标准差降低的干涉对占71.1%，说明利用地面气象信息能更好地削弱对流层延迟误差的影响。

2) 以地面气象信息计算得到的ZTD为真值，NCEP气象再分析资料计算得到的ZTD为观测值。结果表明，最大偏差为15.98 mm，平均均方根误差为4.29 mm。

3) 未进行大气校正的监测结果最差，RMSE为5.62 mm；利用NCEP气象再分析资料进行大气校正后的监测结果略有改善，RMSE为3.86 mm；而利用地面气象信息进行大气校正后的监测结果能更真实地反映地面沉降情况，RMSE为2.78 mm。

4) 2018~2020年鄞州区出现不规则地面沉降，选取5个沉降区域进行分析，结果表明，鄞州区地面沉降主要受地下水过度抽取、铁路隧道开挖、高层建筑物施工及重型载货汽车过往等因素影响。