基于全球大气模型的日本关东平原InSAR对流层延迟改正研究

焦广棋，孙 玉,杨 琰

(1. 福州大学 空间数据挖掘与信息共享教育部重点实验室， 福建 福州 350108； 2. 地理空间信息技术国家地方联合工程研究中心， 福建 福州 350108；3. 卫星空间信息技术综合应用国家地方联合工程研究中心， 福建 福州 350108)

星载合成孔径雷达(synthetic aperture radar，SAR)信号穿透大气层时受到压强、温度和总电子含量等影响，其传播速率及路径发生改变而导致雷达信号相位延迟，称为大气延迟效应[1]。大气延迟包含电离层延迟和对流层延迟：电离层延迟由自由电子密度分布不均引起，可显著影响极地、高纬地区和长波段(L波)信号，但对C波段信号影响较小；对流层中的温度、湿度和水汽等变化引起对流层延迟，该延迟可影响任何波段和地区，水汽含量丰富的地区影响更大[2]。本研究仅讨论对流层延迟对C波段信号的影响，对电离层延迟未做修正。因此，所使用的Sentinel 1A影像干涉原始相位中包含不同程度的对流层延迟，严重干扰了形变监测及其机制的正确解译，是合成孔径雷达干涉测量技术(interferometric SAR, InSAR)数据处理中主要的误差源之一[3]。因此，估计并改正对流层延迟对基于InSAR技术的高精度形变监测至关重要。

对流层延迟可分为由地形起伏引起的垂直分层延迟和具有随机性的大气对流过程引起的湍流混合延迟。目前，国内外学者已针对对流层延迟改正开展了诸多研究，主要分为3类：①时空滤波法。整个大气延迟被假设为具有时空相关的随机分布性质，时空滤波可以有效抑制湍流混合延迟，但选择合适的滤波器和滤波窗口大小成为难点。时空滤波常被用于时序InSAR技术处理流程中去除对流层延迟，比如永久散射体方法(persistent scatterer InSAR，PSInSAR)、短基线子集方法(small baseline subset，SBAS)[4]。②函数模型。以干涉相位与地形高程之间的关联性建立函数模型，包括线性模型和幂律模型，能够有效去除垂直分层延迟[5]。③外部辅助数据集。利用外部数据集提供气象参数计算改正对流层延迟，但光学遥感数据受限于云、雨天气及夜间环境；全球导航定位卫星系统(global navigation satellites system，GNSS)数据受GNSS站点密度的影响导致改正效果不佳；全球大气模型(global atmospheric models，GAM)是近年来发展的改正对流层延迟的主要数据源，以全球覆盖、数据免费、适用性强等优势在单幅、时序干涉图对流层延迟改正中应用广泛，包括现代卫星应用生产再分析资料(Modern-Era retrospective analysis for research and applications, v2，MERRA2)、欧洲中期天气预报中心(European centre for medium-range weather forecasts，ECMWF)发布的再分析资料(ECMWF reanalysis interim/v5，ERA-Interim/5)[6]。InSAR通用性大气改正在线服务(generic atmospheric correction online service，GACOS)融合了ECMWF数据集和连续全球定位系统(global positioning system，GPS)对流层延迟估计数据，采用迭代对流层分解模型处理后向用户免费提供全球性的对流层延迟数据[7]。

针对GAM以及GACOS改正时序InSAR干涉图对流层延迟的可靠性，李鹏等[1]采用线性模型、GACOS及ERA-Interim改正覆盖阿尔金断裂带地区Sentinel 1A和ASAR影像生成干涉图中的对流层延迟，结果表明两类SAR影像生成的干涉图经GACOS改正后的相位标准差削减量分别可达28.8%和14.8%，当去除轨道平面后可达54.5%和68.1%。Wang等[3]以青藏高原地区的宽幅Sentinel 1A影像为例，对比评估ERA-Interim、ERA5、GACOS及MERRA2对流层延迟改正的效果，结果表明ERA5与GACOS具有相近的改正效果，验证了2018年新发布时间分辨率最高的ERA5应用于青藏高原InSAR对流层延迟改正的巨大潜力。ERA5改正对流层延迟的效果在内陆地区已被证明，但在沿海地势平坦地区相比其他GAM的表现尚需进行对比研究。日本关东平原位于沿海地势平坦、地质构造复杂地区，2011年Mw9.0地震引起地表渗透率改变导致关东平原内出现地表隆升现象，更需进一步开展该地区的地壳形变监测研究[8]。

