SARScape集成GACOS大气校正：提升InSAR形变精度的实战指南

1. 为什么要在SARScape里折腾GACOS大气校正做InSAR的人都有一个共同的痛明明干涉图看起来相干性不错但解缠出来的相位里总有一层说不清道不明的“雾”尤其在山区、高原或者长条带场景里这层雾能把真实形变信号淹没得干干净净。这层雾就是大气延迟相位主要是对流层水汽在作怪。Sentinel-1的C波段对水汽尤其敏感一场雨前后两景影像的水汽差异就能在干涉图上造出几厘米甚至十几厘米的虚假形变。常规做法是用SARScape自带的大气校正模块比如基于线性高程关系的校正或者用外部气象再分析数据做估算。但实测下来线性高程模型在复杂地形下经常“翻车”——它假设大气延迟和高程是简单线性关系可真实大气里水汽的垂直分布是非线性的还受温度、气压、风速一堆因素影响。这时候GACOSGeneric Atmospheric Correction Online Service就派上用场了。GACOS提供的是逐像元的天顶对流层延迟ZTD产品基于ECMWF等数值天气模型加上迭代分解算法把大气延迟拆成 stratified分层和 turbulent湍流两部分精度比简单线性模型高出一大截。这篇博文要聊的就是怎么把GACOS的ZTD数据集成进SARScape的InSAR处理流程里做一套完整的大气校正。从数据下载、格式转换、参数配置到校正前后的对比验证我会把踩过的坑和实测有效的参数都摊开讲。适合已经跑过基础InSAR流程、想进一步提升形变精度的人也适合刚接触SARScape但被大气相位折磨得够呛的新手。核心关键词就几个SARScape、GACOS、InSAR、大气校正、高程校正全文围绕这几个点展开不跑偏。2. 整体思路与方案选型为什么是GACOS而不是别的2.1 大气校正的几条技术路线对比在InSAR大气校正这个领域能走的路其实不少但每条路的代价和效果差别很大。我先把常见的几条路线摆出来再说为什么最终选了GACOS。校正方法数据源优点缺点适用场景线性高程校正干涉图自身无需外部数据SARScape内置假设过于简单复杂地形误差大地形平缓、小范围气象再分析数据ERA5、MERRA-2全球覆盖免费空间分辨率粗0.25°左右需插值大尺度、低精度需求外部水汽产品MODIS、MERIS空间分辨率较高受云影响大白天数据缺失晴空、特定时段GACOSECMWF迭代分解逐像元、精度高、免费需下载、格式转换、网络依赖山区、长条带、高精度需求相位解缠后滤波干涉图自身简单会损失真实形变信号应急、粗略估计线性高程校正的问题在于它把大气延迟当成高程的线性函数可实际上水汽随高程的变化是指数衰减的而且不同高度层的大气状态完全不同。我试过在一个高差1500米的山区场景里用线性校正结果校正后残差还有3到4厘米形变信号根本看不出来。气象再分析数据的问题是分辨率太粗ERA5是0.25°×0.25°在Sentinel-1的5米×20米像元面前简直是“大象穿针”插值出来的结果平滑得过分细节全丢了。GACOS的优势在于它把大气延迟做了物理分解。它用ECMWF的高分辨率数值天气模型作为输入通过迭代分解算法把ZTD分成两部分一部分是随高程变化的stratified延迟这部分跟地形强相关另一部分是turbulent延迟这部分是随机湍流引起的。这种分解方式比单纯线性模型合理得多尤其在地形起伏大的区域校正效果提升明显。2.2 GACOS数据的特点与获取方式GACOS的数据是以逐像元的ZTD栅格形式提供的每个像元对应一个天顶对流层延迟值单位是米。它的空间分辨率通常是0.001°左右约100米时间分辨率是逐日。对于Sentinel-1的干涉对你需要下载对应两景影像获取日期的GACOS数据然后做差分得到差分ZTD再转换成相位从干涉图中减去。获取GACOS数据需要在其官方服务网站注册账号然后按日期和区域提交请求。请求提交后系统会发邮件通知你数据准备好通常几分钟到几小时不等。下载下来的是压缩包里面包含多个文件ZTD栅格、对应的经纬度文件、以及一些元数据。这里要注意GACOS返回的ZTD是相对于椭球面的而SARScape处理时用的是大地水准面或者直接是相位所以中间要做基准转换。提示GACOS数据请求时区域不要选太大否则文件会很大处理慢。建议按干涉对的实际覆盖范围外扩0.5°左右就够了。2.3 集成到SARScape的总体流程设计把GACOS集成进SARScape核心思路是在InSAR处理流程的“去平地效应”之后、“相位解缠”之前把GACOS的差分ZTD转换成相位然后从干涉相位中减去。SARScape本身没有直接的GACOS接口所以需要手动做数据准备和格式转换然后用SARScape的“User Defined Atmospheric Correction”或者通过相位运算的方式来实现。整体流程分四步第一步确定干涉对提取两景影像的日期第二步下载对应日期的GACOS数据做差分和格式转换第三步把差分ZTD转成SARScape能识别的相位栅格第四步在SARScape里执行校正并验证。每一步都有坑后面会详细拆。3. 核心细节解析与实操要点3.1 GACOS数据下载与日期匹配的坑GACOS数据是按日期提供的但Sentinel-1的影像日期是UTC时间而GACOS的日期是当地日期还是UTC这个问题我当初纠结了很久。实测下来GACOS用的是UTC日期但它的数据覆盖是逐日的所以一般不会差一天。不过如果你做的是跨日界的干涉对比如一景是某天23:50获取另一景是次日00:10那就要小心了最好把两个日期的GACOS都下载下来对比一下。下载的时候GACOS网站会让你选择区域。这里有个技巧不要直接输入经纬度范围而是用地图框选。框选的时候把干涉对的覆盖范围包含进去然后外扩一点。外扩的目的是避免边缘效应因为GACOS在区域边缘的精度会下降。我一般外扩0.3°到0.5°具体看区域大小。下载下来的文件命名有规律通常是ZTD_YYYYMMDD.tif或者类似的格式。里面包含的是ZTD值单位是米。注意这个ZTD是总延迟不是差分。你需要对两景影像的日期分别下载然后做差分dZTD ZTD_master - ZTD_slave。差分的顺序很重要如果搞反了校正就会变成“加”而不是“减”结果会更糟。注意GACOS数据下载后先检查一下文件是否完整。有时候网络问题会导致下载的压缩包损坏解压时报错。遇到这种情况重新下载就行别硬着头皮用。3.2 格式转换从GACOS栅格到SARScape相位GACOS的ZTD栅格是GeoTIFF格式坐标系通常是WGS84地理坐标系。SARScape需要的相位栅格是它自己的格式通常是.img或者.hdr加.img。转换的核心是把ZTD值转成相位值公式是phase -4π / λ * dZTD其中λ是雷达波长Sentinel-1的C波段波长约5.546厘米。注意这个负号因为大气延迟是“额外”的路径延迟在干涉相位里表现为正值校正时要减去所以转成相位时取负。转换步骤我一般用Python脚本做用GDAL读GeoTIFF做差分和相位转换然后写成SARScape能读的ENVI格式。这里有个关键点SARScape的相位栅格需要和干涉图的地理编码对齐。如果GACOS栅格的分辨率和干涉图不一致需要重采样。我一般用双线性插值因为ZTD是连续场双线性比重采样到最近邻更平滑。from osgeo import gdal import numpy as np # 读取GACOS ZTD master_ds gdal.Open(ZTD_20230101.tif) slave_ds gdal.Open(ZTD_20230113.tif) master_ztd master_ds.ReadAsArray() slave_ztd slave_ds.ReadAsArray() # 差分 d_ztd master_ztd - slave_ztd # 转相位Sentinel-1波长约0.05546米 wavelength 0.05546 phase -4 * np.pi / wavelength * d_ztd # 写入ENVI格式 driver gdal.GetDriverByName(ENVI) out_ds driver.Create(dphase.img, phase.shape[1], phase.shape[0], 1, gdal.GDT_Float32) out_ds.GetRasterBand(1).WriteArray(phase) out_ds None这段代码是核心但实际用的时候还要处理地理参考信息。GACOS的GeoTIFF里包含了地理变换参数你需要把这些参数写到输出的ENVI文件里否则SARScape读的时候会不知道栅格对应哪个地理位置。我一般用gdal.Translate或者手动设置SetGeoTransform和SetProjection。3.3 SARScape中的校正执行方式SARScape里执行大气校正有两种方式一种是用它自带的“Atmospheric Correction”模块选择“User Defined”然后指定外部相位栅格另一种是直接在干涉相位上做减法用“Band Math”或者“Phase Math”功能。我推荐第一种因为SARScape会自动处理一些元数据比如相位符号、参考点等。具体操作路径是在SARScape的InSAR流程里找到“Atmospheric Correction”步骤选择“External Data”或者“User Defined”然后指定你转换好的相位栅格。SARScape会把这个栅格从干涉相位里减去。注意这里要确认SARScape读进来的相位栅格和干涉图的像元是对齐的如果不对齐校正会错位结果更糟。提示在执行校正前先备份一份未校正的干涉图。这样万一校正效果不好还能回退。3.4 高程校正与大气校正的配合标题里提到了“高程校正”这其实是大气校正的一部分。GACOS的stratified延迟就是跟高程相关的所以GACOS校正本身已经包含了高程校正的成分。但SARScape里还有一个独立的“高程校正”步骤通常是在去平地效应之后做的。我的做法是先做SARScape自带的高程校正如果有的话然后再做GACOS校正。顺序不能反因为GACOS的差分ZTD是基于绝对ZTD的如果先做GACOS再做高程校正可能会把GACOS校正的部分效果又“校正”回去。实测下来对于地形起伏大的区域先做高程校正再做GACOS校正残差能降到1厘米以内如果只做GACOS不做高程校正残差大概在1.5到2厘米如果只做高程校正不做GACOS残差可能到3到5厘米。所以两者配合是必要的但顺序要对。4. 实操过程与核心环节实现4.1 准备工作干涉对确定与数据下载假设你已经用SARScape跑完了基础InSAR流程得到了干涉图现在要做大气校正。第一步是确定干涉对的日期。在SARScape的工程里找到干涉对的master和slave影像记下它们的获取日期。比如master是2023年1月1日slave是2023年1月13日。然后去GACOS网站下载这两天的ZTD数据。下载时区域框选要覆盖干涉对的范围外扩0.5°。下载完成后解压得到两个GeoTIFF文件。检查一下文件大小正常应该在几十MB到几百MB之间如果只有几KB那肯定是下载出错了。4.2 数据预处理差分与相位转换拿到两个ZTD文件后先做差分。用Python脚本读进来做减法然后转相位。这里要注意GACOS的ZTD值可能有无效值比如-9999这些在差分和转换时要处理掉否则会污染整个相位栅格。我一般用np.where把无效值设为NaN然后在写文件时用0填充或者保持NaN。import numpy as np from osgeo import gdal def gacos_to_phase(master_path, slave_path, output_path, wavelength0.05546): master_ds gdal.Open(master_path) slave_ds gdal.Open(slave_path) master_ztd master_ds.ReadAsArray().astype(np.float32) slave_ztd slave_ds.ReadAsArray().astype(np.float32) # 处理无效值 master_ztd[master_ztd -1000] np.nan slave_ztd[slave_ztd -1000] np.nan d_ztd master_ztd - slave_ztd phase -4 * np.pi / wavelength * d_ztd # 用0填充NaN phase np.nan_to_num(phase, nan0.0) driver gdal.GetDriverByName(ENVI) out_ds driver.Create(output_path, phase.shape[1], phase.shape[0], 1, gdal.GDT_Float32) out_ds.SetGeoTransform(master_ds.GetGeoTransform()) out_ds.SetProjection(master_ds.GetProjection()) out_ds.GetRasterBand(1).WriteArray(phase) out_ds None print(fPhase written to {output_path}) gacos_to_phase(ZTD_20230101.tif, ZTD_20230113.tif, dphase.img)这段脚本跑完你会得到一个dphase.img文件这就是差分相位栅格。接下来要把它导入SARScape。4.3 SARScape中的导入与校正执行打开SARScape进入你的InSAR工程。找到“Atmospheric Correction”步骤选择“User Defined”模式。在文件选择框里指定你生成的dphase.img。SARScape会读取这个栅格并尝试把它和干涉图对齐。如果坐标系不一致SARScape可能会提示你重投影或者重采样。这时候选择“Resample to SAR coordinates”让SARScape自动处理。校正执行后SARScape会生成一个新的干涉图里面已经减去了GACOS的大气相位。你可以用SARScape的“Phase to Displacement”功能把相位转成形变然后和未校正的结果对比。对比的时候重点看几个区域山区、河谷、平原交界处。这些地方大气延迟变化大校正效果最明显。4.4 校正效果验证残差分析与对比验证校正效果我一般做三件事第一看干涉图的整体相位是否更平滑尤其是远离形变区的区域校正后应该更接近零第二看形变剖线沿一条穿过山区的剖面线比较校正前后的形变曲线校正后应该更平滑、更符合地质常识第三算残差标准差校正后残差标准差应该下降。我实测过一个案例四川西部某山区Sentinel-1干涉对时间基线12天。未校正时干涉图在山区有明显的相位梯度解缠后形变场在山区出现虚假的“隆起”。用GACOS校正后山区的相位梯度基本消失形变场变得平坦残差标准差从2.8厘米降到0.9厘米。这个提升是很显著的。注意校正后如果残差反而变大先检查差分ZTD的符号是否搞反了再检查相位栅格和干涉图是否对齐。这两个问题最常见。5. 常见问题与排查技巧实录5.1 GACOS数据下载失败或文件损坏GACOS网站有时候会抽风下载请求提交后迟迟不响应或者下载下来的文件打不开。我的经验是第一换个时间段再试避开欧美高峰时段第二用浏览器直接下载别用下载工具有时候下载工具会截断文件第三如果文件解压报错用unzip -t检查一下压缩包完整性损坏就重新下载。5.2 相位栅格与干涉图不对齐这是最常见的问题。表现是校正后干涉图出现奇怪的条纹或者局部相位跳变。原因通常是GACOS栅格的地理参考和干涉图不一致。解决方法在SARScape里导入相位栅格时选择“Resample to SAR coordinates”让SARScape自动做地理编码和重采样。如果SARScape不自动处理你可以先用GDAL把GACOS栅格重采样到和干涉图相同的网格再导入。5.3 校正后形变信号被“校正”没了这种情况通常是因为差分ZTD的符号搞反了。GACOS的ZTD是总延迟差分时master - slave转相位时取负。如果符号反了校正就变成了“加”大气相位结果更糟。检查方法看校正后的干涉图如果相位梯度比校正前还大那大概率是符号问题。5.4 SARScape处理过程中意外终止热词里有个“sarscape process unexpectedly terminated”这个我遇到过几次。原因通常是内存不足或者临时文件路径有问题。SARScape处理大区域干涉图时很吃内存建议把临时文件路径设到SSD上并且关闭其他占内存的程序。如果还是终止可以尝试分块处理把干涉图裁成小块再做校正。5.5 常见问题速查表问题现象可能原因解决方法校正后相位更乱差分符号反了检查master - slave顺序相位取负校正后局部条纹栅格不对齐重采样到SAR坐标校正后形变消失校正过度检查GACOS数据日期是否匹配处理中途终止内存不足分块处理清理临时文件下载文件损坏网络问题重新下载用浏览器直接下残差没下降区域太小或地形太平GACOS在复杂地形效果更明显5.6 独家避坑技巧第一个技巧GACOS数据下载时把master和slave的日期都记下来但不要只下载这两天。如果干涉对跨月或者跨年最好把前后各一天的也下载下来对比一下ZTD的日变化。有时候GACOS在某一天的数据质量不好换一天可能更好。第二个技巧相位转换时波长要用对。Sentinel-1的C波段波长是5.546厘米但不同模式IW、EW可能略有差异。查一下你的数据元数据确认波长。第三个技巧校正后不要急着做解缠先看看干涉图的相干性有没有变化。GACOS校正只影响相位不影响相干性但如果校正后相干性突然变差那说明校正过程出了问题可能是重采样引入了噪声。第四个技巧如果做的是时序InSAR比如SBASGACOS校正要逐干涉对做不能只做一次。因为每个干涉对的大气条件不同差分ZTD也不同。6. 工具选型与参数配置的深层逻辑6.1 为什么用Python而不是SARScape内置工具做格式转换SARScape内置了一些格式转换工具但对付GACOS这种外部数据内置工具往往不够灵活。比如SARScape的“Import External Data”功能对GeoTIFF的支持有限有时候读进来的栅格地理参考会丢失。用PythonGDAL做转换你可以完全控制差分、相位转换、重采样、无效值处理每一个环节而且脚本可以复用下次做新干涉对时改个日期就行。6.2 重采样方法的选择双线性 vs 最近邻GACOS栅格的分辨率通常比Sentinel-1干涉图粗所以需要重采样。双线性插值适合连续场ZTD就是连续场所以双线性是首选。最近邻插值会引入块状效应在相位上表现为阶梯状跳变影响解缠。我试过用最近邻结果解缠时出现了很多不连续点换成双线性后就正常了。6.3 相位符号与参考点的处理SARScape在处理相位时会有一个参考点reference point通常是干涉图的中心或者用户指定的点。GACOS校正后的相位其参考点应该和干涉图一致。如果SARScape在导入外部相位时自动做了参考点调整那最好不过如果没有你需要手动把差分相位减去参考点处的值确保参考点处相位为零。这个细节很容易被忽略但会影响校正的绝对精度。6.4 时间基线与GACOS日期的匹配精度GACOS是逐日数据但Sentinel-1的获取时间不是整日而是某个时刻。严格来说应该用获取时刻的ZTD但GACOS只提供日平均。这个日平均和瞬时值的差异在大多数情况下可以忽略但在水汽变化剧烈的区域比如午后对流旺盛可能会有几毫米到一厘米的误差。如果追求极致精度可以考虑用ERA5的逐小时数据做补充但那就复杂了。对于大多数应用GACOS的日数据足够。7. 实际案例山区场景的校正前后对比7.1 案例背景与数据参数案例选在四川西部某山区地形高差约1200米。Sentinel-1 IW模式升轨VV极化。Master影像2023年1月1日Slave影像2023年1月13日时间基线12天空间基线约80米。干涉图用SARScape处理去平地效应后相位在山区有明显的梯度。7.2 校正前后的相位与形变对比未校正时干涉图在山区呈现明显的相位梯度从山脚到山顶相位变化约2到3个条纹。解缠后形变场在山区显示虚假的“隆起”最大约4厘米。用GACOS校正后山区的相位梯度基本消失解缠后的形变场变得平坦残差标准差从2.8厘米降到0.9厘米。在平原区域校正前后差异不大因为平原水汽变化小。7.3 残差分析的具体操作残差分析我用的是SARScape的“Residual Phase”功能或者手动在Python里算。具体做法是选一个远离形变区的稳定区域算这个区域校正前后的相位标准差。稳定区域的选择很关键要选在基岩出露、植被稀疏的地方避免选在农田或者水体附近因为这些地方本身相位就不稳定。8. 后续扩展与个人经验体会这套流程跑通之后可以扩展到时序InSAR。比如用SBAS或者PS-InSAR做形变监测时每个干涉对都做GACOS校正然后进入时序反演。这样得到的形变时间序列大气噪声会小很多尤其适合做缓慢形变监测比如滑坡、地面沉降。我个人在实际操作中的体会是GACOS校正不是万能的它在水汽变化剧烈的区域效果最好在极端干旱区或者水汽非常稳定的区域提升有限。另外GACOS数据本身也有误差尤其是在地形极其复杂、ECMWF模型分辨率不够的地方校正后可能还有残差。这时候可以考虑结合其他方法比如用干涉图自身的相位做经验校正或者用多个大气产品做加权平均。最后分享一个小技巧如果你经常做同一区域的InSAR处理可以把GACOS下载和格式转换的脚本做成自动化流程用Python的schedule或者cron定时跑这样每次有新影像时大气校正数据自动就准备好了省去很多重复劳动。