CN103323028A - 一种基于物方定位一致性的卫星多光谱影像配准方法 - Google Patents

Abstract

一种基于物方定位一致性的卫星多光谱影像配准方法，包括选定多光谱相机中某一谱段作为参考谱段，其余谱段为非参考谱段，针对各非参考谱段与参考谱段之间的相对几何畸变分别进行在轨检校，并将检校结果保存；基于物方定位一致性，利用保存的校验结果建立各谱段的严格几何成像模型，将非参考谱段与参考谱段进行精确配准。本发明提供了一种真正几何意义上的光学卫星多光谱影像高精度自动配准方法，能够在无需影像匹配的情况下对卫星多光谱影像进行高精度自动配准，不仅提高了处理效率，并且配准质量与影像质量、地物类型无关；同时，波段间的几何纠正模型为严格几何成像模型，在理论上具有严密性。实践表明该方法可行、有效，配准质量稳定、精度高。

Description

一种基于物方定位一致性的卫星多光谱影像配准方法

技术领域

本发明属于遥感影像几何处理领域，涉及一种基于物方定位一致性的卫星多光谱影像配准方法。

背景技术

多光谱相机由于能够获取多个波段的影像，通过后续的配准、融合等处理，能够生成各种专题影像产品，极大地丰富了影像信息，提升了影像数据的应用潜力，已成为当前遥感卫星搭载的重要成像载荷。各波段影像之间的配准是多光谱影像处理的重要环节，配准效率、精度以及可靠性直接影响其后续处理和应用的质量。

目前，基于影像匹配的卫星多光谱影像配准方法已经开展了大量的研究与实践工作，这类方法对匹配质量的依赖使其适用性受到一定限制，对于那些纹理特征不明显、辐射差异较大的区域，此类方法的配准质量难以保证；另外，利用该类方法每次进行影像配准前，都必须进行影像匹配，其配准效率大大降低，严重影响后续的处理与应用。基于此，在卫星影像地面预处理中，有必要研究一种不依赖影像匹配的多光谱影像自动配准方法，避免上述问题的出现。

发明内容

本发明所要解决的问题是提供一种无需影像匹配的卫星多光谱影像自动配准方法。

本发明的技术方案为一种基于物方定位一致性的卫星多光谱影像配准方法，包括以下步骤：

步骤1，选定多光谱相机中某一谱段作为参考谱段，其余谱段为非参考谱段，针对各非参考谱段与参考谱段之间的相对几何畸变分别进行在轨检校，并将检校结果保存；

设参考谱段记为B2，任一非参考谱段记为B1，针对非参考谱段B1与参考谱段B2之间的相对几何畸变进行在轨检校实现方式如下，

步骤1.1，设为物方点P在非参考谱段B1与参考谱段B2影像上有同名像点p₁和p₂，

构建参考谱段B2的严格几何成像模型如下，

$\left(\begin{array}{c}{x}_{c2}\\ y_{c2}\\ f\end{array}\right)=m_2{R}_{\mathrm{BS}}{R}_{\mathrm{BJ}2}{R}_{T2}{\left[\begin{array}{c}\mathrm{Xp}-X_{S2}\\ \mathrm{Yp}-Y_{S2}\\ \mathrm{Zp}-Z_{S2}\end{array}\right]}_{\mathrm{WGS}84}$

构建非参考谱段B1的严格几何成像模型如下，

$\left(\begin{array}{c}{x}_{c1}+\mathrm{\Δ}{x}_1\\ y_{c1}+\mathrm{\Δ}{y}_1\\ f\end{array}\right)=m_1{R}_{\mathrm{BS}}{R}_{\mathrm{BJ}1}{R}_{T1}{\left[\begin{array}{c}\mathrm{Xp}-X_{S1}\\ \mathrm{Yp}-Y_{S1}\\ \mathrm{Zp}-Z_{S1}\end{array}\right]}_{\mathrm{WGS}84}$

上式中，(Xp,Yp,Zp)_WGS84为物方点P的WGS84地心直角坐标；(x_C1,y_C1,-f)和(x_C2,y_C2,-f)分别为同名像点p₁和p₂在相机坐标系下的坐标，f代表相机主距；m₁和m₂为摄影比例尺因子；(X_S1,Y_S1,Z_S1)和(X_S2,Y_S2,Z_S2)为投影中心S1和S2的WGS84地心直角坐标；R_BS为相机在卫星本体坐标系下的安装角矩阵；R_BJ1和R_T1分别为p₁成像时卫星本体坐标系与地心惯性坐标系、地心惯性坐标系与WGS84地心直角坐标系之间的旋转矩阵；R_BJ2和R_T2则为p₂成像时相应的旋转矩阵；对于附加参数项Δx₁和Δy₁，采用以探元号为自变量的三次多项式如下，