不同对流层延迟改正数据集的改正效果具有差异性，需要可靠的统计指标去衡量改正效果最优的数据集。目前，最常用的统计指标相位标准差能快速识别最佳对流层延迟改正数据集，但未能从空间尺度的角度去衡量。半变异函数和顾及距离变化的相位标准差不仅可表达对流层延迟改正量的空间关系，还可从空间尺度上评估对流层延迟改正模型的效果[9]。本研究利用覆盖日本关东平原中部2016年12月—2020年8月的57景Sentinel 1A影像获取时序干涉原始相位，采用2种统计指标(顾及距离变化的相位标准差和半变异函数)定性与定量评估2种GAM(ERA5与MERRA2)和GACOS对流层延迟改正的效果。

1 研究区域与数据源

1.1 研究区域概况

日本关东平原西、北部接山地和丘陵，中部包含台地和洼地，东部和南部临太平洋和东京湾，面积约为1.6万 km2，大部分海拔位于100 m以下。该地区属于亚热带季风气候，受海洋作用影响，具有明显的海洋性气候特征，气温年较差较小；沿岸的日本暖流导致年降水偏多，降雨量达1 500 mm以上；夏、秋季多台风，对流层中水汽含量丰富[10]。世界上类似关东平原地形、纬度的地区如中国的华北平原。如图1所示，研究区位于35.03°N～36.67°N，139.09°E～140.40°E(红色矩形所示)，面积约13 300 km2。该区域内中部地形相对开阔，林地分布较少，人口密度大，拥有世界上经济高度发达的东京都市圈。Morishita使用2015—2020年Sentinel 1A影像监测东京周围地表隆升达5 mm/a，在Tochigi地区出现季节性的沉降[11]。为了利用InSAR技术更加精细地反映关东平原地表形变的时空演化规律，需要最大限度减弱或去除对流层延迟对形变结果的干扰。

图1 研究区示意图

1.2 数据简介

“哨兵一号”(Sentinel 1A/B)卫星向公众免费提供高质量的SAR影像，宽幅模式的地面覆盖范围约250 km×160 km，分辨率为5 m×20 m(距离向×方位向)[12]。本研究获取2016年12月8日—2020年8月25日覆盖研究区的Sentinel 1A影像，具体参数信息见表1。欧空局提供与影像日期对应的精密定轨星历数据；美国国家航空航天局(national aeronautics and space administration，NASA)提供30 m分辨率数字高程模型(shuttle radar topography mission)[13]。

表1 Sentinel 1A影像参数信息

本研究采用的3种对流层延迟改正数据集包括：①ERA5是ECMWF于2018年发布的时间分辨率最高的模型，采用了含有卫星数据改正、湿度分析等方法的四维同化技术(4D-VAR)处理后提供的全球气象再分析资料，时间分辨率1 h，空间分辨率0.25°×0.25°，压力层37层[14]；②GACOS是采用迭代对流层分解模型，融合全球连续性GPS的对流层延迟估计数据、地形高程数据与高分辨率的ECMWF产品，提供1979年至今的全球对流层延迟产品，时间分辨率6 h，空间分辨率0.125°×0.125°，压力层137层，是目前最高的模型[7]；③MERRA2是NASA利用大气数据同化系统对卫星、地面观测站、无线电探空及舰船的观测数据进行质量控制、增量分析及同化，提供1980年至今的全球气象再分析资料，时间分辨率6 h，空间分辨率0.5°×0.625°，压力层42层[15]。

2 研究方法

2.1 GAM对流层延迟改正原理

GAM在全球尺度上以均匀分布的网格点形式提供相同时间间隔(1 h或6 h)的气象参数，利用Sentinel 1A卫星在研究区的成像时间提取GAM提供的温度、气压和位势高等气象要素计算对流层延迟。利用大气折射率N来表示对流层延迟，由于云层中水滴引起的液态水延迟对C波段的影响小于0.1～0.4 mm/km，可忽略不计。N由静力延迟和湿延迟组成，可用式(1)表示[16]，利用式(2)对N积分计算得到对流层延迟[7]：

