如果你要用光学或雷达影像测量地表形变,COSI-Corr 这个工具迟早会出现在你面前。不管是地震之后想从卫星影像里抠出同震位移场,还是做冰川流速、滑坡蠕动、沙丘迁移,COSI-Corr 的核心就一件事:通过亚像素级互相关,把两张影像之间每个像元的相对位移算出来,再转换成东西向、南北向(或雷达方位向、距离向)的形变结果。它不靠人工选控制点,也不靠特征匹配,而是基于频率域统计相关,所以在植被稀少、纹理稳定的区域能测到 1/10 甚至 1/20 像素量的位移,这是很多普通配准工具做不到的。
这篇内容适合刚接触 COSI-Corr 的研究生、工程师,也适合已经被参数界面劝退过一次、想搞清楚“窗口大小到底怎么填”的同行。我会把安装、数据准备、核心参数、跑相关、后处理、常见报错完整过一遍,尽量按我实际操作的顺序来写,不绕弯子。
1. COSI-Corr 到底是干什么的
1.1 原理一句话讲清
COSI-Corr 全称很长——“Co-registration of Optically Sensed Images and Correlation”,它是加州理工和 JPL 联合开发的图像互相关软件包,后来以 ENVI 插件的形式集成在软件里。它首先不是用来“拼接影像”的,而是用来“测量两张影像之间的变形量”。
它的思路看起来简单:在同一地区拿到两张重访影像,如果期间地表发生了位移,那么同一地物在两幅影像里的位置就会发生偏移。把影像切成很多小窗口,逐个窗口做统计相关,找到每个窗口的最佳匹配位置,就能得到一个密集的偏移场。关键是,它把处理放在频率域里,相当于对窗口内的图像做相位相关,而不是在原始空间域里逐个像素模板匹配。这种方法在高分辨率光学影像上能获得亚像素级的匹配精度,而且不需要预先布设任何地面控制点。
1.2 它和普通影像配准有什么本质区别
普通影像配准的目标通常是:把两幅影像统一到一个坐标系里,让同名点对齐,然后做融合或变化检测。它求的是“一个整体变换”,比如仿射变换、多项式变换,或者用少量地面控制点拟合一个全局模型。
COSI-Corr 求的是“每个像元局部的位移”,它不假设全局一致。地震断层两侧的位移方向不同,冰川中轴和边缘的速度不同,滑坡体内部也可能有拉裂和挤压,这些都需要高空间分辨率的稠密位移场,而不是一个全局变换参数。COSI-Corr 做的事更像“给每个窗口单独做一次配准,再汇总成一张位移场”。
所以,如果你只是想把两张图严格叠合,不需要用 COSI-Corr;如果你想要的是形变、运动、差异,那它就非常合适。另外,它也常被用来反算两景影像之间是否存在系统性的轨道误差、姿态误差,所以很多人在做 InSAR 之前,也会先用 COSI-Corr 的偏移量结果来校验影像之间的实际偏移情况。
2. 安装和前置条件
2.1 ENVI 版本和 COSI-Corr 版本怎么选
COSI-Corr 不是一个独立运行的软件,它必须挂在 ENVI 上。我接触过的版本里,比较常见的是 COSI-Corr 4.0 及后续补丁版,对应的 ENVI 主要是 ENVI 5.1、5.3、5.5 这些。新版 ENVI(比如 6.0 之后)我建议先看看官方扩展兼容表,不要盲目装最新版,否则菜单可能不显示。
这里说一个我踩过的坑:把 COSI-Corr 目录直接解压到 ENVI 安装目录后,在新版 ENVI 的 Extension Manager 里可能找不到它,原因往往是插件没有注册到正确的位置。最稳妥的办法是,打开 ENVI 后,在菜单栏右侧的“Extensions”里面选“COSI-Corr”,如果没有,就手动把插件路径添加进 ENVI 的扩展搜索路径里。
2.2 装完之后怎么验证
安装完成后,你至少能看到一个独立的 COSI-Corr 菜单,里面会有相关性、正射校正、变形后处理等几个模块。如果菜单看不到,进入 ENVI 经典模式(Classic)里找,很多老版本只在经典模式下加载完整 UI。
验证时不要急着跑大影像。拿一小块 1000×1000 像素左右的测试区域,手动平移几个像素,生成一张模拟“形变图”,跑一次相关。如果输出偏移量的平均值和预设平移基本一致,说明安装正常。这一步能帮你把“软件坏了”和“参数填错”区分开,后面排查问题会省很多时间。
2.3 数据层面的硬门槛
COSI-Corr 对输入数据有三个基本要求,漏掉任何一个,后面都会出问题:
- 两幅影像必须已经做了地理参考或者至少具有相同的投影坐标系,像素大小最好一致。如果不一致,先预重采样。
- 影像必须保留足够的纹理信息。纯水面、大片均匀雪地、浓密云层覆盖区,相关性会直接失败。
- 输入文件最好是 16 位整型或浮点型的单波段灰度图像,多光谱影像可以先选择波段或合成灰度。
这里特别提醒,别把“两幅影像投影一样”理解成“文件名里有坐标系就行”。最好在 ENVI 里同时打开两幅图,右键查看 Map Info,确认投影、像元大小、影像四个角坐标都对齐。像元大小一个是 10 米一个是 10.0001 米的场景我见过好几次,表面看都是 UTM,但重采样网格错开以后,位移结果会沿轨道方向出现系统性条带。
3. 数据准备:决定成败的第一关
3.1 光学影像选择的取舍
光学影像相关适合地表纹理稳定的地区,比如岩石裸地、城市建成区、干旱区的冲积扇。对冰川来说,夏季积雪消融后露出的冰碛和冰面纹理也能提相关性,但要避开云和新生积雪。
选影像时我最看重三点:时间间隔不能太长、太阳高度角差异不能太大、地表覆盖变化尽量小。同样是地震形变研究,同震对用震前一期和震后一期就够了;如果做滑坡年际位移,最好选同一季节的数据,否则植被和阴影变化会淹没真实位移信号。
分辨率方面,COSI-Corr 常用影像包括 SPOT-5/6/7、WorldView 全色、GeoEye、Planet 等。简单经验是:想测厘米级到分米级的位移,原始分辨率最好优于 5 米;如果目标位移是几十米量级,中分辨率影像也能用,但精度会明显打折。
3.2 SAR 影像到底能不能用
很多人以为 COSI-Corr 只支持光学影像,其实它也支持雷达影像的幅度图。SAR 数据的优势是能昼夜工作、不依赖太阳光照,而且重访周期稳定。你可以用单视复数影像的幅度分量,也可以用法哨兵一号的 GRD 强度数据。
对 SAR 幅度图做相关,思路和光学影像一致,但要注意斑点噪声:窗口如果太小,相关结果会被噪声主导。我一般会把窗口尺寸调大,比如初始窗口 256、最终窗口 64 起步,同时多做一些中值滤波。另一个问题是,雷达侧视成像造成叠掩和阴影区域相关值会很低,这些区域应该在后期做掩膜剔除。
3.3 DEM 和坐标系:现场常见坑
COSI-Corr 里有一个正射校正模块,作用是利用 DEM 消除地形引起的投影差。如果你拿到的影像本身已经做过严格正射校正(比如第三方提供的正射产品),可以跳过这一步;但如果是 Level 1 级产品,并且地形起伏大于几百米,不做 DEM 校正会引入虚假位移。
DEM 我用过 SRTM 1 秒、3 秒和 AW3D30,都能用。注意 DEM 的投影最好和影像保持一致,尤其在高纬度地区,投影不匹配会带来明显的南北向误差。做正射校正之前,先把 DEM 裁剪到影像范围,再重采样到影像像元大小。
3.4 转出 COSI-Corr 的标准输入
数据进入 COSI-Corr 之前,通常通过工具菜单转成它自己的格式。转换过程主要做三件事:统一波段、把浮点影像转为 16 位整型并记录最小最大值、给背景值打上特殊标记。这样做是为了后续相关计算能够通过固定零值或背景值来排除无效像元。
转换时建议把两幅影像都处理成同一像元深度和像素大小。如果你有两景影像像元大小不一样,先重采样到较小像元,再转换,这样能保留更多空间细节。注意背景值一般设为一个极端负数,比如 -32768,这样相关窗口碰到背景时会自动判为无效区域,不会参与匹配。很多“位移场边缘出现巨大错误值”的问题,都是因为背景值没有正确设置。
4. 完整操作流程详解
4.1 第一步:定义图像对
打开 COSI-Corr 模块后,第一步通常是创建图像对。界面里会让你选主影像和从影像:
- 主影像:作为参考基准,通常是形变发生前的影像。
- 从影像:需要被匹配到主影像上的影像,通常是形变后的影像。
这里的命名很容易让人混淆,因为“主”不等于“重要”,只是一个坐标参考。最终输出的位移场含义是“从影像相对主影像移动了多少”。如果主从选反了,位移符号就会反向。比如地震研究里,把震前当主影像、震后当从影像,正东西位移表示朝东运动;反过来的话结果符号全部对调。
定义图像对的时候,界面还会要求你确认两幅影像的范围和重叠区。如果重叠区域太小,边缘区域会因为没有有效对应窗口而输出大片空白,这是正常的,关键是确保研究区在重叠范围内。
4.2 第二步:粗配准
正式做亚像素相关之前,最好先做一次粗配准,把两幅影像之间的整体位移粗略估计出来。这个整体位移可能来自轨道差异、基线差异,也可能是真实形变里的“大分量”。如果不做粗配准,后续固定窗口大小的相关很容易在搜索范围内找不到正确匹配,尤其当位移超过窗口尺寸的 1/3 时,结果会直接崩掉。
粗配准可以用人工选点的方式,也可以让软件自动估算一个全局偏移。我的经验是:先让自动方式算一遍,再看一下全局偏移值是否在合理范围。例如两景 WorldView 影像之间常有几个像素到几十个像素的整体偏差,如果算出来是几百上千像元,先检查输入影像是不是出现了反转或错位。
4.3 第三步:相关参数设置全解析
这是 COSI-Corr 最核心、也最容易把人绕晕的界面。以下参数我挨个说清楚:
- 初始窗口(Initial Window):用于估算低频大势,捕捉整体位移和较大变形。一般设 256 或 128。
- 最终窗口(Final Window):用于输出最终亚像素位移场,决定了测量的空间光滑程度。常用 32 或 64,如果研究区纹理细腻,可用 16。
- 步长(Step):位移场采样间距。步长越小,输出像元越多,空间分辨率越高,但计算时间成倍上升,随机噪声也更大。
- 过采样(Oversampling):频率域互相关时为了获得亚像素精度而做的插值倍数,一般取 8 或 16。
- 掩膜阈值或信噪比阈值:用于剔除相关性差的窗口。
还有一个参数叫“搜索半径”或“最大位移”,这直接约束了相关窗口在从影像里允许移动的最大范围。如果粗配准后还有几十个像素的残余偏移,那么这个值可以给大点;如果你有准确的轨道参数和粗配准结果,就尽量收紧,这样能减少误匹配。
4.4 第四步:跑相关,得到东西向与南北向位移
设置完成,点 Run。软件会逐窗口计算两幅影像之间的相位相关峰值,得到每个窗口的偏移量。输出结果通常是一组多波段文件:
- 东西向偏移(East-West offset)
- 南北向偏移(North-South offset)
- 信噪比(SNR)
- 有效像元计数(number of valid pixels)
如果跑的是雷达影像,可能还会输出距离向和方位向偏移。这些偏移量单位是“像元”,不是“米”。要换算成米的位移,需要乘上影像像元大小。比如 UTM 投影下 10 米分辨率的 WorldView 全色影像,某像元东西偏移 0.3,代表实际位移 3 米。
我第一次跑出来的结果往往会很花,像噪点一样,这时候千万别急着做分析,先把信噪比波段打开,看看低信噪比区域是否和植被、水体、云层、阴影重合。如果是,说明相关失败是地表变化造成的,而不是参数问题。
4.5 第五步:信噪比滤波
COSI-Corr 对每个窗口除了输出偏移量,还会输出一个信噪比值,表示相关峰值的尖锐程度和匹配可靠性。信噪比越低,越可能是误匹配或不可靠匹配。
我自己的习惯是,先按信噪比大于 0.85 或 0.9 做一个初步掩膜,把低置信度像元直接置空。然后针对保留的像元,再做一个 3×3 或 5×5 的中值滤波,用来消除个别跳变像元。注意不要先用平滑滤波器再设掩膜,否则无效像元会向有效区域扩散,形成一圈伪影。
5. 参数怎么调才能不瞎跑
5.1 窗口大小的白刃战
窗口大小是 COSI-Corr 里最关键的参数,没有之一。它决定了“在多大一个范围内寻找匹配”。窗口太小,参与统计的像元少,抗噪声能力差;窗口太大,窗口内部可能混入不同方向的位移,结果等于被平均掉了,形变细节丢失。典型例子:断层两侧的位移符号截然相反,如果窗口跨过断层,相关结果会落在“中间值”上,断层边界被明显模糊。
所以实际工作中很少只用同一个窗口,而是用“从大窗口到小窗口”的逐步细化策略。先在大窗口 256 下算一个整体趋势,再用 128、64 逐级向下重新计算,每级以上一级结果作为初始偏移,这样既利用了窗口大的稳健性,又保住了窗口小的分辨率。
5.2 步长的意义
步长控制输出网格的密度。步长等于 1 时,每个原始像元都会输出一个位移值,看起来爽,但计算量巨大;而且由于窗口本身就是滑动相关,相邻像元结果高度相关,并不会带来真正的信息增益。步长等于 8 或 16 时,输出间隔大,后续插值会平滑掉一些细节。
我通常先用步长 4 跑一版整体结果,看形变信号的规模;如果研究区范围小、位移梯度大,再改用步长 2。只有在最终出图或做矢量分析时,再考虑插值到更高密度。这样能大幅缩短试错时间。
5.3 过采样次数为什么不能省
过采样这个参数,说人话就是“在频域里把相关峰多插几个点”。过采样 8 意味着在峰值附近细分到 1/8 像元精度,过采样 16 就是 1/16 像元精度。理论上过采样倍数越高精度越高,但有两个问题:一是计算量快速上升,二是其他误差源(如大气折光、地形校正残余)会限制真实精度。当 SNR 本身很低时,过采样再高也很难挽回。
经验值:中低分辨率光学和雷达影像用 8,高分辨率全色影像想追求极限精度可以用 16。更高反而容易出现数值振荡。
5.4 不同场景推荐参数表
| 场景 | 初始窗口 | 最终窗口 | 步长 | 过采样 | 说明 |
|---|---|---|---|---|---|
| 地震同震位移场(高分辨率光学) | 256 | 32~64 | 2~4 | 8 | 需要刷新断层细节时窗口取小值 |
| 冰川流速(中分辨率光学/SAR) | 256 | 64~128 | 8 | 8 | 冰川表面纹理弱,窗口需大些 |
| 滑坡年际位移 | 128 | 32 | 2 | 16 | 位移梯度大,窗口不宜过大 |
| 沙丘或冲积扇迁移 | 128 | 64 | 4 | 8 | 纹理较稳定,兼顾计算效率 |
| SAR 幅度影像 | 256 | 64~128 | 8~16 | 8 | 斑点噪声大,需要大窗口抗噪 |
表格只是一个起点。每个数据集的特性都不一样,拿到新数据后最好先用一个小块区域做参数敏感性测试。如果两个窗口参数都得到基本相同的大尺度位移分布,说明参数稳健;如果结果差异明显,就该去找是哪个环节引入了不确定性。
6. 后处理:把位移场变成能用的图
6.1 扣除整体偏移
即使做了粗配准,最终相关结果里依然可能残留一个全局偏移。原因是影像重采样网格、轨道参数、正射校正模型存在系统误差,这个误差会叠加在真实形变上。处理办法是:在研究区内找一块肯定没发生变形的稳定区域,统计这块区域的位移中值,然后把它当作系统偏移,从整个位移场中减掉。
这里有个细节:稳定区域必须避开断层、滑坡边界、冰川两侧的侧碛、山体阴影等容易产生误匹配的位置。最好是平坦基岩、城市道路等纹理清晰且明确稳定的区域。如果场区没有一个区域是稳定的,那就要依赖外部控制点来做系统偏移估计,否则结果没法信。
6.2 坏值和低相干区处理
COSI-Corr 输出的信噪比波段,是处理坏值的首选依据。除了信噪比阈值,你还要检查位移值的物理合理性。比如一个低海拔河流阶地,位移突然出现几十米,这在大多数地质场景下都不合理,肯定有问题。
我常用的流程是:
- 根据信噪比(比如小于 0.85)生成掩膜,置为 NaN。
- 对位移场做 3×3 中值滤波,孤立跳变点基本能被压掉。
- 看上一步滤波前后的差值,差值大的像元就是突变点,把突变点也掩掉。
- 最后用平滑插值把空白区域补上,但要在成图时用透明度把插值区域标示出来。
这套流程看起来保守,但能避免为了好看而把错误值真的当成有效信号。
6.3 位移矢量导出
位移场本身是东西方向和南北方向两个分量,你可以把这两个分量合成矢量和幅度。在 ENVI 里生成灰度图后,再导出到 GIS 或 MATLAB 里做箭头图。
导出时注意保留地理坐标信息,不要把数据裸存成一个没有投影的二进制文件。我一般会把结果先输出为 GeoTIFF,每个波段带一个坐标文件,然后在 QGIS 或 ArcGIS 里按坐标叠加到影像上,这样期次对比、剖面分析都方便很多。如果要做剖面,建议在垂直于断层的方向等间距取几条线,求沿线的 EW、NS 分量并画成曲线,比逐像元看图直觉得多。
7. 常见问题与排查实录
7.1 位移场全为 0 或者全是 NaN
几乎所有新手都会遇到。先别怀疑数据,大部分原因就三个:输入影像没有地理参考信息、背景值没有正确设置、影像重叠区太小。还有一个很低级但常见的原因:跑相关的时候把主影像和从影像的波段顺序选反了,选了全色以外的某个全零波段。
排查方法:先在 ENVI 里分别查看两幅影像的显示,确认能正常输出灰度;再查看 Map Info 中的投影参数;最后打印一个中间像元的原始值,确保不是背景值。
7.2 位移场出现规律性条纹
条纹通常出现在高分辨率影像里,方向和影像行或列方向平行。主要原因是两幅影像之间有一个亚像素级的整体偏移,但你在输出位移场后没有扣除这个“全局偏移”,于是地表真实形变叠加上一个周期性的插值误差,看起来就像条纹。
这时先不要调窗口参数,而是检查稳定区的位移直方图。如果整体偏移只有 0.3 像元到 0.5 像元,就把它直接从全局位移图里扣掉,再观察条纹是否消失。如果条纹还留在断层两侧,再考虑是不是窗口跨断层造成的振铃效应。
7.3 地形起伏区边缘大面积出错
山区的叠掩、阴影、太阳方位角差异会造成相关失败。严格来说这不是 COSI-Corr 的问题,而是影像几何和光照的问题。处理方法:
- 先用正射校正模块配合 DEM 对影像做地形纠正;
- 把太阳高度角和方位角差异大的影像对换掉;
- 对阴影叠掩区做掩膜,不参与后续反演;
- 如果只有一个影像对,没有替代数据,那么宁可保留空白也不要强行内插。
7.4 内存占用失控或者中断
大范围高分辨率影像,比如两景 WorldView 全色拼接起来达到几万乘几万像素,COSI-Corr 跑起来会非常吃内存。我的建议是分块处理:把研究区按重叠关系切成两到三块,每一块单独跑相关,然后再拼接。分块之间保留至少一个缓冲窗口大小的重叠区,这样拼接后不会出现明显接缝。
此外,关闭其他占用内存的程序,把输出历史参数关掉,能显著降低崩溃概率。如果软件在 Windows 下经常崩溃,试试在 ENVI 经典模式里跑,有些版本打开现代界面就会触发 UI 内存泄漏。
7.5 版本不兼容导致菜单消失
版本问题很烦人。换电脑或升级 ENVI 后,COSI-Corr 菜单不见是重新安装后的常见问题。解决思路是:先卸载插件,重启 ENVI,再关闭 ENVI 重新安装插件目录;如果插件是 dll 或 sav 文件,检查杀毒软件是否隔离了它;另外,有些旧插件不支持新版 ENVI 的“Visualized bands”机制,只能从经典界面进。
8. 我踩过的最深的坑
说一个我印象最深的事,不是参数问题,而是“数据源不一致”问题。有一次我拿到两景光学影像,官网上都写着同一投影,但实际查看时,一景的角点坐标已经做过“轨道精炼”微调,另一景没有。跑完 COSI-Corr 后,形变场上出现了一个横跨整幅图的斜坡状伪影,让我误以为是区域性构造运动,差点写进报告里。后来我用稳定区中位数组进行拟合并减掉了一个一次项平面,才把真实信号从伪影里剥出来。
这件事教给我的习惯是:每次拿到新数据,先同时打开两景影像,在统一显示窗口里叠合,目视检查它们在稳定区域是不是真能对上;然后在跑相关之前就预留几个稳定区作为校验点。做完相关后,先把稳定区位移分布画出来,再判断“哪些位移是信号、哪些是系统误差”。
如果你刚上手 COSI-Corr,我建议你从 64 最终窗口、4 步长、2 倍过采样、信噪比 0.85 阈值这个组合开始,不要一上来追求最大分辨率。跑通之后,再用你的研究区特性去照样板参数多试几组。愿你的位移场清晰、信噪比饱满、误匹配为零。