刚处理完一组某试验区的C波段数据,干涉图出来之后,条纹怎么看怎么不对劲——不是那种干净的地表形变环,而是整个图面上铺了一层细密的“水波纹”。当时第一个怀疑是大气延迟,差分之后又重新滤波,折腾了两个多小时,最后把配准结果单独导出来一看,主副影像之间的偏移量在方位向上差了快0.3个像素,而粗配准报告的误差只有0.1个像素。问题恰恰就出在这个看似不起眼的偏差上。
这就是InSAR处理里最磨人、也最容易被低估的一环:图像配准。很多人觉得配准嘛,就是让两张SAR影像“大致对齐”就行了,反正后面有滤波和相位解缠兜底。但实际做过几轮高分数据就会明白,配准精度如果到不了亚像素级,后面所有步骤都是在浪费计算资源,甚至会把本来有效的形变信号彻底淹没在相位噪声里。
这篇文章就把我对InSAR图像配准从像素级到亚像素级这条技术路线的一些经验和理解整理出来,重点放在为什么配准精度如此敏感、粗配准和精配准各自要做到什么程度、几种主流的亚像素配准方法怎么选、以及实际处理中那些文档里很少写的细节。适用对象是刚进门的学生、想系统踩一遍流程的工程师,也包括那些已经跑了无数单数据但偶尔还是会栽在配准上的老手。
1. 影像没对齐,干涉图就是一团噪声——配准误差如何侵蚀相位
1.1 干涉测量不是在比强度,而是在比相位差
要理解为什么配准精度要求这么变态,首先得搞清楚InSAR干涉图到底在算什么。SAR影像每个像素存的是一个复数,既有幅度信息也有相位信息。干涉处理拿主影像和副影像逐像素共轭相乘,得到的就是每个像素上的相位差。如果同一地面目标在主副影像里没有落在同一个像素上,那这个相位差里混进去的就不只是地表形变、地形、大气这些“正经信号”,还有一大块由错位引入的随机相位噪声。
这里的关键在于,SAR影像是相干成像系统,同一个目标在不同轨道的视角下,散射特性会发生一定变化,但只要我们能把影像对齐到足够高的精度,干涉相位就能保持高相干。反过来,只要错开那么零点几个像素,相位噪声就会迅速抬升。早期文献里有个常被引用的结论:要想把干涉相位噪声压到可接受范围,配准精度通常要优于1/10到1/20个像素。在实际工程中处理5米分辨率的数据时,这个要求换算出来不过零点几米的空间误差,听起来很苛刻,但确实是硬指标。
1.2 相位噪声与配准误差的定量关系
我习惯用一个简单的类比去解释这个误差传递过程。假设两幅影像之间在方位向存在一个配准误差Δ,若该方向上的局部相位变化梯度为k(单位是弧度每像素),那么干涉相位里就会叠加一个约等于k乘以Δ的误差项。局部相位变化越剧烈的地方,配准误差的影响越明显,这就是为什么高分辨率数据里城市边缘、陡峭地形区对配准极度敏感。
定量地看,如果距离向和方位向的相位梯度不同,配准误差就会在两个方向分别产生相位误差项。通常在平地情况下,由于几何关系,距离向的相位梯度相对稳定,而方位向在存在多普勒质心差异时会产生明显的系统性偏移。如果这片区域里恰好有一些垂直向的陡峭结构,比如建筑物密集区,局部相位梯度会变得很大,轻微的配准失配就会直接让相干性从0.8掉到0.5以下。这种局部去相干后续是没法通过滤波弥补的,因为相位已经被污染了,滤波只会把噪声“抹匀”,但并不能还原真实信号。
1.3 像素级、亚像素级和“过分追求”的边界
配准精度是有等级的。像素级配准就是把偏移量定位到整像素,这在早期分辨率较低、主要用于生成相干图或者幅度影像应用时勉强够用。亚像素级配准则要把偏移量定到零点几甚至零点零几个像素,这是干涉测量真正可用的前提。雷达遥感圈子里常说的“1/8像元”“1/16像元”标准,就是围绕亚像素配准质量提出的经验阈值。
但也不代表配准精度“越高越好”到无脑追求。配准精度本身受限于影像信噪比、地物散射稳定性、窗口内是否存在空间变化的形变梯度。在一块均匀裸地或者低相干水体上,哪怕你算法的理论精度能做到1/100像素,实际估计出来的偏移量也会因为测不准而抖动。所以成熟的流程不是盲目堆精度,而是设置合理的质量控制门限,把低置信度区域的偏移量拿掉或者平滑掉。
2. 粗配准阶段:把误差异常控制住的三个环节
亚像素精配准看起来很高级,但前提是粗配准得先给出一个足够好的初值。粗配准如果离真实偏移差了三五个像素,精配准的搜索窗口稍微设小一点就直接跳进局部极值出不来。所以粗配准阶段反而特别考验对全局几何关系的理解。
2.1 基于轨道参数的初值和基于强度的互相关修正
粗配准的第一步通常是用两景影像的轨道状态矢量计算研究区内某个中心位置的偏移量初始估计。星载SAR数据的轨道参数一般给得比较准,这个初值通常能到几十米甚至几米的精度,换算成像素可能还有几十个像素的偏差,不能直接用,但方向已经给对了。
接下来就要做强度互相关粗配准。做法很直观:在主影像上开若干个匹配窗口,到副影像里找一块“看起来最像”的区域,以强度图像的互相关最大作为匹配准则。为了照顾计算效率,第一步通常在大范围内粗搜索,比如每隔一个整像素去算一次相关,找出最大值位置,然后再在最大值附近做小范围精搜索,比如0.2像素步长的插值重采样。这样走下来,粗配准能把两幅影像对齐到一两像素以内,个别区域可能误差更大些,但整体趋势已经可以交给精配准去收拾了。
2.2 匹配窗口的布置远比想象中讲究
粗配准的窗口不是随便开的。如果研究区里有大面积水域、阴影区或者叠掩区,这些地方的强度信号不稳定,互相关峰值往往不明显,匹配结果自然不可靠。实操里我一般会把影像均匀划分成若干子块,每个子块大小为256乘256或512乘512像素,然后在每个子块内部再选取若干小块做匹配。事后要统计每个匹配点的信噪比,也就是相关峰与旁瓣的比值,信噪比过低的点直接剔除,不参与后续偏移场拟合。
还有一点容易被忽略:窗口跨度过大时,如果研究区内存在较大的轨道不平行或者地形起伏,单一偏移量根本描述不了全局变形。粗配准阶段就要能识别这种情况,最好分块处理,每一块单独估计偏移量,再用一个二维多项式去拟合偏移量场,而不是拿一个全局偏移量糊弄过去。
2.3 用多窗口投票机制抑制坏匹配点
粗配准偶尔会碰到个别窗口匹配结果明显异常的情况,比如在郊区种满规则农田的区域,强度图案具有周期性纹理,互相关会出现多个并列峰值,算法一旦选了错误的那个峰值,偏移量就会错开一个或几个周期。这种错误在粗配准里非常坑,因为局部偏差很大,但全局拟合时又容易被当噪声抹掉。
我的习惯做法是:把每个子块的匹配偏移量画成散点图,先看有没有空间突变点,然后以中值滤波为基准,对偏离中值超过一定门限(比如三个像素)的匹配点做重匹配或者剔除。这样处理之后,粗配准的可靠性会高很多,给精配准的初值也就不容易带偏。这个环节多花十分钟,后面精配准能少很多麻烦。
3. 亚像素精配准的三条技术路线与选型逻辑
粗配准完成之后,主副影像之间的误差已经控制在一两个像素范围内,接下来才是真正决定成败的亚像素级精配准,也是整个InSAR流程里算法最密集的部分。工程上用得最多的有三条技术路线:过采样峰值内插、复数互相关最大化、以及相位相关法。三条路各有利弊,实际流程里往往还会组合使用。
3.1 过采样峰值内插:最简单也最稳妥的粗到细策略
过采样峰值内插的思路很朴素:互相关曲面在整数像素位置取到最大值,但如果真实偏移量位于两个整数像素之间,那么相关峰周围的形状会在亚像素位置上有一个略微偏移的峰值。通过对互相关计算结果进行频域补零(也就是过采样),等效于把相关函数在连续域上近似出来,然后再从过采样后的相关表面里找最大值位置。
具体操作时,我通常先把粗配准搜索范围缩小到3乘3或者5乘5像素,在原始分辨率下计算这一小块区域的互相关矩阵。然后对这个矩阵做FFT,补零到原来的16倍长度,再反变换回去。这时相关表面的网格间距已经相当于原像素的1/16,直接取最大值位置就能得到亚像素偏移量。如果想再精细一点,还可以在最大值附近做一个抛物线或者高斯拟合,拟合的峰值位置能进一步推到1/32像素量级。
这条路最大的优点是稳定。只要粗配准结果离真实值不超过一两个像素,互相关峰值形态基本是单峰的,过采样内插几乎不会出什么幺蛾子。缺点则是计算量偏大,如果全图开几百个窗口,每个窗口都做16倍过采样,处理时间会明显上升。工程上一般只对选取的稀疏窗口做高倍数过采样,估计出偏移场后再插值到全影像上的每个像素,而不是对全图每个像素都做这种操作。
3.2 复数互相关最大化:直接对干涉质量负责
如果辅以干涉图的最终质量来评价配准成败,那么复数互相关最大化其实是更“对症”的方法。它不再只看幅度图案的相似程度,而是直接对主影像与副影像的复数值做互相关,计算复数相关系数(也常被称为相干性)。当两幅影像达到最佳配准时,复数互相关的模最大,相位也最集中。
复数互相关和强度互相关有一个本质区别:复数信号的相位包含了地物散射的细微信息,对亚像素偏移的响应要比幅度信号灵敏得多。因此用复数互相关做配准时,即便地形特征不明显,只要地物在两次观测期间的相干性足够高,也能通过相位一致性锁定亚像素偏移。这一点对于草地、裸岩等纹理较弱的区域尤为有价值。
实操里我会在粗配准给出的初值附近,按0.1像素或者0.05像素的步长生成一系列候选偏移,对副影像做对应的相位平移(不需要完整重采样,用频域相位斜坡就能模拟)后,与主影像计算复数相关系数。相关系数最大的偏移即作为该窗口的最佳偏移量。由于每次只做频域相位乘积,计算效率比重复做完整重采样快得多,精度也足以满足常规InSAR流程需求。
3.3 相位相关法:用频域相位差异直接解偏移
相位相关法是图像配准领域的老牌方法,在光学遥感里广泛使用,在SAR复数据上同样适用。它的原理是利用傅里叶变换的位移性质:两幅影像如果存在纯平移,它们在频域里的相位差就是一个线性相位斜坡,斜坡的斜率直接对应空间偏移量。
实际操作中,把主副影像的复数窗口分别做二维FFT,取互功率谱并归一化(也就是只用相位信息,丢弃幅度),再做逆FFT,得到的脉冲函数在偏移位置出现一个尖峰。亚像素精度可以通过对尖峰周围做最小二乘平面拟合来估计相位斜坡斜率,从而解算出精确偏移量,也可以在互功率谱上做加权拟合,抑制噪声频率成分的干扰。
相位相关法的抗噪性在低信噪比区域表现很不错,因为它剥离了幅度信息,只利用相位一致性。但对SAR数据来说有个坑:SAR影像存在固有斑点噪声,虽然相位相关的统计特性仍能保持无偏,但方差会增大。所以实际使用时,我通常会把相位相关法限制在强散射体比较集中的窗口内,比如城区或人造地物密集区,效果会更好。而在均匀地表区域,复数互相关最大化往往更可靠。
3.4 三条路线怎么选:组合比选边更重要
在我自己做过的不少实验里,这三种方法在良好相干区得到的亚像素偏移量通常能互相验证到0.05像素以内,差距不大。真正拉开差距的是在恶劣区域:低相干水体、阴影区域、植被茂密区。这时候复数互相关最大化因为直接使用相位信息,通常表现最好;过采样峰值内插次之;相位相关法在斑点噪声严重时略微吃亏。
但“好”也是要付出代价的。复数互相关要反复生成候选偏移并计算相干系数,计算量比频域相位相关要大不少。所以在实际流程设计里,我喜欢用这样一套组合:先在全图范围用过采样峰值内插快速估计偏移场的大致形态,再在相干性较高、且形变梯度较大的重点区域用复数互相关做局部精化,最后用相位相关法做交叉验证。这样既保证效率,也能守住质量。
4. 决定性参数:窗口尺寸、迭代次数与偏移场的平滑逻辑
方法本身各花入各眼,但在实际工程中,我发现真正影响配准结果的往往是那几个看起来不那么炫酷的参数:匹配窗口多大、要不要迭代、偏移场怎么平滑。这些参数选得合理,简单方法也能用出高端效果;选得不当,再先进的算法也救不回来。
4.1 匹配窗口尺寸:要扛得住噪声,也要留得住细节
窗口的尺寸直接决定偏移估计的分辨率和精度之间的平衡。窗口选太大(比如512乘512像素),内部如果本身存在空间变化的偏移量,会被平均掉,得到的只是一个“平均偏移”,在形变梯度大的区域会造成局部配准误差。窗口选太小(比如32乘32像素),参与统计的独立样本数量不足,相关性估计方差很大,偏移量结果容易剧烈跳动。
我的经验是:在高分辨率星载数据(比如3米分辨率)上,常用的匹配窗口是128乘128到256乘256像素;如果地形平坦、相干性高,可以适当缩小到64乘64像素以获得更细的偏移场空间变化特征;如果是在低相干区,则要增大窗口来换取统计稳定性。但这里有一个平衡技巧:在最终的全图配准重采样阶段,不必每个像素都用独立估计的偏移量,而是用一个平滑后的偏移场插值到全图。这样可以兼顾局部分辨率和整体平稳性。
4.2 由粗到精的迭代设计
亚像素配准极少一步到位,几乎所有商用或开源的InSAR软件都会做多层迭代。第一次迭代用比较大的窗口和比较粗的搜索步长,估计出整体偏移场,然后将副影像按这个偏移场重采样,得到初次校正后的副影像。第二次迭代时,在校正后的副影像与主影像之间重新估计残余偏移量。由于此时的偏移量已经很小,可以用更小的窗口和更细的步长去提取残余的高阶空间变化分量。
这个“由粗到精”的设计,本质上是在做多尺度优化:粗尺度保证算法不会陷入局部极值,细尺度保证最终精度能推高到亚像素。我通常做两到三轮迭代就够用,第一轮窗口256乘256,步长1像素;第二轮窗口128乘128,步长0.2像素;第三轮窗口64乘64,步长0.05像素。这样下来,大多数数据都能稳定收敛到0.1像素以内。
4.3 偏移场空间变化与曲面拟合
配准偏移量在实际数据里往往不是一个常数。轨道不平行、地形起伏、电离层扰动等,都会让偏移量在研究区范围内呈现出一个缓慢变化的空间分布。尤其是地形起伏大的地区,斜距几何会引入明显的距离向偏移梯度,如果只用单一全局偏移量,山区的干涉条纹会出现系统性的残差。
处理办法是用二维多项式对离散窗口估计到的偏移量场做拟合。常见做法是使用二阶或三阶多项式,将距离向和方位向偏移量分别建模为坐标的函数,然后用最小二乘拟合求解多项式系数。拟合之后,不仅可以去除个别异常点的影响,还能得到一个平滑的、适合逐像素插值的偏移场模型。
不过拟合也有风险:如果某个区域存在局部的大气扰动或者电离层扰动导致的非多项式型偏移异常,多项式拟合会把这种局部信号“抹匀”到邻域里,造成局部配准质量下降。这种情况下,我会先做一次残差分析,把拟合残差超过门限的窗口点列出,人为判断是坏点需要剔除,还是真实的空间异常不能动。这种“人机结合”的判断,往往是保证配准质量的最后一道防线。
5. 精度检验远不止看相关系数——用干涉图本身反查配准质量
配准做没做好的直观判断,当然可以看相干性是否提升,干涉条纹是否干净,但这些指标过于整体化,真正出了问题定位不到原因。我习惯用三个维度同时验证配准结果:定量指标、干涉图案形态、以及残差空间分布。
5.1 定量指标:相关系数、相位残差和相干性模型
相关系数是最直观的指标。粗配准前后主副影像的强度互相关从0.3提升到0.8以上,说明大体对齐了。但强度相关系数达到0.9以上就几乎没有区分度了——肉眼看不出来差别,实际亚像素偏差仍然可能存在。这时候需要看复数相干性。
复数相干系数(相干性估计)通常用一定的窗口(如5乘5、9乘9像素)内复数样本计算。在无噪理想情况下,配准精度越高,复相干性越接近理论上限。但实际场景中,时间去相干、体散射去相干等因素会拉低相干性,所以不能只看绝对值,而是要看配准前后相干性是否系统性提升。我更常用的办法是,在低相干区域检查配准后是否出现明显的相位残差条纹沿距离向或者方位向排列。如果有这种系统性条纹,说明残余偏移量在空间上存在常数梯度,直接通过残差条纹的周期可以反推出亚像素配准残余量。
5.2 条纹异常的快速目视检测法
干涉图生成之后,第一步我会先看两件事:一是整体条纹是否平滑连续,二是是否存在沿着距离向或方位向的“平行条纹状”系统性噪声。如果平行条纹出现,高概率是配准残余偏移量未被完全消除。这个特征很典型——真实的形变条纹通常与地表构造或形变场形态相关,很少是纯直线、平行且均匀分布在整个图幅内的。
遇到这种情况,可以在干涉图上沿距离向或者方位向取一条剖面线,看相位是否呈周期性波动。相位剖面的周期如果换算成像素,正好对应残余偏移量的大小,就能验证配准误差来源。我最夸张的一次经历里,剖面线显示出一个约10个像素周期的相位抖动,算下来配准残余偏移量为0.2像素,和后来复数互相关重新配准后的结果对得上。
5.3 残余偏移量评估的数学工具
如果想做更严格的定量评估,可以利用干涉相位与配准误差之间的敏感关系。在失配状态下,干涉相位可以表示为真实相位加上一个由偏移量梯度引起的附加相位。利用这一关系,可以从干涉相位图中估计残余偏移量,但这个过程需要对方位向和距离向的相位梯度做稳健估计,同时去掉地形相位和形变相位的贡献,操作起来比较繁琐。
工程上我更推荐另一种方法:利用二次配准的残差直方图来判断。就是把最终配准后的副影像再次与主影像做一次复数互相关,但这次把窗口设得更小更密,得到的偏移量分布理论上应当围绕零值窄幅分布。如果分布的中心偏离零点,说明当前配准存在系统性残余偏移,其均值就代表未校正的偏差;如果分布的方差过大,说明局部区域仍然存在配准噪声。把这两项指标控制在0.05像素以内,通常就能保证后续干涉解缠结果稳定可靠。
6. 实战中容易栽跟头的几件事
配准流程跑通容易,跑好则需要大量试错。我把这几年积累的几个高频“坑”集中列出来,每一个都曾经让我多熬过好几个晚上。
6.1 忽略频谱偏移预滤波
星载SAR影像在距离向和方位向都存在一定的频谱偏移。距离向的频谱偏移主要由视角差异造成,方位向主要由多普勒质心差异造成。如果两幅影像的公共频谱带宽很小,直接做配准和干涉,相干性会很低,但问题并不在于配准算法本身,而在于没有做公共频谱滤波。
预滤波通常分两步:距离向用基线参数计算频谱偏移量,对方位向先估计多普勒质心差,然后设计带通滤波器让主副影像只保留公共频带。没做这一步就强行配准,偏移量估计结果会偏向同一个方向,而且相干性很难提升到0.5以上。这算是我第一次处理高分辨率星载数据时踩过最深的坑,当时一度怀疑是软件安装出了毛病。
6.2 大形变梯度区域的局部失配
在地震破裂带、滑坡体或者冰川快速运动区,局部形变梯度很大,地表在一个雷达数据获取周期内已经被“拉扯”得面目全非。这种情况下,即便全局配准做得很完美,局部区域还是会失配。这类区域的干涉图通常表现为边界分明的高相位梯度条纹或完全去相干区域。
强行用更大的窗口去“压住”噪声解决不了问题,因为偏移量在窗口内就不是一个常数。我常用的变通手段是,在形变梯度区域缩小匹配窗口并加密网格,让偏移场能表达局部快速变化。如果仍然失配,就只能把这块区域作为低相干区掩膜掉,不参与解缠。强行解缠只会把误差传递到周边区域,导致大范围相位错误连接。
6.3 重采样核的选择直接影响相位精度
确定偏移场之后,需要将副影像重采样到主影像网格上。这里有一个很多人容易忽略的细节:重采样核的函数会对干涉相位引入系统性误差。简单最近邻插值或线性插值会显著降低信号质量,尤其在亚像素偏移量大、内部频率成分复杂的区域,相位误差甚至可以达到0.1弧度量级。
实操中我坚持使用sinc插值核做重采样。标准的sinc核理论上最优,但实际实现时都会加窗截断,常用的有8点或16点sinc核。窗口越长,重采样精度越高,但计算量也越大,对边缘区域容易产生振铃效应。经过多次测试,我认为8点sinc核在精度和计算效率之间是最平衡的选择。如果数据处理量允许,也可以考虑使用更高阶次的插值核,但收益已经很小了。
6.4 不要迷信单一质量指标
最后想说一个心态上的坎。不少流程把“相干性大于某个阈值”作为配准是否合格的唯一标准,这是危险的。相干性会受到时间去相干、地表含水量变化、植被生长等多重因素影响,一个0.5的相干性既可能来自完美的配准加中等时间去相干,也可能来自0.3像素配准误差加极低时间去相干。同样数值背后的质量含义完全不同。
所以我的习惯是同时看三样东西:相干性数值、干涉相位剖面形态、以及配准残差直方图。只有当三者都符合预期时,才敢放心进行下一步相位解缠。如果其中一个指标报警,宁可停下来重新精配准,也不要抱着“先跑通再说”的侥幸心理。数据处理的教训往往就是:你在某一个环节省下的时间,最后会以十倍的成本在后面还回来。
最后分享一个我最近常用的检查习惯:配准完成但不急于生成正式干涉图,先快速生成一张低分辨率的干涉图,把条纹图放大到300%以上,沿着距离向和方位向各扫两三条剖面线。这个操作只需要两分钟,但能非常直观地捕捉到配准残余的系统性条纹。熟练之后,光看干涉图上的条纹形态就能判断出配准方向和大小的问题,比反复看一堆数值指标高效得多。这套方法用了很久,一直是我的最后一道“人工把关”。