(1)

(2)

式中：Nhydro和Nwet指静力延迟、湿延迟；P指大气压，hPa；T指温度，K；e指水汽分压，hPa；k1、k2、k3系数值分别是0.776、0.716和3.75×103K2Pa-1；θ指雷达信号入射角，°；λ指雷达波长，m；z指地面高程，m；zref指参考高度值15 km，认为超过15 km以上的高度，大气折射率不随高程改变。

根据影像获取时间t，提取GAM中网格点上的气象参数，利用相对湿度计算水汽分压，再用式(3)～(5)计算GAM网格点上的天顶绝对对流层延迟[16]：

(3)

(4)

(5)

式中：P=Pd+e，为总气压；Rd和Rv分别为静力延迟气体和湿延迟气体常数，287.05、461.495 Jkg-1K-1；gm为重力加速度。

由于GAM的分辨率较低，在垂直向上采用三次样条插值及水平向上采用双线性插值得到每景SAR影像中像素点的天顶对流层延迟。由于SAR卫星系统侧视成像，InSAR技术只能监测视线向的形变信息，因此需要利用每个像素的入射角信息将天顶方向的对流层延迟投影到视线向，以主、辅影像中对应像素的对流层延迟相减得到干涉图的对流层延迟，并将对流层延迟改正数从原始干涉图中扣除来完成对流层延迟改正。Sentinel 1A卫星为C波段卫星，波长为5.6 cm，利用相位和距离之间的转换关系，采用式(6)将对流层延迟转换成对流层延迟相位：

(6)

式中，φ为对流层延迟相位，d为对流层延迟距离。

2.2 统计指标

采用不同的数据集改正对流层延迟可能会引入数据集本身的误差，严重影响形变解译的可靠性和精度，合适的评估指标对于识别最佳改正数据集尤为重要。统计指标相位标准差基于全局或局部范围内的干涉原始相位中无较大形变信号的准则来评估对流层延迟改正前后的效果[1]。因此，筛选未进行对流层延迟改正而解算的地表形变速率±0.5 mm/a之内的数据点为评估数据集。本研究采用以下2种评价指标衡量对流层延迟改正前后的效果：

1) 考虑距离变化的相位标准差(standard deviation,STD)。随机选择100个半径从1～100 km的窗口，计算每个相同大小的窗口内未进行对流层延迟改正和每种数据集进行对流层延迟改正后所有数据点的相位标准差，如式(7)所示。随着窗口距离增大会包含更多的数据点，然后计算100个大小相同窗口的相位标准差的平均值，作为识别时序InSAR对流层延迟改正最佳数据集的指标：

(7)

式中，干涉图数量M，干涉原始相位、3种数据集改正对流层延迟后的相位均表示为xi。

2) 半变异函数。该指标表示设定固定的间隔滞后距离h，计算变量在x和x+h处的样本方差值。距离相近的数据点的性质较为相似，半方差值随空间距离的增大而增大，如式(8)所示[9]。为提高计算效率，对前期筛选的数据点进行降采样(1∶100)，设置滞后距离h为1 km，从0～100 km，计算数据点对的半方差值：

(8)

式中，干涉原始相位、3种数据集改正对流层延迟后相位在点x处的值为Z(x)，h为滞后距离，Var[·]为计算采样点方差函数。

2.3 数据处理流程

本研究时序InSAR干涉图的对流层延迟估计与改正的具体技术路线见图2。①基于SNAP软件进行时序差分干涉：选择时间2018年8月12日的影像为主影像，其余56景辅影像分别与主影像进行配准，将配准后的影像进行差分干涉处理得到时序差分干涉图；②识别永久散射体点(persistent scatters，PS)：采用StaMPS软件根据振幅离差指数(0.4)和相干点的相位空间选择长时间内保持稳定雷达散射特性的PS点[17]，得到时序PS点集；③估计与改正对流层延迟：以TRAIN平台计算GACOS、MERRA2和ERA5估计的对流层延迟[18]；④评估最优的对流层延迟改正数据集：采用顾及距离变化的相位标准差和半变异函数评估对流层延迟改正前后的效果，得到最佳改正数据集。

