CN103323028B

CN103323028B - One locates conforming satellite multispectral image method for registering based on object space

Info

Publication number: CN103323028B
Application number: CN201310236939.3A
Authority: CN
Inventors: 王密; 杨博; 金淑英; 李德仁
Original assignee: Wuhan University WHU
Current assignee: Land Sea Space Yantai Information Technology Co ltd
Priority date: 2013-06-14
Filing date: 2013-06-14
Publication date: 2015-10-21
Anticipated expiration: 2033-06-14
Also published as: CN103323028A

Abstract

A satellite multi-spectral image registration method based on object-space positioning consistency, including selecting a certain spectral segment in the multi-spectral camera as a reference spectral segment, and the remaining spectral segments as non-reference spectral segments. The relative geometric distortion between the spectrum segments is calibrated separately on the track, and the calibration results are saved; based on the consistency of object positioning, the strict geometric imaging model of each spectrum segment is established by using the saved calibration results, and the non-reference spectrum segment Accurate registration with reference spectrum. The present invention provides a high-precision automatic registration method for optical satellite multi-spectral images in a true geometric sense, which can perform high-precision automatic registration of satellite multi-spectral images without image matching, not only improving processing efficiency, but also The registration quality has nothing to do with the image quality and the type of ground objects; at the same time, the geometric correction model between bands is a strict geometric imaging model, which is rigorous in theory. Practice shows that the method is feasible and effective, with stable registration quality and high precision.

Description

A satellite multi-spectral image registration method based on object-space positioning consistency

技术领域technical field

本发明属于遥感影像几何处理领域，涉及一种基于物方定位一致性的卫星多光谱影像配准方法。The invention belongs to the field of geometric processing of remote sensing images, and relates to a satellite multispectral image registration method based on the consistency of object positioning.

背景技术Background technique

多光谱相机由于能够获取多个波段的影像，通过后续的配准、融合等处理，能够生成各种专题影像产品，极大地丰富了影像信息，提升了影像数据的应用潜力，已成为当前遥感卫星搭载的重要成像载荷。各波段影像之间的配准是多光谱影像处理的重要环节，配准效率、精度以及可靠性直接影响其后续处理和应用的质量。Because multispectral cameras can acquire images of multiple bands, through subsequent processing such as registration and fusion, they can generate various thematic image products, which greatly enriches image information and improves the application potential of image data. It has become the current remote sensing satellite An important imaging payload carried. The registration of images in various bands is an important part of multispectral image processing, and the registration efficiency, accuracy and reliability directly affect the quality of its subsequent processing and application.

目前，基于影像匹配的卫星多光谱影像配准方法已经开展了大量的研究与实践工作，这类方法对匹配质量的依赖使其适用性受到一定限制，对于那些纹理特征不明显、辐射差异较大的区域，此类方法的配准质量难以保证；另外，利用该类方法每次进行影像配准前，都必须进行影像匹配，其配准效率大大降低，严重影响后续的处理与应用。基于此，在卫星影像地面预处理中，有必要研究一种不依赖影像匹配的多光谱影像自动配准方法，避免上述问题的出现。At present, a lot of research and practical work have been carried out on the satellite multispectral image registration method based on image matching. The dependence of this kind of method on the quality of matching makes its applicability limited. The registration quality of this method is difficult to guarantee; in addition, image matching must be performed before each image registration using this method, and the registration efficiency is greatly reduced, which seriously affects subsequent processing and application. Based on this, in the ground preprocessing of satellite images, it is necessary to study an automatic multispectral image registration method that does not rely on image matching to avoid the above problems.

发明内容Contents of the invention

本发明所要解决的问题是提供一种无需影像匹配的卫星多光谱影像自动配准方法。The problem to be solved by the present invention is to provide an automatic registration method of satellite multispectral images without image matching.

本发明的技术方案为一种基于物方定位一致性的卫星多光谱影像配准方法，包括以下步骤：The technical solution of the present invention is a satellite multi-spectral image registration method based on object positioning consistency, comprising the following steps:

步骤1，选定多光谱相机中某一谱段作为参考谱段，其余谱段为非参考谱段，针对各非参考谱段与参考谱段之间的相对几何畸变分别进行在轨检校，并将检校结果保存；Step 1. Select a certain spectral segment in the multispectral camera as the reference spectral segment, and the rest of the spectral segments are non-reference spectral segments, and carry out on-orbit calibration for the relative geometric distortion between each non-reference spectral segment and the reference spectral segment. and save the calibration results;

设参考谱段记为B2，任一非参考谱段记为B1，针对非参考谱段B1与参考谱段B2之间的相对几何畸变进行在轨检校实现方式如下，Let the reference spectrum be denoted as B2, and any non-reference spectrum be denoted as B1. The implementation method of on-orbit calibration for the relative geometric distortion between the non-reference spectrum B1 and the reference spectrum B2 is as follows,

步骤1.1，设为物方点P在非参考谱段B1与参考谱段B2影像上有同名像点p₁和p₂，Step 1.1, assuming that the object space point P has image points p ₁ and p ₂ with the same name on the images of the non-reference spectrum segment B1 and the reference spectrum segment B2,

构建参考谱段B2的严格几何成像模型如下，The strict geometric imaging model for constructing the reference spectrum B2 is as follows,

$(\begin{matrix} {x x}_{c c 22} \\ {y the y}_{c c 22} \\ f f \end{matrix}) = = {m m}_{22} {R R}_{BS BS} {R R}_{BJ BJ 22} {R R}_{T T 22} {[\begin{matrix} Xp XP - - {X x}_{S S 22} \\ Yp Yp - - {Y Y}_{S S 22} \\ Zp z - - {Z Z}_{S S 22} \end{matrix}]}_{WGS WGS 8484}$

构建非参考谱段B1的严格几何成像模型如下，The strict geometric imaging model for constructing the non-reference band B1 is as follows,

$(\begin{matrix} {x x}_{c c 11} + + Δ Δ {x x}_{11} \\ {y the y}_{c c 11} + + Δ Δ {y the y}_{11} \\ f f \end{matrix}) = = {m m}_{11} {R R}_{BS BS} {R R}_{BJ BJ 11} {R R}_{T T 11} {[\begin{matrix} Xp XP - - {X x}_{S S 11} \\ Yp Yp - - {Y Y}_{S S 11} \\ Zp z - - {Z Z}_{S S 11} \end{matrix}]}_{WGS WGS 8484}$

