CN105701830A - 基于几何模型的lasis波段影像配准方法及系统 - Google Patents

Abstract

本发明提供一种基于几何模型的LASIS波段影像配准方法及系统，包括首先将LASIS成像模型等效为多线阵推扫模型，先进行相机标定，利用不同波段影像重叠区域构建平行观测值，利用两个波段的平行观测量检测平台颤振，并根据利用检测到的颤振特征补偿其他波段颤振造成的像素偏移；更新成像模型的姿态信息，利用精化所得新的成像模型与无畸变的成像模型对准，对原始影像进行重采样，实现波段影像的自动配准。传统方法两两波段匹配，其构建的干涉序列中存在累积误差从而精度较差，本发明规避了以上缺陷，能够实现干涉序列影像的整体高精度配准，改进了各波段自我补偿的传统方法，能够极大提升处理效率。

Description

基于几何模型的LASIS波段影像配准方法及系统

技术领域

本发明属于摄影测量与遥感的影像处理领域，涉及一种基于严密成像几何模型以及平行观测法消除颤振的LASIS波段影像配准方法及系统。

背景技术

大孔径静态干涉成像光谱仪(LASIS)是基于推扫的傅里叶变换型成像光谱仪，该类型光谱仪具有稳定性高，结构紧凑等特点，在轨运行性能较好。LASIS型干涉仪光谱性能依赖于不同波段影像间的配准精度，然而由于其时空联合调制的特点，成像过程中，同一地物在不同波段影像成像时刻不一致，一旦卫星平台产生扰动，将对不同波段影像像元位置产生不同影响。而姿态抖动现象是卫星运行过程中的常见现象，因此，如何消除由于姿态抖动以及光学畸变等因素导致的不同波段影像中的各异性误差，从而实现干涉序列的精确配准是本发明研究的重点。

LASIS高光谱影像具有时空联合调制、星上结构稳定、地面处理复杂、波段众多等特点。传统的高光谱干涉序列提取技术，是针对两两波段，使用相位相关算法或者大型平差矩阵等方法等纠正像元位置，实现像素间配准，这样获取的干涉序列会由于LASIS波段较多存在累积误差而导致整体序列质量下降，同时LASIS影像波段众多，两两配准效率较低。

国内外学者为了实现LASIS干涉图的精确获取做了大量研究，主要从三个方面进行了努力：一是配准算法的完善，例如：CASTRO等提出了基于亚像素精度的相位配准算法进行配准；李湘滨等分析姿态角随机变化造成的像面畸变规律，提出用相位相关等算法校正配准影像，从而校正LASIS影像旋转和偏移；曾琪明和叶沅鑫等分别提出了基于谱相似和局部自相似的影像配准方法；二是在配准算法的基础上，利用平差方法求解像元精确位置，例如:Gengzhang提出了通过亚像素配准算法获取多波段影像配准点，建立大规模误差方程求解像点的最或是值作为像元最优位置解；三是从成像链路方面分析LASIS影像的误差源，并讨论影响LASIS配准的几何原因。例如YannFerrec等系统地分析了降低LASIS光谱成像质量的因子，并指出姿态颤振将直接影响干涉图的质量，但是并未给出具体的解决方案。

发明内容

针对现有LASIS干涉序列配准方法存在的缺点，本发明的目的是为了提供一种能够整体配准LASIS干涉序列的亚像素级配准技术方案。

本发明技术方案提供一种基于几何模型的LASIS波段影像配准方法，包括以下步骤，

步骤a，将LASIS面阵探测器等效为多线阵探测器，建立多线阵推扫成像模型，所述多线阵推扫成像模型中的每条垂轨向CCD模型符合单线阵推扫模型，线阵之间无缝排列，每个线阵成像对应LASIS一个单波段影像，其中各个波段p的严密成像模型如下，

$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\left[\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right]+{\mathrm{mR}}_{BW}{R}_{CB}{R}_U\left[\begin{array}{c}tan{\left({\mathrm{\ψ}}_x\right)}_p\\ \tan {\left({\mathrm{\ψ}}_y\right)}_p\\ -1\end{array}\right]---\left(1\right)$

式中，[XYZ]^T为像点对应的地面点的WGS84空间直角坐标系的坐标值，m为比例系数，[X_SY_SZ_S]^T为成像时卫星在WGS84坐标系下的位置矢量，R_BW为本体坐标系相对WGS84坐标系的旋转矩阵，R_CB为相机坐标系相对本体坐标系的旋转矩阵，R_U代表偏置矩阵，(ψ_x)_p,(ψ_y)_p分别为第p行CCD任一探元对应的沿轨向指向角和垂轨向指向角，计算如下，

$\begin{array}{c}tan{\left({\mathrm{\ψ}}_x\right)}_p=\frac{x-x_0-\mathrm{\Δ}x}{-f}\\ tan{\left({\mathrm{\ψ}}_y\right)}_p=\frac{y-y_0-\mathrm{\Δ}y}{-f}\end{array}---\left(2\right)$

式中，(x,y)为像点影像坐标、(x₀,y₀,f)表示内方位元素，(Δx,Δy)为相机畸变引起的像点偏移；

步骤b，根据控制数据进行相机标定，包括利用Harris算法提取基准波段影像角点，求解物方信息；利用最小二乘配准算法在非基准波段影像获取同名点作为对应的像方信息，基于相应波段的严密成像模型利用同名点求解指向角，并把指向角和对应的垂轨向坐标s代入如下公式，利用最小二乘计算拟合的系数a_t,b_t，0≤t≤5，