图2 评估3种数据集改正对流层延迟技术路线图

3 结果与讨论

3.1 对流层延迟改正结果

对预处理得到的56幅时序InSAR干涉图进行对流层延迟改正，采用3种数据集改正后给出不同的改正表现，选择2期单幅干涉图进行说明。单幅2018年7月19日—2018年8月12日干涉图时间基线24 d，空间基线-49.4 m，图3(a)给出了3种数据集估计的对流层延迟差分相位与干涉图改正前后的结果对比，对流层延迟差分相位呈明显的空间变化特征，从相模湾、东京湾沿岸向内陆地区(从东南向西北)递增。GACOS模型估计的对流层延迟相位在空间上表现出更平滑且有利于捕捉小尺度变化，但GACOS改正后的结果显示对流层延迟相位仍残留在西北山区地带。图4(a)、4(c)显示在25～100 km的空间尺度上，3种数据集均有明显的对流层延迟改善效果。如表2所示，计算了整幅干涉图GACOS、ERA5及MERRA2改正后的相位标准差分别减少了33.5%、27.9%及25.8%。

单幅2018年8月12日—2019年4月9日干涉图时间基线216 d，空间基线-15.7 m，采用3种数据集估计的对流层延迟差分相位具有明显的差异，其中图3(b)显示MERRA2在相模湾沿岸表现异常。图4(b)、4(d)显示MERRA2在13～100 km的空间尺度上，相比另外2种数据集改正结果出现了较大的偏离，可能与MERRA2的空间分辨率有关。在50～100 km的空间尺度上，GACOS与ERA5表现出较为明显的改善效果。如表2所示，整幅干涉图GACOS、ERA5改正后的相位标准差分别减少了41.7%和12.9%，而MERRA2改正之后相位标准差增加了65.9%。

表2 干涉图对流层延迟改正后结果对比

图3 3种数据集对InSAR干涉图对流层延迟的估计与改正结果

图4 统计指标评估2幅干涉图对流层延迟改正的效果

采用上述2种统计指标评估了56幅时序InSAR干涉图对流层延迟改正前后的效果，并以评估结果的平均值作为总体评价指标。由图5可知2种统计指标的结果基本一致，不同的颜色带表示57景影像生成56幅时序干涉图，在相同的空间尺度上，每幅干涉图经过对流层延迟改正后会对应56个评估值，然后计算56个评估值的标准差，表明每种数据集改正时序干涉图对流层延迟后评估值的离散程度。相位标准差是子区域内对流层改正后相位值变化的总体指标，顾及距离变化的相位标准差是在100个子区域内各自计算相位标准差后再取所有子区域标准差的平均值。半变异函数只判断计算滞后距离内目标点的相位值，该指标能够在离散的空间尺度上有效捕捉改正后相位值的空间变化。在大于30 km的空间尺度上，半变异函数值明显表现出ERA5比GACOS能够更有效改正干涉图中对流层延迟的优势。56幅干涉图对流层延迟改正前后的评估结果平均值显示0～30 km的空间尺度上，3种数据集均不能有效去除对流层延迟的影响，其中在大于30 km的空间尺度上MERRA2比GACOS和ERA5出现了较大的偏离。在30～100 km的空间尺度上，ERA5能够有效改正对流层延迟，但GACOS与原始干涉相位的曲线相近也取得了微弱的改正效果。

图5 统计指标评估56幅InSAR干涉图对流层延迟的平均结果

另外，从定量角度比较了3种数据集改正对流层延迟的效果，表3给出了统计整幅干涉图相位标准差的结果，可看出GACOS、ERA5及MERRA2对流层延迟改正后干涉图的相位标准差数量分别降低了38、36和12幅，对应的干涉图平均降低了21.6%、27.3%及17.6%。ERA5改正后的相位标准差平均值为1.18 cm，相比原始相位标准差平均值(1.43 cm)降低了17.5%，GACOS改正后相比原始值降低了13.3%。因此，相比GACOS及MERRA2，采用ERA5进行估计与改正时序干涉图对流层延迟取得了较好的效果。

表 3 3种数据集时序InSAR对流层延迟改正结果

3.2 地表形变解算结果