$\left\{\begin{array}{c}\mathrm{\Δ}{x}_1={\mathrm{ax}}_0^1+{\mathrm{ax}}_1^1\×s+{\mathrm{ax}}_2^1\×s^2+{\mathrm{ax}}_3^1\×s^3\\ \mathrm{\Δ}{y}_1={\mathrm{ay}}_0^1+{\mathrm{ay}}_1^1\×s+{\mathrm{ay}}_2^1\×s^2+{\mathrm{ay}}_3^1\×s^3\end{array}\right.$

其中，

为三次多项式的系数，s代表探元号；

将参考谱段B2相对于非参考谱段B1的几何畸变参数记为，则

$X_I^1=({\mathrm{ax}}_0^1,{\mathrm{ax}}_1^{1,}{\mathrm{ax}}_2^1,{\mathrm{ax}}_3^1,{\mathrm{ay}}_0^1,{\mathrm{ay}}_1^1,{\mathrm{ay}}_2^1,{\mathrm{ay}}_3^1)$

采用符号f₁表示坐标正投影换算，建立将像点(x,y)正投影至物方获取其对应的物方点的WGS84地心直角坐标(X,Y,Z)_WGS84的公式如下，

采用符号f₂表示坐标反投影换算，建立将物方点的WGS84地心直角坐标(X,Y,Z)_WGS84反投影至像方所得到的像点坐标(x,y)的公式如下，

步骤1.2，在非参考谱段B1与参考谱段B2影像上量测n对同名像点

同名像点

对应物方点Pⁱ，n为同名像点对的总数，i＝1,…,n；对每个同名点对根据参考谱段B2影像上的像点

的坐标，利用参考谱段B2的严格几何成像模型及物方高程信息，执行坐标正投影换算f₁，将像点

投影至物方，获取相应物方点Pⁱ的WGS84地心直角坐标 ${(X_{P^i},Y_{P^i},Z_{P^i})}_{\mathrm{WGS}84};$

步骤1.3，利用步骤1.2所得物方点Pⁱ的坐标

基于非参考谱段B1的严格几何成像模型，利用空间后方交会的原理解算参考谱段B2相对于非参考谱段B1的几何畸变参数

消除非参考谱段B1与参考谱段B2之间的相对几何畸变；

步骤2，基于物方定位一致性，利用各谱段的严格几何成像模型，将非参考谱段与参考谱段进行精确配准；针对非参考谱段B1与参考谱段B2进行精确配准实现方式如下，

步骤2.1，根据参考谱段影像像点坐标和物方高程信息获取物方点坐标，包括对参考谱段B2影像上的每个像元p₂(x_c2,y_c2)执行坐标正投影换算f₁，获取其物方点P的WGS84地心直角坐标(X_P,Y_P,Z_P)_WGS84；

步骤2.2，获取非参考谱段影像像点坐标，包括利用步骤1所得参考谱段B2相对于非参考谱段B1的几何畸变参数

构建非参考谱段B1的严格几何成像模型；对步骤2.1所得物方点P的坐标(X_P,Y_P,Z_P)_WGS84执行坐标反投影计算f₂，获得非参考谱段B1影像上对应的像点p₁的坐标(x_c1,y_c1)；

步骤2.3，根据步骤2.2所得非参考谱段B1影像上对应的像点p₁的坐标(x_c1,y_c1)进行灰度重采样，完成非参考谱段B1与参考谱段B2的精确配准。

而且，坐标正投影换算和坐标反投影换算的实现方式如下，

设参考谱段B2影像上像点p₂的坐标为(x_c2,y_c2)，对应的物方点P的WGS84地心直角坐标为(Xp,Yp,Zp)_WGS84，令

$R_{\mathrm{BS}}{R}_{\mathrm{BJ}2}{R}_{T2}=\left[\begin{array}{ccc}{a}_1& b_1& c_1\\ a_2& b_2& c_2\\ a_3& b_3& c_3\end{array}\right]$

a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃为矩阵的元素；

根据物方点P的WGS84地心直角坐标(Xp,Yp,Zp)_WGS84与其大地坐标(Bp,Lp,Hp)之间的如下关系，

${\left[\begin{array}{c}\mathrm{Xp}\\ \mathrm{Yp}\\ \mathrm{Zp}\end{array}\right]}_{\mathrm{WGS}84}=\left[\begin{array}{c}(N+\mathrm{Hp})\·\cos \mathrm{Bp}\·\cos \mathrm{Lp}\\ (N+\mathrm{Hp})\·\cos \mathrm{Bp}\·\sin \mathrm{Lp}\\ (N\·(1-e^2)+\mathrm{Hp})\·\sin \mathrm{Bp}\end{array}\right]$

其中，e代表地球椭球扁率，变量

a代表地球椭球长半轴；

得到像点p₂的坐标(x_c2,y_c2)与大地坐标(Bp,Lp,Hp)关系式如下，

$\left\{\begin{array}{c}(N+\mathrm{Hp})\·\cos \mathrm{Bp}\·\cos \mathrm{Lp}=\frac{a_1{x}_{c2}+a_2{x}_{c2}+a_{3}f}{c_1{x}_{c2}+c_2{x}_{c2}+c_{3}f}\·((N\·(1-e^2)+\mathrm{Hp})\·\sin \mathrm{Bp}-Z_S)+X_S\\ (N+\mathrm{Hp})\·\cos \mathrm{Bp}\·\sin \mathrm{Lp}=\frac{b_1{x}_{c2}+b_1{x}_{c2}+b_{3}f}{c_1{x}_{c2}+c_2{x}_{c2}+c_{3}f}\·((N\·(1-e^2)+\mathrm{Hp})\·\sin \mathrm{Bp}-Z_S)+Y_S\end{array}\right.$

坐标正投影换算时，由像点p₂的坐标(x_c2,y_c2)以及给定的物方高程值Hp，利用上式解算物方点P的大地经纬度(Bp,Lp)；再根据物方点P的大地坐标(Bp,Lp,Hp)获取相应WGS84地心直角坐标(Xp,Yp,Zp)_WGS84；

坐标反投影换算时，基于参考谱段B2的严格几何成像模型，由物方点P的WGS84地心直角坐标(Xp,Yp,Zp)_WGS84，利用下式解算，

$\left\{\begin{array}{c}{x}_{c2}=\frac{a_1(\mathrm{Xp}-X_{S2})+b_1(\mathrm{Yp}-Y_{S2})+c_1(\mathrm{Zp}-Z_{S2})}{a_3(\mathrm{Xp}-X_{S2})+b_3(\mathrm{Yp}-Y_{S2})+c_3(\mathrm{Zp}-Z_{S2})}\·f\\ y_{c2}=\frac{a_2(\mathrm{Xp}-X_{S2})+b_2(\mathrm{Yp}-Y_{S2})+c_2(\mathrm{Zp}-Z_{S2})}{a_3(\mathrm{Xp}-X_{S2})+b_3(\mathrm{Yp}-Y_{S2})+c_3(\mathrm{Zp}-Z_{S2})}\·f\end{array}\right.$

得到对应的像点p₂的坐标(x_c2,y_c2)。

而且，步骤1.3解算参考谱段B2相对于非参考谱段B1的几何畸变参数

的实现方式如下，

令

$\left[\begin{array}{c}\mathrm{Ux}\\ \mathrm{Uy}\\ \mathrm{Uz}\end{array}\right]=R_{\mathrm{BS}}{R}_{\mathrm{BJ}1}{R}_{T1}{\left[\begin{array}{c}\mathrm{Xp}-X_{S1}\\ \mathrm{Yp}-Y_{S1}\\ \mathrm{Zp}-Z_{S1}\end{array}\right]}_{\mathrm{WGS}84}$

上式中，矢量 $\left[\begin{array}{c}\mathrm{Ux}\\ \mathrm{Uy}\\ \mathrm{Uz}\end{array}\right]$ 代表从相机投影中心到物方点的矢量在相机坐标系下的坐标；

根据非参考谱段B1的严格几何成像模型，得到下式

$\left\{\begin{array}{c}(x_{c1}+{\mathrm{\Δx}}_1)-\frac{\mathrm{Ux}\·f}{\mathrm{Uz}}=0\\ (y_{c1}+{\mathrm{\Δy}}_1)-\frac{\mathrm{Uy}\·f}{\mathrm{Uz}}=0\end{array}\right.$

令

$\left\{\begin{array}{c}{v}_{\mathrm{xi}}=(x_{c1}+{\mathrm{\Δx}}_1)-\frac{\mathrm{Ux}\·f}{\mathrm{Uz}}\\ v_{\mathrm{yi}}=(y_{c1}+\mathrm{\Δ}{y}_1)-\frac{\mathrm{Uy}\·f}{\mathrm{Uz}}\end{array}\right.$

上式中，v_xi和v_yi分别代表沿轨和垂轨方向的像方残差；

将步骤1.2中所得物方点Pⁱ坐标

代入上式中，对每个物方点Pⁱ构建如如下误差方程，

V_i＝A_iX-L_i，W_i，

其中，

$V_i=\left[\begin{array}{c}{v}_{\mathrm{xi}}\\ v_{\mathrm{yi}}\end{array}\right]$ $A_i=\left[\begin{array}{cccccccc}1& s& s^2& s^3& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& s& s^2& s^3\end{array}\right]$

$L_i=\left[\begin{array}{c}{(\frac{\mathrm{Ux}\·f}{\mathrm{Uz}}-x_{c1})}_i\\ {(\frac{\mathrm{Uy}\·f}{\mathrm{Uz}}-y_{c1})}_i\end{array}\right]$

$X={\left(X_I^1\right)}^T={({\mathrm{ax}}_0^1,{\mathrm{ax}}_1^1,{\mathrm{ax}}_2^1,{\mathrm{ax}}_3^1,{\mathrm{ay}}_0^1,{\mathrm{ay}}_1^1,{\mathrm{ay}}_2^1,{\mathrm{ay}}_3^1)}^T$

上式中，V_i、A_i、L_i分别是利用物方点Pⁱ构建的误差方程，的残差向量、待解参数的系数矩阵以及常向量；X代表参考谱段B2相对于非参考谱段B1的几何畸变参数 ${\left(X_I^1\right)}^T={({\mathrm{ax}}_0^1,{\mathrm{ax}}_1^1,{\mathrm{ax}}_2^1,{\mathrm{ax}}_3^1,{\mathrm{ay}}_0^1,{\mathrm{ay}}_1^1,{\mathrm{ay}}_2^1,{\mathrm{ay}}_3^1)}^T;$ W_i是非参考B1谱段影像上的像点的量测精度对应的权；

基于最小二乘平差原理，利用下式计算参考谱段B2相对于非参考谱段B1的几何畸变参数，

$X={\left(X_I^1\right)}^T={\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{{A_i}^{T}W}_i{A}_i\right)}^{-1}\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{{A_i}^{T}W}_i{L}_i\right)$

记录计算所得几何畸变参数。

本发明选定多光谱相机中某一谱段作为参考谱段，其余谱段为非参考谱段。利用一景质量较优的影像，基于各谱段严格几何成像模型，仅利用谱段间的同名像点，对各非参考谱段其与参考谱段之间的相对几何畸变进行在轨检校；对非参考波段其与参考波段间的相对几何畸变进行在轨检校后，基于同名像元物方定位一致性的约束条件，利用各谱段严格几何成像模型，在无需影像匹配的情况下实现子像素级的多光谱影像自动配准，波段间的几何纠正模型为严格几何成像模型，理论上具有严密性。本发明的优点在于：1.影像配准之前无需进行影像匹配，一方面大大提升了数据处理效率，同时，配准质量与影像质量、地物类型无关，对于水域、沙漠以及山地等纹理特征不丰富、匹配质量难以保证的区域，本发明的配准质量也能得到保证。2.一种真正几何意义上的配准，各波段影像之间的几何纠正模型为严格几何成像模型，在理论上具有严密性，能够进一步保证配准精度。

附图说明

图1为本发明实施例的配准流程示意图；

图2为本发明基于物方定位一致性的卫星多光谱影像配准原理示意图。

具体实施方式

以下结合附图和实施例详细说明本发明具体实施方式。实施例的流程可以分为两个步骤，每个步骤实施的具体方法、公式以及流程如下：

步骤1，选定多光谱相机中某一谱段作为参考谱段，其余谱段为非参考谱段，利用一景质量较优的影像，针对各非参考谱段，对其与参考谱段之间的相对几何畸变进行在轨检校，并将检校结果保存，用于后续的影像配准。设参考谱段记为B2，任一非参考谱段记为B1，针对非参考谱段B1，其与参考谱段B2之间的相对几何畸变在轨检校的具体步骤及公式如下：

步骤1.1，利用姿态、轨道、时间等辅助数据以及相机参数，构建参考谱段B2的严格几何成像模型；如式（1）；

$\left(\begin{array}{c}{x}_{c2}\\ y_{c2}\\ f\end{array}\right)=m_2{R}_{\mathrm{BS}}{R}_{\mathrm{BJ}2}{R}_{T2}{\left[\begin{array}{c}\mathrm{Xp}-X_{S2}\\ \mathrm{Yp}-Y_{S2}\\ \mathrm{Zp}-Z_{S2}\end{array}\right]}_{\mathrm{WGS}84}---\left(1\right)$

同理，构建B1谱段的严格几何成像模型，并在其内定向参数模型中引入附加参数项Δx₁和Δy₁，用于补偿B1谱段相对于B2谱段的几何畸变，实现构建基于扩展共线条件方程的自检校平差模型，如式（2）。

$\left(\begin{array}{c}{x}_{c1}+\mathrm{\Δ}{x}_1\\ y_{c1}+\mathrm{\Δ}{y}_1\\ f\end{array}\right)=m_1{R}_{\mathrm{BS}}{R}_{\mathrm{BJ}1}{R}_{T1}{\left[\begin{array}{c}\mathrm{Xp}-X_{S1}\\ \mathrm{Yp}-Y_{S1}\\ \mathrm{Zp}-Z_{S1}\end{array}\right]}_{\mathrm{WGS}84}---\left(2\right)$

上式中，(Xp,Yp,Zp)_WGS84为物方点P的WGS84地心直角坐标；(x_C1,y_C1,-f)和(x_C2,y_C2,-f)分别为物方点P在B1谱段与B2谱段影像上同名像点p₁和p₂在相机坐标系下的坐标，f代表相机主距；m₁和m₂为摄影比例尺因子；(X_S1,Y_S1,Z_S1)和(X_S2,Y_S2,Z_S2)为投影中心S1和S2的WGS84地心直角坐标；R_BS为相机在卫星本体坐标系下的安装角矩阵；R_BJ1和R_T1分别为p₁成像时卫星本体坐标系与地心惯性坐标系、地心惯性坐标系与WGS84地心直角坐标系之间的旋转矩阵，地心惯性坐标系为J2000坐标系等；R_BJ2和R_T2则为p₂成像时相应的旋转矩阵。对于B1谱段成像模型中的附加参数项Δx₁和Δy₁，采用以探元号为自变量的三次多项式，如式（3）所示：

$\left\{\begin{array}{c}\mathrm{\Δ}{x}_1={\mathrm{ax}}_0^1+{\mathrm{ax}}_1^1\×s+{\mathrm{ax}}_2^1\×s^2+{\mathrm{ax}}_3^1\×s^3\\ \mathrm{\Δ}{y}_1={\mathrm{ay}}_0^1+{\mathrm{ay}}_1^1\×s+{\mathrm{ay}}_2^1\×s^2+{\mathrm{ay}}_3^1\×s^3\end{array}\right.---\left(3\right)$

其中，

为三次多项式的系数，s代表探元号。

将B1谱段相对于B2谱段的几何畸变参数记为

（上标1代表谱段B1），则：

$X_I^1=({\mathrm{ax}}_0^1,{\mathrm{ax}}_1^1,{\mathrm{ax}}_2^1,{\mathrm{ax}}_3^1,{\mathrm{ay}}_0^1,{\mathrm{ay}}_1^1,{\mathrm{ay}}_2^1,{\mathrm{ay}}_3^1)$

利用卫星获取的姿态、轨道、时间、相机参数以及物方高程信息，基于严格几何成像模型，能够实现像点坐标与物方点坐标之间的换算。以谱段B2为例，对具体的换算过程及公式进行阐述。

假设谱段B2影像上某像点p₂的坐标为(x_c2,y_c2)，其对应的物方点P的WGS84地心直角坐标为(Xp,Yp,Zp)_WGS84，令公式（1）中：

$R_{\mathrm{BS}}{R}_{\mathrm{BJ}2}{R}_{T2}=\left[\begin{array}{ccc}{a}_1& b_1& c_1\\ a_2& b_2& c_2\\ a_3& b_3& c_3\end{array}\right]$

a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃为矩阵的元素。

则利用公式（4）及物方点P的WGS84地心直角坐标(Xp,Yp,Zp)_WGS84可解算像点p₂的坐标(x_c2,y_c2)：

$\left\{\begin{array}{c}{x}_{c2}=\frac{a_1(\mathrm{Xp}-X_{S2})+b_1(\mathrm{Yp}-Y_{S2})+c_1(\mathrm{Zp}-Z_{S2})}{a_3(\mathrm{Xp}-X_{S2})+b_3(\mathrm{Yp}-Y_{S2})+c_3(\mathrm{Zp}-Z_{S2})}\·f\\ y_{c2}=\frac{a_2(\mathrm{Xp}-X_{S2})+b_2(\mathrm{Yp}-Y_{S2})+c_2(\mathrm{Zp}-Z_{S2})}{a_3(\mathrm{Xp}-X_{S2})+b_3(\mathrm{Yp}-Y_{S2})+c_3(\mathrm{Zp}-Z_{S2})}\·f\end{array}\right.---\left(4\right)$

根据大地测量学的基本原理可知，物方点P的WGS84地心直角坐标(Xp,Yp,Zp)_WGS84与其大地坐标(Bp,Lp,Hp)（Bp、Lp分别为物方点P的纬度和经度，Hp为其椭球高）之间有如下关系：

${\left[\begin{array}{c}\mathrm{Xp}\\ \mathrm{Yp}\\ \mathrm{Zp}\end{array}\right]}_{\mathrm{WGS}84}=\left[\begin{array}{c}(N+\mathrm{Hp})\·\cos \mathrm{Bp}\·\cos \mathrm{Lp}\\ (N+\mathrm{Hp})\·\cos \mathrm{Bp}\·\sin \mathrm{Lp}\\ (N\·(1-e^2)+\mathrm{Hp})\·\sin \mathrm{Bp}\end{array}\right]---\left(5\right)$

其中，e代表地球椭球扁率，变量a代表地球椭球长半轴。将上式（5）代入公式（4）中有公式（6）：

$\left\{\begin{array}{c}(N+\mathrm{Hp})\·\cos \mathrm{Bp}\·\cos \mathrm{Lp}=\frac{a_1{x}_{c2}+a_2{x}_{c2}+a_{3}f}{c_1{x}_{c2}+c_2{x}_{c2}+c_{3}f}\·((N\·(1-e^2)+\mathrm{Hp})\·\sin \mathrm{Bp}-Z_S)+X_S\\ (N+\mathrm{Hp})\·\cos \mathrm{Bp}\·\sin \mathrm{Lp}=\frac{b_1{x}_{c2}+b_1{x}_{c2}+b_{3}f}{c_1{x}_{c2}+c_2{x}_{c2}+c_{3}f}\·((N\·(1-e^2)+\mathrm{Hp})\·\sin \mathrm{Bp}-Z_S)+Y_S\end{array}\right.---\left(6\right)$

坐标正投影换算时，由像点p₂的坐标(x_c2,y_c2)以及给定的物方高程值Hp，利用公式（6）即可解算物方点P的大地经纬度(Bp,Lp)；将物方点P的大地坐标(Bp,Lp,Hp)代入公式（5）即可获取其WGS84地心直角坐标(Xp,Yp,Zp)_WGS84。坐标反投影换算时，基于参考谱段B2的严格几何成像模型，由物方点P的WGS84地心直角坐标(Xp,Yp,Zp)_WGS84，利用式（4）解算对应的像点p₂的坐标(x_c2,y_c2)

为了便于描述，本文利用公式（7）和（8）简要表达上述的像点坐标与物方点坐标正反换算。其中，公式（7）中的符号f₁表示将像点(x,y)正投影至物方获取其对应的物方点的WGS84地心直角坐标(X,Y,Z)_WGS84，即坐标正投影换算；公式（8）中的符号f₂则表示将物方点的WGS84地心直角坐标(X,Y,Z)_WGS84反投影至像方所得到的像点坐标(x,y)，即坐标反投影换算。

步骤1.2，在B1和B2两个谱段影像上量测n对同名像点

表示B1和B2谱段影像上的一对同名像点），同名像点对应物方点Pⁱ。n为同名像点对的总数，可由本领域技术人员根据具体情况设定。具体量测实现为现有技术，为了提高检校结果的精度以及可靠性，建议在实施时对像点

在非参考谱段影像上的分布作如下要求（对应的同名像点

在参考谱段影像上的分布形状与

在非参考谱段上是一致的）：（1）在影像行方向（即推扫方向）上尽量分布在较短的一段区域内；（2）在影像列方向（即沿CCD方向）上均匀覆盖整个线阵CCD。其中，第一点要求是为了降低外方位元素误差对解算结果的影响；第二点则是为了保证解算结果对线阵CCD所有探元均适用。对每个同名点对(i＝1,…,n)，将B2谱段影像上的像点

的坐标代入公式（7）中的(x,y)，利用式（1）所示B2谱段的严格几何成像模型及物方高程信息，执行坐标正投影换算f₁，将像点

投影至物方，获取其物方点Pⁱ的WGS84地心直角坐标即公式（7）中的(X,Y,Z)_WGS84；

步骤1.3，利用前述步骤1.2得到的物方点Pⁱ的坐标(i＝1,…,n)，基于B1谱段的严格几何成像模型（式2），利用空间后方交会的原理解算，消除B1谱段与B2谱段影像之间的相对几何畸变。具体解算公式如下。

1）令

$\left[\begin{array}{c}\mathrm{Ux}\\ \mathrm{Uy}\\ \mathrm{Uz}\end{array}\right]=R_{\mathrm{BS}}{R}_{\mathrm{BJ}1}{R}_{T1}{\left[\begin{array}{c}\mathrm{Xp}-X_{S1}\\ \mathrm{Yp}-Y_{S1}\\ \mathrm{Zp}-Z_{S1}\end{array}\right]}_{\mathrm{WGS}84}---\left(9\right)$

上式中，矢量 $\left[\begin{array}{c}\mathrm{Ux}\\ \mathrm{Uy}\\ \mathrm{Uz}\end{array}\right]$ 称为物方矢量U，代表从相机投影中心到物方点的矢量在相机坐标系下的坐标；

式（2）可转化为式（10）：

$\left\{\begin{array}{c}(x_{c1}+{\mathrm{\Δx}}_1)-\frac{\mathrm{Ux}\·f}{\mathrm{Uz}}=0\\ (y_{c1}+\mathrm{\Δ}{y}_1)-\frac{\mathrm{Uy}\·f}{\mathrm{Uz}}=0\end{array}\right.---\left(10\right)$

令

$\left\{\begin{array}{c}{v}_{\mathrm{xi}}=(x_{c1}+{\mathrm{\Δx}}_1)-\frac{\mathrm{Ux}\·f}{\mathrm{Uz}}\\ v_{\mathrm{yi}}=(y_{c1}+\mathrm{\Δ}{y}_1)-\frac{\mathrm{Uy}\·f}{\mathrm{Uz}}\end{array}\right.---\left(11\right)$

上式中，v_xi和v_yi分别代表沿轨和垂轨方向的像方残差。

2）将步骤1.2中解算的物方点Pⁱ坐标

(i＝1,…,n)代入式（11）中，对每个物方点Pⁱ均可构建如式（12）所示的误差方程（下标i代表利用物方点Pⁱ建立的误差方程）：

V_i＝A_iX-L_i，W_i(i＝1,…,n)         （12）

其中，

$V_i=\left[\begin{array}{c}{v}_{\mathrm{xi}}\\ v_{\mathrm{yi}}\end{array}\right]$ $A_i=\left[\begin{array}{cccccccc}1& s& s^2& s^3& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& s& s^2& s^3\end{array}\right]$

$L_i=\left[\begin{array}{c}{(\frac{\mathrm{Ux}\·f}{\mathrm{Uz}}-x_{c1})}_i\\ {(\frac{\mathrm{Uy}\·f}{\mathrm{Uz}}-y_{c1})}_i\end{array}\right]$

$X={\left(X_I^1\right)}^T={({\mathrm{ax}}_0^1,{\mathrm{ax}}_1^1,{\mathrm{ax}}_2^1,{\mathrm{ax}}_3^1,{\mathrm{ay}}_0^1,{\mathrm{ay}}_1^1,{\mathrm{ay}}_2^1,{\mathrm{ay}}_3^1)}^T$

上式中，V_i、A_i、L_i分别是利用物方点Pⁱ构建的误差方程式的残差向量、待解参数的系数矩阵以及常向量；X代表B1谱段相对于B2谱段的几何畸变参数 ${\left(X_I^1\right)}^T={({\mathrm{ax}}_0^1,{\mathrm{ax}}_1^1,{\mathrm{ax}}_2^1,{\mathrm{ax}}_3^1,{\mathrm{ay}}_0^1,{\mathrm{ay}}_1^1,{\mathrm{ay}}_2^1,{\mathrm{ay}}_3^1)}^T;$ W_i是B1谱段影像上的像点的量测精度对应的权。

3）基于最小二乘平差原理，利用式（13）计算B1谱段相对于B2谱段的几何畸变参数，实现对B1谱段相对于B2谱段之间的相对几何畸变进行在轨检校，并将几何畸变参数用文件记录下来，用于后续多光谱影像波段配准。

$X={\left(X_I^1\right)}^T={\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{{A_i}^{T}W}_i{A}_i\right)}^{-1}\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{A_i}^T{W}_i{L}_i\right)---\left(13\right)$

步骤2，通过步骤1预先获取了各非参考谱段其相对于参考谱段的几何畸变参数后，基于物方定位一致性，利用各谱段影像的严格几何成像模型，对任意一景多光谱影像，将其各非参考谱段与参考谱段分别进行精确配准。仍设参考谱段记为B2，任一非参考谱段记为B1，结合图1和图2对本步骤中的具体原理和流程进行阐述如下。其中，图1为流程示意图，图2为原理示意图。图2中，p₁和p₂分别为B1和B2谱段影像上的一对同名像点，P则代表其对应的物方点；S1和S2分别为像点p₁和p₂成像时的相机投影中心；f₁和f₂则代表公式7和8中的像点坐标与物方点坐标之间的正反换算。

步骤2.1，根据参考谱段影像像点坐标和物方高程信息获取物方点坐标：对B2谱段影像上的每个像元p₂(x_c2,y_c2)，将其像点坐标(x_c2,y_c2)代入公式（7）中的(x,y)，即利用B2谱段的严格几何成像模型（式1）及物方高程信息，执行坐标正投影换算f₁，将像点p₂(x_c2,y_c2)投影至物方，获取其物方点P的WGS84地心直角坐标(X_P,Y_P,Z_P)_WGS84（即公式（7）中的(X,Y,Z)_WGS84）；

步骤2.2，获取非参考谱段影像像点坐标：利用步骤1获取的B1谱段其相对于B2谱段的几何畸变参数，构建式（2）所示的严格几何成像模型。将获取的物方点P的坐标(X_P,Y_P,Z_P)_WGS84代入式（8）中的(X,Y,Z)_WGS84，执行坐标反投影计算f₂，获得B1谱段影像上对应的像点p₁的坐标(x_c1,y_c1)，即公式（8）中的(x,y)；由于式（2）已经对B1和B2两谱段的相对几何畸变进行了补偿，基于同名像点物方定位一致性的约束关系，所得到的像点p₁(s₁,l₁)即为B2谱段影像上像点p₂(s,l)在B1谱段影像上对应的同名像点，由此可建立B1、B2两波段影像上同名像点之间的映射关系，实现多光谱影像的配准，这就是基于物方定位一致性的多光谱影像配准原理，如图2所示。

步骤2.3，灰度重采样：在步骤2.1获取B2谱段影像的每个像元p₂(s,l)在B1谱段影像上对应的同名像点坐标p₁(s₁,l₁)后，通过灰度重采样，即可实现B1与B2两谱段影像的配准。其中，灰度重采样可采用现有技术中的双线性内插算法实现，本发明不予赘述。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (3)

1.一种基于物方定位一致性的卫星多光谱影像配准方法，包括以下步骤：

步骤1，选定多光谱相机中某一谱段作为参考谱段，其余谱段为非参考谱段，针对各非参考谱段与参考谱段之间的相对几何畸变分别进行在轨检校，并将检校结果保存；

设参考谱段记为B2，任一非参考谱段记为B1，针对非参考谱段B1与参考谱段B2之间的相对几何畸变进行在轨检校实现方式如下，

步骤1.1，设为物方点P在非参考谱段B1与参考谱段B2影像上有同名像点p₁和p₂，

构建参考谱段B2的严格几何成像模型如下，

$\left(\begin{array}{c}{x}_{c2}\\ y_{c2}\\ f\end{array}\right)=m_2{R}_{\mathrm{BS}}{R}_{\mathrm{BJ}2}{R}_{T2}{\left[\begin{array}{c}\mathrm{Xp}-X_{S2}\\ \mathrm{Yp}-Y_{S2}\\ \mathrm{Zp}-Z_{S2}\end{array}\right]}_{\mathrm{WGS}84}$

构建非参考谱段B1的严格几何成像模型如下，

$\left(\begin{array}{c}{x}_{c1}+\mathrm{\Δ}{x}_1\\ y_{c1}+\mathrm{\Δ}{y}_1\\ f\end{array}\right)=m_1{R}_{\mathrm{BS}}{R}_{\mathrm{BJ}1}{R}_{T1}{\left[\begin{array}{c}\mathrm{Xp}-X_{S1}\\ \mathrm{Yp}-Y_{S1}\\ \mathrm{Zp}-Z_{S1}\end{array}\right]}_{\mathrm{WGS}84}$

上式中，(Xp,Yp,Zp)_WGS84为物方点P的WGS84地心直角坐标；(x_C1,y_C1,-f)和(x_C2,y_C2,-f)分别为同名像点p₁和p₂在相机坐标系下的坐标，f代表相机主距；m₁和m₂为摄影比例尺因子；(X_S1,Y_S1,Z_S1)和(X_S2,Y_S2,Z_S2)为投影中心S1和S2的WGS84地心直角坐标；R_BS为相机在卫星本体坐标系下的安装角矩阵；R_BJ1和R_T1分别为p₁成像时卫星本体坐标系与地心惯性坐标系、地心惯性坐标系与WGS84地心直角坐标系之间的旋转矩阵；R_BJ2和R_T2则为p₂成像时相应的旋转矩阵；对于附加参数项Δx₁和Δy₁，采用以探元号为自变量的三次多项式如下，

$\left\{\begin{array}{c}\mathrm{\Δ}{x}_1={\mathrm{ax}}_0^1+{\mathrm{ax}}_1^1\×s+{\mathrm{ax}}_2^1\×s^2+{\mathrm{ax}}_3^1\×s^3\\ \mathrm{\Δ}{y}_1={\mathrm{ay}}_0^1+{\mathrm{ay}}_1^1\×s+{\mathrm{ay}}_2^1\×s^2+{\mathrm{ay}}_3^1\×s^3\end{array}\right.$

其中，

为三次多项式的系数，s代表探元号；

将参考谱段B2相对于非参考谱段B1的几何畸变参数记为

则

$X_I^1=({\mathrm{ax}}_0^1,{\mathrm{ax}}_1^1,{\mathrm{ax}}_2^1,{\mathrm{ax}}_3^1,{\mathrm{ay}}_0^1,{\mathrm{ay}}_1^1,{\mathrm{ay}}_2^1,{\mathrm{ay}}_3^1)$

采用符号f₁表示坐标正投影换算，建立将像点(x,y)正投影至物方获取其对应的物方点的WGS84地心直角坐标(X,Y,Z)_WGS84的公式如下，

采用符号f₂表示坐标反投影换算，建立将物方点的WGS84地心直角坐标(X,Y,Z)_WGS84反投影至像方所得到的像点坐标(x,y)的公式如下，

步骤1.2，在非参考谱段B1与参考谱段B2影像上量测n对同名像点同名像点

对应物方点Pⁱ，n为同名像点对的总数，i＝1,…,n；对每个同名点对

根据参考谱段B2影像上的像点

的坐标，利用参考谱段B2的严格几何成像模型及物方高程信息，执行坐标正投影换算f₁，将像点

投影至物方，获取相应物方点Pⁱ的WGS84地心直角坐标 $(X_{P^i},Y_{I^i}{Z}_{P^i}{)}_{\mathrm{WGS}84};$

步骤1.3，利用步骤1.2所得物方点Pⁱ的坐标

基于非参考谱段B1的严格几何成像模型，利用空间后方交会的原理解算参考谱段B2相对于非参考谱段B1的几何畸变参数

消除非参考谱段B1与参考谱段B2之间的相对几何畸变；

步骤2，基于物方定位一致性，利用各谱段的严格几何成像模型，将非参考谱段与参考谱段进行精确配准；针对非参考谱段B1与参考谱段B2进行精确配准实现方式如下，

步骤2.1，根据参考谱段影像像点坐标和物方高程信息获取物方点坐标，包括对参考谱段B2影像上的每个像元p₂(x_c2,y_c2)执行坐标正投影换算f₁，获取其物方点P的WGS84地心直角坐标(X_P,Y_P,Z_P)_WGS84；

步骤2.2，获取非参考谱段影像像点坐标，包括利用步骤1所得参考谱段B2相对于非参考谱段B1的几何畸变参数构建非参考谱段B1的严格几何成像模型；对步骤2.1所得物方点P的坐标(X_P,Y_P,Z_P)_WGS84执行坐标反投影计算f₂，获得非参考谱段B1影像上对应的像点p₁的坐标(x_c1,y_c1)；

步骤2.3，根据步骤2.2所得非参考谱段B1影像上对应的像点p₁的坐标(x_c1,y_c1)进行灰度重采样，完成非参考谱段B1与参考谱段B2的精确配准。

2.如权利要求1所述基于物方定位一致性的卫星多光谱影像配准方法，其特征在于：坐标正投影换算和坐标反投影换算的实现方式如下，

设参考谱段B2影像上像点p₂的坐标为(x_c2,y_c2)，对应的物方点P的WGS84地心直角坐标为(Xp,Yp,Zp)_WGS84，令

$R_{\mathrm{BS}}{R}_{\mathrm{BJ}2}{R}_{T2}=\left[\begin{array}{ccc}{a}_1& b_1& c_1\\ a_2& b_2& c_2\\ a_3& b_3& c_3\end{array}\right]$

a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃为矩阵的元素；

根据物方点P的WGS84地心直角坐标(Xp,Yp,Zp)_WGS84与其大地坐标(Bp,Lp,Hp)之间的如下关系，

${\left[\begin{array}{c}\mathrm{Xp}\\ \mathrm{Yp}\\ \mathrm{Zp}\end{array}\right]}_{\mathrm{WGS}84}=\left[\begin{array}{c}(N+\mathrm{Hp})\·\cos \mathrm{Bp}\·\cos \mathrm{Lp}\\ (N+\mathrm{Hp})\·\cos \mathrm{Bp}\·\sin \mathrm{Lp}\\ (N\·(1-e^2)+\mathrm{Hp})\·\sin \mathrm{Bp}\end{array}\right]$

其中，e代表地球椭球扁率，变量

a代表地球椭球长半轴；

得到像点p₂的坐标(x_c2,y_c2)与大地坐标(Bp,Lp,Hp)关系式如下，

$\left\{\begin{array}{c}(N+\mathrm{Hp})\·\cos \mathrm{Bp}\·\cos \mathrm{Lp}=\frac{a_1{x}_{c2}+a_2{x}_{c2}+a_{3}f}{c_1{x}_{c2}+c_2{x}_{c2}+c_{3}f}\·((N\·(1-e^2)+\mathrm{Hp})\·\sin \mathrm{Bp}-Z_S)+X_S\\ (N+\mathrm{Hp})\·\cos \mathrm{Bp}\·\sin \mathrm{Lp}=\frac{b_1{x}_{c2}+b_1{x}_{c2}+b_{3}f}{c_1{x}_{c2}+c_2{x}_{c2}+c_{3}f}\·((N\·(1-e^2)+\mathrm{Hp})\·\sin \mathrm{Bp}-Z_S)+Y_S\end{array}\right.$

坐标正投影换算时，由像点p₂的坐标(x_c2,y_c2)以及给定的物方高程值Hp，利用上式解算物方点P的大地经纬度(Bp,Lp)；再根据物方点P的大地坐标(Bp,Lp,Hp)获取相应WGS84地心直角坐标(Xp,Yp,Zp)_WGS84；

坐标反投影换算时，基于参考谱段B2的严格几何成像模型，由物方点P的WGS84地心直角坐标(Xp,Yp,Zp)_WGS84，利用下式解算，

$\left\{\begin{array}{c}{x}_{c2}=\frac{a_1(\mathrm{Xp}-X_{S2})+b_1(\mathrm{Yp}-Y_{S2})+c_1(\mathrm{Zp}-Z_{S2})}{a_3(\mathrm{Xp}-X_{S2})+b_3(\mathrm{Yp}-Y_{S2})+c_3(\mathrm{Zp}-Z_{S2})}\·f\\ y_{c2}=\frac{a_2(\mathrm{Xp}-X_{S2})+b_2(\mathrm{Yp}-Y_{S2})+c_2(\mathrm{Zp}-Z_{S2})}{a_3(\mathrm{Xp}-X_{S2})+b_3(\mathrm{Yp}-Y_{S2})+c_3(\mathrm{Zp}-Z_{S2})}\·f\end{array}\right.$

得到对应的像点p₂的坐标(x_c2,y_c2)。

3.如权利要求2所述基于物方定位一致性的卫星多光谱影像配准方法，其特征在于：步骤1.3解算参考谱段B2相对于非参考谱段B1的几何畸变参数

的实现方式如下，

令

$\left[\begin{array}{c}\mathrm{Ux}\\ \mathrm{Uy}\\ \mathrm{Uz}\end{array}\right]=R_{\mathrm{BS}}{R}_{\mathrm{BJ}1}{R}_{T1}{\left[\begin{array}{c}\mathrm{Xp}-X_{S1}\\ \mathrm{Yp}-Y_{S1}\\ \mathrm{Zp}-Z_{S1}\end{array}\right]}_{\mathrm{WGS}84}$

上式中，矢量 $\left[\begin{array}{c}\mathrm{Ux}\\ \mathrm{Uy}\\ \mathrm{Uz}\end{array}\right]$ 代表从相机投影中心到物方点的矢量在相机坐标系下的坐标；

根据非参考谱段B1的严格几何成像模型，得到下式

$\left\{\begin{array}{c}(x_{c1}+{\mathrm{\Δx}}_1)-\frac{\mathrm{Ux}\·f}{\mathrm{Uz}}=0\\ (y_{c1}+\mathrm{\Δ}{y}_1)-\frac{\mathrm{Uy}\·f}{\mathrm{Uz}}=0\end{array}\right.$

令

$\left\{\begin{array}{c}{v}_{\mathrm{xi}}=(x_{c1}+{\mathrm{\Δx}}_1)-\frac{\mathrm{Ux}\·f}{\mathrm{Uz}}\\ v_{\mathrm{yi}}=(y_{c1}+\mathrm{\Δ}{y}_1)-\frac{\mathrm{Uy}\·f}{\mathrm{Uz}}\end{array}\right.$

上式中，v_xi和v_yi分别代表沿轨和垂轨方向的像方残差；

将步骤1.2中所得物方点Pⁱ坐标

(i＝1,…,n)代入上式中，对每个物方点Pⁱ构建如如下误差方程，

_Vi＝A_iX-L_i，W_i，

其中，

$V_i=\left[\begin{array}{c}{v}_{\mathrm{xi}}\\ v_{\mathrm{yi}}\end{array}\right]$ $A_i=\left[\begin{array}{cccccccc}1& s& s^2& s^3& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& s& s^2& s^3\end{array}\right]$

$L_i=\left[\begin{array}{c}{(\frac{\mathrm{Ux}\·f}{\mathrm{Uz}}-x_{c1})}_i\\ {(\frac{\mathrm{Uy}\·f}{\mathrm{Uz}}-y_{c1})}_i\end{array}\right]$

$X={\left(X_I^1\right)}^T={({\mathrm{ax}}_0^1,{\mathrm{ax}}_1^1,{\mathrm{ax}}_2^1,{\mathrm{ax}}_3^1,{\mathrm{ay}}_0^1,{\mathrm{ay}}_1^1,{\mathrm{ay}}_2^1,{\mathrm{ay}}_3^1)}^T$

上式中，V_i、A_i、L_i分别是利用物方点Pⁱ构建的误差方程，的残差向量、待解参数的系数矩阵以及常向量；X代表参考谱段B2相对于非参考谱段B1的几何畸变参数 ${\left(X_I^1\right)}^T={({\mathrm{ax}}_0^1,{\mathrm{ax}}_1^1,{\mathrm{ax}}_2^1,{\mathrm{ax}}_3^1,{\mathrm{ay}}_0^1,{\mathrm{ay}}_1^1,{\mathrm{ay}}_2^1,{\mathrm{ay}}_3^1)}^T;$ W_i是非参考B1谱段影像上的像点

的量测精度对应的权；

基于最小二乘平差原理，利用下式计算参考谱段B2相对于非参考谱段B1的几何畸变参数，

$X={\left(X_I^1\right)}^T={\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{{A_i}^{T}W}_i{A}_i\right)}^{-1}\left(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{A_i}^T{W}_i{L}_i\right)$

记录计算所得几何畸变参数。

Images