以下文章来源于经纬石旁话遥测 ，作者张永军等
本文改编自学术论文
《高分七号卫星影像融合中的全色-多光谱配准误差补偿模型》
已刊载于《武汉大学学报（信息科学版）》2023年第7期
智能导航与位置服务专刊
张永军 1 王梦欣 1 万一 1 周碧莲 2
丰洁 2 曹炳霞 2 姚永祥 1
1.武汉大学遥感信息工程学院，湖北 武汉，430079
2.广东省国土资源技术中心，广东 广州，510062
张永军
博士，教授，主要从事航空航天摄影测量、实景三维重建、遥感解译与空间知识图谱相关研究。zhangyj@whu.edu.cn
摘 要
高分七号是中国第一颗面向1∶10 000立体测图的遥感卫星，其后视全色影像与多光谱影像通过融合可以得到最高0.65 m分辨率的多光谱融合影像。但在实际生产过程中，发现融合影像存在局部晕边现象。针对这一问题，分析了高分七号卫星影像的误差来源，并提出了全色-多光谱配准误差补偿模型，包含用于补偿线性误差的线性项、用于补偿子线阵电荷耦合器件安置误差和镜头畸变的分段四次多项式以及用于补偿周期性误差的三角级数项。根据该模型设计了求解模式和复用模式两种融合方案。利用中国广东地区和青海地区的高分七号影像数据进行实验验证，结果表明所提模型可以使融合配准精度由最大0.6多光谱像素提升至约0.07多光谱像素，平面误差小于0.25倍多光谱像素的区域占比由不足50%提升到约98%，完全满足融合影像生产要求。
引 用
张永军,王梦欣,万一，等.高分七号卫星影像融合中的全色-多光谱配准误差补偿模型[J].武汉大学学报（信息科学版）,2023,48(7):1029-1038.DOI:10.13203/j.whugis20220714
影像融合是光学遥感影像处理中的重要过程。光学遥感卫星一般以高分辨率全色影像与低分辨率多光谱影像相组合的形式获取数据，在地面处理中进行影像融合，既获得了高分辨率的多波段图像，又节约了宝贵的卫星通信带宽、星上存储空间和信号处理功耗。在影像融合之前，一般需要先进行全色影像与多光谱影像的配准，称为融合配准。这是因为受全色光学传感器与多光谱光学传感器的电荷耦合器件（charged coupled device， CCD）安置精度和镜头畸变、外方位元素误差[6]、卫星平台颤振等因素的影响，全色影像与多光谱影像之间往往存在非线性、不连续的几何错位。如不能进行有效的几何补偿，融合影像往往会出现晕边现象，导致图像锐度下降，严重影响融合影像的可检测性、可分辨性和可量测性。
高分七号是中国首颗民用亚米级高分辨率光学传输型立体测绘卫星，其后视传感器可以获得约0.65 m分辨率的全色影像和约2.6m分辨率的多光谱影像。在实验和生产中，发现部分高分七号卫星影像产品存在全色/多光谱影像配准误差，导致了融合产品的晕边现象。影像融合是遥感领域的重要研究方向，目前常用的融合方法包括主成分变换、高通滤波、小波变换等。这些方法都是在全色影像和多光谱影像准确配准的条件下使用的，因此实际生产中需要首先解决配准问题。目前，常用的商业软件在对一级影像产品（即通过有理多项式系数（rational polynomial coefficient，RPC）表达地理信息的卫星影像）进行融合误差补偿时，一般采用小面元方法。但基于小面元校正的融合配准误差补偿模型参数过多，需要依赖均匀分布的匹配点才能稳定求解，不仅限制了运算效率，还难以应对高分七号后视全色和多光谱相机子线阵CCD错位等情况。
本文以高分七号卫星影像为研究对象，基于多光谱影像和全色影像的成像几何规律，分析融合配准误差的产生机理，提出了适用于整幅影像的全色-多光谱同名像点映射误差补偿，以此进行融合配准误差的定量计算与精确补偿，最大限度地提高融合影像产品质量。利用中国广东地区和青海地区的高分七号一级数据验证了融合配准误差检测与补偿的精度，同时提出了分段项复用方法，为高效率、高精度地生产融合影像提供参考。
1 高分七号融合映射误差分析和补偿模型
为了设计合理的融合映射误差补偿模型，从成像机理的角度分析高分七号卫星成像误差源及其对融合配准精度的影响。
1.1
全色-多光谱影像的双像几何关系
线阵推扫式高分辨率卫星影像的每个扫描行在同一时刻拍摄，物方点与像方点的严密共线关系为：
式中，P是物方点在WGS-84坐标系下的坐标向量；PS是扫描行成像时刻的卫星平台中心在WGS-84坐标系下的坐标向量；λ是尺度参数；
是成像时刻J2000天球坐标系相对WGS-84地球坐标系的转换矩阵；
是卫星本体坐标系与J2000坐标系之间的转换矩阵，一般通过定姿设备（星敏感器和陀螺仪）获取；d是卫星平台中心与相机成像中心之间的偏移线元素；
是卫星相机姿态与卫星本体坐标系之间的转换矩阵，一般由安装角确定，通过实验室标定获取，并通过星上在轨标定优化；p为物方点P对应的像方坐标。
线阵卫星影像融合时需要确定全色影像与多光谱影像的位置映射关系，在产品生产中，一般采用RPC模型进行映射，但为了分析映射误差来源，本文给出了基于严格成像模型的全色-多光谱影像同名像点映射关系：
式中，PS、λ、R1、R2、R3和d为全色影像成像参数；p为全色影像像点的无畸变坐标；
和
为多光谱影像成像参数；p'为多光谱影像上的同名像点的无畸变坐标。调整式（2），获得p与p'的函数关系为：
式中，第二项中的(PS-
)是全色影像成像时刻和多光谱影像成像时刻的卫星位移，误差约1×10－2m水平；λ'是多光谱影像的尺度参数，根据卫星航高，其数量级在1×106水平；(PS-
)/λ' 的误差水平大约在1×10－8 m水平，相对于一般为1×10－6~1×10－5m水平的卫星相机CCD像元大小，(PS-
)/λ' 可以忽略不计，因此可以将第二项近似看作常量；第三项均由卫星传感器安置参数组成，也可以视为常量。因此，误差分析主要针对第一项进行。定义式（3）中第一项的比例因子和旋转因子为：
式中，K是全色-多光谱影像成像尺度参数之比，为常数。经测试，高分七号全色-多光谱相机的主光轴夹角约为0.5°，扫描线夹角约为0.02°，由于对同一位置的成像时间差异不到0.5 s，因此旋转因子R中，地球自转差异极小，
近似为单位阵，卫星姿态变化极小，
近似为单位阵。而由于主光轴夹角和扫描线夹角极小，安置矩阵
近似为单位阵。
综上所述，旋转因子R≈I，近似为一个单位矩阵，在融合配准误差分析过程中，可以忽略R本身，仅考虑其误差的影响。综上所述，式（3）可简化为：
式中，D是式（3）中第二项和第三项之和，为常量。由此，可以得到理想条件下全色像点与多光谱同名像点的空间关系，但当卫星传感器的内外方位元素误差不可忽视时，式（3）中的等式则不再成立，需要通过几何补偿，才能实现全色-多光谱影像的精确配准。
1.2
全色-多光谱同名像点映射误差补偿模型
根据式（5），全色-多光谱同名像点映射误差的来源主要有（见图1）：（1）角元素误差，包括卫星姿态矩阵R2和
的误差，传感器安置矩阵R3和
的误差。（2）内定向误差，即无畸变像坐标p和p'与其对应的有畸变坐标（图像实际量测坐标）的差异。
图 1 全色-多光谱影像同名点映射误差源示意图
在卫星影像融合配准中，构建全色-多光谱同名像点映射误差补偿模型的目标是求解其中一个传感器（一般为多光谱传感器）影像像点坐标改正值，从而在姿态误差和内定向误差存在的情况下配平等式（3），保证融合解算过程中获得准确的全色-多光谱同名像点。定义多光谱影像的像点坐标改正值为δ，定义全色影像和多光谱影像的有畸变坐标p0和p0'，内定向误差∇p和∇p'，定义旋转因子的误差中的非颤振项为∇(R)1，颤振项为∇(R)2，根据式（5），误差改正值为：
要实现全色-多光谱同名像点映射误差的补偿，则需要构建δ的拟合函数模型，在不依赖外部控制数据的情况下求解。
旋转因子的误差项∇(R)1包含常量误差、随时间缓慢变化的卫星定姿误差，误差项∇(R)2包含随时间周期性变化的颤振误差，因此在拟合过程中，需要将时间作为自变量。但在卫星影像融合产品生产中，一般采用RPC模型进行定向，已无法获取每个扫描行对应的时间戳。针对这一问题，考虑到卫星影像扫描周期的稳定性，采用影像的行号（l）替代时间项，作为旋转因子误差的拟合自变量。
卫星影像的内定向误差主要包括内方位元素误差、镜头畸变、子线阵CCD安置错位3种误差源。其中内方位元素误差引起的像点偏移与探元指向角呈线性关系，镜头畸变引起的像点偏移与探元指向角呈多项式关系，而子线阵CCD安置错位引起的像点偏移与探元指向角构成分段线性函数的关系。由于卫星传感器线阵CCD的几何特性，探元指向角与影像列号（s）近似呈线性关系，因此卫星影像的内定向误差使用影像的列号s作为自变量进行拟合，拟合函数为一个线性函数（内方位元素误差）、一个多项式（镜头畸变）和一个分段线性函数（子线阵CCD安置错位）之和，在模型设计中，将其合并为分段多项式进行拟合，分段个数为子线阵CCD个数。
由于内定向误差一般仅为数个像元大小， 为1×10－6~1×10－5 m水平，因此忽略内定向误差与旋转因子误差的乘积项，考虑到旋转因子R≈I，式（6）可以转化为：
内定向误差项(K∙∇p-∇p')由全色影像内定向误差乘以一个尺度比值，与多光谱影像内定向误差相减得到。但由于高分七号卫星的后视全色传感器和后视多光谱传感器均采用3个等长子线阵CCD构造，且扫描线方向近似平行，全色影像和多光谱影像同名像点的列号近似呈线性关系，因此内定向误差项在拟合中可以合并为一个分段多项式，并采用多光谱影像列号s作为唯一的自变量。定义内定向误差(K∙∇p-∇p')引起的高分七号多光谱影像像点坐标改正值为
，其拟合方式为：
由于推扫式成像周期的稳定性，高分七号全色影像与多光谱影像同名像点行号也近似呈线性关系。旋转因子误差的非颤振项在一景影像范围内引起的像点偏移可以通过像方线性模型近似进行拟合。定义误差项K∙∇(R)1∙(p+R3^(-1) d)与式（7）中它的常数项引起的像点偏移为
，其拟合模型为线性模型，以多光谱影像行列号s和l为自变量：
式中，A1，A2，…，A6为线性模型参数；As(s,l)、Al (s,l)分别为列方向和行方向线性模型；
和
分别为列方向和行方向线性补偿值。
定义旋转因子的颤振项K∙(p+R3^(-1) d)∙∇(R)2引起的像点偏移为
，其拟合模型为以多光谱影像行号l为自变量的三角级数与多光谱影像坐标列号s线性项的乘积：
式中，B1、B2为影像列号s的系数；Ti、Si为第i个三角函数模型的幅值；pi、qi为第i个三角函数模型的频率值；φi、ψi为第i个三角函数模型的相位值；Js(s,l)、Jl (s,l)为三角级数补偿模型；
为颤振项补偿值。
结合上述误差项，高分七号卫星影像的全色-多光谱同名像点映射误差拟合模型为：
在实际应用中，采用四次分段多项式拟合内定向误差的影响，而颤振项的拟合则需要根据不同时期高分七号卫星影像产品的具体误差特性。
1.3
误差补偿模型的求解
高质量融合影像的生产需要保证全色-多光谱同名点映射精度达到1个全色像素以内，即0.25多光谱像素，因此需要获得高精度的稠密匹配点进行配准误差补偿模型的解算。定义全色、多光谱影像之间的匹配点为{(p,p')}，将全色影像上的像点p=(s,l)通过“反投影-重投影”过程（即在航天飞机雷达地形测绘任务（shuttle radar topography mission， SRTM）等高程数据的辅助下利用全色影像RPC参数模型反投影到物方空间，再通过多光谱影像RPC参数模型投影到其多光谱像方空间，见式（3）投影到多光谱影像像平面得到像点p''=(s'',l'')，从而得到映射偏差的观测值：
在这一过程中，SRTM数据的高程误差影响可以忽略，因为高分七号全色相机和多光谱相机的光线夹角仅为约0.094°，而SRTM的标称精度为16 m，对反投影-重投影过程的影响约为0.01多光谱像素。
根据式（12）即可得到观测方程，对融合配准模型进行求解。但在模型中，线性项部分和分段多项式部分存在一定的相关性，而三角级数的求解又需要在排除其他误差的干扰后，结合频率域分析等手段进行，因此本文采用分别求解的方式获得配准模型参数的估值，其步骤如下所示。
1） 求解线性项，将p'=(s',l')和p''=(s'',l'')的坐标差异作为误差进行求解：
2）定义步骤1）中求解得到的线性项模型为
，从误差中减去线性项模型部分，剩余部分用于求解分段项参数：
3）定义步骤2）中求解的分段项模型为
，再从上述误差中减去分段项模型部分，剩余部分用于进行周期性误差的分析求解。
其中，三角级数项的求解方法参考文献［20］。
在映射误差补偿模型（12）中，分段项和颤振项具有较多的参数，需要在同名点均匀覆盖整张影像的情况下进行求解。当影像存在大量水面或云层覆盖时，无法获得均匀的同名点，则无法稳定求解分段项和颤振项。但是，分段项对应的误差源是卫星相机内参数误差、镜头畸变、子线阵CCD安置误差，这一类误差可以长时间保持不变，因此获取时间相近的影像可以复用同一个分段项。在复用分段项的情况下，只需要通过匹配稀疏同名点，求解线性项即可，求解时的观测值为坐标差异与分段项预测值之差：
2 高分七号卫星影像融合配准补偿实验
为了验证本文融合配准模型的有效性，采用了中国广东地区6景和青海地区3景高分七号一级影像进行实验分析。广东地区影像拍摄时间为2021年1月、2月和4月，位置分布如图2（a）所示，影像文件的前缀为“GF701”。青海地区影像拍摄时间均为2022年1月，位置分布如图2（b）所示，影像文件的前缀均为“GF7_DLC”。影像详细信息如表1所示。
图 2 实验影像分布
表 1 实验影像信息说明
2.1
融合配准精度及评价指标
首先对高分七号全色影像与多光谱影像进行匹配，采用§1.3所述的过程对补偿模型参数进行求解。本文在实验中采用了局部图像特征匹配（local feature transformer，LoFTR）算法进行匹配，LoFTR算法首先在粗粒度上建立图像特征的检测、描述和匹配，然后在精粒度上细化亚像素级别的密集匹配，采用预训练的室外场景匹配模型参数，具有极强的弱纹理影像匹配能力。在高分七号全色/多光谱影像进行匹配任务中，匹配点密度可以达到每100像素至少一个匹配点，且分布均匀。
融合配准的精度目标一般是达到全色影像的1像素水平（针对高分七号，相当于0.25多光谱像素），从而消除可见的晕边，而匹配算法往往只能达到0.3~0.5多光谱像素的精度水平。因此实验通过将影像划分为规则格网，对其中匹配点对应的配准误差用鲁棒性递增求解方法进行加权平均，得到更高精度的配准误差值。为了验证格网配准误差的估算精度，本文将多光谱影像进行亚像素平移（如沿行方向和列方向平移0.2~3.8像素）和重采样后，考察平移对格网配准误差估算的影响，发现格网误差的变化（观测值）与实际平移的大小（真值）差异保持在0.1像素以内，这一结果间接证明格网配准误差的量测精度达到了0.1像素水平，能满足融合配准要求。
定义格网中心和配准偏差组成的一组观测值为(s',l',∇s,∇l)，其中s'、l'为格网中心坐标，∇s、∇l为格网中匹配点配准误差平均值。为了评价配准模型中3个子项对融合配准的影响，§2.2统计了对式（14）和式（15）对应的线性项和分段项先后进行拟合以后，垂轨方向和沿轨方向的配准残差中误差和P(∇xy≤0.25)（格网配准误差小于0.25倍多光谱像素的比例），同时为了定量评价融合配准操作对融合影像质量的影响，还统计了求解两组误差后得到的融合影像的无参考图像质量参数（quality with no reference，QNR）评价指标。QNR指标是影像融合中的一种无参考数据的质量评价指标，它包含两项，一项是光谱偏差参考值，一项是几何偏差参考值，通过计算和整合两种参考值，得到融合质量的客观评价指标。§2.3分析了两组影像映射偏差在l方向上的分布情况，发现广州组影像可能存在颤振误差，对这些影像经过线性项和分段项拟合后的残差进行了三角级数拟合（参考式（16）），并统计分析。§2.4测试了利用求解出来的分段项校正同一轨道其他影像映射误差的效果（参考式（17））。
2.2
基于线性项和分段项的融合配准补偿
选取广东地区1~3号影像和青海地区1~2号影像进行融合配准模型求解实验。上述影像的匹配点数均达到1 000万个以上。表2给出了两组影像的格网偏移值的初值、线性项补偿残差、分段项补偿残差的统计结果。图3和图4分别给出了广东组影像1和青海组影像1的补偿残差图。
表 2 实验影像不同模型残差补偿精度
图 3 广东组影像1残差图
图 4 青海组影像1残差图
广东组的3景影像在线性项补偿之前，普遍具有列方向约0.60 像素、行方向约0.30 像素的融合配准中误差，P(∇xy≤0.25)低于1%。经过线性项补偿后，精度提升明显，列方向中误差均低于0.12 像素，行方向中误差均达到约0.20 像素，一定程度上消除了系统性偏移，使得P(∇xy≤0.25)大幅提升到约70%水平，但是，结合图3（b）、3（e）给出的线性项补偿后残差分布，可以发现仍然存在明显的非线性断裂误差。经过分段项补偿后，有效消除了非线性误差，中误差均达到0.06~0.08像素，P(∇xy≤0.25)提升到97%以上水平，满足高质量融合的要求。青海组的影像1在通过线性项和分段项进行补偿后，有效消除了非线性误差，列方向和行方向的中误差均达到0.1像素以内，P(∇xy≤0.25)提升到97%以上，满足高质量融合的要求。而影像2含有一定的云覆盖，校正后P(∇xy≤0.25)仅提升到83.8%，这一结果并不能反映实际的精度，要得到准确的评价结果需要先排除含云区域。经检查，去除云雪覆盖率为50%以上的区块后，P(∇xy≤0.25)达到90%以上。
线性项和分段项补偿前后的融合影像的部分截图如图5所示，补偿前的融合影像图5（a）~5（d）中，在地物颜色突变的位置，往往存在可见的晕边现象，其中图5（a）、5（b）为补偿前融合影像建筑物区域显示效果，补偿前影像的建筑物边缘呈现多彩颜色，以致无法准确判断地物边缘；图5（c）、5（d）为在补偿前融合影像水域的显示效果，补偿前，影像的水体边缘不明确，产生了明显的位移。而在补偿后的融合影像5（e）~5（h）中，晕边现象得到显著的缓解，建筑物、水域不再存在模糊的多重边缘，提高了影像可解释性和可量测性。
图 5 广东组影像2补偿前后的假彩色融合影像细节对比
2.3
基于三角级数项的融合配准补偿
实验结果表明，补偿线性项和分段项后，高分七号影像的配准精度已基本满足不出现融合晕边这一要求，但是由图6和图7中分段项补偿后残差与行方向的误差分布曲线可知，广东组影像的融合配准误差呈现一定的周期性。通过8个三角函数组成的三角级数进行拟合，发现颤振引起的周期性误差振幅均不到0.1像素，拟合后的配准精度有一定提升，列方向的残差中误差由0.07 像素下降为0.06像素，行方向的残差中误差由0.08像素下降为0.07像素。而青海组影像的残差没有沿行方向呈现出周期性。实验结果表明高分七号的颤振误差基本没有对融合质量造成明显影响，可以在生产中忽略不计。
图 6 广东组影像1三角级数项补偿残差图
图7 青海组影像1三角级数项补偿残差图
2.4
分段项的复用性验证
融合补偿模型中分段项的主要作用是补偿由子线阵CCD安置误差和镜头畸变引起的配准误差，由于CCD安置误差和镜头畸变误差可以在较长时期内保持不变，因此分段项参数理论上在一段时间内具有复用性。为验证补偿模型分段项参数的复用性，将广东组影像一的分段项参数直接对影像4、影像5、影像6进行补偿，将青海组影像1的分段参数对影像3进行补偿。补偿的方法是先用分段项参数进行补偿，然后利用分段项的残差求解线性项，最终得到补偿模型。分段项参数复用模式下的补偿精度如表3所示，补偿后的广东组影像中误差基本达到0.07像素水平，青海组影像3中误差达到0.1像素以内，4组影像的P(∇xy≤0.25)均提升到98%水平，能满足融合要求。这一实验结果表明分段项参数在一定时间内具有较好的复用性，可以在影像含云、水等不利于匹配的情况下进行复用。
表 3 复用模式下分段项前后的配准精度
3 结 语
本文从高分七号卫星影像的成像机理出发，对影像融合过程中内外方位元素误差和影像内部畸变对配准精度的影响进行分析和建模，提出了由线性模型、四次分段多项式模型、三角级数模型组成的融合映射误差补偿模型。本文对两个地区的高分七号影像进行融合配准模型拟合，实验结果表明：线性模型可以有效消除内外方位元素误差引起的线性像方误差，将垂轨方向的中误差降至0.1像素左右，沿轨方向中误差降至0.2 像素以下；四次分段多项式模型可明显消除子线阵CCD安置误差引起的像方配准残差不连续性，将垂轨方向中误差降至0.1像素以下，沿轨方向的中误差从约0.2像素降至0.1 像素以内。通过线性项和分段项的结合，使得P(∇xy≤0.25)提升至约98%水平，基本达到“不出现融合晕边”的要求。而颤振模型的求解和补偿有效抑制了周期性的配准误差，达到采用高分七号一级数据生产融合影像的理论极限精度。
本文测试了分段多项式参数的复用，即采用从其他（同轨或异轨）影像中求解的分段多项式模型参数和本影像中求解的线性项参数进行补偿。实验结果表明，复用模式下融合配准的精度和P(∇xy≤0.25)基本达到求解模式的同一水平，证明了分段多项式主要补偿内部镜头畸变和子线阵CCD安置误差，在不改变零级数据生产相关参数的情况下，基本保持不变。因此，可考虑通过如下流程进行高分七号影像融合生产：（1）从整轨或同一批次获取的所有高分七号影像中，选出含云、水较少且纹理清晰（即沙漠、冰雪等不利于匹配的区域较少）的影像，通过LoFTR等先进的匹配算法进行全影像稠密匹配，然后求解本文提出的融合配准模型；（2）对其他影像，复用已求解的融合配准模型的分段项参数，然后仅匹配少量同名点求解线性项，通过复用模式进行补偿。这一方案可以灵活调整，兼顾精度和效率。
转自：“科研圈内人”微信公众号
如有侵权，请联系本站删除！