根据上述3种数据集进行对流层延迟改正后的评估结果可知ERA5改正表现较好，利用ERA5去除对流层延迟后，采用PSInSAR方法获取了研究区雷达视线向形变速率，如图6(a)所示。研究区内的地表形变机制可分为自然因素和人类活动2个方面：复杂的地质环境主要影响隆升区，发生在东京周围((图6(b))。 由于2011年Mw9.0东北大地震可能导致该地区的地表渗透率增强和孔隙压力改变， 最大的隆起区东南部边界清晰，可能受地下地质边界的影响，这与前人研究结果一致[8]。人类活动的影响主要体现在沉降区域，在古贺-野木町周围(图6(c))和相模湾沿岸(图6(d))发生了地表下沉，古贺-野木町周围拥有易收缩的土壤，可能由于夏季农业地下水使用引起；而相模湾沿岸沉降主要由独特的地质构造背景和围填海工程引起[11]。

(a)形变速率场;(b)隆升区域分布在东京周围;(c)沉降1区域分布在古贺-野木周围;(d)沉降2分布在相模湾沿岸

利用卫星平台的方位角和每个PS点的入射角信息，将GPS的NEU方向的观测值投影到SAR卫星系统的雷达视线向并计算GPS站点的形变速率以验证InSAR监测结果的可靠性。以每个GPS站点为中心，半径100 m内的InSAR观测值进行反距离加权插值得到该站点的InSAR观测值并计算形变速率。以G171站点为参考点计算其他GPS站点的相对形变速率，发现InSAR与GPS的较差小于2.5 mm/yr，具有较高的一致性(见表4)。

表4 InSAR与GPS结果比较

3.3 讨论

日本关东平原地区的时序干涉图中对流层延迟明显，受海洋的影响导致该地区大气对流活动频繁，更容易产生空间小尺度且无规则的湍流混合延迟。研究区内大部分地表高程变化较小，仅在地形起伏剧烈的西北山区地带少量存在由大气厚度变化引起的垂直分层延迟。对流层延迟改正前后的相位标准差均值和削减量反映了不同改正数据集对InSAR对流层延迟的敏感性具有差异，可能是由于不同改正数据集与Sentinel 1A影像的采集时间存在不同的时间差。GAM依赖于地面观测气象资料的采集，关东平原地区毗邻太平洋，易受台风等事件影响引起水汽的快速变化，难以被现有精度的GAM所捕捉。不同时空分辨率、气压分层数的数据集提供的气象参数也存在差异性，也导致了改正效果不同，图7表明了ERA5和MERRA2提供的温度、位势和相对水汽含量在第二压力层的显著差异性。MERRA2的时空分辨率比GACOS、ERA5更低，提供的大气参数难以计算研究区内对流层的真实延迟。GACOS融合了全球连续GPS天顶延迟数据进行增强处理，提高了空间与时间上的改正精度，研究区内GACOS略低于ERA5的改正效果。ERA5在时间上更接近Sentinel 1A影像的获取时间，进而在时序InSAR对流层延迟改正中降低了对流层估计的不确定性。

图7 ERA5与MERRA2在第2压力层的温度、水汽及位势参数比较(白色区域为空值)

4 结论

本研究以覆盖日本关东平原的57景Sentinel 1A影像为例，采用考虑距离变化的相位标准差和半变异函数评估分析了GACOS、ERA5及MERRA2这3种数据集改正时序InSAR干涉图对流层延迟的效果。结果表明，考虑距离变化的相位标准差和半变异函数可用于不同空间尺度对流层延迟改正效果的评估，二者评估结果一致。对于单幅干涉图来说，3种数据集的改正效果具有差异；对于56幅时序InSAR干涉图改正效果综合评估来看，在0～30 km空间尺度上，3种数据集在该研究区内的对流层延迟改正效果均不明显；在30～100 km空间尺度上，ERA5与GACOS具有一定改正效果，MERRA2效果较差。

以日本关东平原作为典型的沿海平原地区，采用InSAR技术在海岸带地区开展大范围的地表形变监测受对流层延迟误差的严重影响，结果显示，时间分辨率最高的ERA5数据集在该地区的对流层延迟改正中具有较好的效果，能够为InSAR技术在国内外沿海平原地区的应用提供参考。