1. 为什么我要折腾GMTSAR处理陆探一号数据陆探一号LT-1上天之后做InSAR的人手里终于多了一套自主可控的L波段数据源。L波段波长大约24厘米比C波段哨兵一号约5.6厘米长出一大截穿透植被的能力强在山区、林地、矿区这些地方做形变监测相干性保持得明显更好。但问题也随之而来主流的开源InSAR处理链比如GMTSAR、ISCE、SNAP早期都是围绕哨兵一号、ALOS这些数据设计的LT-1的元数据格式、轨道文件、参数文件跟它们并不完全兼容直接拿现成流程去跑十有八九卡在第一步读数据上。我自己前前后后花了大概两周时间把GMTSAR处理LT-1的整条链路跑通中间踩的坑不算少——从元数据解析报错到轨道改正对不上再到解缠相位大面积跳变基本每个环节都交过学费。这篇东西就是把这套流程完整记录下来包括每一步为什么这么做、参数怎么定、报错怎么查。适合已经对InSAR有基本概念、手里有LT-1数据、想用GMTSAR做时序或差分处理的同行参考。如果你连干涉相位、去平地效应这些词都还没搞明白建议先补一下基础不然中间很多操作你会知其然不知其所以然。GMTSAR这套工具链的特点是“脚本驱动GMT绘图”核心处理靠C程序流程编排靠shell脚本输出结果用GMT画图。它不像SNAP那样有图形界面点点点但胜在透明、可改、可批量尤其适合做长时序、大批量的自动化处理。LT-1数据本身是国产格式官方提供的数据包里有SAR影像、精密轨道、参数文件我们要做的就是把这些“原料”翻译成GMTSAR能认的“语言”然后按标准流程走一遍。下面我按“整体设计思路—核心细节—实操过程—问题排查”四块来讲尽量把每个参数背后的道理说清楚而不是只丢一堆命令让你抄。2. 整体方案设计与思路拆解2.1 为什么选GMTSAR而不是SNAP或ISCE先说选型。做LT-1的InSAR能用的开源工具其实就那几个SNAP、ISCE2、GMTSAR。SNAP图形化好上手但它的LT-1支持主要靠社区补丁批量处理时内存占用大跑长时序容易崩ISCE2功能全但依赖环境复杂conda装一遍能折腾半天而且它的LT-1 reader更新不算及时。GMTSAR相对轻量核心依赖就是GMT、NETCDF、FFTW这几个编译一次能用很久脚本改起来也直观。更关键的一点GMTSAR的处理逻辑是“先把所有数据统一成它自己的中间格式再走标准流程”。这意味着只要我们把LT-1的元数据正确翻译成GMTSAR的格式后面的干涉、滤波、解缠、地理编码就完全复用了它成熟的代码不用为每种新数据重写算法。这个设计思路对我们这种“数据源新、算法想复用”的场景特别友好。提示GMTSAR对LT-1没有官方支持本文所有适配都是基于数据格式的手动对接属于社区常见做法不是官方推荐流程。跑之前建议先用一小景数据验证别一上来就批量。2.2 LT-1数据长什么样和哨兵一号差在哪要对接先得搞清楚LT-1数据包的结构。一景LT-1的SLC数据通常包含单视复数影像SLC一般是TIFF或GeoTIFF格式分实部虚部或者直接复数存储参数文件XML或文本记录成像时间、中心频率、采样率、脉冲重复频率、轨道状态矢量等精密轨道文件POD用于轨道改正辅助的定标、几何参数文件。和哨兵一号比差异主要在三个地方。第一元数据格式不同哨兵是标准的XML注解LT-1是另一套字段命名GMTSAR自带的make_s1a_tops之类脚本读不了。第二轨道文件格式不同哨兵用EOF格式的精密轨道LT-1用自己的POD格式需要转换。第三LT-1是条带模式Stripmap为主而哨兵一号是TOPS模式两者的成像几何和相位特性不一样GMTSAR里对应的处理分支也不同——LT-1走的是make_raw加条带处理那条线不是TOPS那条线。这个区别很关键。很多人拿哨兵的流程直接套LT-1第一步就错了因为TOPS的burst结构、方位向频谱处理跟条带完全不是一回事。LT-1按条带处理流程反而更接近早期的ALOS PALSAR简单一些。2.3 整体处理链路的设计我把整条链路拆成六个阶段每个阶段解决一个明确问题数据准备与格式转换把LT-1的SLC和元数据翻译成GMTSAR的中间格式主要是.PRM参数文件和原始影像影像配准主从影像对齐估计方位向和距离向偏移干涉图生成共轭相乘得到干涉相位去平地效应、去地形相位滤波与相干性估计压制噪声计算相干图判断质量相位解缠把缠绕相位还原成连续相位地理编码与形变解算把结果投影到地理坐标转成形变量。这个顺序是GMTSAR的标准流程我们做的所有适配工作80%集中在第一阶段。第一阶段打通了后面基本就是调参数的事。注意LT-1的轨道改正必须在配准之前做或者至少在干涉之前做。轨道误差会直接体现为相位斜坡如果放到后面处理解缠会大面积出错返工成本很高。3. 核心细节解析与实操要点3.1 元数据解析把LT-1参数翻译成PRM文件GMTSAR每个影像对应一个.PRM文件里面是一堆键值对记录影像的几何和雷达参数。LT-1的元数据字段和GMTSAR要求的字段不是一一对应的需要手动映射。核心字段包括GMTSAR字段含义LT-1对应来源SC_vel卫星速度轨道状态矢量计算SC_height卫星高度轨道状态矢量计算radar_wavelength雷达波长中心频率换算PRF脉冲重复频率参数文件rng_samp_rate距离向采样率参数文件num_rng_bins距离向像元数影像宽度num_lines方位向行数影像高度near_range近距延迟参数文件rng_pixel_size距离向像元间距采样率换算这里最容易出错的是radar_wavelength和near_range。波长用光速除以中心频率得到LT-1的中心频率在参数文件里有注意单位是Hz算出来波长单位是米。near_range是近距斜距单位是米有些LT-1产品给的是采样延迟需要乘以光速再除以2换算成斜距。我踩过的坑一开始直接把参数文件里的“近距时间”填进near_range结果干涉图整个是错的相位完全对不上。后来才反应过来GMTSAR要的是斜距米不是时间秒。换算公式是near_range near_range_time * c / 2其中c是光速约299792458 m/s。这个换算不做后面全盘皆输。3.2 轨道文件转换精密轨道怎么对接LT-1的精密轨道文件给的是卫星在若干时刻的位置和速度矢量格式是时间XYZ位置XYZ速度。GMTSAR需要的是它自己的轨道格式通常是每景影像对应一个轨道文件记录方位向每一行的卫星位置。转换思路读取LT-1的POD文件按成像时间插值出每一方位向行的卫星位置写成GMTSAR的轨道格式。插值用拉格朗日或者样条都行我一般用三次样条精度够用。这里要注意时间系统LT-1用的是UTC还是GPS时要和影像的成像时间对齐差一秒轨道位置就差好几公里相位斜坡会很明显。提示如果暂时拿不到精密轨道用广播星历也能凑合但形变精度会下降。做科研或者精度要求高的项目一定要用精密轨道。3.3 影像配准的参数选择配准是干涉成败的关键。GMTSAR的配准分两步先粗配准估计大偏移再精配准细化到亚像元。LT-1条带数据的配准相对TOPS简单因为没有burst之间的相位跳变问题。粗配准用esarp或者xcorr我一般用xcorr它对条带数据更稳。关键参数是搜索窗口大小LT-1的轨道精度如果够好初始偏移不会太大窗口设64到128像元就够。精配准用fitoffset多项式阶数一般取2阶数太高会过拟合太低拟合不够。实测下来LT-1的配准精度做到0.1像元以内干涉条纹就很干净了。如果配准残差超过0.5像元条纹会明显模糊相干性掉得厉害。3.4 去平地效应与去地形相位去平地效应flatten是把参考椭球面引起的相位去掉这一步GMTSAR用phasefilt配合轨道自动完成。去地形相位topographic phase removal需要外部DEM常用的是SRTM或AW3D30。这里有个细节LT-1是L波段对地形相位更敏感DEM的精度要求比C波段高。如果DEM有空洞或者误差大去地形之后会残留明显的地形条纹。我的做法是先用SRTM 1弧秒数据如果残留大再换AW3D30。DEM的坐标系要和影像对齐GMTSAR用dem2topo_ra把DEM转成雷达坐标。注意去地形相位之前一定要确认轨道改正已经做好否则地形相位和轨道误差混在一起根本分不清哪个是哪个。4. 实操过程与核心环节实现4.1 环境准备与依赖安装先把环境搭起来。GMTSAR依赖GMT、NETCDF、FFTW、GMTSAR本体。我用的系统是Ubuntu 20.04GMT版本6.3GMTSAR是GitHub上的最新版。# 安装依赖 sudo apt-get install -y gmt gmt-gshhg netcdf-bin libnetcdf-dev libfftw3-dev # 编译GMTSAR git clone https://github.com/gmtsar/gmtsar.git cd gmtsar autoconf ./configure --prefix/usr/local make sudo make install编译完把GMTSAR的bin目录加到PATH里。验证一下which esarp which phasefilt能找到就说明装好了。这里提醒一句GMT版本和GMTSAR版本要匹配GMT 6配新版GMTSARGMT 5配老版混用会报一堆函数找不到的错。4.2 数据准备与PRM文件生成假设你手里有两景LT-1的SLC一景主影像一景从影像还有对应的参数文件和轨道文件。目录结构我一般这样组织project/ raw/ master/ slc.tiff param.xml orbit.pod slave/ slc.tiff param.xml orbit.pod proc/先生成PRM文件。我写了个Python脚本解析LT-1的XML参数文件输出GMTSAR格式的PRM。核心代码逻辑import xml.etree.ElementTree as ET def parse_lt1_param(xml_path): tree ET.parse(xml_path) root tree.getroot() params {} # 读取中心频率、PRF、采样率等 params[radar_wavelength] 299792458.0 / float(root.find(.//center_frequency).text) params[PRF] float(root.find(.//prf).text) params[rng_samp_rate] float(root.find(.//range_sample_rate).text) # ... 其他字段 return params def write_prm(params, out_path): with open(out_path, w) as f: for k, v in params.items(): f.write(f{k} {v}\n)这个脚本要根据你实际拿到的LT-1参数文件字段名调整不同批次的产品字段命名可能略有差异。生成PRM之后把SLC影像转成GMTSAR能读的格式一般是把TIFF转成二进制或者直接用GMT支持的格式。4.3 轨道改正与配准轨道改正我用自己写的脚本读POD文件插值出每行的卫星位置写成GMTSAR轨道格式。然后跑配准# 粗配准 xcorr master.PRM slave.PRM # 精配准 fitoffset.csh 2 2 master.PRM slave.PRMfitoffset.csh的两个参数分别是距离向和方位向的多项式阶数我一般都用2。跑完看offset.out如果残差RMS小于0.2像元说明配准没问题。4.4 干涉图生成与滤波配准好了就可以做干涉# 生成干涉图 esarp master.PRM slave.PRM -topo topo_ra.grd # 滤波 phasefilt intf.grd corr.grd 4 4esarp的-topo选项是去地形相位需要提前用dem2topo_ra把DEM转成雷达坐标。phasefilt的滤波窗口我一般用4x4或者8x8窗口越大越平滑但细节损失也越多。L波段相干性好窗口可以小一点保留更多细节。4.5 相位解缠解缠用snaphuGMTSAR有封装脚本snaphu_interp.csh 0.05 20 corr.grd intf.grd第一个参数是相干性阈值低于这个值的像元不参与解缠第二个参数是解缠的平滑约束。L波段相干性高阈值可以设0.05到0.1比C波段低一些。解缠完检查一下有没有大面积跳变如果有多半是相干性太低或者轨道残差太大。4.6 地理编码与形变转换最后把解缠相位转成形变再地理编码# 相位转形变 phase2d.sh intf.unw master.PRM # 地理编码 geocode.csh master.PRMphase2d.sh会把相位乘以波长除以4π得到视线向形变。地理编码用geocode.csh需要DEM和影像的几何参数。输出是经纬度网格上的形变场可以直接用GMT画图或者导出成GeoTIFF。5. 常见问题与排查技巧实录5.1 元数据解析报错速查报错信息原因解决方法radar_wavelength not foundPRM缺字段检查中心频率换算num_rng_bins mismatch影像宽高填错用gdalinfo确认实际尺寸near_range out of range单位没换算时间乘光速除2PRF invalid字段读错核对XML字段名5.2 干涉图质量差的排查思路干涉图如果全是噪声或者条纹混乱按这个顺序查配准看offset.out残差超过0.5像元就是配准问题轨道看有没有系统性斜坡有的话是轨道改正没做好时间去相干两景间隔太久地表变化大相干性自然低DEM去地形后残留条纹换更高精度DEM。我遇到过一次干涉图半边好半边坏查了半天发现是轨道文件时间戳差了几秒插值出来的轨道位置偏移导致相位斜坡。改时间对齐之后立马正常。5.3 解缠跳变的处理解缠跳变最常见的原因是相干性太低。解决办法提高相干性阈值或者先做多视降低噪声。L波段虽然相干性好但在水体、裸地这些地方还是会失相干。如果跳变集中在特定区域可以在解缠前把这些区域掩掉。另一个原因是解缠参数不合适。snaphu的平滑约束设太小噪声会被当成真实相位设太大真实形变会被抹平。我的经验是先从中间值试看结果再调。提示解缠结果一定要和干涉图、相干图叠着看光看解缠相位很难判断对错。5.4 地理编码错位的排查地理编码错位一般是DEM和影像没对齐。检查dem2topo_ra的输出看地形条纹和干涉条纹是否吻合。如果不吻合多半是DEM坐标系或者分辨率设错了。LT-1的L波段对DEM分辨率要求高建议用30米或更高精度的DEM。6. 我在实操中总结的几条经验第一条LT-1处理最花时间的不是算法是数据格式对接。把元数据解析和轨道转换这两个脚本写扎实后面就是流水线。我建议把这两个脚本做成可配置的换一批数据只改路径不改逻辑。第二条L波段虽然穿透好但对轨道误差和DEM误差也更敏感。同样的轨道残差在C波段可能看不出来在L波段就是明显的斜坡。所以轨道改正和DEM选择要格外上心。第三条GMTSAR的脚本很多是csh写的调试起来不如bash方便。遇到报错先看脚本里调用的C程序输出往往错误信息在C程序那层csh只是把它吞了。第四条做时序分析的话LT-1的数据量比哨兵大单景覆盖范围也大磁盘和内存要提前规划。我一般先做小区域裁剪验证流程通了再上全图。最后分享一个小技巧GMTSAR处理完的中间结果很多建议每步都保留尤其是配准后的偏移文件和轨道文件。一旦后面出问题可以快速定位是哪一步引入的不用从头重跑。这个习惯帮我省了至少一半的返工时间。