上式中，(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₁，采用以探元号为自变量的三次多项式如下，In the above formula, (Xp, Yp, Zp) _WGS84 is the WGS84 geocentric Cartesian coordinates of object space point P; (x _C1 , y _C1 , -f) and (x _C2 , y _C2 , -f) are image points with the same name The coordinates of p ₁ and p ₂ in the camera coordinate system, f represents the main distance of the camera; m ₁ and m ₂ are the photographic scale factors; (X _S1 , Y _S1 , Z _S1 ) and (X _S2 , Y _S2 , Z _S2 ) is the WGS84 earth-centered Cartesian coordinates of the projection centers S1 and S2; R _BS is the installation angle matrix _of the camera in the satellite body coordinate system; R _BJ1 and R _T1 are the satellite body coordinate system and the earth-centered inertial coordinate system, The rotation matrix between the earth-centered inertial coordinate system and the WGS84 earth-centered rectangular coordinate system; R _BJ2 and R _T2 are the corresponding rotation matrices when p ₂ is imaged; for the additional parameter items Δx ₁ and Δy ₁ , the probe number is used as The cubic polynomial of the independent variable is as follows,

$\{\begin{matrix} Δ Δ {x x}_{11} = = {ax ax}_{00}^{11} + + {ax ax}_{11}^{11} \times \times s the s + + {ax ax}_{22}^{11} \times \times {s the s}^{22} + + {ax ax}_{33}^{11} \times \times {s the s}^{33} \\ Δ Δ {y the y}_{11} = = {ay ay}_{00}^{11} + + {ay ay}_{11}^{11} \times \times s the s + + {ay ay}_{22}^{11} \times \times {s the s}^{22} + + {ay ay}_{33}^{11} \times \times {s the s}^{33} \end{matrix}$

其中，为三次多项式的系数，s代表探元号；in, is the coefficient of the cubic polynomial, and s represents the probe number;

将参考谱段B2相对于非参考谱段B1的几何畸变参数记为，则The geometric distortion parameters of the reference spectrum B2 relative to the non-reference spectrum B1 are recorded as ,but

${X x}_{I I}^{11} = = (({ax ax}_{00}^{11},, {ax ax}_{11}^{11,,} {ax ax}_{22}^{11},, {ax ax}_{33}^{11},, {ay ay}_{00}^{11},, {ay ay}_{11}^{11},, {ay ay}_{22}^{11},, {ay ay}_{33}^{11}))$

采用符号f₁表示坐标正投影换算，建立将像点(x,y)正投影至物方获取其对应的物方点的WGS84地心直角坐标(X,Y,Z)_WGS84的公式如下，The symbol f ₁ is used to represent the coordinate orthographic conversion, and the WGS84 geocentric Cartesian coordinates (X, Y, Z) _WGS84 formula for establishing the orthographic projection of the image point (x, y) to the object space to obtain its corresponding object space point is as follows,

采用符号f₂表示坐标反投影换算，建立将物方点的WGS84地心直角坐标(X,Y,Z)_WGS84反投影至像方所得到的像点坐标(x,y)的公式如下，The symbol f ₂ is used to represent the coordinate back-projection conversion, and the formula for establishing the coordinates (x, y) of the image point obtained by back-projecting the _WGS84 geocentric Cartesian coordinates (X, Y, Z) of the object space point to the image space is as follows,

步骤1.2，在非参考谱段B1与参考谱段B2影像上量测n对同名像点同名像点对应物方点Pⁱ，n为同名像点对的总数，i＝1,…,n；对每个同名点对根据参考谱段B2影像上的像点的坐标，利用参考谱段B2的严格几何成像模型及物方高程信息，执行坐标正投影换算f₁，将像点投影至物方，获取相应物方点Pⁱ的WGS84地心直角坐标 ${(X_{P^{i}}, Y_{P^{i}}, Z_{P^{i}})}_{WGS 84};$ Step 1.2, measure n pairs of image points with the same name on the non-reference spectrum B1 and reference spectrum B2 images The same name as the dot Corresponding object space point P ⁱ , n is the total number of image point pairs with the same name, i=1,...,n; for each point pair with the same name According to the image points on the reference spectrum B2 image coordinates of , using the strict geometric imaging model of the reference spectrum segment B2 and the elevation information of the object space, the coordinate orthographic conversion f ₁ is performed, and the image point Project to the object space, and obtain the WGS84 geocentric Cartesian coordinates of the corresponding object point P ⁱ ${(x_{P^{i}}, Y_{P^{i}}, Z_{P^{i}})}_{WGS 84};$

步骤1.3，利用步骤1.2所得物方点Pⁱ的坐标基于非参考谱段B1的严格几何成像模型，利用空间后方交会的原理解算参考谱段B2相对于非参考谱段B1的几何畸变参数消除非参考谱段B1与参考谱段B2之间的相对几何畸变；Step 1.3, use the coordinates of the object space point P ⁱ obtained in step 1.2 Based on the strict geometric imaging model of the non-reference spectrum B1, use the principle of spatial resection to calculate the geometric distortion parameters of the reference spectrum B2 relative to the non-reference spectrum B1 Eliminate the relative geometric distortion between the non-reference spectrum B1 and the reference spectrum B2;

步骤2，基于物方定位一致性，利用各谱段的严格几何成像模型，将非参考谱段与参考谱段进行精确配准；针对非参考谱段B1与参考谱段B2进行精确配准实现方式如下，Step 2, based on the consistency of object space positioning, using the strict geometric imaging model of each spectrum segment, the non-reference spectrum segment and the reference spectrum segment are accurately registered; the non-reference spectrum segment B1 and the reference spectrum segment B2 are accurately registered to achieve The way is as follows,

步骤2.1，根据参考谱段影像像点坐标和物方高程信息获取物方点坐标，包括对参考谱段B2影像上的每个像元p₂(x_c2,y_c2)执行坐标正投影换算f₁，获取其物方点P的WGS84地心直角坐标(X_P,Y_P,Z_P)_WGS84；Step 2.1, obtain the object space point coordinates according to the image point coordinates of the reference spectrum segment image and the object space elevation information, including performing coordinate orthographic conversion f for each pixel p ₂ (x _c2 ,y _c2 ) on the reference spectrum segment B2 image _1. Obtain the WGS84 geocentric Cartesian coordinates (X _P , Y _P , Z _P ) _WGS84 of its object space point P;

步骤2.2，获取非参考谱段影像像点坐标，包括利用步骤1所得参考谱段B2相对于非参考谱段B1的几何畸变参数构建非参考谱段B1的严格几何成像模型；对步骤2.1所得物方点P的坐标(X_P,Y_P,Z_P)_WGS84执行坐标反投影计算f₂，获得非参考谱段B1影像上对应的像点p₁的坐标(x_c1,y_c1)；Step 2.2, obtain the image point coordinates of the non-reference spectrum segment, including using the geometric distortion parameters of the reference spectrum segment B2 obtained in step 1 relative to the non-reference spectrum segment B1 Construct a strict geometric imaging model of the non-reference spectrum segment B1; perform coordinate back-projection calculation f ₂ on the coordinates (X _P , Y _P , Z _P ) _WGS84 of the object space point P obtained in step 2.1, and obtain the correspondence on the image of the non-reference spectrum segment B1 The coordinates of the image point p ₁ (x _c1 , y _c1 );

步骤2.3，根据步骤2.2所得非参考谱段B1影像上对应的像点p₁的坐标(x_c1,y_c1)进行灰度重采样，完成非参考谱段B1与参考谱段B2的精确配准。Step 2.3, according to the coordinates (x _c1 , y _c1 ) of the corresponding image point p ₁ on the non-reference spectrum B1 image obtained in step 2.2, perform grayscale resampling to complete the precise registration of the non-reference spectrum B1 and the reference spectrum B2 .

而且，坐标正投影换算和坐标反投影换算的实现方式如下，Moreover, the implementation of coordinate forward projection conversion and coordinate back-projection conversion is as follows,

设参考谱段B2影像上像点p₂的坐标为(x_c2,y_c2)，对应的物方点P的WGS84地心直角坐标为(Xp,Yp,Zp)_WGS84，令Let the coordinates of image point p ₂ on the reference spectrum B2 image be (x _c2 , y _c2 ), and the WGS84 geocentric Cartesian coordinates of the corresponding object space point P be (Xp, Yp, Zp) _WGS84 , let

${R R}_{BS BS} {R R}_{BJ BJ 22} {R R}_{T T 22} = = [\begin{matrix} {a a}_{11} & {b b}_{11} & {c c}_{11} \\ {a a}_{22} & {b b}_{22} & {c c}_{22} \\ {a a}_{33} & {b b}_{33} & {c c}_{33} \end{matrix}]$

a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃为矩阵的元素；a ₁ , a ₂ , a ₃ , b ₁ , b ₂ , b ₃ , c ₁ , c ₂ , and c ₃ are elements of the matrix;

根据物方点P的WGS84地心直角坐标(Xp,Yp,Zp)_WGS84与其大地坐标(Bp,Lp,Hp)之间的如下关系，According to the following relationship between the _WGS84 geocentric Cartesian coordinates (Xp, Yp, Zp) of the object space point P and its geocentric coordinates (Bp, Lp, Hp),

${[\begin{matrix} Xp XP \\ Yp Yp \\ Zp z \end{matrix}]}_{WGS WGS 8484} = = [\begin{matrix} ((N N + + Hp HP)) \cdot \cdot cos cos Bp bp \cdot \cdot cos cos Lp LP \\ ((N N + + Hp HP)) \cdot \cdot cos cos Bp bp \cdot \cdot sin sin Lp LP \\ ((N N \cdot \cdot ((11 - - {e e}^{22})) + + Hp HP)) \cdot \cdot sin sin Bp bp \end{matrix}]$

其中，e代表地球椭球扁率，变量a代表地球椭球长半轴；Among them, e represents the oblateness of the earth ellipsoid, and the variable a represents the semi-major axis of the earth ellipsoid;

得到像点p₂的坐标(x_c2,y_c2)与大地坐标(Bp,Lp,Hp)关系式如下，The relationship between the coordinates (x _c2 , y _c2 ) of the image point p ₂ and the earth coordinates (Bp, Lp, Hp) is obtained as follows,

$\{\begin{matrix} ((N N + + Hp HP)) \cdot &Center Dot; cos cos Bp bp \cdot &Center Dot; cos cos Lp LP = = \frac{{a a}_{11} {x x}_{c c 22} + + {a a}_{22} {x x}_{c c 22} + + {a a}_{33} f f}{{c c}_{11} {x x}_{c c 22} + + {c c}_{22} {x x}_{c c 22} + + {c c}_{33} f f} \cdot &Center Dot; ((((N N \cdot &Center Dot; ((11 - - {e e}^{22})) + + Hp HP)) \cdot &Center Dot; sin sin Bp bp - - {Z Z}_{S S})) + + {X x}_{S S} \\ ((N N + + Hp HP)) \cdot &Center Dot; cos cos Bp bp \cdot &Center Dot; sin sin Lp LP = = \frac{{b b}_{11} {x x}_{c c 22} + + {b b}_{11} {x x}_{c c 22} + + {b b}_{33} f f}{{c c}_{11} {x x}_{c c 22} + + {c c}_{22} {x x}_{c c 22} + + {c c}_{33} f f} \cdot &Center Dot; ((((N N \cdot &Center Dot; ((11 - - {e e}^{22})) + + Hp HP)) \cdot &Center Dot; sin sin Bp bp - - {Z Z}_{S S})) + + {Y Y}_{S S} \end{matrix}$

坐标正投影换算时，由像点p₂的坐标(x_c2,y_c2)以及给定的物方高程值Hp，利用上式解算物方点P的大地经纬度(Bp,Lp)；再根据物方点P的大地坐标(Bp,Lp,Hp)获取相应WGS84地心直角坐标(Xp,Yp,Zp)_WGS84；In coordinate orthographic conversion, the coordinates (x _c2 , y _c2 ) of the image point p ₂ and the given object elevation value Hp are used to calculate the geodetic latitude and longitude (Bp, Lp) of the object-space point P; then according to The geodetic coordinates (Bp, Lp, Hp) of the object space point P obtain the corresponding WGS84 geocentric Cartesian coordinates (Xp, Yp, Zp) _WGS84 ;

坐标反投影换算时，基于参考谱段B2的严格几何成像模型，由物方点P的WGS84地心直角坐标(Xp,Yp,Zp)_WGS84，利用下式解算，In coordinate back-projection conversion, based on the strict geometric imaging model of the reference spectrum segment B2, the WGS84 geocentric Cartesian coordinates (Xp, Yp, Zp) _WGS84 of the object space point P are calculated using the following formula,

$\{\begin{matrix} {x x}_{c c 22} = = \frac{{a a}_{11} ((Xp XP - - {X x}_{S S 22})) + + {b b}_{11} ((Yp Yp - - {Y Y}_{S S 22})) + + {c c}_{11} ((Zp z - - {Z Z}_{S S 22}))}{{a a}_{33} ((Xp XP - - {X x}_{S S 22})) + + {b b}_{33} ((Yp Yp - - {Y Y}_{S S 22})) + + {c c}_{33} ((Zp z - - {Z Z}_{S S 22}))} \cdot \cdot f f \\ {y the y}_{c c 22} = = \frac{{a a}_{22} ((Xp XP - - {X x}_{S S 22})) + + {b b}_{22} ((Yp Yp - - {Y Y}_{S S 22})) + + {c c}_{22} ((Zp z - - {Z Z}_{S S 22}))}{{a a}_{33} ((Xp XP - - {X x}_{S S 22})) + + {b b}_{33} ((Yp Yp - - {Y Y}_{S S 22})) + + {c c}_{33} ((Zp z - - {Z Z}_{S S 22}))} \cdot &Center Dot; f f \end{matrix}$

得到对应的像点p₂的坐标(x_c2,y_c2)。The coordinates (x _c2 , y _c2 ) of the corresponding image point p ₂ are obtained.

而且，步骤1.3解算参考谱段B2相对于非参考谱段B1的几何畸变参数的实现方式如下，Moreover, step 1.3 calculates the geometric distortion parameters of the reference spectrum B2 relative to the non-reference spectrum B1 is implemented as follows,

令make

$[\begin{matrix} Ux Ux \\ Uy uy \\ Uz Uz \end{matrix}] = = {R R}_{BS BS} {R R}_{BJ BJ 11} {R R}_{T T 11} {[\begin{matrix} Xp XP - - {X x}_{S S 11} \\ Yp Yp - - {Y Y}_{S S 11} \\ Zp z - - {Z Z}_{S S 11} \end{matrix}]}_{WGS WGS 8484}$

上式中，矢量 $[\begin{matrix} Ux \\ Uy \\ Uz \end{matrix}]$ 代表从相机投影中心到物方点的矢量在相机坐标系下的坐标；In the above formula, the vector $[\begin{matrix} Ux \\ uy \\ Uz \end{matrix}]$ Represents the coordinates of the vector from the camera projection center to the object space point in the camera coordinate system;

根据非参考谱段B1的严格几何成像模型，得到下式According to the strict geometric imaging model of the non-reference spectrum B1, the following formula is obtained

$\{\begin{matrix} (({x x}_{c c 11} + + {Δx Δx}_{11})) - - \frac{Ux Ux \cdot &Center Dot; f f}{Uz Uz} = = 00 \\ (({y the y}_{c c 11} + + {Δy Δy}_{11})) - - \frac{Uy uy \cdot \cdot f f}{Uz Uz} = = 00 \end{matrix}$

令make

$\{\begin{matrix} {v v}_{xi xi} = = (({x x}_{c c 11} + + {Δx Δx}_{11})) - - \frac{Ux Ux \cdot &Center Dot; f f}{Uz Uz} \\ {v v}_{yi yi} = = (({y the y}_{c c 11} + + Δ Δ {y the y}_{11})) - - \frac{Uy uy \cdot &Center Dot; f f}{Uz Uz} \end{matrix}$

上式中，v_xi和v_yi分别代表沿轨和垂轨方向的像方残差；In the above formula, v _xi and v _yi represent the image square residuals in the along-track and vertical-track directions, respectively;

将步骤1.2中所得物方点Pⁱ坐标代入上式中，对每个物方点Pⁱ构建如如下误差方程，The coordinates of the object space point P ⁱ obtained in step 1.2 Substituting into the above formula, construct the following error equation for each object space point P ⁱ ,

V_i＝A_iX-L_i，W_i，V _i =A _i XL _i , W _i ,

其中，in,

${V V}_{i i} = = [\begin{matrix} {v v}_{xi xi} \\ {v v}_{yi yi} \end{matrix}]$ ${A A}_{i i} = = [\begin{matrix} 11 & s the s & {s the s}^{22} & {s the s}^{33} & 00 & 00 & 00 & 00 \\ 00 & 00 & 00 & 00 & 11 & s the s & {s the s}^{22} & {s the s}^{33} \end{matrix}]$

${L L}_{i i} = = [\begin{matrix} {((\frac{Ux Ux \cdot \cdot f f}{Uz Uz} - - {x x}_{c c 11}))}_{i i} \\ {((\frac{Uy uy \cdot &Center Dot; f f}{Uz Uz} - - {y the y}_{c c 11}))}_{i i} \end{matrix}]$

$X x = = {(({X x}_{I I}^{11}))}^{T T} = = {(({ax ax}_{00}^{11},, {ax ax}_{11}^{11},, {ax ax}_{22}^{11},, {ax ax}_{33}^{11},, {ay ay}_{00}^{11},, {ay ay}_{11}^{11},, {ay ay}_{22}^{11},, {ay ay}_{33}^{11}))}^{T T}$

上式中，V_i、A_i、L_i分别是利用物方点Pⁱ构建的误差方程，的残差向量、待解参数的系数矩阵以及常向量；X代表参考谱段B2相对于非参考谱段B1的几何畸变参数 ${(X_{I}^{1})}^{T} = {({ax}_{0}^{1}, {ax}_{1}^{1}, {ax}_{2}^{1}, {ax}_{3}^{1}, {ay}_{0}^{1}, {ay}_{1}^{1}, {ay}_{2}^{1}, {ay}_{3}^{1})}^{T};$ W_i是非参考B1谱段影像上的像点的量测精度对应的权；In the above formula, V _i , A _i , and _Li are respectively the residual vector of the error equation constructed using the object space point P ⁱ , the coefficient matrix of the parameters to be solved, and the constant vector; X represents the relative Geometric Distortion Parameters of Spectrum B1 ${(x_{I}^{1})}^{T} = {({ax}_{0}^{1}, {ax}_{1}^{1}, {ax}_{2}^{1}, {ax}_{3}^{1}, {ay}_{0}^{1}, {ay}_{1}^{1}, {ay}_{2}^{1}, {ay}_{3}^{1})}^{T};$ W _i is the image point on the non-reference B1 spectrum image The weight corresponding to the measurement accuracy;

基于最小二乘平差原理，利用下式计算参考谱段B2相对于非参考谱段B1的几何畸变参数，Based on the principle of least squares adjustment, use the following formula to calculate the geometric distortion parameters of the reference spectrum B2 relative to the non-reference spectrum B1,

$X x = = {(({X x}_{I I}^{11}))}^{T T} = = {(({Σ Σ}_{i i = = 11}^{n no} {{A A}_{i i}^{T T} W W}_{i i} {A A}_{i i}))}^{- - 11} (({Σ Σ}_{i i = = 11}^{n no} {{A A}_{i i}^{T T} W W}_{i i} {L L}_{i i}))$

记录计算所得几何畸变参数。Record the calculated geometric distortion parameters.

本发明选定多光谱相机中某一谱段作为参考谱段，其余谱段为非参考谱段。利用一景质量较优的影像，基于各谱段严格几何成像模型，仅利用谱段间的同名像点，对各非参考谱段其与参考谱段之间的相对几何畸变进行在轨检校；对非参考波段其与参考波段间的相对几何畸变进行在轨检校后，基于同名像元物方定位一致性的约束条件，利用各谱段严格几何成像模型，在无需影像匹配的情况下实现子像素级的多光谱影像自动配准，波段间的几何纠正模型为严格几何成像模型，理论上具有严密性。本发明的优点在于：1.影像配准之前无需进行影像匹配，一方面大大提升了数据处理效率，同时，配准质量与影像质量、地物类型无关，对于水域、沙漠以及山地等纹理特征不丰富、匹配质量难以保证的区域，本发明的配准质量也能得到保证。2.一种真正几何意义上的配准，各波段影像之间的几何纠正模型为严格几何成像模型，在理论上具有严密性，能够进一步保证配准精度。The present invention selects a certain spectral segment in the multispectral camera as a reference spectral segment, and the rest of the spectral segments are non-reference spectral segments. Using an image with better quality, based on the strict geometric imaging model of each spectrum, only the same-named image points between the spectrums are used to perform on-orbit calibration of the relative geometric distortion between each non-reference spectrum and the reference spectrum ; After the on-orbit calibration of the relative geometric distortion between the non-reference band and the reference band, based on the constraints of the same-name pixel object space positioning consistency, using the strict geometric imaging model of each band, without image matching Realize the automatic registration of multi-spectral images at the sub-pixel level, and the geometric correction model between the bands is a strict geometric imaging model, which is rigorous in theory. The advantages of the present invention are as follows: 1. Image matching is not required before image registration, on the one hand, the data processing efficiency is greatly improved, and at the same time, registration quality is independent of image quality and feature type, and is not suitable for texture features such as waters, deserts, and mountains. The registration quality of the present invention can also be guaranteed for areas that are abundant and whose matching quality is difficult to guarantee. 2. A true geometric registration. The geometric correction model between images in each band is a strict geometric imaging model, which is rigorous in theory and can further guarantee the registration accuracy.

附图说明Description of drawings

图1为本发明实施例的配准流程示意图；FIG. 1 is a schematic diagram of a registration process in an embodiment of the present invention;

图2为本发明基于物方定位一致性的卫星多光谱影像配准原理示意图。Fig. 2 is a schematic diagram of the principle of satellite multi-spectral image registration based on object space positioning consistency in the present invention.

具体实施方式Detailed ways

以下结合附图和实施例详细说明本发明具体实施方式。实施例的流程可以分为两个步骤，每个步骤实施的具体方法、公式以及流程如下：The specific implementation of the present invention will be described in detail below in conjunction with the accompanying drawings and examples. The process of the embodiment can be divided into two steps, and the specific method, formula and process implemented in each step are as follows:

步骤1，选定多光谱相机中某一谱段作为参考谱段，其余谱段为非参考谱段，利用一景质量较优的影像，针对各非参考谱段，对其与参考谱段之间的相对几何畸变进行在轨检校，并将检校结果保存，用于后续的影像配准。设参考谱段记为B2，任一非参考谱段记为B1，针对非参考谱段B1，其与参考谱段B2之间的相对几何畸变在轨检校的具体步骤及公式如下：Step 1. Select a certain spectral segment in the multi-spectral camera as the reference spectral segment, and the rest of the spectral segments are non-reference spectral segments. Using an image with better quality, for each non-reference spectral segment, the difference between the non-reference spectral segment and the reference spectral segment The relative geometric distortion between them is checked on the track, and the check results are saved for subsequent image registration. Assuming that the reference spectrum is denoted as B2, and any non-reference spectrum is denoted as B1, for the non-reference spectrum B1, the specific steps and formulas for on-orbit calibration of the relative geometric distortion between it and the reference spectrum B2 are as follows:

步骤1.1，利用姿态、轨道、时间等辅助数据以及相机参数，构建参考谱段B2的严格几何成像模型；如式（1）；Step 1.1, using auxiliary data such as attitude, orbit, time, and camera parameters to construct a strict geometric imaging model of the reference spectrum B2; as in formula (1);

$(\begin{matrix} {x x}_{c c 22} \\ {y the y}_{c c 22} \\ f f \end{matrix}) = = {m m}_{22} {R R}_{BS BS} {R R}_{BJ BJ 22} {R R}_{T T 22} {[\begin{matrix} Xp XP - - {X x}_{S S 22} \\ Yp Yp - - {Y Y}_{S S 22} \\ Zp z - - {Z Z}_{S S 22} \end{matrix}]}_{WGS WGS 8484} - - - - - - ((11))$

同理，构建B1谱段的严格几何成像模型，并在其内定向参数模型中引入附加参数项Δx₁和Δy₁，用于补偿B1谱段相对于B2谱段的几何畸变，实现构建基于扩展共线条件方程的自检校平差模型，如式（2）。In the same way, a strict geometric imaging model of the B1 spectrum is constructed, and additional parameter items Δx ₁ and Δy ₁ are introduced into its internal orientation parameter model to compensate for the geometric distortion of the B1 spectrum relative to the B2 spectrum, so as to realize the construction based on the extended The self-checking adjustment difference model of the collinear conditional equation, such as formula (2).

$(\begin{matrix} {x x}_{c c 11} + + Δ Δ {x x}_{11} \\ {y the y}_{c c 11} + + Δ Δ {y the y}_{11} \\ f f \end{matrix}) = = {m m}_{11} {R R}_{BS BS} {R R}_{BJ BJ 11} {R R}_{T T 11} {[\begin{matrix} Xp XP - - {X x}_{S S 11} \\ Yp Yp - - {Y Y}_{S S 11} \\ Zp z - - {Z Z}_{S S 11} \end{matrix}]}_{WGS WGS 8484} - - - - - - ((22))$

上式中，(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）所示：In the above formula, (Xp, Yp, Zp) _WGS84 is the WGS84 geocentric Cartesian coordinates of the object space point P; (x _C1 , y _C1 , -f) and (x _C2 , y _C2 , -f) are the object space points The coordinates of p 1 and p ₂ in the camera coordinate system of the image point p ₁ and p 2 with the same name on the B1 spectrum and B2 spectrum images, f represents the main distance of the camera; m ₁ and m ₂ are the photographic scale factors; (X _S1 , Y _S1 , Z _S1 ) and (X _S2 , Y _S2 , Z _S2 ) are the WGS84 geocentric Cartesian coordinates of projection centers S1 and S2; R _BS is the installation angle matrix of the camera in the satellite body coordinate system; R _BJ1 and R _T1 are p ₁ The rotation matrix between the satellite body coordinate system and the earth-centered inertial coordinate system, the earth-centered inertial coordinate system and the WGS84 earth-centered Cartesian coordinate system during imaging, the earth-centered inertial coordinate system is the J2000 coordinate system, etc.; R _BJ2 and R _T2 are The corresponding rotation matrix when p ₂ is imaged. For the additional parameter items Δx ₁ and Δy ₁ in the B1 spectrum imaging model, a cubic polynomial with the probe number as an independent variable is used, as shown in formula (3):

$\{\begin{matrix} Δ Δ {x x}_{11} = = {ax ax}_{00}^{11} + + {ax ax}_{11}^{11} \times \times s the s + + {ax ax}_{22}^{11} \times \times {s the s}^{22} + + {ax ax}_{33}^{11} \times \times {s the s}^{33} \\ Δ Δ {y the y}_{11} = = {ay ay}_{00}^{11} + + {ay ay}_{11}^{11} \times \times s the s + + {ay ay}_{22}^{11} \times \times {s the s}^{22} + + {ay ay}_{33}^{11} \times \times {s the s}^{33} \end{matrix} - - - - - - ((33))$

其中，为三次多项式的系数，s代表探元号。in, is the coefficient of the cubic polynomial, and s represents the probe number.

将B1谱段相对于B2谱段的几何畸变参数记为（上标1代表谱段B1），则：The geometric distortion parameter of the B1 spectrum relative to the B2 spectrum is recorded as (The superscript 1 represents the spectrum B1), then:

${X x}_{I I}^{11} = = (({ax ax}_{00}^{11},, {ax ax}_{11}^{11},, {ax ax}_{22}^{11},, {ax ax}_{33}^{11},, {ay ay}_{00}^{11},, {ay ay}_{11}^{11},, {ay ay}_{22}^{11},, {ay ay}_{33}^{11}))$

利用卫星获取的姿态、轨道、时间、相机参数以及物方高程信息，基于严格几何成像模型，能够实现像点坐标与物方点坐标之间的换算。以谱段B2为例，对具体的换算过程及公式进行阐述。Using the attitude, orbit, time, camera parameters and object space elevation information obtained by the satellite, based on the strict geometric imaging model, the conversion between the image point coordinates and the object space point coordinates can be realized. Taking spectrum B2 as an example, the specific conversion process and formula are described.

假设谱段B2影像上某像点p₂的坐标为(x_c2,y_c2)，其对应的物方点P的WGS84地心直角坐标为(Xp,Yp,Zp)_WGS84，令公式（1）中：Assuming that the coordinates of a certain image point p ₂ on the image of spectrum segment B2 are (x _c2 , y _c2 ), and the WGS84 geocentric Cartesian coordinates of the corresponding object space point P are (Xp,Yp,Zp) _WGS84 , let the formula (1) middle:

a₁、a₂、a₃、b₁、b₂、b₃、c₁、c₂、c₃为矩阵的元素。a ₁ , a ₂ , a ₃ , b ₁ , b ₂ , b ₃ , c ₁ , c ₂ , and c ₃ are elements of the matrix.

则利用公式（4）及物方点P的WGS84地心直角坐标(Xp,Yp,Zp)_WGS84可解算像点p₂的坐标(x_c2,y_c2)：Then, the coordinates (x _c2 , y _c2 ) of the image point p ₂ can be calculated by using the formula (4) and the _WGS84 geocentric Cartesian coordinates (Xp, Yp, Zp) of the object space point P:

$\{\begin{matrix} {x x}_{c c 22} = = \frac{{a a}_{11} ((Xp XP - - {X x}_{S S 22})) + + {b b}_{11} ((Yp Yp - - {Y Y}_{S S 22})) + + {c c}_{11} ((Zp z - - {Z Z}_{S S 22}))}{{a a}_{33} ((Xp XP - - {X x}_{S S 22})) + + {b b}_{33} ((Yp Yp - - {Y Y}_{S S 22})) + + {c c}_{33} ((Zp z - - {Z Z}_{S S 22}))} \cdot \cdot f f \\ {y the y}_{c c 22} = = \frac{{a a}_{22} ((Xp XP - - {X x}_{S S 22})) + + {b b}_{22} ((Yp Yp - - {Y Y}_{S S 22})) + + {c c}_{22} ((Zp z - - {Z Z}_{S S 22}))}{{a a}_{33} ((Xp XP - - {X x}_{S S 22})) + + {b b}_{33} ((Yp Yp - - {Y Y}_{S S 22})) + + {c c}_{33} ((Zp z - - {Z Z}_{S S 22}))} \cdot &Center Dot; f f \end{matrix} - - - - - - ((44))$

根据大地测量学的基本原理可知，物方点P的WGS84地心直角坐标(Xp,Yp,Zp)_WGS84与其大地坐标(Bp,Lp,Hp)（Bp、Lp分别为物方点P的纬度和经度，Hp为其椭球高）之间有如下关系：According to the basic principles of geodesy, it can be known that the _WGS84 geocentric Cartesian coordinates (Xp, Yp, Zp) of object space point P and its geodetic coordinates (Bp, Lp, Hp) (Bp, Lp are the latitude and Longitude, Hp is its ellipsoid height) has the following relationship:

${[\begin{matrix} Xp XP \\ Yp Yp \\ Zp z \end{matrix}]}_{WGS WGS 8484} = = [\begin{matrix} ((N N + + Hp HP)) \cdot &Center Dot; cos cos Bp bp \cdot &Center Dot; cos cos Lp LP \\ ((N N + + Hp HP)) \cdot \cdot cos cos Bp bp \cdot \cdot sin sin Lp LP \\ ((N N \cdot &Center Dot; ((11 - - {e e}^{22})) + + Hp HP)) \cdot &Center Dot; sin sin Bp bp \end{matrix}] - - - - - - ((55))$

其中，e代表地球椭球扁率，变量a代表地球椭球长半轴。将上式（5）代入公式（4）中有公式（6）：Among them, e represents the oblateness of the earth ellipsoid, and the variable a represents the semi-major axis of the Earth's ellipsoid. Substituting the above formula (5) into formula (4) gives formula (6):

$\{\begin{matrix} ((N N + + Hp HP)) \cdot &Center Dot; cos cos Bp bp \cdot &Center Dot; cos cos Lp LP = = \frac{{a a}_{11} {x x}_{c c 22} + + {a a}_{22} {x x}_{c c 22} + + {a a}_{33} f f}{{c c}_{11} {x x}_{c c 22} + + {c c}_{22} {x x}_{c c 22} + + {c c}_{33} f f} \cdot &Center Dot; ((((N N \cdot &Center Dot; ((11 - - {e e}^{22})) + + Hp HP)) \cdot &Center Dot; sin sin Bp bp - - {Z Z}_{S S})) + + {X x}_{S S} \\ ((N N + + Hp HP)) \cdot &Center Dot; cos cos Bp bp \cdot &Center Dot; sin sin Lp LP = = \frac{{b b}_{11} {x x}_{c c 22} + + {b b}_{11} {x x}_{c c 22} + + {b b}_{33} f f}{{c c}_{11} {x x}_{c c 22} + + {c c}_{22} {x x}_{c c 22} + + {c c}_{33} f f} \cdot &Center Dot; ((((N N \cdot &Center Dot; ((11 - - {e e}^{22})) + + Hp HP)) \cdot &Center Dot; sin sin Bp bp - - {Z Z}_{S S})) + + {Y Y}_{S S} \end{matrix} - - - - - - ((66))$

坐标正投影换算时，由像点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) _In coordinate orthographic conversion, the geodetic latitude and longitude ( _Bp , _Lp ); substituting the geodetic coordinates (Bp, Lp, Hp) of the object space point P into formula (5) to obtain its WGS84 geocentric Cartesian coordinates (Xp, Yp, Zp) _WGS84 . In coordinate back-projection conversion, based on the strict geometric imaging model of the reference spectrum segment B2, from the WGS84 geocentric Cartesian coordinates (Xp, Yp, Zp) _WGS84 of the object space point P, use formula (4) to solve the corresponding image point p ₂ The coordinates of (x _c2 , y _c2 )

为了便于描述，本文利用公式（7）和（8）简要表达上述的像点坐标与物方点坐标正反换算。其中，公式（7）中的符号f₁表示将像点(x,y)正投影至物方获取其对应的物方点的WGS84地心直角坐标(X,Y,Z)_WGS84，即坐标正投影换算；公式（8）中的符号f₂则表示将物方点的WGS84地心直角坐标(X,Y,Z)_WGS84反投影至像方所得到的像点坐标(x,y)，即坐标反投影换算。For the convenience of description, this paper uses formulas (7) and (8) to briefly express the positive and negative conversion between the above-mentioned image point coordinates and object space point coordinates. Among them, the symbol f ₁ in the formula (7) means that the WGS84 geocentric Cartesian coordinates (X, Y, Z) _WGS84 of the corresponding object space point obtained by orthographic projection of the image point (x, y) onto the object space, that is, the coordinates are positive Projection conversion; the symbol f ₂ in the formula (8) represents the image point coordinates (x, y) obtained by back-projecting the _WGS84 geocentric Cartesian coordinates (X, Y, Z) of the object space point to the image space, namely Coordinate back projection conversion.

步骤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；Step 1.2, measure n pairs of image points with the same name on the two spectral images of B1 and B2 Indicates a pair of image points with the same name on the B1 and B2 spectrum images), and the image points with the same name Corresponding to object space point P ⁱ . n is the total number of image point pairs with the same name, which can be set by those skilled in the art according to specific situations. The specific measurement implementation is an existing technology. In order to improve the accuracy and reliability of the calibration results, it is recommended to The distribution on the non-reference spectrum image is required as follows (corresponding image points with the same name The distribution shape and It is consistent in the non-reference spectrum): (1) try to distribute in a short area in the image line direction (ie, push-broom direction); (2) uniform in the image column direction (ie, along the CCD direction) Cover the entire linear array CCD. Among them, the first requirement is to reduce the influence of the outer azimuth element error on the calculation result; the second requirement is to ensure that the calculation result is applicable to all detectors of the linear array CCD. For each dot pair with the same name (i=1,...,n), the image points on the B2 spectral segment image Substituting the coordinates of the coordinates into (x, y) in formula (7), using the strict geometric imaging model of the B2 spectrum shown in formula (1) and the elevation information of the object space, the coordinate orthographic conversion f ₁ is performed, and the image point Project to the object space, and obtain the WGS84 geocentric Cartesian coordinates of its object point P ⁱ That is (X,Y,Z) _WGS84 in formula (7);

步骤1.3，利用前述步骤1.2得到的物方点Pⁱ的坐标(i＝1,…,n)，基于B1谱段的严格几何成像模型（式2），利用空间后方交会的原理解算，消除B1谱段与B2谱段影像之间的相对几何畸变。具体解算公式如下。Step 1.3, using the coordinates of the object space point P ⁱ obtained in the previous step 1.2 (i=1,...,n), based on the strict geometric imaging model of the B1 spectrum (Equation 2), using the principle of spatial resection to calculate , to eliminate the relative geometric distortion between the images of the B1 spectrum and the B2 spectrum. The specific solution formula is as follows.

1）令1) order

$[\begin{matrix} Ux Ux \\ Uy uy \\ Uz Uz \end{matrix}] = = {R R}_{BS BS} {R R}_{BJ BJ 11} {R R}_{T T 11} {[\begin{matrix} Xp XP - - {X x}_{S S 11} \\ Yp Yp - - {Y Y}_{S S 11} \\ Zp z - - {Z Z}_{S S 11} \end{matrix}]}_{WGS WGS 8484} - - - - - - ((99))$

上式中，矢量 $[\begin{matrix} Ux \\ Uy \\ Uz \end{matrix}]$ 称为物方矢量U，代表从相机投影中心到物方点的矢量在相机坐标系下的坐标；In the above formula, the vector $[\begin{matrix} Ux \\ uy \\ Uz \end{matrix}]$ It is called the object space vector U, which represents the coordinates of the vector from the camera projection center to the object space point in the camera coordinate system;

式（2）可转化为式（10）：Formula (2) can be transformed into formula (10):

$\{\begin{matrix} (({x x}_{c c 11} + + {Δx Δx}_{11})) - - \frac{Ux Ux \cdot &Center Dot; f f}{Uz Uz} = = 00 \\ (({y the y}_{c c 11} + + Δ Δ {y the y}_{11})) - - \frac{Uy uy \cdot \cdot f f}{Uz Uz} = = 00 \end{matrix} - - - - - - ((1010))$

令make

$\{\begin{matrix} {v v}_{xi xi} = = (({x x}_{c c 11} + + {Δx Δx}_{11})) - - \frac{Ux Ux \cdot \cdot f f}{Uz Uz} \\ {v v}_{yi yi} = = (({y the y}_{c c 11} + + Δ Δ {y the y}_{11})) - - \frac{Uy uy \cdot \cdot f f}{Uz Uz} \end{matrix} - - - - - - ((1111))$

上式中，v_xi和v_yi分别代表沿轨和垂轨方向的像方残差。In the above formula, v _xi and v _yi represent the image square residuals in the along-track and vertical-track directions, respectively.

2）将步骤1.2中解算的物方点Pⁱ坐标(i＝1,…,n)代入式（11）中，对每个物方点Pⁱ均可构建如式（12）所示的误差方程（下标i代表利用物方点Pⁱ建立的误差方程）：2) The coordinates of the object space point P ⁱ calculated in step 1.2 (i=1,...,n) are substituted into formula (11), and the error equation shown in formula (12) can be constructed for each object space point P ⁱ (the subscript ⁱ represents the error equation):

V_i＝A_iX-L_i，W_i(i＝1,…,n) （12）V _i ＝A _i XL _i ，W _i (i=1,...,n) (12)

其中，in,

${L L}_{i i} = = [\begin{matrix} {((\frac{Ux Ux \cdot &Center Dot; f f}{Uz Uz} - - {x x}_{c c 11}))}_{i i} \\ {((\frac{Uy uy \cdot &Center Dot; f f}{Uz Uz} - - {y the y}_{c c 11}))}_{i i} \end{matrix}]$

上式中，V_i、A_i、L_i分别是利用物方点Pⁱ构建的误差方程式的残差向量、待解参数的系数矩阵以及常向量；X代表B1谱段相对于B2谱段的几何畸变参数 ${(X_{I}^{1})}^{T} = {({ax}_{0}^{1}, {ax}_{1}^{1}, {ax}_{2}^{1}, {ax}_{3}^{1}, {ay}_{0}^{1}, {ay}_{1}^{1}, {ay}_{2}^{1}, {ay}_{3}^{1})}^{T};$ W_i是B1谱段影像上的像点的量测精度对应的权。In the above formula, V _i , A _i , and _Li are respectively the residual vector of the error equation constructed using the object space point P ⁱ , the coefficient matrix of the parameters to be solved, and the constant vector; X represents the difference between the B1 spectrum and the B2 spectrum Geometric Distortion Parameters ${(x_{I}^{1})}^{T} = {({ax}_{0}^{1}, {ax}_{1}^{1}, {ax}_{2}^{1}, {ax}_{3}^{1}, {ay}_{0}^{1}, {ay}_{1}^{1}, {ay}_{2}^{1}, {ay}_{3}^{1})}^{T};$ W _i is the image point on the B1 spectrum image The weight corresponding to the measurement accuracy.

3）基于最小二乘平差原理，利用式（13）计算B1谱段相对于B2谱段的几何畸变参数，实现对B1谱段相对于B2谱段之间的相对几何畸变进行在轨检校，并将几何畸变参数用文件记录下来，用于后续多光谱影像波段配准。3) Based on the principle of least squares adjustment, use formula (13) to calculate the geometric distortion parameters of the B1 spectrum relative to the B2 spectrum, and realize the on-orbit calibration of the relative geometric distortion between the B1 spectrum and the B2 spectrum , and record the geometric distortion parameters in a file for subsequent multispectral image band registration.

$X x = = {(({X x}_{I I}^{11}))}^{T T} = = {(({Σ Σ}_{i i = = 11}^{n no} {{A A}_{i i}^{T T} W W}_{i i} {A A}_{i i}))}^{- - 11} (({Σ Σ}_{i i = = 11}^{n no} {A A}_{i i}^{T T} {W W}_{i i} {L L}_{i i})) - - - - - - ((1313))$

步骤2，通过步骤1预先获取了各非参考谱段其相对于参考谱段的几何畸变参数后，基于物方定位一致性，利用各谱段影像的严格几何成像模型，对任意一景多光谱影像，将其各非参考谱段与参考谱段分别进行精确配准。仍设参考谱段记为B2，任一非参考谱段记为B1，结合图1和图2对本步骤中的具体原理和流程进行阐述如下。其中，图1为流程示意图，图2为原理示意图。图2中，p₁和p₂分别为B1和B2谱段影像上的一对同名像点，P则代表其对应的物方点；S1和S2分别为像点p₁和p₂成像时的相机投影中心；f₁和f₂则代表公式7和8中的像点坐标与物方点坐标之间的正反换算。Step 2. After pre-obtaining the geometric distortion parameters of each non-reference spectrum segment relative to the reference spectrum segment through step 1, based on the consistency of object space positioning, using the strict geometric imaging model of each spectrum segment image, the multispectral Image, each non-reference spectrum segment and the reference spectrum segment are precisely registered. It is still assumed that the reference spectrum segment is marked as B2, and any non-reference spectrum segment is marked as B1. The specific principles and processes in this step are described below in conjunction with Figure 1 and Figure 2 . Wherein, FIG. 1 is a schematic flow chart, and FIG. 2 is a schematic schematic diagram of the principle. In Fig. 2, p ₁ and p ₂ are a pair of image points with the same name on the B1 and B2 spectrum images respectively, and P represents the corresponding object space point; S1 and S2 are the imaging points of image points p ₁ and p ₂ respectively The camera projection center; f ₁ and f ₂ represent the positive and negative conversion between the image point coordinates and the object space point coordinates in formulas 7 and 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）；Step 2.1, obtain object space point coordinates according to the image point coordinates and object space elevation information of the reference spectral segment image: For each pixel p ₂ (x _c2 ,y _c2 ) on the B2 spectral segment image, its image point coordinates (x _c2 , y _c2 ) into (x, y) in formula (7), that is, using the strict geometric imaging model (Formula 1) of the B2 spectrum and the elevation information of the object space, the coordinate orthographic projection conversion f ₁ is performed, and the image point p ₂ (x _c2 ,y _c2 ) is projected to the object space, and the WGS84 geocentric Cartesian coordinates (X _P ,Y _P ,Z _P ) of the object space point P are obtained _WGS84 (that is, (X,Y,Z in formula (7) ) _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所示。Step 2.2, obtain image point coordinates of the non-reference spectral segment: use the geometric distortion parameters of the B1 spectral segment obtained in step 1 relative to the B2 spectral segment, and construct the strict geometric imaging model shown in formula (2). Substitute the obtained coordinates (X _P , Y _P , Z _P ) _WGS84 of the object space point P into (X, Y, Z) _WGS84 in formula (8), perform coordinate back-projection calculation f ₂ , and obtain the B1 spectral segment image The coordinates (x _c1 , y _c1 ) of the corresponding image point p ₁ are (x, y) in formula (8); since formula (2) has compensated the relative geometric distortion of the two spectral bands B1 and B2, Based on the constrained relationship of the consistency of object space positioning of the image point with the same name, the obtained image point p ₁ (s ₁ ,l ₁ ) is the corresponding image point p ₂ (s,l) on the B2 spectrum image to the B1 spectrum image The same-named image points, so the mapping relationship between the same-named image points on the B1 and B2 band images can be established to realize the registration of multi-spectral images. This is the principle of multi-spectral image registration based on the consistency of object space positioning, such as Figure 2 shows.

步骤2.3，灰度重采样：在步骤2.1获取B2谱段影像的每个像元p₂(s,l)在B1谱段影像上对应的同名像点坐标p₁(s₁,l₁)后，通过灰度重采样，即可实现B1与B2两谱段影像的配准。其中，灰度重采样可采用现有技术中的双线性内插算法实现，本发明不予赘述。Step 2.3, grayscale resampling: After obtaining the pixel coordinates p ₁ (s ₁ ,l ₁ ) corresponding to each pixel p ₂ (s,l) of the B2 spectrum image on the B1 spectrum image in step 2.1 , through grayscale resampling, the registration of the B1 and B2 spectral band images can be realized. Wherein, the gray level resampling can be realized by using the bilinear interpolation algorithm in the prior art, which will not be described in detail in the present invention.

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。The specific embodiments described herein are merely illustrative of the spirit of the invention. Those skilled in the art to which the present invention belongs can make various modifications or supplements to the described specific embodiments or adopt similar methods to replace them, but they will not deviate from the spirit of the present invention or go beyond the definition of the appended claims range.

Claims

1. A method for registration of satellite multispectral images based on object positioning consistency, comprising the following steps:

Step 1. Select a certain spectral segment in the multispectral camera as the reference spectral segment, and the rest of the spectral segments are non-reference spectral segments, and carry out on-orbit calibration for the relative geometric distortion between each non-reference spectral segment and the reference spectral segment. and save the calibration results;

Let the reference spectrum be denoted as B2, and any non-reference spectrum be denoted as B1, and the implementation method of on-orbit calibration for the relative geometric distortion between the non-reference spectrum B1 and the reference spectrum B2 is as follows,

Step 1.1, assuming that the object space point P has image points p ₁ and p ₂ with the same name on the images of the non-reference spectrum segment B1 and the reference spectrum segment B2,

The strict geometric imaging model for constructing the reference spectrum B2 is as follows,

(\begin{matrix} {x x}_{c c 22} \\ {y the y}_{c c 22} \\ f f \end{matrix}) = = {m m}_{22} {R R}_{B B S S} {R R}_{B B J J 22} {R R}_{T T 22} {[\begin{matrix} X x p p - - {X x}_{S S 22} \\ Y Y p p - - {Y Y}_{S S 22} \\ Z Z p p - - {Z Z}_{S S 22} \end{matrix}]}_{W W G G S S 8484}

The strict geometric imaging model for constructing the non-reference band B1 is as follows,

(\begin{matrix} {x x}_{c c 11} + + Δ Δ {x x}_{11} \\ {y the y}_{c c 11} + + {Δy Δy}_{11} \\ f f \end{matrix}) = = {m m}_{11} {R R}_{B B S S} {R R}_{B B J J 11} {R R}_{T T 11} {[\begin{matrix} X x p p - - {X x}_{S S 11} \\ Y Y p p - - {Y Y}_{S S 11} \\ Z Z p p - - {Z Z}_{S S 11} \end{matrix}]}_{W W G G S S 8484}

In the above formula, (Xp, Yp, Zp) _WGS84 is the WGS84 geocentric Cartesian coordinates of object space point P; (x _C1 , y _C1 , -f) and (x _C2 , y _C2 , -f) are image points with the same name The coordinates of p ₁ and p ₂ in the camera coordinate system, f represents the main distance of the camera; m ₁ and m ₂ are the photographic scale factors; (X _S1 , Y _S1 , Z _S1 ) and (X _S2 , Y _S2 , Z _S2 ) is the WGS84 earth-centered Cartesian coordinates of the projection centers S1 and S2; R _BS is the installation angle matrix _of the camera in the satellite body coordinate system; R _BJ1 and R _T1 are the satellite body coordinate system and the earth-centered inertial coordinate system, The rotation matrix between the earth-centered inertial coordinate system and the WGS84 earth-centered rectangular coordinate system; R _BJ2 and R _T2 are the corresponding rotation matrices when p ₂ is imaged; for the additional parameter items Δx ₁ and Δy ₁ , the probe number is used as The cubic polynomial of the independent variable is as follows,

\{\begin{matrix} {Δx Δx}_{11} = = {ax ax}_{00}^{11} + + {ax ax}_{11}^{11} \times \times s the s + + {ax ax}_{22}^{11} \times \times {s the s}^{22} + + {ax ax}_{33}^{11} \times \times {s the s}^{33} \\ {Δy Δy}_{11} = = {ay ay}_{00}^{11} + + {ay ay}_{11}^{11} \times \times s the s + + {ay ay}_{22}^{11} \times \times {s the s}^{22} + + {ay ay}_{33}^{11} \times \times {s the s}^{33} \end{matrix}

in, is the coefficient of the cubic polynomial, and s represents the probe number;

The geometric distortion parameters of the reference spectrum B2 relative to the non-reference spectrum B1 are recorded as but

{X x}_{I I}^{11} = = (({ax ax}_{00}^{11},, {ax ax}_{11}^{11},, {ax ax}_{22}^{11},, {ax ax}_{33}^{11},, {ay ay}_{00}^{11},, {ay ay}_{11}^{11},, {ay ay}_{22}^{11},, {ay ay}_{33}^{11}))

The symbol f ₁ is used to represent the coordinate orthographic conversion, and the WGS84 geocentric Cartesian coordinates (X, Y, Z) _WGS84 formula for establishing the orthographic projection of the image point (x, y) to the object space to obtain its corresponding object space point is as follows,

((x x,, y the y)) \overset{{f f}_{11}}{&RightArrow; &Right Arrow;} {((X x,, Y Y,, Z Z))}_{W W G G S S 8484}

The symbol f ₂ is used to represent the coordinate back-projection conversion, and the formula for establishing the coordinates (x, y) of the image point obtained by back-projecting the _WGS84 geocentric Cartesian coordinates (X, Y, Z) of the object space point to the image space is as follows,

{((X x,, Y Y,, Z Z))}_{W W G G S S 8484} \overset{{f f}_{22}}{&RightArrow; &Right Arrow;} ((x x,, y the y))

Step 1.2, measure n pairs of image points with the same name on the non-reference spectrum B1 and reference spectrum B2 images The same name as the dot Corresponding object space point P ⁱ , n is the total number of image point pairs with the same name, i=1,...,n; for each point pair with the same name According to the image points on the reference spectrum B2 image coordinates of , using the strict geometric imaging model of the reference spectrum segment B2 and the elevation information of the object space, the coordinate orthographic conversion f ₁ is performed, and the image point Project to the object space, and obtain the WGS84 geocentric Cartesian coordinates of the corresponding object point P ⁱ

Step 1.3, use the coordinates of the object space point P ⁱ obtained in step 1.2 Based on the strict geometric imaging model of the non-reference spectrum B1, use the principle of spatial resection to calculate the geometric distortion parameters of the reference spectrum B2 relative to the non-reference spectrum B1 Eliminate the relative geometric distortion between the non-reference spectrum B1 and the reference spectrum B2; step 1.3 solve the geometric distortion parameters of the reference spectrum B2 relative to the non-reference spectrum B1 is implemented as follows,

make

[\begin{matrix} U u x x \\ U u y the y \\ U u z z \end{matrix}] = = {R R}_{B B S S} {R R}_{B B J J 11} {R R}_{T T 11} {[\begin{matrix} X x p p - - {X x}_{S S 11} \\ Y Y p p - - {Y Y}_{S S 11} \\ Z Z p p - - {Z Z}_{S S 11} \end{matrix}]}_{W W G G S S 8484}

In the above formula, the vector

[\begin{matrix} u x \\ u the y \\ u z \end{matrix}]

Represents the coordinates of the vector from the camera projection center to the object space point in the camera coordinate system;

According to the strict geometric imaging model of the non-reference spectrum B1, the following formula is obtained

\{\begin{matrix} (({x x}_{c c 11} + + Δ Δ {x x}_{11})) - - \frac{U u x x \cdot &Center Dot; f f}{U u z z} = = 00 \\ (({y the y}_{c c 11} + + Δ Δ {y the y}_{11})) - - \frac{U u y the y \cdot &Center Dot; f f}{U u z z} = = 00 \end{matrix}

make

\{\begin{matrix} {v v}_{x x i i} = = (({x x}_{c c 11} + + Δ Δ {x x}_{11})) - - \frac{U u x x \cdot &Center Dot; f f}{U u z z} \\ {v v}_{y the y i i} = = (({y the y}_{c c 11} + + Δ Δ {y the y}_{11})) - - \frac{U u y the y \cdot &Center Dot; f f}{U u z z} \end{matrix}

In the above formula, v _xi and v _yi represent the image square residuals in the along-track and vertical-track directions, respectively;

The coordinates of the object space point P ⁱ obtained in step 1.2 Substituting into the above formula, construct the following error equation for each object space point P ⁱ ,

V _i =A _i XL _i , W _i ,

in,

\begin{matrix} {V V}_{i i} = = [\begin{matrix} {v v}_{x x i i} \\ {v v}_{y the y i i} \end{matrix}] & {A A}_{i i} = = [\begin{matrix} 11 & s the s & {s the s}^{22} & {s the s}^{33} & 00 & 00 & 00 & 00 \\ 00 & 00 & 00 & 00 & 11 & s the s & {s the s}^{22} & {s the s}^{33} \end{matrix}] \end{matrix}

{L L}_{i i} [\begin{matrix} {((\frac{U u x x \cdot &Center Dot; f f}{U u z z} - - {x x}_{c c 11}))}_{i i} \\ {((\frac{U u y the y \cdot &Center Dot; f f}{U u z z} - - {y the y}_{c c 11}))}_{i i} \end{matrix}]

X x = = {(({X x}_{I I}^{11}))}^{T T} = = {(({ax ax}_{00}^{11},, {ax ax}_{11}^{11},, {ax ax}_{22}^{11},, {ax ax}_{33}^{11},, {ay ay}_{00}^{11},, {ay ay}_{11}^{11},, {ay ay}_{22}^{11},, {ay ay}_{33}^{11}))}^{T T}

In the above formula, V _i , A _i , and _Li are respectively the residual vector of the error equation constructed using the object space point P ⁱ , the coefficient matrix of the parameters to be solved, and the constant vector; X represents the reference spectrum segment B2 relative to the non-reference spectrum Geometric distortion parameters of segment B1 W _i is the image point on the non-reference B1 spectrum image The weight corresponding to the measurement accuracy;

Based on the principle of least squares adjustment, use the following formula to calculate the geometric distortion parameters of the reference spectrum B2 relative to the non-reference spectrum B1,

X x = = {(({X x}_{I I}^{11}))}^{T T} = = (({Σ Σ}_{i i = = 11}^{n no} {A A}_{i i}^{T T} {W W}_{i i} {A A}_{i i} {))}^{- - 11} (({Σ Σ}_{i i = = 11}^{n no} {A A}_{i i}^{T T} {W W}_{i i} {L L}_{i i}))

Record the calculated geometric distortion parameters;

Step 2, based on the consistency of object space positioning, using the strict geometric imaging model of each spectrum segment, the non-reference spectrum segment and the reference spectrum segment are accurately registered; the non-reference spectrum segment B1 and the reference spectrum segment B2 are accurately registered to achieve The way is as follows,

Step 2.1, obtain the object space point coordinates according to the image point coordinates of the reference spectrum segment image and the object space elevation information, including performing coordinate orthographic conversion f for each pixel p ₂ (x _c2 ,y _c2 ) on the reference spectrum segment B2 image _1. Obtain the WGS84 geocentric Cartesian coordinates (X _P , Y _P , Z _P ) _WGS84 of its object space point P;

Step 2.2, obtain the image point coordinates of the non-reference spectrum segment, including using the geometric distortion parameters of the reference spectrum segment B2 obtained in step 1 relative to the non-reference spectrum segment B1 Construct a strict geometric imaging model of the non-reference spectrum segment B1; perform coordinate back-projection calculation f ₂ on the coordinates (X _P , Y _P , Z _P ) _WGS84 of the object space point P obtained in step 2.1, and obtain the correspondence on the image of the non-reference spectrum segment B1 The coordinates of the image point p ₁ (x _c1 , y _c1 );

Step 2.3, according to the coordinates (x _c1 , y _c1 ) of the corresponding image point p ₁ on the non-reference spectrum B1 image obtained in step 2.2, perform grayscale resampling to complete the precise registration of the non-reference spectrum B1 and the reference spectrum B2 .

2. as claimed in claim 1, based on the satellite multi-spectral image registration method of object side positioning consistency, it is characterized in that: the realization mode of coordinate front projection conversion and coordinate back projection conversion is as follows,

Let the coordinates of image point p ₂ on the reference spectrum B2 image be (x _c2 , y _c2 ), and the WGS84 geocentric Cartesian coordinates of the corresponding object space point P be (Xp, Yp, Zp) _WGS84 , let

{R R}_{B B S S} {R R}_{B B J J 22} {R R}_{T T 22} = = [\begin{matrix} {a a}_{11} & {b b}_{11} & {c c}_{11} \\ {a a}_{22} & {b b}_{22} & {c c}_{22} \\ {a a}_{33} & {b b}_{33} & {c c}_{33} \end{matrix}]

a ₁ , a ₂ , a ₃ , b ₁ , b ₂ , b ₃ , c ₁ , c ₂ , and c ₃ are elements of the matrix;

According to the following relationship between the _WGS84 geocentric Cartesian coordinates (Xp, Yp, Zp) of the object space point P and its geocentric coordinates (Bp, Lp, Hp),

{[\begin{matrix} X x p p \\ Y Y p p \\ Z Z p p \end{matrix}]}_{W W G G S S 8484} = = [\begin{matrix} ((N N + + H h p p)) \cdot &Center Dot; cos cos B B p p \cdot \cdot cos cos L L p p \\ ((N N + + H h p p)) \cdot \cdot c c o o s the s B B p p \cdot \cdot sin sin L L p p \\ ((N N \cdot \cdot ((11 - - {e e}^{22})) + + H h p p)) \cdot \cdot sin sin B B p p \end{matrix}]

Among them, e represents the oblateness of the earth ellipsoid, and the variable a represents the semi-major axis of the earth ellipsoid;

The relationship between the coordinates (x _c2 , y _c2 ) of the image point p ₂ and the earth coordinates (Bp, Lp, Hp) is obtained as follows,

\{\begin{matrix} ((N N + + H h p p)) \cdot \cdot cos cos B B p p \cdot &Center Dot; cos cos L L p p = = \frac{{a a}_{11} {x x}_{c c 22} + + {a a}_{22} {x x}_{c c 22} + + {a a}_{33} f f}{{c c}_{11} {x x}_{c c 22} + + {c c}_{22} {x x}_{c c 22} + + {c c}_{33} f f} \cdot &Center Dot; ((((N N \cdot \cdot ((11 - - {e e}^{22})) + + H h p p)) \cdot \cdot sin sin B B p p - - {Z Z}_{S S})) + + {X x}_{S S} \\ ((N N + + H h p p)) \cdot \cdot cos cos B B p p \cdot \cdot sin sin L L p p = = \frac{{b b}_{11} {x x}_{c c 22} + + {b b}_{22} {x x}_{c c 22} + + {b b}_{33} f f}{{c c}_{11} {x x}_{c c 22} + + {c c}_{22} {x x}_{c c 22} + + {c c}_{33} f f} \cdot \cdot ((((N N \cdot \cdot ((11 - - {e e}^{22})) + + H h p p)) \cdot \cdot S S i i n no B B p p - - {Z Z}_{S S})) + + {Y Y}_{S S} \end{matrix}

In coordinate orthographic conversion, the coordinates (x _c2 , y _c2 ) of the image point p ₂ and the given object elevation value Hp are used to calculate the geodetic latitude and longitude (Bp, Lp) of the object-space point P; then according to The geodetic coordinates (Bp, Lp, Hp) of the object space point P obtain the corresponding WGS84 geocentric Cartesian coordinates (Xp, Yp, Zp) _WGS84 ;

In coordinate back-projection conversion, based on the strict geometric imaging model of the reference spectrum segment B2, the WGS84 geocentric Cartesian coordinates (Xp, Yp, Zp) _WGS84 of the object space point P are calculated using the following formula,

\{\begin{matrix} {x x}_{c c 22} = = \frac{{a a}_{11} ((X x p p - - {X x}_{S S 22})) + + {b b}_{11} ((Y Y p p - - {Y Y}_{S S 22})) + + {c c}_{11} ((Z Z p p - - {Z Z}_{S S 22}))}{{a a}_{33} ((X x p p - - {X x}_{S S 22})) + + {b b}_{33} ((Y Y p p - - {Y Y}_{S S 22})) + + {c c}_{33} ((Z Z p p - - {Z Z}_{S S 22}))} \cdot \cdot f f \\ {y the y}_{c c 22} = = \frac{{a a}_{22} ((X x p p - - {X x}_{S S 22})) + + {b b}_{22} ((Y Y p p - - {Y Y}_{S S 22})) + + {c c}_{22} ((Z Z p p - - {Z Z}_{S S 22}))}{{a a}_{33} ((X x p p - - {X x}_{S S 22})) + + {b b}_{33} ((Y Y p p - - {Y Y}_{S S 22})) + + {c c}_{33} ((Z Z p p - - {Z Z}_{S S 22}))} \cdot \cdot f f \end{matrix}

The coordinates (x _c2 , y _c2 ) of the corresponding image point p ₂ are obtained.