$\left\{\begin{array}{c}\tan \left({\mathrm{\ψ}}_x\right)=a_0+a_{1}s+a_2{s}^2+a_3{s}^3+a_4{s}^4+a_5{s}^5\\ \tan \left({\mathrm{\ψ}}_y\right)=b_0+b_{1}s+b_2{s}^2+b_3{s}^3+b_4{s}^4+b_5{s}^5\end{array}\right.---\left(3\right)$

对姿轨模型进行平滑，并按照平滑后的姿轨模型对全部波段影像重采样，从而实现影像成像范围的一致性；

步骤c，利用不同波段影像重叠区域构建平行观测值，利用两个波段的平行观测量检测平台颤振，并根据利用检测到的颤振特征补偿其他波段颤振造成的像素偏移，获取精确像点坐标和对应的控制点坐标；实现方式如下，

设基准波段的颤振模型为

其中，V_x、V_y分别代表沿轨向像差与垂轨像像差在沿轨向的投影；ε_x、ε_y代表沿轨向像差与垂轨向像差中的随机成分；代表某个沿轨向颤振分量的振幅，频率和相位；(B_jω_jψ_j)代表某个垂轨向颤振分量的振幅，频率和相位；i代表某个沿轨向颤振成分的序号，j代表某个垂轨向颤振成分的序号，m代表沿轨向颤振成分的个数，n代表垂轨向颤振成分的个数；

构建基准波段影像的平行观测量，其他波段的颤振模型为

其中，分别代表沿轨向像差与垂轨像像差在沿轨向的投影；代表沿轨向像差与垂轨向像差中的随机成分；代表某个沿轨向颤振分量的振幅，频率和相位；代表某个垂轨向颤振分量的振幅，频率和相位；i代表某个沿轨向颤振成分的序号，j代表某个沿轨向颤振成分的序号，m代表沿轨向颤振成分的个数，n代表垂轨向颤振成分的个数；

得到两波段同名点配准残差中颤振影响对应的公式如下，

其中，ΔV_x、ΔV_y代表某一波段与基准波段配准残差的沿轨向坐标和垂轨向坐标；

设为存在以下等价条件，

则

利用焦面排列特点和成像模型根据波段成像时间差，改变参量t，得到其他波段与参考波段由于颤振造成的像差量；

步骤d，利用步骤c所得结果重新计算控制点对应的姿态信息，即更新其对应的相机坐标系相对本体坐标系的旋转矩阵R_CB；

步骤e，利用步骤d更新结果进一步精化式(1)所示成像模型，利用精化所得新的成像模型与无畸变的成像模型对准，对原始影像进行重采样，实现波段影像的自动配准。

而且，步骤e中，利用步骤d更新结果进一步精化成像模型实现方式如下，

根据下式，

其中，分别代表俯仰角和侧滚角由于平台振荡引起的角度变化，a是沿轨向归化系数，是沿轨向平移量；b是垂轨向归化系数，是垂轨向平移量；

利用姿态文件中各姿态点对应的时间t参量计算并让姿态文件中的姿态点中的俯仰角和侧滚角分别减去

本发明还相应提供一种基于几何模型的LASIS波段影像配准系统，包括以下模块，

第一模块，用于将LASIS面阵探测器等效为多线阵探测器，建立多线阵推扫成像模型，所述多线阵推扫成像模型中的每条垂轨向CCD模型符合单线阵推扫模型，线阵之间无缝排列，每个线阵成像对应LASIS一个单波段影像，其中各个波段p的严密成像模型如下，

$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\left[\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right]+{\mathrm{mR}}_{BW}{R}_{CB}{R}_U\left[\begin{array}{c}tan{\left({\mathrm{\ψ}}_x\right)}_p\\ \tan {\left({\mathrm{\ψ}}_y\right)}_p\\ -1\end{array}\right]---\left(1\right)$

式中，[XYZ]^T为像点对应的地面点的WGS84空间直角坐标系的坐标值，m为比例系数，[X_SY_SZ_S]^T为成像时卫星在WGS84坐标系下的位置矢量，R_BW为本体坐标系相对WGS84坐标系的旋转矩阵，R_CB为相机坐标系相对本体坐标系的旋转矩阵，R_U代表偏置矩阵，(ψ_x)_p,(ψ_y)_p分别为第p行CCD任一探元对应的沿轨向指向角和垂轨向指向角，计算如下，

$\begin{array}{c}tan{\left({\mathrm{\ψ}}_x\right)}_p=\frac{x-x_0-\mathrm{\Δ}x}{-f}\\ tan{\left({\mathrm{\ψ}}_y\right)}_p=\frac{y-y_0-\mathrm{\Δ}y}{-f}\end{array}---\left(2\right)$

式中，(x,y)为像点影像坐标、(x₀,y₀,f)表示内方位元素，(Δx,Δy)为相机畸变引起的像点偏移；

第二模块，用于根据控制数据进行相机标定，包括利用Harris算法提取基准波段影像角点，求解物方信息；利用最小二乘配准算法在非基准波段影像获取同名点作为对应的像方信息，基于相应波段的严密成像模型利用同名点求解指向角，并把指向角和对应的垂轨向坐标s代入如下公式，利用最小二乘计算拟合的系数a_t,b_t，0≤t≤5，

$\left\{\begin{array}{c}\tan \left({\mathrm{\ψ}}_x\right)=a_0+a_{1}s+a_2{s}^2+a_3{s}^3+a_4{s}^4+a_5{s}^5\\ \tan \left({\mathrm{\ψ}}_y\right)=b_0+b_{1}s+b_2{s}^2+b_3{s}^3+b_4{s}^4+b_5{s}^5\end{array}\right.---\left(3\right)$

对姿轨模型进行平滑，并按照平滑后的姿轨模型对全部波段影像重采样，从而实现影像成像范围的一致性；

第三模块，用于利用不同波段影像重叠区域构建平行观测值，利用两个波段的平行观测量检测平台颤振，并根据利用检测到的颤振特征补偿其他波段颤振造成的像素偏移，获取精确像点坐标和对应的控制点坐标；实现方式如下，

设基准波段的颤振模型为

其中，V_x、V_y分别代表沿轨向像差与垂轨像像差在沿轨向的投影；ε_x、ε_y代表沿轨向像差与垂轨向像差中的随机成分；代表某个沿轨向颤振分量的振幅，频率和相位；(B_jω_jψ_j)代表某个垂轨向颤振分量的振幅，频率和相位；i代表某个沿轨向颤振成分的序号，j代表某个垂轨向颤振成分的序号，m代表沿轨向颤振成分的个数，n代表垂轨向颤振成分的个数；

构建基准波段影像的平行观测量，其他波段的颤振模型为

其中，分别代表沿轨向像差与垂轨像像差在沿轨向的投影；代表沿轨向像差与垂轨向像差中的随机成分；代表某个沿轨向颤振分量的振幅，频率和相位；代表某个垂轨向颤振分量的振幅，频率和相位；i代表某个沿轨向颤振成分的序号，j代表某个沿轨向颤振成分的序号，m代表沿轨向颤振成分的个数，n代表垂轨向颤振成分的个数；

得到两波段同名点配准残差中颤振影响对应的公式如下，

其中，ΔV_x、ΔV_y代表某一波段与基准波段配准残差的沿轨向坐标和垂轨向坐标；

设为存在以下等价条件，

则

利用焦面排列特点和成像模型根据波段成像时间差，改变参量t，得到其他波段与参考波段由于颤振造成的像差量；

第四模块，用于利用第三模块所得结果重新计算控制点对应的姿态信息,更新其对应的相机坐标系相对本体坐标系的旋转矩阵R_CB；

第五模块，用于利用第四模块更新结果进一步精化式(1)所示成像模型，利用精化所得新的成像模型与无畸变的成像模型对准，对原始影像进行重采样，实现波段影像的自动配准。

而且，第五模块中，利用第四模块更新结果进一步精化成像模型实现方式如下，

根据下式，

其中，分别代表俯仰角和侧滚角由于平台振荡引起的角度变化，a是沿轨向归化系数，是沿轨向平移量；b是垂轨向归化系数，是垂轨向平移量；

利用姿态文件中各姿态点对应的时间t参量计算并让姿态文件中的姿态点中的俯仰角和侧滚角分别减去

综上所述，本发明具有如下特点：

(1)针对实际成像过程中的误差源进行消除，效果理想。

(2)针对LASIS成像特点与焦面特点，提出利用多线阵推扫模型代替实际面阵模型，充分利用各线阵之间的成像相关性，利用两个波段即可消除所有波段的颤振影响，可操作性强。

与其他LASIS干涉序列获取方式相比，本发明方法针对成像过程中影响各波段成像一致性的误差源进行处理，能够实现各波段的整体配准，方法简便，适用性强，性能稳定，能够达到亚像素级的配准精度。

附图说明

图1为本发明实施例的流程图。

图2为本发明实施例的A景第1波段直接定位残差分布图，其中图2a为垂轨残差沿列向变化规律示意图，图2b为沿轨残差沿列向变化规律示意图。

图3为本发明实施例的A景第1波段内检校后定位残差分布图，其中图3a为垂轨残差沿列向变化规律示意图，图3b为沿轨残差沿列向变化规律示意图。

图4为本发明实施例的B景影像波段配准精度示意图。

图5为本发明实施例的C景影像波段配准精度示意图。

图6为本发明实施例的D景影像波段配准精度示意图。

具体实施方式

下面结合附图和实施例对本发明做进一步说明。

为了实现上述目的，本发明的技术方案是一种基于严密成像几何模型以消除光学畸变以及姿态颤振的摄影测量方法，主要包含两方面的创新，一是将面阵模型等效成多线阵几何模型进行建模，二是使用两个波段检测颤振特征，补偿所有波段颤振导致的像差，改进了以前各个波段自我补偿颤振影响的方案。

本发明的实施例是以天津地区影像(下称A)为检校数据，对日本三沢市(下称B)以及澳大利亚昂佩利附近两景影像(下称C、D)进行实验，预先输入相应辅助文件，一般包括姿态、轨道、行时以及相机参数文件。参照图1，本发明实施例的具体步骤如下：

步骤a，将LASIS面阵探测器等效为多线阵探测器，建立多线阵推扫成像模型，其中线阵之间无缝排列，每个线阵成像对应LASIS一个单波段影像；

步骤a的具体实现如下：

多线阵推扫成像模型中的每条垂轨向CCD模型符合单线阵推扫模型，每个线阵成像对应LASIS一个单波段影像，其中各个波段p的严密成像模型如下：

$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\left[\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right]+{\mathrm{mR}}_{BW}{R}_{CB}{R}_U\left[\begin{array}{c}tan{\left({\mathrm{\ψ}}_x\right)}_p\\ \tan {\left({\mathrm{\ψ}}_y\right)}_p\\ -1\end{array}\right]---\left(1\right)$

式中，[XYZ]^T为像点对应的地面点的WGS84空间直角坐标系的坐标值，m为比例系数，[X_SY_SZ_S]^T为成像时卫星在WGS84坐标系下的位置矢量，R_BW为本体坐标系相对WGS84坐标系的旋转矩阵，R_CB为相机坐标系相对本体坐标系的旋转矩阵；R_U代表偏置矩阵，(ψ_x)_p,(ψ_y)_p分别为第p行CCD(对应第p波段影像)任一探元对应的沿轨向指向角和垂轨向指向角；

其中，指向角计算公式如下，式中，(x,y)为像点影像坐标、(x₀,y₀,f)表示内方位元素，(Δx,Δy)为相机畸变引起的像点偏移，(ψ_x)_p,(ψ_y)_p表示对应的探元指向角。

$\begin{array}{c}tan{\left({\mathrm{\ψ}}_x\right)}_p=\frac{x-x_0-\mathrm{\Δ}x}{-f}\\ tan{\left({\mathrm{\ψ}}_y\right)}_p=\frac{y-y_0-\mathrm{\Δ}y}{-f}\end{array}---\left(2\right)$

实施例中，建立各个波段的严密成像模型(针对影像A、B、C、D)。

步骤b，根据控制数据进行相机标定：根据多线阵推扫成像模型和控制数据，消除光学畸变，纠正CCD排列畸变和由于星上姿轨不平滑导致的成像范围不一致问题。

步骤b的具体实现如下：

利用Harris算法提取基准波段(LASIS共有180个波段，一般可预设为中间的波段，例如第90波段)影像角点，加辅助信息(DEM,DOM)正算求解物方信息；利用最小二乘配准算法在非基准影像获取同名点作为对应的像方信息，从而得到控制信息。建立指向角标定模型。假定某非基准波段p与基准波段配准获得了同名点，根据公式(1)基于波段p的严密成像模型计算这些点的指向角；如公式(3)所示，其中，a_t,b_t(0≤t≤5)分别代表沿轨向拟合系数；s代表像点的垂轨像坐标。即如公式(1)所示，基于成像模型利用配准点(同名点)求解指向角，并把指向角和其对应的垂轨向坐标s代入如下的公式(3)，利用最小二乘计算拟合的系数a_t,b_t(0≤t≤5)。

$\left\{\begin{array}{c}\tan \left({\mathrm{\ψ}}_x\right)=a_0+a_{1}s+a_2{s}^2+a_3{s}^3+a_4{s}^4+a_5{s}^5\\ \tan \left({\mathrm{\ψ}}_y\right)=b_0+b_{1}s+b_2{s}^2+b_3{s}^3+b_4{s}^4+b_5{s}^5\end{array}\right.---\left(3\right)$

并将求解后的系数a_t,b_t(0≤t≤5)重新代入(3)，利用该式纠正各个整数像元的指向角。

为实现成像均匀，飞行平台的状态(姿态、轨道状态等)应均匀维持。可以对姿轨模型进行平滑，并按照平滑后的姿轨模型对全部波段影像重采样，从而实现影像成像范围的一致性(现有技术)。

实施例中，利用A景影像进行相机标定，探测实际成像过程中成像系统的畸变情况，提取定标参数(偏置矩阵R_U以及指向角拟合系数a_t,b_t(0≤t≤5))，对B、C、D景影像进行处理。对B、C、D景影像进行纠正之后，并将像差在垂轨和沿轨方向显示，如果存在时变性误差，则可推断存在姿态抖动现象。将重叠区域多个同名像点的像差在垂轨向展开，观察像面畸变情况。如图2所示，其中图3a为垂轨残差沿列向变化规律示意图，图3b为沿轨残差沿列向变化规律示意图，垂轨向残差在列向的投影，像差的标准方差在0.1像素左右，说明第一行CCD的探元整体存在向左偏移的规律；沿轨向残差在垂轨向的投影显示一个非线性曲线，说明第一行CCD在沿轨方向也存在非线性畸变。同一行CCD上的探元，如果焦面畸变不纠正，则这行探元投影到沿轨向时，是焦面畸变与颤振畸变的综合效果，焦面畸变有可能淹没颤振效果，从而无法探测颤振规律，利用指向角标定法进行检校后，残差如图3所示，其中图3a为垂轨残差沿列向变化规律示意图，图3b为沿轨残差沿列向变化规律示意图。

步骤c，利用不同波段影像重叠区域构建平行观测值，利用两个波段的平行观测量检测平台颤振，并根据成像模型(即波段间的成像模型的关联)补偿其他波段颤振造成的像素偏移。即利用两个波段进行颤振特征检测，利用检测到的颤振特征实现补偿所有波段颤振导致的像差，获取精确像点坐标和对应的控制点坐标。

步骤c的具体实现如下：

将经过步骤b中预处理的同名点对残差展示在沿轨方向上，观察是否存在规律性的像差抖动成分。因为各种造成像差的原因在步骤b中进行了分析，只有姿态颤振跟时间相关，因此规律性三角函数成分可以认为是由姿态颤振引起的。

对不同波段重叠区域进行LSM匹配获取同名点以及消除CCD排列畸变是准备工作。在此基础上，对规律性的像差抖动成分进行频谱分析:将像差曲线采样作为离散傅里叶样本点，利用傅里叶变换进行频谱分析，用得到的主频成分模拟合成的波形与残差的形状比对，确定波形的一致性。

由公式1可知，假设基准波段的颤振模型为：

其中，V_xV_y分别代表沿轨向像差与垂轨像像差在沿轨向的投影。ε_x,ε_y代表沿轨向像差与垂轨向像差中的随机成分；代表某个沿轨向颤振分量的振幅，频率和相位；(B_jω_jψ_j)代表某个垂轨向颤振分量的振幅，频率和相位；i代表某个沿轨向颤振成分的序号，j代表某个垂轨向颤振成分的序号，m代表沿轨向颤振成分的个数，n代表垂轨向颤振成分的个数，t为时间。

其他波段(基准波段影像的平行观测量)的颤振模型为：

其中，分别代表沿轨向像差与垂轨像像差在沿轨向的投影。代表沿轨向像差与垂轨向像差中的随机成分；代表某个沿轨向颤振分量的振幅，频率和相位；代表某个垂轨向颤振分量的振幅，频率和相位；i代表某个沿轨向颤振成分的序号，j代表某个沿轨向颤振成分的序号，m代表沿轨向颤振成分的个数，n代表垂轨向颤振成分的个数，t'为时间。

则(4)与(5)相减，能够得到两波段同名点配准残差中颤振影响(已经消除光学畸变以及姿轨不平滑影响)对应的公式：

其中，ΔV_x、ΔV_y代表某一波段与基准波段配准残差的沿轨向坐标和垂轨向坐标；其他符号含义与公式(4)、(5)中的对应。

由于一景影像的成像时间很短，在该段时间内，认为平台颤振是规律性颤振，则(4)与(5)中的颤振成分是相等的，即存在以下等价条件:

因此，(6)式又可以表示为：

实际的处理过程中，消除光学畸变并经过姿轨平滑后的不同波段同名点像差可以经过由控制点反算像点坐标减去配准点像面坐标获得，此时的像差提取的颤振的振幅频率和相位 如(8)中所示，利用焦面排列特点和成像模型根据波段成像时间差，改变参量t，例如辅助数据中会给出一景影像的起始成像时间t₀和截止成像时间t₁，成像行数为z,根据推扫成像原理，则第f个波段的第k行地物的成像时间为即可得到其他波段与参考波段由于颤振造成的像差量。

实施例中，姿态颤振的检测与补偿：对第1波段与基准波段(第90波段)的像面残差进行谐波分析，并且设定阈值条件提取主频率(振幅值超过其余成分的10倍以上)，分别提取B、C、D三景验证景影像的频谱特征。以第1波段与基准波段提取的颤振特征为基础参量，由于LASIS焦面的密集连接与其时空联合调制型特点，波段间颤振导致的像差是相关的。根据第1波段模型与基准模型颤振导致的像差能够推导出其他波段模型与基准模型的颤振导致的像差，并消除此像差,利用原始控制点的像面坐标减去颤振引起的像差，得到更精确的像点坐标。

步骤d，利用步骤c获取的精确像点坐标和对应的控制点坐标重新计算更新的像点对应的姿态值，即根据公式(1)，利用控制点与像点计算相机的姿态R_CB，具体计算实现为现有技术)。R_CB是姿态角的旋转矩阵表示，还可以将姿态表示为三个方向上的角度，即欧拉角表示，本发明考虑欧拉角中侧滚和翻转方向的变化。

实施例中，利用新的像点坐标值作为控制点，根据线阵推扫严密几何模型重算欧拉角。研究表明侧滚角(俯仰角)抖动和像点的垂轨(沿轨)像差是线性关系，对两者进行线性回归分析，获取回归系数。从而获取侧滚角(俯仰角)的振荡模型。分别代表俯仰角和侧滚角由于平台振荡引起的角度变化，a是沿轨向归化系数，是沿轨向平移量；b是垂轨向归化系数，是垂轨向平移量；其余参数含义与公式(4)同。

步骤e，利用步骤d更新的姿态数据进一步精化公式(1)中的成像模型(即根据公式(9)，利用姿态文件中各姿态点对应的时间t参量计算并让姿态文件中的姿态点中的俯仰角(侧滚角)减去理论上，精化后的模型对应真实的成像模型，通过建立一个理想的成像模型(假定焦面无畸变，即模型焦面完全符合面阵设计)，将真实的影像重采样，使之符合无畸变模型(通过真实模型和理想模型的之间的正算、反算即可实现，本发明不予赘述)，即可实现波段影像的自动配准。

具体实施时，利用新的成像模型与无畸变的成像模型对准，对原始影像进行重采样，即可实现波段影像的自动配准。

用精化后的成像模型代替真实成像模型，理想的成像模型是指代一种姿态稳定变化，轨道参数稳定变化，并且焦面无畸变的一种成像系统，这种系统是被人为构造出来的，为的是达到各个波段影像自动配准的目的，由于它的参数也是已知的，真实成像模型的参数也是已知的，真实影像也是有的，通过这三者就可以求解波段自动配准的影像，采用正算加反算的方法即可实现求解。实施例中，根据以上的焦面畸变情况与颤振特征更新原始模型，同时建立无畸变焦面的理想模型，对原始影像(符合真实模型)重采样，使其满足无畸变焦面模型的成像条件，即可实现各个波段影像的自动配准。

对重采样后的影像可以进行精度评定：

B、C、D景影像经处理后的结果分别如图4、图5、图6所示，从影像配准精度结果上分析，三景影像B、C、D均利用第一波段与参考波段提取颤振特征，并补偿所有波段的颤振残差导致的影响，经验证，三景的波段配准精度均达到了0.1像素以上，实验结果非常稳定，适用性很好。

按照上述实验结果主观评价如下：

传统的LASIS干涉序列提取是针对波段影像变形利用配准算法校正，实现波段配准，实际处理过程中，计算量大，实际效果不理想。LASIS波段影像的变形是由于相机内部焦面畸变以及卫星平台颤振对各波段影像影响不一致导致的，通过基于严密成像模型，消除光学畸变，并针对其焦面特点及时空联合调制特点，利用两个波段提取整体平台颤振特并补偿整体波段，能够有效消除平台以及相机对波段成像的影响，从而实现LASIS各波段影像的整体配准，有效服务于干涉序列提取。

具体实施时，本发明所提供方法可基于软件技术实现自动运行流程，也可采用模块化方式实现相应系统。本发明实施例还相应提供一种基于几何模型的LASIS波段影像配准系统，包括以下模块，

第一模块，用于将LASIS面阵探测器等效为多线阵探测器，建立多线阵推扫成像模型，所述多线阵推扫成像模型中的每条垂轨向CCD模型符合单线阵推扫模型，线阵之间无缝排列，每个线阵成像对应LASIS一个单波段影像，其中各个波段p的严密成像模型如下，

$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\left[\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right]+{\mathrm{mR}}_{BW}{R}_{CB}{R}_U\left[\begin{array}{c}tan{\left({\mathrm{\ψ}}_x\right)}_p\\ \tan {\left({\mathrm{\ψ}}_y\right)}_p\\ -1\end{array}\right]---\left(1\right)$

式中，[XYZ]^T为像点对应的地面点的WGS84空间直角坐标系的坐标值，m为比例系数，[X_SY_SZ_S]^T为成像时卫星在WGS84坐标系下的位置矢量，R_BW为本体坐标系相对WGS84坐标系的旋转矩阵，R_CB为相机坐标系相对本体坐标系的旋转矩阵，R_U代表偏置矩阵，(ψ_x)_p,(ψ_y)_p分别为第p行CCD任一探元对应的沿轨向指向角和垂轨向指向角，计算如下，

$\begin{array}{c}tan{\left({\mathrm{\ψ}}_x\right)}_p=\frac{x-x_0-\mathrm{\Δ}x}{-f}\\ tan{\left({\mathrm{\ψ}}_y\right)}_p=\frac{y-y_0-\mathrm{\Δ}y}{-f}\end{array}---\left(2\right)$

式中，(x,y)为像点影像坐标、(x₀,y₀,f)表示内方位元素，(Δx,Δy)为相机畸变引起的像点偏移；

第二模块，用于根据控制数据进行相机标定，包括利用Harris算法提取基准波段影像角点，求解物方信息；利用最小二乘配准算法在非基准波段影像获取同名点作为对应的像方信息，基于相应波段的严密成像模型利用同名点求解指向角，并把指向角和对应的垂轨向坐标s代入如下公式，利用最小二乘计算拟合的系数a_t,b_t，0≤t≤5，

$\left\{\begin{array}{c}\tan \left({\mathrm{\ψ}}_x\right)=a_0+a_{1}s+a_2{s}^2+a_3{s}^3+a_4{s}^4+a_5{s}^5\\ \tan \left({\mathrm{\ψ}}_y\right)=b_0+b_{1}s+b_2{s}^2+b_3{s}^3+b_4{s}^4+b_5{s}^5\end{array}\right.---\left(3\right)$

对姿轨模型进行平滑，并按照平滑后的姿轨模型对全部波段影像重采样，从而实现影像成像范围的一致性；

第三模块，用于利用不同波段影像重叠区域构建平行观测值，利用两个波段的平行观测量检测平台颤振，并根据利用检测到的颤振特征补偿其他波段颤振造成的像素偏移，获取精确像点坐标和对应的控制点坐标；实现方式如下，

设基准波段的颤振模型为

其中，V_x、V_y分别代表沿轨向像差与垂轨像像差在沿轨向的投影；ε_x、ε_y代表沿轨向像差与垂轨向像差中的随机成分；代表某个沿轨向颤振分量的振幅，频率和相位；(B_jω_jψ_j)代表某个垂轨向颤振分量的振幅，频率和相位；i代表某个沿轨向颤振成分的序号，j代表某个垂轨向颤振成分的序号，m代表沿轨向颤振成分的个数，n代表垂轨向颤振成分的个数；

构建基准波段影像的平行观测量，其他波段的颤振模型为

其中，分别代表沿轨向像差与垂轨像像差在沿轨向的投影；代表沿轨向像差与垂轨向像差中的随机成分；代表某个沿轨向颤振分量的振幅，频率和相位；代表某个垂轨向颤振分量的振幅，频率和相位；i代表某个沿轨向颤振成分的序号，j代表某个沿轨向颤振成分的序号，m代表沿轨向颤振成分的个数，n代表垂轨向颤振成分的个数；

得到两波段同名点配准残差中颤振影响对应的公式如下，

其中，ΔV_x、ΔV_y代表某一波段与基准波段配准残差的沿轨向坐标和垂轨向坐标；

设为存在以下等价条件，

则

利用焦面排列特点和成像模型根据波段成像时间差，改变参量t，得到其他波段与参考波段由于颤振造成的像差量；

第四模块，用于利用第三模块所得结果重新计算更新相机坐标系相对本体坐标系的旋转矩阵R_CB；

第五模块，用于利用第四模块更新结果进一步精化式(1)所示成像模型，利用精化所得新的成像模型与无畸变的成像模型对准，对原始影像进行重采样，实现波段影像的自动配准。

各模块具体实现可参见相应步骤，本发明不予赘述。

Claims (4)

1.一种基于几何模型的LASIS波段影像配准方法，其特征在于：包括以下步骤，

步骤a，将LASIS面阵探测器等效为多线阵探测器，建立多线阵推扫成像模型，所述多线阵推扫成像模型中的每条垂轨向CCD模型符合单线阵推扫模型，线阵之间无缝排列，每个线阵成像对应LASIS一个单波段影像，其中各个波段p的严密成像模型如下，

$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\left[\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right]+{\mathrm{mR}}_{BW}{R}_{CB}{R}_U\left[\begin{array}{c}tan{\left({\mathrm{\ψ}}_x\right)}_p\\ tan{\left({\mathrm{\ψ}}_y\right)}_p\\ -1\end{array}\right]---\left(1\right)$

式中，[XYZ]^T为像点对应的地面点的WGS84空间直角坐标系的坐标值，m为比例系数，[X_SY_SZ_S]^T为成像时卫星在WGS84坐标系下的位置矢量，R_BW为本体坐标系相对WGS84坐标系的旋转矩阵，R_CB为相机坐标系相对本体坐标系的旋转矩阵，R_U代表偏置矩阵，(ψ_x)_p,(ψ_y)_p分别为第p行CCD任一探元对应的沿轨向指向角和垂轨向指向角，计算如下，

$\begin{array}{c}\tan {\left({\mathrm{\ψ}}_x\right)}_p=\frac{x-x_0-\mathrm{\Δ}x}{-f}\\ \tan {\left({\mathrm{\ψ}}_y\right)}_p=\frac{y-y_0-\mathrm{\Δ}y}{-f}\end{array}---\left(2\right)$

式中，(x,y)为像点影像坐标、(x₀,y₀,f)表示内方位元素，(Δx,Δy)为相机畸变引起的像点偏移；

步骤b，根据控制数据进行相机标定，包括利用Harris算法提取基准波段影像角点，求解物方信息；利用最小二乘配准算法在非基准波段影像获取同名点作为对应的像方信息，基于相应波段的严密成像模型利用同名点求解指向角，并把指向角和对应的垂轨向坐标s代入如下公式，利用最小二乘计算拟合的系数a_t,b_t，0≤t≤5，

$\left\{\begin{array}{c}tan\left({\mathrm{\ψ}}_x\right)=a_0+a_{1}s+a_2{s}^2+a_3{s}^3+a_4{s}^4+a_5{s}^5\\ tan\left({\mathrm{\ψ}}_y\right)=b_0+b_{1}s+b_2{s}^2+b_3{s}^3+b_4{s}^4+b_5{s}^5\end{array}\right.---\left(3\right)$

对姿轨模型进行平滑，并按照平滑后的姿轨模型对全部波段影像重采样，从而实现影像成像范围的一致性；

步骤c，利用不同波段影像重叠区域构建平行观测值，利用两个波段的平行观测量检测平台颤振，并根据利用检测到的颤振特征补偿其他波段颤振造成的像素偏移，获取精确像点坐标和对应的控制点坐标；实现方式如下，

设基准波段的颤振模型为

其中，V_x、V_y分别代表沿轨向像差与垂轨像像差在沿轨向的投影；ε_x、ε_y代表沿轨向像差与垂轨向像差中的随机成分；代表某个沿轨向颤振分量的振幅，频率和相位；(B_jω_jψj)代表某个垂轨向颤振分量的振幅，频率和相位；i代表某个沿轨向颤振成分的序号，j代表某个垂轨向颤振成分的序号，m代表沿轨向颤振成分的个数，n代表垂轨向颤振成分的个数；

构建基准波段影像的平行观测量，其他波段的颤振模型为

其中，V’_x、V’_y分别代表沿轨向像差与垂轨像像差在沿轨向的投影；ε’_x、ε’_y代表沿轨向像差与垂轨向像差中的随机成分；代表某个沿轨向颤振分量的振幅，频率和相位；(B'_jω'_jψ'_j)代表某个垂轨向颤振分量的振幅，频率和相位；i代表某个沿轨向颤振成分的序号，j代表某个沿轨向颤振成分的序号，m代表沿轨向颤振成分的个数，n代表垂轨向颤振成分的个数；

得到两波段同名点配准残差中颤振影响对应的公式如下，

其中，ΔV_x、ΔV_y代表某一波段与基准波段配准残差的沿轨向坐标和垂轨向坐标；

设为存在以下等价条件，

则

利用焦面排列特点和成像模型根据波段成像时间差，改变参量t，得到其他波段与参考波段由于颤振造成的像差量；

步骤d，利用步骤c所得结果重新计算控制点对应的姿态信息，更新其对应的相机坐标系相对本体坐标系的旋转矩阵R_CB；

步骤e，利用步骤d更新结果进一步精化式(1)所示成像模型，利用精化所得新的成像模型与无畸变的成像模型对准，对原始影像进行重采样，实现波段影像的自动配准。

2.根据权利要求1基于几何模型的LASIS波段影像配准方法，其特征在于：步骤e中，利用步骤d更新结果进一步精化成像模型实现方式如下，

根据下式，

其中，分别代表俯仰角和侧滚角由于平台振荡引起的角度变化，a是沿轨向归化系数，是沿轨向平移量；b是垂轨向归化系数，是垂轨向平移量；

利用姿态文件中各姿态点对应的时间t参量计算并让姿态文件中的姿态点中的俯仰角和侧滚角分别减去

3.一种基于几何模型的LASIS波段影像配准系统，其特征在于：包括以下模块，

第一模块，用于将LASIS面阵探测器等效为多线阵探测器，建立多线阵推扫成像模型，所述多线阵推扫成像模型中的每条垂轨向CCD模型符合单线阵推扫模型，线阵之间无缝排列，每个线阵成像对应LASIS一个单波段影像，其中各个波段p的严密成像模型如下，

$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\left[\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right]+{\mathrm{mR}}_{BW}{R}_{CB}{R}_U\left[\begin{array}{c}tan{\left({\mathrm{\ψ}}_x\right)}_p\\ tan{\left({\mathrm{\ψ}}_y\right)}_p\\ -1\end{array}\right]---\left(1\right)$

式中，[XYZ]^T为像点对应的地面点的WGS84空间直角坐标系的坐标值，m为比例系数，[X_SY_SZ_S]^T为成像时卫星在WGS84坐标系下的位置矢量，R_BW为本体坐标系相对WGS84坐标系的旋转矩阵，R_CB为相机坐标系相对本体坐标系的旋转矩阵，R_U代表偏置矩阵，(ψ_x)_p,(ψ_y)_p分别为第p行CCD任一探元对应的沿轨向指向角和垂轨向指向角，计算如下，

$\begin{array}{c}\tan {\left({\mathrm{\ψ}}_x\right)}_p=\frac{x-x_0-\mathrm{\Δ}x}{-f}\\ \tan {\left({\mathrm{\ψ}}_y\right)}_p=\frac{y-y_0-\mathrm{\Δ}y}{-f}\end{array}---\left(2\right)$

式中，(x,y)为像点影像坐标、(x₀,y₀,f)表示内方位元素，(Δx,Δy)为相机畸变引起的像点偏移；

第二模块，用于根据控制数据进行相机标定，包括利用Harris算法提取基准波段影像角点，求解物方信息；利用最小二乘配准算法在非基准波段影像获取同名点作为对应的像方信息，基于相应波段的严密成像模型利用同名点求解指向角，并把指向角和对应的垂轨向坐标s代入如下公式，利用最小二乘计算拟合的系数a_t,b_t，0≤t≤5，

$\left\{\begin{array}{c}tan\left({\mathrm{\ψ}}_x\right)=a_0+a_{1}s+a_2{s}^2+a_3{s}^3+a_4{s}^4+a_5{s}^5\\ tan\left({\mathrm{\ψ}}_y\right)=b_0+b_{1}s+b_2{s}^2+b_3{s}^3+b_4{s}^4+b_5{s}^5\end{array}\right.---\left(3\right)$

对姿轨模型进行平滑，并按照平滑后的姿轨模型对全部波段影像重采样，从而实现影像成像范围的一致性；

第三模块，用于利用不同波段影像重叠区域构建平行观测值，利用两个波段的平行观测量检测平台颤振，并根据利用检测到的颤振特征补偿其他波段颤振造成的像素偏移，获取精确像点坐标和对应的控制点坐标；实现方式如下，

设基准波段的颤振模型为

其中，V_x、V_y分别代表沿轨向像差与垂轨像像差在沿轨向的投影；ε_x、ε_y代表沿轨向像差与垂轨向像差中的随机成分；代表某个沿轨向颤振分量的振幅，频率和相位；(B_jω_jψ_j)代表某个垂轨向颤振分量的振幅，频率和相位；i代表某个沿轨向颤振成分的序号，j代表某个垂轨向颤振成分的序号，m代表沿轨向颤振成分的个数，n代表垂轨向颤振成分的个数；

构建基准波段影像的平行观测量，其他波段的颤振模型为

其中，V’_x、V’_y分别代表沿轨向像差与垂轨像像差在沿轨向的投影；ε’_x、ε’_y代表沿轨向像差与垂轨向像差中的随机成分；代表某个沿轨向颤振分量的振幅，频率和相位；(B'_jω'_jψ'_j)代表某个垂轨向颤振分量的振幅，频率和相位；i代表某个沿轨向颤振成分的序号，j代表某个沿轨向颤振成分的序号，m代表沿轨向颤振成分的个数，n代表垂轨向颤振成分的个数；

得到两波段同名点配准残差中颤振影响对应的公式如下，

其中，ΔV_x、ΔV_y代表某一波段与基准波段配准残差的沿轨向坐标和垂轨向坐标；

设为存在以下等价条件，

利用焦面排列特点和成像模型根据波段成像时间差，改变参量t，得到其他波段与参考波段由于颤振造成的像差量；

第四模块，用于利用第三模块所得结果重新计算控制点对应的姿态信息，更新其对应的相机坐标系相对本体坐标系的旋转矩阵R_CB；

第五模块，用于利用第四模块更新结果进一步精化式(1)所示成像模型，利用精化所得新的成像模型与无畸变的成像模型对准，对原始影像进行重采样，实现波段影像的自动配准。

4.根据权利要求3基于几何模型的LASIS波段影像配准系统，其特征在于：第五模块中，利用第四模块更新结果进一步精化成像模型实现方式如下，

根据下式，

其中，分别代表俯仰角和侧滚角由于平台振荡引起的角度变化，a是沿轨向归化系数，是沿轨向平移量；b是垂轨向归化系数，是垂轨向平移量；

利用姿态文件中各姿态点对应的时间t参量计算并让姿态文件中的姿态点中的俯仰角和侧滚角分别减去