CN103679711A - Method for calibrating in-orbit exterior orientation parameters of push-broom optical cameras of remote sensing satellite linear arrays - Google Patents

Abstract

Translated fromChinese

一种遥感卫星线阵推扫光学相机在轨外方位参数标定方法，步骤为：（1）采集高精度控制点信息；（2）获取控制点处相机内方位元素，以及轨道等辅助数据，建立各控制点处严格成像几何模型；（3）计算各控制点处的定位误差，剔除粗差控制点，得到初步检校控制点信息；（4）建立相机外方位元素几何检校数学模型，计算各控制点处理论指向矢量与实际指向矢量；（5）建立理论指向矢量与实际指向矢量之间的三轴旋转误差方程，对方程进行线性化处理；（6）最小二乘法解算误差方程参数，得到标定外参数。本发明利用标定外参数进行系统误差补偿后，相机无控定位精度稳定。

A calibration method for off-orbit azimuth parameters of linear array push-broom optical cameras of remote sensing satellites, the steps of which are: (1) collecting high-precision control point information; Strictly image the geometric model at each control point; (3) Calculate the positioning error at each control point, eliminate the gross error control points, and obtain the preliminary calibration control point information; Each control point processes the theoretical pointing vector and the actual pointing vector; (5) Establish the three-axis rotation error equation between the theoretical pointing vector and the actual pointing vector, and linearize the equation; (6) Solve the parameters of the error equation by the least square method , to get the calibrated extrinsic parameters. After the system error compensation is performed by using the calibration external parameters in the present invention, the uncontrolled positioning accuracy of the camera is stable.

Description

Translated fromChinese

一种遥感卫星线阵推扫光学相机在轨外方位参数标定方法A calibration method for off-orbit azimuth parameters of a linear array push-broom optical camera for remote sensing satellites

技术领域technical field

本发明涉及一种遥感卫星线阵推扫式光学相机在轨外方位参数标定方法，属于遥感卫星图像处理技术领域，用于遥感卫星搭载的线性推扫式相机的在轨定位精度优化提升。The invention relates to an off-orbit azimuth parameter calibration method of a linear array push-broom optical camera of a remote sensing satellite, belongs to the technical field of remote sensing satellite image processing, and is used for optimizing and improving the on-orbit positioning accuracy of a linear push-broom camera mounted on a remote sensing satellite.

背景技术Background technique

随着国内外遥感卫星的不断发射，遥感卫星的分辨率不断提高，卫星遥感的大面积、短周期、低成本等特点使其成为遥感发展的主导方向。截止2012年底，我国已经发射了环境、资源、遥感等多系列多颗光学遥感卫星，但是大量的国产遥感数据没有得到充分利用，究其原因，几何定位精度低是一个制约遥感数据的后续应用效果的重要因素。国外商业光学遥感卫星几何定位精度以及达到米级的精度，如Geoeye、WorldView等，而我国已发射的遥感卫星中资源一号02星定位精度为7Km，资源二号03星定位精度200m。国内遥感卫星无控定位精度差的一个重要原因是发射后地面没有对相机进行外方位元素的精确标定。With the continuous launch of remote sensing satellites at home and abroad, the resolution of remote sensing satellites has been continuously improved. The characteristics of large area, short period, and low cost of satellite remote sensing make it the dominant direction of remote sensing development. As of the end of 2012, my country has launched multiple series of optical remote sensing satellites for environment, resources, remote sensing, etc., but a large amount of domestic remote sensing data has not been fully utilized. The reason is that the low accuracy of geometric positioning is a follow-up application effect that restricts remote sensing data important factor. The geometric positioning accuracy of foreign commercial optical remote sensing satellites reaches meter-level accuracy, such as Geoeye and WorldView. Among the remote sensing satellites launched in my country, the positioning accuracy of the 02 star of Ziyuan No. 1 is 7Km, and the positioning accuracy of 03 satellite of Ziyuan No. 2 is 200m. An important reason for the poor positioning accuracy of domestic remote sensing satellites without control is that the ground has not accurately calibrated the camera's outer azimuth elements after launch.

线阵推扫式光学相机是目前光学遥感卫星传感器的主流，线阵推扫式光学相机随卫星沿着预先定义好的轨道向前推进动态成像，逐条扫描后形成一幅二维影像。在某摄影瞬间，影像上对应的行与所摄地面存在严格的中心投影关系。受地面测量误差、发射后环境冲力等因素的影响，遥感卫星发射后的外方位参数发生改变，使得其无控几何定位精度严重下降，地面处理中缺少对在轨卫星外方位元素的直接标定手段。The linear array push-broom optical camera is the mainstream of optical remote sensing satellite sensors at present. The linear array push-broom optical camera advances with the satellite along the predefined orbit for dynamic imaging, and scans one by one to form a two-dimensional image. At a moment of photography, there is a strict central projection relationship between the corresponding line on the image and the ground being photographed. Affected by factors such as ground measurement errors and post-launch environmental momentum, the outer azimuth parameters of remote sensing satellites have changed after launch, causing a serious drop in the accuracy of their uncontrolled geometric positioning. In ground processing, there is a lack of direct calibration means for the outer azimuth elements of satellites in orbit .

发明内容Contents of the invention

本发明解决的技术问题是：克服现有技术的不足，提供一种遥感卫星线阵推扫光学相机在轨外方位参数标定方法，利用一景影像采集精确地面控制点，可对相机的严格成像模型进行修正，解算精确外方位参数，可用于提升相机获取影像的无控几何定位精度，支持遥感卫星线阵推扫式光学相机的几何检校。The technical problem solved by the present invention is: to overcome the deficiencies of the prior art, to provide a remote sensing satellite linear array push-broom optical camera off-orbit azimuth parameter calibration method, using a scene image to collect accurate ground control points, which can strictly image the camera The model is corrected and the precise outer orientation parameters are calculated, which can be used to improve the uncontrolled geometric positioning accuracy of the image acquired by the camera, and support the geometric calibration of the linear array push-broom optical camera of remote sensing satellites.

本发明的技术方案是：一种遥感卫星线阵推扫光学相机在轨外方位参数标定方法，包括下列步骤：The technical solution of the present invention is: a method for calibrating the off-orbit azimuth parameters of a linear array push-broom optical camera of a remote sensing satellite, comprising the following steps:

1）在待检校影像I₁、控制影像I₂上利用控制点自动匹配算法，采集控制点信息，所述的控制点信息为T个匹配控制点对{PixI_1i，PixI_2i}，其中i=1,2,3,…,T）；记录每一个控制点对{PixI_1i，PixI_2i}的待检校影像I₁图像点坐标（i,j）以及对应控制影像特征点在WGS84坐标系下的位置坐标（lat,lon,height）；1) Use the control point automatic matching algorithm on the image I₁ to be checked and the control image I₂ to collect control point information. The control point information is T matching control point pairs {PixI_1i , PixI_2i }, where i =1,2,3,...,T); record the image point coordinates (i, j) of each control point pair {PixI_1i , PixI_2i } in the image to be checked I₁ and the corresponding control image feature points in the WGS84 coordinate system The position coordinates under (lat, lon, height);

2）获取各控制点处的相机内方位元素以及辅助数据，所述的辅助数据包括卫星的轨道、姿态、行时；建立各控制点处的严格成像模型$\left[\begin{array}{c}{X}_G\\ Y_G\\ Z_G\end{array}\right]=\begin{array}{c}\left[\begin{array}{c}{X}_s\\ Y_s\\ Z_s\end{array}\right]\end{array}+{\mathrm{uM}}_3{M}_2{M}_1{M}_0\left[\begin{array}{c}\tan \mathrm{psiX}\\ -\tan \mathrm{psiY}\\ 1\end{array}\right];$其中[X_S,Y_S,Z_S]为各控制点对应的行时时刻下卫星在协议地心坐标系中的位置，即控制点处轨道位置（PX,PY,PZ）；[X_G,Y_G,Z_G]为各控制点对应像元所对应的地面目标点在协议地心坐标系中的坐标；psiX、psiY分别为各控制点对应图像的像元主光轴单位矢量与卫星本体坐标系X轴、Y轴的夹角；u为比例因子；M₀为卫星本体坐标系相对于相机的安装矩阵，在卫星发射前由地面测量获取；M₁为各控制点对应的行时时刻下卫星至轨道坐标系旋转矩阵，由星上测量的姿态角构成；M₂为各控制点对应的行时时刻下轨道至J2000.0坐标系旋转矩阵；M₃为该时刻下J2000.0至WGS84坐标系旋转矩阵；2) Acquire the in-camera orientation elements and auxiliary data at each control point, the auxiliary data includes satellite orbit, attitude, and travel time; establish a strict imaging model at each control point$\left[\begin{array}{c}{x}_G\\ Y_G\\ Z_G\end{array}\right]=\begin{array}{c}\left[\begin{array}{c}{x}_{\mathrm{the\; s}}\\ Y_{\mathrm{the\; s}}\\ Z_{\mathrm{the\; s}}\end{array}\right]\end{array}+{\mathrm{uM}}_3{m}_2{m}_1{m}_0\left[\begin{array}{c}\mathrm{the\; tan}\mathrm{psiX}\\ -\mathrm{the\; tan}\mathrm{psiY}\\ 1\end{array}\right];$ Where [X_S , Y_S , Z_S ] is the position of the satellite in the agreed geocentric coordinate system at the travel time corresponding to each control point, that is, the orbital position (PX, PY, PZ) at the control point; [X_G , Y_G , Z_G ] are the coordinates of the ground target point corresponding to the pixel corresponding to each control point in the protocol geocentric coordinate system; The angle between the X-axis and the Y-axis of the coordinate system; u is the scale factor; M₀ is the installation matrix of the satellite body coordinate system relative to the camera, which is obtained by ground measurement before the satellite launch; M₁ is the travel time corresponding to each control point The rotation matrix from the satellite to the orbit coordinate system is composed of the attitude angle measured on the satellite; M₂ is the rotation matrix from the orbit to the J2000.0 coordinate system at the time corresponding to each control point; M₃ is the rotation matrix from J2000.0 to WGS84 coordinate system rotation matrix;

3）根据步骤2）得到的严格成像模型，计算获得各控制点处的定位误差大小以及方向，初步剔除粗差控制点，得到初始控制点对；3) According to the strict imaging model obtained in step 2), the magnitude and direction of the positioning error at each control point are calculated and obtained, and the gross error control points are preliminarily eliminated to obtain the initial control point pair;

4）建立相机外方位元素几何检校数学模型，根据步骤3）得到的初始控制点对，计算获得各控制点处理论指向矢量与实际指向矢量；4) Establish a mathematical model for geometric calibration of the camera's outer orientation elements, and calculate the theoretical pointing vector and actual pointing vector of each control point according to the initial control point pair obtained in step 3);

41）在严格成像模型中引入误差矩阵M，建立相机外方位元素几何检校数学模型41) Introduce the error matrix M into the strict imaging model, and establish a mathematical model for geometric calibration of the camera's external orientation elements

$\left[\begin{array}{c}{\mathrm{<span><span>X</span>x</span>}}_{\mathrm{<span><span>G</span>G</span>}}\\ {\mathrm{<span><span>Y</span>Y</span>}}_{\mathrm{<span><span>G</span>G</span>}}\\ {\mathrm{<span><span>Z</span>Z</span>}}_{\mathrm{<span><span>G</span>G</span>}}\end{array}\right]<span><span>=</span>=</span>\left[\begin{array}{c}{\mathrm{<span><span>X</span>x</span>}}_{\mathrm{<span><span>s</span>the\; s</span>}}\\ {\mathrm{<span><span>Y</span>Y</span>}}_{\mathrm{<span><span>s</span>the\; s</span>}}\\ {\mathrm{<span><span>Z</span>Z</span>}}_{\mathrm{<span><span>s</span>the\; s</span>}}\end{array}\right]<span><span>+</span>+</span>{\mathrm{<span><span>uM</span>uM</span>}}_{\mathrm{<span><span>3</span>3</span>}}{\mathrm{<span><span>M</span>m</span>}}_{\mathrm{<span><span>2</span>2</span>}}{\mathrm{<span><span>M</span>m</span>}}_{\mathrm{<span><span>1</span>1</span>}}{\mathrm{<span><span>M</span>m</span>}}_{\mathrm{<span><span>0</span>0</span>}}\mathrm{<span><span>M</span>m</span>}\left[\begin{array}{c}\mathrm{<span><span>tan</span>the\; tan</span>}\mathrm{<span><span>psiX</span>psiX</span>}\\ <span><span>-</span>-</span>\mathrm{<span><span>tan</span>the\; tan</span>}\mathrm{<span><span>psiY</span>psiY</span>}\\ \mathrm{<span><span>1</span>1</span>}\end{array}\right]$

42）对任一组控制点，利用控制点处的内方位元素（psiX，psiY）计算控制点处光轴在相机坐标系中指向矢量V，将其进行归一化处理后作为该控制点处的实际指向矢量V_实际；42) For any set of control points, use the internal orientation elements (psiX, psiY) at the control point to calculate the optical axis at the control point pointing to the vector V in the camera coordinate system, and normalize it as the control point The actual pointing vector_Vactual ;

43）对任一组控制点，利用控制点地面点位置坐标（lat,lon,height）计算控制点处光轴在相机坐标系中指向矢量V，将其进行归一化处理后作为该控制点处的理论指向矢量V_理论；43) For any set of control points, use the position coordinates (lat, lon, height) of the ground point of the control point to calculate the pointing vector V of the optical axis at the control point in the camera coordinate system, and normalize it as the control point The theory at pointing vector V_theory ;

其中$\left[\begin{array}{c}{X}_G\\ Y_G\\ Z_G\end{array}\right]=\left[\begin{array}{c}(a+\mathrm{height})\&times;\cos \left(\mathrm{lon}\right)\&times;\cos \left(\mathrm{lat}\right)\\ (a+\mathrm{height})\&times;\sin \left(\mathrm{lon}\right)\&times;\cos \left(\mathrm{lat}\right)\\ (b+\mathrm{height})\&times;\sin \left(\mathrm{lat}\right)\end{array}\right],$a为地球长半轴，b为地球短半轴；in$\left[\begin{array}{c}{x}_G\\ Y_G\\ Z_G\end{array}\right]=\left[\begin{array}{c}(a+\mathrm{height})\&times;\cos \left(\mathrm{the\; lon}\right)\&times;\cos \left(\mathrm{lat}\right)\\ (a+\mathrm{height})\&times;\sin \left(\mathrm{the\; lon}\right)\&times;\cos \left(\mathrm{lat}\right)\\ (b+\mathrm{height})\&times;\sin \left(\mathrm{lat}\right)\end{array}\right],$ a is the semi-major axis of the earth, b is the semi-minor axis of the earth;

5）建立理论指向矢量与实际指向矢量之间的三轴旋转误差方程，对该三轴旋转误差方程进行线性化处理，得到线性方程组；5) Establish a three-axis rotation error equation between the theoretical pointing vector and the actual pointing vector, and linearize the three-axis rotation error equation to obtain a linear equation set;

51）设误差矩阵M为三轴旋转正交矩阵，建立理论指向矢量与实际指向矢量之间的矩阵转换模型：51) Let the error matrix M be a three-axis rotation orthogonal matrix, and establish a matrix transformation model between the theoretical pointing vector and the actual pointing vector:

其中$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]$为理论指向矢量，$\left[\begin{array}{c}x\\ y\\ z\end{array}\right]$为实际指向矢量；

in$\left[\begin{array}{c}x\\ Y\\ Z\end{array}\right]$ is the theoretical pointing vector,$\left[\begin{array}{c}x\\ \mathrm{the\; y}\\ z\end{array}\right]$ is the actual pointing vector;

52）将模型进行转换，得到理论指向矢量各参数与实际指向矢量各参数的非线性方程；52) Convert the model to obtain the nonlinear equations of the parameters of the theoretical pointing vector and the parameters of the actual pointing vector;

53）对非线性方程在原点（0，0，0）处进行泰勒展开，得到非线性方程的线性近似方程；53) Perform Taylor expansion on the nonlinear equation at the origin (0, 0, 0) to obtain the linear approximation equation of the nonlinear equation;

54）针对每一组理论指向矢量和实际指向矢量，计算线性方程参数，得到线性方程组；54) For each set of theoretical pointing vectors and actual pointing vectors, calculate the parameters of the linear equation to obtain a set of linear equations;

6）采用最小二乘法迭代解算步骤5）得到的线性方程组，获得标定外参数。6) Use the least square method to iteratively solve the linear equations obtained in step 5) to obtain the calibration external parameters.

步骤1）中采集控制点信息的具体方法为：The specific method of collecting control point information in step 1) is as follows:

11）基于SIFT算法对待检校原始图像I₁进行特征点提取，得到M个特征点PixI_1i（i=1,2,3,…,M），M为正整数；记录各特征点的SIFT特征向量；11) Based on the SIFT algorithm, extract the feature points of the original image I₁ to be corrected, and obtain M feature points PixI_1i (i=1,2,3,...,M), M is a positive integer; record the SIFT features of each feature point vector;

12）基于SIFT算法对高精度控制影像I₂进行特征点提取，得到N个特征点PixI_2i（i=1,2,3,…,N），N为正整数；记录各特征点处的SIFT特征向量；12) Based on the SIFT algorithm, perform feature point extraction on the high-precision control image I₂ to obtain N feature points PixI_2i (i=1,2,3,…,N), where N is a positive integer; record the SIFT at each feature point Feature vector;

13）采用欧式距离作为相似性度量准则对两幅影像的特征点进行匹配，得到T个匹配控制点对{PixI_1i，PixI_2i}（i=1,2,3,…,T）。13) Use the Euclidean distance as the similarity measure criterion to match the feature points of the two images, and get T matching control point pairs {PixI_1i , PixI_2i } (i=1,2,3,...,T).

步骤2）中建立各控制点处的严格成像模型的具体方法为：The specific method of establishing the strict imaging model at each control point in step 2) is as follows:

21）根据控制点对{PixI_1i，PixI_2i}的图像坐标（i,j），获取其对应的行时；21) According to the image coordinates (i, j) of the control point pair {PixI_1i , PixI_2i }, obtain its corresponding row time;

22）利用拉格朗日插值算法插值计算控制点对{PixI_1i，PixI_2i}对应行时的轨道位置（PX,PY,PZ,VX,VY,VZ）；22) Use the Lagrange interpolation algorithm to interpolate and calculate the orbital position (PX, PY, PZ, VX, VY, VZ) of the control point pair {PixI_1i , PixI_2i } corresponding to the row time;

23）对卫星下传的姿态四元数进行坐标系转换处理，利用拉格朗日算法插值计算控制点对{PixI_1i，PixI_2i}对应行时时刻处相机相对于轨道坐标系的三轴姿态角（Roll，Pitch，Yaw）；23) Perform coordinate system conversion processing on the attitude quaternion transmitted by the satellite, and use the Lagrangian algorithm to interpolate and calculate the three-axis attitude of the camera relative to the orbital coordinate system at the corresponding travel time of the control point pair {PixI_1i , PixI_2i } angle(Roll, Pitch, Yaw);

24）根据相机实验室测量数据或内方位元素几何检校结果，计算控制点对{PixI_1i，PixI_2i}对应的CCD探元光轴指向角（psiX,psiY）；24) According to the camera laboratory measurement data or the geometric calibration results of the internal orientation elements, calculate the CCD detector optical axis pointing angle (psiX, psiY) corresponding to the control point pair {PixI_1i , PixI_2i };

25）根据步骤21）-步骤24）得到的结果，建立控制点处的严格成像模型，针对线阵推扫成像相机某一时刻获取的图像上的某一点，构建遥感影像的严格成像模型；25) According to the results obtained in step 21)-step 24), a strict imaging model at the control point is established, and a strict imaging model of remote sensing images is constructed for a certain point on the image acquired by the linear push-broom imaging camera at a certain moment;

步骤3）中初始控制点对的具体获取方法为：The specific method of obtaining the initial control point pair in step 3) is:

31）根据步骤2）得到的严格成像模型，计算各控制点对{PixI_1i，PixI_2i}中待检校影像点PixI_1i的初始地面点位置坐标{lon’,lat’}；31) According to the strict imaging model obtained in step 2), calculate the initial ground point position coordinate {lon',lat'} of the image point PixI_1i to be checked in each control point pair {PixI_1i , PixI_2i };

32）根据各控制点对{PixI_1i，PixI_2i}中控制影像点PixI_2i位置坐标信息，计算控制点PixI_1i的初始定位误差{Δxi,Δyi}；32) Calculate the initial positioning error {_Δxi ,Δyi} of the control point PixI_{1i according to the position coordinate information of the control image point PixI 2iin} each control point pair {PixI 1i , PixI_2i };

33）计算控制点对{PixI_1i，PixI_2i}定位误差平均值{Δx,Δy}；33) Calculate the average positioning error {Δx,Δy} of the control point pair {PixI_1i , PixI_2i };

34）对每一个控制点对，判断其初始定位误差{Δxi,Δyi}是否在{Δx±20%Δx,Δy±20%Δy}范围内，将不在{Δx±20%Δx,Δy±20%Δy}范围内的控制点对剔除；34) For each control point pair, judge whether its initial positioning error {Δxi, Δyi} is within the range of {Δx±20%Δx, Δy±20%Δy}, and will not be within {Δx±20%Δx, Δy±20% The control point pairs within the range of Δy} are eliminated;

35）剔除误差控制点对后，得到K组初始控制点对{PixI_1i，PixI_2i}（i=1,2,3,…,K）。35) After eliminating error control point pairs, K groups of initial control point pairs {PixI_1i , PixI_2i } (i=1,2,3,...,K) are obtained.

本发明与现有技术相比的有益效果在于：The beneficial effect of the present invention compared with prior art is:

（1）本发明对遥感卫星线阵推扫式光学相机严格成像模型以及辅助数据特性进行了深入分析，设计了包含影像高精度地面控制点信息获取、相机轨道、姿态、行时等辅助数据优化处理、严格几何检校数学模型建立、理想指向矢量与实际指向矢量计算、三轴误差方程建立以及方程线性化处理、最小二乘解算外方位参数等步骤的高精度外参数检校流程，满足卫星在轨几何检校需求。(1) The present invention conducts an in-depth analysis of the strict imaging model and auxiliary data characteristics of the remote sensing satellite linear array push-broom optical camera, and designs auxiliary data optimization including image high-precision ground control point information acquisition, camera orbit, attitude, travel time, etc. Processing, strict geometric calibration mathematical model establishment, ideal pointing vector and actual pointing vector calculation, three-axis error equation establishment and equation linearization processing, least squares calculation of external azimuth parameters and other steps of high-precision external parameter calibration process, to meet Satellite in-orbit geometric calibration requirements.

（2）在高精度地面控制点信息提取方面，针对基本影像与高精度控制点影像的旋转、尺度等特点，采用了基于SIFT算法的控制点提取以及匹配，提高了控制点提取的准确度。(2) In the aspect of high-precision ground control point information extraction, according to the characteristics of rotation and scale of the basic image and high-precision control point image, the control point extraction and matching based on the SIFT algorithm is used to improve the accuracy of control point extraction.

（3）在整个检校过程中，考虑了基于严格成像模型的初始控制点容错，利用影像各点定位误差基本一致的原则对匹配控制点对进行筛选，保证采集控制点对的正确性。(3) During the entire calibration process, the fault tolerance of the initial control points based on the strict imaging model is considered, and the matching control point pairs are screened using the principle that the positioning errors of each point in the image are basically consistent to ensure the correctness of the collected control point pairs.

（4）在建立几何检校数学模型时，引入三轴旋转正交矩阵作为误差矩阵，将误差方程线性化后，在降低误差方程解算的难度的基础上很好的保障了线性化误差方程对真实误差的拟合相似度。(4) When establishing the geometric calibration mathematical model, the three-axis rotation orthogonal matrix is introduced as the error matrix. After the error equation is linearized, the linearized error equation is well guaranteed on the basis of reducing the difficulty of solving the error equation. Fit similarity to true error.

（5）在解算外方位参数时采用了基于最小二乘的迭代误差剔除设计，可以进一步剔除粗差控制点，提高了参数检校的精度。(5) The iterative error elimination design based on least squares is adopted when calculating the external orientation parameters, which can further eliminate gross error control points and improve the accuracy of parameter calibration.

附图说明Description of drawings

图1为本发明方法流程图。Fig. 1 is a flow chart of the method of the present invention.

具体实施方式Detailed ways

下面结合附图1对本发明的具体实施方式进行进一步的详细描述：Below in conjunction with accompanying drawing 1, the specific embodiment of the present invention is described in further detail:

1.在待检校影像I₁、高精度控制影像I₂上利用控制点自动匹配算法采集高精度控制点对信息，1. Use the control point automatic matching algorithm to collect high-precision control point pair information on the image I₁ to be checked and the high-precision control image I₂ ,

（1）基于SIFT算法对待检校原始图像I₁进行特征点提取，得到m个特征点PixI_1i（i=1,2…,m），记录各特征点的SIFT特征向量；(1) Extract the feature points of the original image I₁ to be corrected based on the SIFT algorithm to obtain m feature points PixI_1i (i=1,2...,m), and record the SIFT feature vectors of each feature point;

SIFT特征匹配算法表征如下：The SIFT feature matching algorithm is characterized as follows:

（1.1）确定特征点位置坐标与所在尺度。建立图像高斯金字塔，在金字塔尺度空间中的26个邻域中检测极值，若某点（x，y）在金字塔尺度空间本层以及上下两层的26个邻域中是最大或最小值时，定义该点是图像在该尺度L下的一个特征点。(1.1) Determine the position coordinates and scale of the feature points. Build an image Gaussian pyramid, and detect extreme values in the 26 neighborhoods in the pyramid scale space. If a point (x, y) is the maximum or minimum value in the pyramid scale space and the 26 neighborhoods in the upper and lower layers , defining this point as a feature point of the image at this scale L.

（1.2）特征点方向参数计算。在以特征点（x，y）为中心的邻域窗口内采样，用直方图统计领域像素的梯度方向。将梯度直方图的范围划分为36个方向，每10度代表一个方向。定义直方图峰值所在的方向为该特征点的方向参数。(1.2) Calculation of feature point direction parameters. Sampling in the neighborhood window centered on the feature point (x, y), and using the histogram to count the gradient direction of the domain pixels. Divide the range of the gradient histogram into 36 directions, and every 10 degrees represents a direction. Define the direction where the histogram peak is located as the direction parameter of the feature point.

（1.3）计算SIFT特征向量描述。在特征点中心周围取8*8的窗口，将窗口切成2*2的子窗口。在每个子窗口上计算8个方向的梯度方向直方图，统计每个方向的累加值，作为该子窗口的方向信息。4个子窗口统计完成后最终生成该特征点处的32维度特征向量。(1.3) Calculate the SIFT feature vector description. Take an 8*8 window around the center of the feature point, and cut the window into 2*2 sub-windows. Calculate the gradient direction histogram of 8 directions on each sub-window, and count the cumulative value of each direction as the direction information of the sub-window. After the statistics of the four sub-windows are completed, the 32-dimensional feature vector at the feature point is finally generated.

（2）基于SIFT算法对高精度控制影像I₂进行特征点提取，得到n个特征点PixI_2i（i=1,2…,n），记录各特征点处的SIFT特征向量；(2) Based on the SIFT algorithm, extract the feature points of the high-precision control image I₂ to obtain n feature points PixI_2i (i=1,2...,n), and record the SIFT feature vectors at each feature point;

（3）采用欧式距离作为相似性度量准则对两幅影像的特征点进行匹配，得到M个匹配控制点对{PixI_1i，PixI_2i}(i=1,2…,M）；(3) Use the Euclidean distance as the similarity measure criterion to match the feature points of the two images, and obtain M matching control point pairs {PixI_1i , PixI_2i }(i=1,2...,M);

（3.1）对于待检校原始图像I₁的任一控制点PixI_1i，计算其SIFT特征向量与高精度控制影像I₂上提取得到的m个特征向量之间的欧式距离。对于两个向量l₁(x₁,x₂,…,x_n)、l₂(y₁,y₂,…,y_n)，其欧式距离表征如下：(3.1) For any control point PixI_1i of the original image I₁ to be checked, calculate the Euclidean distance between its SIFT feature vector and the m feature vectors extracted from the high-precision control image I₂ . For two vectors l₁ (x₁ ,x₂ ,…,x_n ), l₂ (y₁ ,y₂ ,…,y_n ), the Euclidean distance is characterized as follows:

$\mathrm{<span><span>\&Delta;l</span>\&Delta;l</span>}<span><span>=</span>=</span>\sqrt{\underset{\mathrm{<span><span>i</span>i</span>}<span><span>=</span>=</span>\mathrm{<span><span>0</span>0</span>}}{\overset{\mathrm{<span><span>n</span>no</span>}}{\mathrm{<span><span>\&Sigma;</span>\&Sigma;</span>}}}{<span><span>(</span>(</span>{\mathrm{<span><span>x</span>x</span>}}_{\mathrm{<span><span>i</span>i</span>}}<span><span>-</span>-</span>{\mathrm{<span><span>y</span>the\; y</span>}}_{\mathrm{<span><span>i</span>i</span>}}<span><span>)</span>)</span>}^{\mathrm{<span><span>2</span>2</span>}}}$

（3.2）计算m个欧式距离中的最小值，最小值所在的特征点PixI_2j即为待检校原始图像I₁特征点PixI_1i的匹配特征点。(3.2) Calculate the minimum value among the m Euclidean distances, and the feature point PixI_2j where the minimum value is located is the matching feature point of the feature point PixI_1i of the original image I₁ to be checked.

（4）对每一个匹配特征点对{PixI_1i，PixI_2j}，记录待检校原始图像I₁特征点PixI_1i处像点平面坐标（sample,line），以及高精度控制影像I₂特征点PixI_2j处三维位置坐标（lat,lon,height），作为采集到的高精度控制点信息。(4) For each matching feature point pair {PixI_1i , PixI_2j }, record the image point plane coordinates (sample, line) of the original image I₁ feature point PixI_1i to be checked, and the high-precision control image I₂ feature point The three-dimensional position coordinates (lat, lon, height) at PixI_2j are used as the collected high-precision control point information.

2.获取各控制点处内方位元素、轨道、姿态、行时等辅助数据，建立各控制点处的严格成像模型。2. Obtain auxiliary data such as internal orientation elements, orbits, attitudes, and travel times at each control point, and establish a strict imaging model at each control point.

（1）从卫星下传的原始数据中解析待检校影像成像时间范围内的轨道、姿态、行时等数据；(1) Analyze the orbit, attitude, travel time and other data within the imaging time range of the image to be checked from the original data downloaded by the satellite;

（2）对任一控制点对Points_i{sample,line,lat,lon,height}，根据控制点对的图像坐标（sample,line），获取其对应的摄影时间scanTime；(2) For any control point pair Points_i {sample, line, lat, lon, height}, according to the image coordinates (sample, line) of the control point pair, obtain its corresponding photography time scanTime;

控制点对应的摄影时间可以从卫星下传的辅助数据中直接解析，第line行辅助数据对应的成像时间即为该控制点对应的摄影时间。The photographing time corresponding to the control point can be directly analyzed from the auxiliary data downloaded from the satellite, and the imaging time corresponding to the auxiliary data in the first row is the photographing time corresponding to the control point.

（3）利用拉格朗日插值算法插值计算摄影时间scanTime的卫星轨道位置（PX,PY,PZ,VX,VY,VZ）；卫星按照一定频次下传轨道数据，因此，摄影时间scanTime对应的卫星轨道位置需要利用摄影时刻前后几组的轨道数据进行插值计算。本发明采用拉格朗日插值算法，利用摄影时间前后三组轨道数据计算摄影时间的卫星轨道位置。(3) Use the Lagrangian interpolation algorithm to interpolate and calculate the satellite orbit position (PX, PY, PZ, VX, VY, VZ) of the photography time scanTime; satellites download orbit data according to a certain frequency, therefore, the satellite corresponding to the photography time scanTime The orbital position needs to be interpolated using several sets of orbital data before and after the shooting time. The invention adopts the Lagrangian interpolation algorithm, and uses three sets of orbit data before and after the photographing time to calculate the satellite orbit position at the photographing time.

（3.1）从第一组轨道数据开始，判断该组轨道数据的生成时间、下一组轨道数据的生成时间与scanTime之间的关系，若scanTime大于第i组轨道数据的生成时间，同时小于第i+1组轨道数据的生成时间，则记录i为距离摄影时间最近的轨道数据序号。(3.1) Starting from the first group of orbit data, judge the relationship between the generation time of this group of orbit data, the generation time of the next group of orbit data and scanTime, if scanTime is greater than the generation time of the i-th group of orbit data, and less than the generation time of the i-th group of orbit data The generation time of the i+1 group of track data, record i as the track data sequence number closest to the shooting time.

（3.2）利用第i-1，i,i+1组轨道数据，基于拉格朗日算法计算摄影时刻的轨道位置以及速度。拉格朗日插值算法表述如下：(3.2) Use the i-1, i, and i+1 sets of orbital data to calculate the orbital position and velocity at the time of photography based on the Lagrangian algorithm. The Lagrangian interpolation algorithm is expressed as follows:

对于已知y=f(x)的函数表(x_i,f(x_i))（i=0,1,…,n），对于在[x_o,x_n]范围内任一x，有：For the function table (x_i ,f(x_i )) (i=0,1,…,n) of known y=f(x), for any x in the range of [x_o ,x_n ], we have :

${\mathrm{<span><span>P</span>P</span>}}_{\mathrm{<span><span>n</span>no</span>}}<span><span>(</span>(</span>\mathrm{<span><span>x</span>x</span>}<span><span>)</span>)</span><span><span>=</span>=</span>\underset{\mathrm{<span><span>i</span>i</span>}<span><span>=</span>=</span>\mathrm{<span><span>0</span>0</span>}}{\overset{\mathrm{<span><span>n</span>no</span>}}{\mathrm{<span><span>\&Sigma;</span>\&Sigma;</span>}}}<span><span>(</span>(</span>{\mathrm{<span><span>y</span>the\; y</span>}}_{\mathrm{<span><span>i</span>i</span>}}<span><span>*</span>*</span>\underset{\mathrm{<span><span>j</span>j</span>}<span><span>=</span>=</span>\mathrm{<span><span>0</span>0</span>}<span><span>,</span>,</span>\mathrm{<span><span>j</span>j</span>}<span><span>\&NotEqual;</span>\&NotEqual;</span>\mathrm{<span><span>i</span>i</span>}}{\overset{\mathrm{<span><span>n</span>no</span>}}{\mathrm{<span><span>\&Pi;</span>\&Pi;</span>}}}\frac{\mathrm{<span><span>x</span>x</span>}<span><span>-</span>-</span>{\mathrm{<span><span>x</span>x</span>}}_{\mathrm{<span><span>j</span>j</span>}}}{{\mathrm{<span><span>x</span>x</span>}}_{\mathrm{<span><span>i</span>i</span>}}<span><span>-</span>-</span>{\mathrm{<span><span>x</span>x</span>}}_{\mathrm{<span><span>j</span>j</span>}}}<span><span>)</span>)</span>$

（4）利用拉格朗日算法插值计算摄影时间scanTime相机相对于轨道坐标系的三轴姿态角（Roll,Pitch,Yaw）；(4) Use the Lagrange algorithm to interpolate to calculate the three-axis attitude angle (Roll, Pitch, Yaw) of the scanTime camera relative to the orbital coordinate system;

（5）根据相机实验室测量数据，读取控制点控制点对的图像坐标（sample,line）对应的第sample个CCD探元的光轴指向角（psiX,psiY）。(5) According to the camera laboratory measurement data, read the optical axis pointing angle (psiX, psiY) of the first sample CCD detector corresponding to the image coordinates (sample, line) of the control point control point pair.

（6）建立控制点处的严格成像模型，针对任一控制点，利用其内方位、轨道、姿态、行时等数据构建遥感影像的严格几何成像模型Model(6) Establish a strict imaging model at the control point. For any control point, use its internal orientation, orbit, attitude, travel time and other data to construct a strict geometric imaging model Model of remote sensing images

线阵推扫相机的严格成像几何模型如下所述。The strict imaging geometry model of the linear push-broom camera is described below.

$\left[\begin{array}{c}{\mathrm{<span><span>X</span>x</span>}}_{\mathrm{<span><span>G</span>G</span>}}\\ {\mathrm{<span><span>Y</span>Y</span>}}_{\mathrm{<span><span>G</span>G</span>}}\\ {\mathrm{<span><span>Z</span>Z</span>}}_{\mathrm{<span><span>G</span>G</span>}}\end{array}\right]<span><span>=</span>=</span>\left[\begin{array}{c}{\mathrm{<span><span>X</span>x</span>}}_{\mathrm{<span><span>s</span>the\; s</span>}}\\ {\mathrm{<span><span>Y</span>Y</span>}}_{\mathrm{<span><span>s</span>the\; s</span>}}\\ {\mathrm{<span><span>Z</span>Z</span>}}_{\mathrm{<span><span>s</span>the\; s</span>}}\end{array}\right]<span><span>+</span>+</span>{\mathrm{<span><span>uM</span>uM</span>}}_{\mathrm{<span><span>3</span>3</span>}}{\mathrm{<span><span>M</span>m</span>}}_{\mathrm{<span><span>2</span>2</span>}}{\mathrm{<span><span>M</span>m</span>}}_{\mathrm{<span><span>1</span>1</span>}}{\mathrm{<span><span>M</span>m</span>}}_{\mathrm{<span><span>0</span>0</span>}}\left[\begin{array}{c}\mathrm{<span><span>tan</span>the\; tan</span>}\mathrm{<span><span>psiX</span>psiX</span>}\\ <span><span>-</span>-</span>\mathrm{<span><span>tan</span>the\; tan</span>}\mathrm{<span><span>psiY</span>psiY</span>}\\ \mathrm{<span><span>1</span>1</span>}\end{array}\right]$

其中：in:

[X_S,Y_S,Z_S]为该时刻下卫星在协议地心坐标系中的位置，即控制点处轨道位置（PX,PY,PZ）；[X_S , Y_S , Z_S ] is the position of the satellite in the agreed geocentric coordinate system at this moment, that is, the orbit position at the control point (PX, PY, PZ);

[X_G,Y_G,Z_G]为该像元对应的地面目标点在协议地心坐标系中的坐标。[X_G , Y_G , Z_G ] are the coordinates of the ground target point corresponding to the pixel in the agreed geocentric coordinate system.

psiX、psiY分别为图像对应的相机像元主光轴单位矢量与卫星本体坐标系X轴、Y轴的夹角psiX and psiY are the angles between the main optical axis unit vector of the camera pixel corresponding to the image and the X-axis and Y-axis of the satellite body coordinate system

u—比例因子。u—scale factor.

M₀为卫星本体坐标系相对于相机的安装矩阵，在卫星发射前由地面测量获取；M₀ is the installation matrix of the satellite body coordinate system relative to the camera, which is obtained by ground measurement before the satellite is launched;

M₁为该时刻下卫星至轨道坐标系旋转矩阵，由星上测量的姿态角构成。M₁ is the rotation matrix from the satellite to the orbit coordinate system at this moment, which is composed of the attitude angle measured on the satellite.

M₂为该时刻下轨道至J2000.0坐标系旋转矩阵，由升交点赤经、轨道倾角、幅角等构成。M₂ is the rotation matrix from the lower orbit to the J2000.0 coordinate system at this moment, which is composed of right ascension of ascending node, orbital inclination, argument, etc.

M₃为该时刻下J2000.0至WGS84坐标系旋转矩阵，需进行岁差改正、章动改正、格林尼治恒星时改正和极移改正。M₃ is the rotation matrix of the J2000.0 to WGS84 coordinate system at this moment, which needs to be corrected for precession, nutation, Greenwich mean sidereal time and pole shift.

3.计算各控制点处的定位误差大小以及方向，初步剔除粗差控制点，得到初始控制点对3. Calculate the magnitude and direction of the positioning error at each control point, preliminarily eliminate the gross error control points, and obtain the initial control point pair

（1）根据建立的严格成像模型，基于各控制点处内方位元素、轨道、姿态、行时等辅助数据，计算各控制点中原始图像点（sample,line）的初始地面点位置坐标（lon’,lat’）(1) According to the strict imaging model established, based on the auxiliary data such as orientation elements, orbits, attitudes, and travel times at each control point, calculate the initial ground point position coordinates (lon ',lat')

（2）根据各控制点高精度控制影像点位置(lon,lat,height)坐标信息，计算图像点的初始定位误差（Δxi,Δyi）。对于每一个控制点，(2) Calculate the initial positioning error (Δxi, Δyi) of the image point according to the coordinate information of each control point to control the position (lon, lat, height) of the image point with high precision. For each control point,

（2.1）将初始地面点位置坐标（lon’,lat’）经UTM投影后，得到其平面坐标（x’,y’）(2.1) After the initial ground point position coordinates (lon’, lat’) are projected by UTM, their plane coordinates (x’, y’) are obtained

（2.2）将控制点影像点坐标(lon,lat,height)经UTM投影后，得到其平面坐标(x,y)(2.2) After projecting the control point image point coordinates (lon, lat, height) through UTM, its plane coordinates (x, y) are obtained

（2.3）定义初始定位误差（Δxi,Δyi）=（x-x’,y-y’）(2.3) Define the initial positioning error (Δxi, Δyi) = (x-x’, y-y’)

（3）计算控制点对定位误差平均值（Δx,Δy）(3) Calculate the average value of the positioning error of the control point pair (Δx, Δy)

（4）对每一个控制点对，判断其初始定位误差{Δxi,Δyi}是否在[Δx±20%Δx,Δy±20%Δy]范围内，若不在该范围内，将该控制点剔除。剔除误差控制点后，得到初始控制点对，用于建立外方位元素几何检校数学模型(4) For each control point pair, judge whether its initial positioning error {Δxi, Δyi} is within the range of [Δx±20%Δx, Δy±20%Δy], and if not, remove the control point. After the error control points are eliminated, the initial control point pairs are obtained, which are used to establish the geometric calibration mathematical model of the outer orientation elements

4.建立相机外方位元素几何检校数学模型，计算各控制点处理论指向矢量与实际指向矢量。4. Establish a mathematical model for geometric calibration of the camera's external orientation elements, and calculate the theoretical pointing vector and actual pointing vector of each control point.

（1）建立相机外方位元素几何检校数学模型，在严格成像几何模型中引入误差矩阵M，如下所示(1) Establish a mathematical model for geometric calibration of the camera's external orientation elements, and introduce the error matrix M into the strict imaging geometric model, as shown below

$\left[\begin{array}{c}{\mathrm{<span><span>X</span>x</span>}}_{\mathrm{<span><span>G</span>G</span>}}\\ {\mathrm{<span><span>Y</span>Y</span>}}_{\mathrm{<span><span>G</span>G</span>}}\\ {\mathrm{<span><span>Z</span>Z</span>}}_{\mathrm{<span><span>G</span>G</span>}}\end{array}\right]<span><span>=</span>=</span>\left[\begin{array}{c}{\mathrm{<span><span>X</span>x</span>}}_{\mathrm{<span><span>s</span>the\; s</span>}}\\ {\mathrm{<span><span>Y</span>Y</span>}}_{\mathrm{<span><span>s</span>the\; s</span>}}\\ {\mathrm{<span><span>Z</span>Z</span>}}_{\mathrm{<span><span>s</span>the\; s</span>}}\end{array}\right]<span><span>+</span>+</span>{\mathrm{<span><span>uM</span>uM</span>}}_{\mathrm{<span><span>3</span>3</span>}}{\mathrm{<span><span>M</span>m</span>}}_{\mathrm{<span><span>2</span>2</span>}}{\mathrm{<span><span>M</span>m</span>}}_{\mathrm{<span><span>1</span>1</span>}}{\mathrm{<span><span>M</span>m</span>}}_{\mathrm{<span><span>0</span>0</span>}}\mathrm{<span><span>M</span>m</span>}\left[\begin{array}{c}\mathrm{<span><span>tan</span>the\; tan</span>}\mathrm{<span><span>psiX</span>psiX</span>}\\ <span><span>-</span>-</span>\mathrm{<span><span>tan</span>the\; tan</span>}\mathrm{<span><span>psiY</span>psiY</span>}\\ \mathrm{<span><span>1</span>1</span>}\end{array}\right]$

（2）对任一组控制点，计算原始图像像点（sample,line）处的实际指向矢量V_实际：(2) For any set of control points, calculate the actual pointing vector_Vactual at the original image point (sample, line):

利用控制点处的内方位元素（psiX，psiY）计算控制点处光轴在相机坐标系中指向矢量V，将其进行归一化处理后作为该控制点处的实际指向矢量V_实际；Use the internal orientation elements (psiX, psiY) at the control point to calculate the pointing vector V of the optical axis at the control point in the camera coordinate system, and normalize it as the actual pointing vector Vactual_at the control point;

（3）对任一组控制点，计算控制点处的理论指向矢量V_理论：(3) For any set of control points, calculate the theoretical pointing vector V_theory at the control point:

（3.1）利用控制点地面点位置坐标（lat,lon,height）计算地面点在协议地心坐标系中的三轴位置（XG,YG,ZG），计算公式如下：(3.1) Use the ground point position coordinates (lat, lon, height) of the control point to calculate the three-axis position (XG, YG, ZG) of the ground point in the protocol geocentric coordinate system. The calculation formula is as follows:

$\left[\begin{array}{c}{\mathrm{<span><span>X</span>x</span>}}_{\mathrm{<span><span>G</span>G</span>}}\\ {\mathrm{<span><span>Y</span>Y</span>}}_{\mathrm{<span><span>G</span>G</span>}}\\ {\mathrm{<span><span>Z</span>Z</span>}}_{\mathrm{<span><span>G</span>G</span>}}\end{array}\right]<span><span>=</span>=</span>\left[\begin{array}{c}<span><span>(</span>(</span>\mathrm{<span><span>a</span>a</span>}<span><span>+</span>+</span>\mathrm{<span><span>height</span>height</span>}<span><span>)</span>)</span><span><span>\&times;</span>\&times;</span>\mathrm{<span><span>cos</span>cos</span>}<span><span>(</span>(</span>\mathrm{<span><span>lon</span>the\; lon</span>}<span><span>)</span>)</span><span><span>\&times;</span>\&times;</span>\mathrm{<span><span>cos</span>cos</span>}<span><span>(</span>(</span>\mathrm{<span><span>lat</span>lat</span>}<span><span>)</span>)</span>\\ <span><span>(</span>(</span>\mathrm{<span><span>a</span>a</span>}<span><span>+</span>+</span>\mathrm{<span><span>height</span>height</span>}<span><span>)</span>)</span><span><span>\&times;</span>\&times;</span>\mathrm{<span><span>sin</span>sin</span>}<span><span>(</span>(</span>\mathrm{<span><span>lon</span>the\; lon</span>}<span><span>)</span>)</span><span><span>\&times;</span>\&times;</span>\mathrm{<span><span>cos</span>cos</span>}<span><span>(</span>(</span>\mathrm{<span><span>lat</span>lat</span>}<span><span>)</span>)</span>\\ <span><span>(</span>(</span>\mathrm{<span><span>b</span>b</span>}<span><span>+</span>+</span>\mathrm{<span><span>height</span>height</span>}<span><span>)</span>)</span><span><span>\&times;</span>\&times;</span>\mathrm{<span><span>sin</span>sin</span>}<span><span>(</span>(</span>\mathrm{<span><span>lat</span>lat</span>}<span><span>)</span>)</span>\end{array}\right]$

（3.2）利用控制点处的三轴姿态角（Roll,Pitch，Yaw）计算由相机坐标系至轨道坐标系的旋转矩阵M1(3.2) Use the three-axis attitude angle (Roll, Pitch, Yaw) at the control point to calculate the rotation matrix M1 from the camera coordinate system to the orbit coordinate system

（3.3）利用控制点处的卫星轨道位置（PX,PY,PZ,VX,VY,VZ）计算卫星在J2000.0坐标系下的轨道六根数，由此计算轨道坐标系至J2000.0坐标系的旋转矩阵M2(3.3) Use the orbital position of the satellite at the control point (PX, PY, PZ, VX, VY, VZ) to calculate the six orbital numbers of the satellite in the J2000.0 coordinate system, and then calculate the orbital coordinate system to the J2000.0 coordinate system The rotation matrix M2

（3.4）利用控制点处的成像时间，计算J2000.0坐标系至WGS84坐标系的旋转矩阵M3(3.4) Using the imaging time at the control point, calculate the rotation matrix M3 from the J2000.0 coordinate system to the WGS84 coordinate system

（3.5）计算控制点处的理论指向矢量V_理论，计算公式为：(3.5) Calculate the theoretical pointing vector V_theory at the control point, the calculation formula is:

5.建立理论指向矢量与实际指向矢量之间的三轴旋转误差方程，对方程进行线性化处理5. Establish the three-axis rotation error equation between the theoretical pointing vector and the actual pointing vector, and linearize the equation

（1）设误差矩阵M为三轴旋转正交矩阵，建立理论指向矢量与实际指向矢量之间的矩阵转换模型；本发明中误差矩阵M表征为绕X,Y,Z坐标轴旋转的角度构成的正交矩阵，如下所示：(1) Let the error matrix M be a three-axis rotation orthogonal matrix, and establish a matrix transformation model between the theoretical pointing vector and the actual pointing vector; in the present invention, the error matrix M is represented by the angle of rotation around the X, Y, and Z coordinate axes Orthogonal matrix of , as follows:

其中：in:

$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]$为理论指向矢量$\left[\begin{array}{c}x\\ Y\\ Z\end{array}\right]$ is the theoretical pointing vector

$\left[\begin{array}{c}x\\ y\\ z\end{array}\right]$为实际指向矢量$\left[\begin{array}{c}x\\ \mathrm{the\; y}\\ z\end{array}\right]$ is the actual pointing vector

（2）将模型进行转换，得到理论指向矢量各参数与实际指向矢量各参数的非线性方程，如下所示：(2) Convert the model to obtain the nonlinear equations of the parameters of the theoretical pointing vector and the parameters of the actual pointing vector, as follows:

（3）对非线性方程在原点（0，0，0）处进行泰勒展开，得到非线性方程的线性近似方程。(3) Perform Taylor expansion on the nonlinear equation at the origin (0, 0, 0) to obtain the linear approximation equation of the nonlinear equation.

其中：in:

${\mathrm{<span><span>a</span>a</span>}}_{\mathrm{<span><span>1</span>1</span>}}<span><span>=</span>=</span><span><span>-</span>-</span>\frac{\mathrm{<span><span>XY</span>X\; Y</span>}}{{\mathrm{<span><span>Z</span>Z</span>}}^{\mathrm{<span><span>2</span>2</span>}}}$

${\mathrm{<span><span>a</span>a</span>}}_{\mathrm{<span><span>2</span>2</span>}}<span><span>=</span>=</span>\mathrm{<span><span>1</span>1</span>}<span><span>+</span>+</span>\frac{{\mathrm{<span><span>X</span>x</span>}}^{\mathrm{<span><span>2</span>2</span>}}}{{\mathrm{<span><span>Z</span>Z</span>}}^{\mathrm{<span><span>2</span>2</span>}}}$

${\mathrm{<span><span>a</span>a</span>}}_{\mathrm{<span><span>3</span>3</span>}}<span><span>=</span>=</span><span><span>-</span>-</span>\frac{\mathrm{<span><span>Y</span>Y</span>}}{\mathrm{<span><span>Z</span>Z</span>}}$

${\mathrm{<span><span>c</span>c</span>}}_{\mathrm{<span><span>1</span>1</span>}}<span><span>=</span>=</span>\mathrm{<span><span>x</span>x</span>}<span><span>-</span>-</span>\frac{\mathrm{<span><span>X</span>x</span>}}{\mathrm{<span><span>Z</span>Z</span>}}$

${\mathrm{<span><span>b</span>b</span>}}_{\mathrm{<span><span>1</span>1</span>}}<span><span>=</span>=</span><span><span>-</span>-</span><span><span>(</span>(</span>\mathrm{<span><span>1</span>1</span>}<span><span>+</span>+</span>\frac{{\mathrm{<span><span>Y</span>Y</span>}}^{\mathrm{<span><span>2</span>2</span>}}}{{\mathrm{<span><span>Z</span>Z</span>}}^{\mathrm{<span><span>2</span>2</span>}}}<span><span>)</span>)</span>$

${\mathrm{<span><span>b</span>b</span>}}_{\mathrm{<span><span>2</span>2</span>}}<span><span>=</span>=</span>\frac{\mathrm{<span><span>XY</span>X\; Y</span>}}{{\mathrm{<span><span>Z</span>Z</span>}}^{\mathrm{<span><span>2</span>2</span>}}}$

${\mathrm{<span><span>b</span>b</span>}}_{\mathrm{<span><span>3</span>3</span>}}<span><span>=</span>=</span>\frac{\mathrm{<span><span>X</span>x</span>}}{\mathrm{<span><span>Z</span>Z</span>}}$

${\mathrm{<span><span>c</span>c</span>}}_{\mathrm{<span><span>2</span>2</span>}}<span><span>=</span>=</span>\mathrm{<span><span>y</span>the\; y</span>}<span><span>-</span>-</span>\frac{\mathrm{<span><span>Y</span>Y</span>}}{\mathrm{<span><span>Z</span>Z</span>}}$

（4）针对每一组理论指向矢量和实际指向矢量，计算线性方程参数，得到线性方程组(4) For each set of theoretical pointing vectors and actual pointing vectors, calculate the parameters of the linear equations to obtain a set of linear equations

利用三轴旋转误差方程建立其理论指向矢量与实际指向矢量之间的线性方程组，针对N个控制点，可得到线性方程组如下：Using the three-axis rotation error equation to establish the linear equations between the theoretical pointing vector and the actual pointing vector, for N control points, the linear equations can be obtained as follows:

6.最小二乘法迭代解算误差方程参数，得到标定外参数6. The least square method iteratively solves the parameters of the error equation to obtain the calibration external parameters

（6.1）针对N个控制点建立的误差方程M1，最小二乘迭代解算方程参数(6.1) For the error equation M1 established for N control points, the least squares iterative solution equation parameters

对于误差方程演算出的大型稀疏线性方程组Ax=b，采用迭代法求解方程近似解。对于方程Ax=b，建立其递推公式：For the large sparse linear equations Ax=b calculated from the error equation, the approximate solution of the equations is solved by an iterative method. For the equation Ax=b, establish its recursive formula:

x_i+1=(A^TA+E)^-1A^Tb+(A^tA+E)^-1x_ix_i+1 =(A^T A+E)^-1 A^T b+(A^t A+E)^-1 x_i

建立最小二乘解误差方程：Set up the least squares solution to the error equation:

Δ=x_i+1-x_i=(A^TA+E)^-1A^Tb+(A^tA+E)^-1x_i-x_iΔ=x_i+1 -x_i =(A^T A+E)^-1 A^T b+(A^t A+E)^-1 x_i -x_i

正交误差矩阵M各参数都近似等于0，因此为x赋初始值（0,0,0）。Each parameter of the orthogonal error matrix M is approximately equal to 0, so the initial value (0,0,0) is assigned to x.

按照递推公式解算向量序列x₀,x₁,x₂,...,x_n，利用最小二乘误差方程解算Δ₁,Δ₂,...,Δ_n，直至Δ_k≤min，此时即为该方程的最小二乘解。Solve the vector sequence x₀ ,x₁ ,x₂ ,...,x_n according to the recursion formula, and use the least squares error equation to solve Δ₁ ,Δ₂ ,...,Δ_n until Δ_k ≤min ,at this time is the least squares solution of this equation.

（6.2）将方程M1的最小二乘解作为初步外方位参数标定结果，代入各控制点处的严格成像几何模型，计算各控制点处的检校残差，剔除残差大的控制点。(6.2) The least squares solution of equation M1 is used as the preliminary calibration result of the external orientation parameters, which is substituted into the strict imaging geometric model at each control point, the calibration residual at each control point is calculated, and the control points with large residuals are eliminated.

（6.3）利用剔除过大残差后的M个控制点，建立误差方程M2，最小二乘迭代解算方程参数。解算出的最小二乘解作为外方位参数标定结果。将标定结果应用于遥感图像几何校正算法中，对严格成像几何模型参数进行修正，可以大幅提高该相机遥感图像产品的无控几何定位精度。(6.3) Using the M control points after removing excessive residuals, establish the error equation M2, and solve the equation parameters iteratively by least squares. The least squares solution calculated by the solution is used as the calibration result of the outer orientation parameters. Applying the calibration results to the remote sensing image geometric correction algorithm to correct the parameters of the strict imaging geometric model can greatly improve the uncontrolled geometric positioning accuracy of the remote sensing image product of the camera.

如表1-表3所示，表1为本发明解算出的外方位参数标定结果，表2为未使用本发明前VRSS-1卫星几何定位精度测量结果，表3为利用本发明进行在轨检校后VRSS-1卫星的几何定位精度测量结果。从表1可以看出VRSS-1卫星外方位参数在滚动和俯仰两个方面均发生了比较大的改变，达到了0.06度量级，根据VRSS-1卫星的平均轨道高度（官方标称值为640公里）测算，滚动方向-0.06239度的外参数可导致地面无控定位精度约为640000×tan(-0.06239)=696.1米；俯仰方向0.06804度外参数可导致的地面无控定位精度约为640000×tan(0.06804)=760.07米；滚动俯仰方向可导致的综合误差为

米。从表2可以看出在采用本发明前卫星的实测无控定位精度为1073.18米（平均值）量级，可见外方位参数导致的无控定位误差与实测值基本相符。从表3可以看出，采用本发明对卫星进行外方位参数标定后卫星的实测无控定位精度为36.31米（平均值）量级，本发明将该卫星的无控定位精度提高了29倍。As shown in Table 1-Table 3, Table 1 is the calibration result of the outer azimuth parameters calculated by the present invention, Table 2 is the measurement result of geometric positioning accuracy of the VRSS-1 satellite before the present invention is not used, and Table 3 is the on-orbit using the present invention Measurement results of geometric positioning accuracy of VRSS-1 satellite after calibration. It can be seen from Table 1 that the outer azimuth parameters of the VRSS-1 satellite have undergone relatively large changes in both roll and pitch, reaching a magnitude of 0.06. According to the average orbital altitude of the VRSS-1 satellite (the official nominal value is 640 km) calculation, the external parameter of -0.06239 degrees in the rolling direction can lead to an uncontrolled positioning accuracy of about 640000×tan(-0.06239)=696.1 meters; the external parameter of 0.06804 degrees in the pitch direction can result in an uncontrolled positioning accuracy of about 640000× tan(0.06804)=760.07 meters; the comprehensive error caused by the rolling and pitching direction is

rice. It can be seen from Table 2 that the measured uncontrolled positioning accuracy of the satellite before the present invention is adopted is on the order of 1073.18 meters (average value), and it can be seen that the uncontrolled positioning error caused by the outer azimuth parameters is basically consistent with the measured value. It can be seen from Table 3 that the actual measured uncontrolled positioning accuracy of the satellite after the outer azimuth parameter calibration of the satellite is on the order of 36.31 meters (average value), and the invention improves the uncontrolled positioning accuracy of the satellite by 29 times.

表1Table 1

误差矩阵参数Error Matrix Parameters检校结果（单位：度）Calibration result (unit: degree)Rollroll-0.062391107351712-0.062391107351712Pitchpitch0.0680455338945130.068045533894513YawYaw0.0008432357594850.000843235759485

表2Table 2

表3table 3

本发明未详细阐述的部分属于本领域公知技术。The parts not described in detail in the present invention belong to the well-known technology in the art.

Claims (4)

1. remote sensing satellite linear array push is swept an optical camera outer orientation parameter calibration method in-orbit, it is characterized in that comprising the following steps:

1) treating calibration image I₁, control image I₂on utilize reference mark Auto-matching algorithm, gather reference mark information, described reference mark information is that T registering control points is to { PixI_1i, PixI_2i, i=1 wherein, 2,3 ..., T); Record each dominating pair of vertices { PixI_1i, PixI_2itreat calibration image I₁picture point coordinate (i, j) and the corresponding position coordinates (lat, lon, height) of image feature point under WGS84 coordinate system of controlling;

2) obtain camera internal position element and the auxiliary data at each place, reference mark, when described auxiliary data comprises the track of satellite, attitude, row; Set up the strict imaging model at each place, reference mark[X wherein_s, Y_s, Z_s] for row corresponding to each reference mark, inscribe constantly the position of satellite in agreement geocentric coordinate system, orbital position (PX, PY, PZ) is located at reference mark; [X_g, Y_g, Z_g] be the coordinate of each reference mark corresponding terrain object point of corresponding pixel in agreement geocentric coordinate system; PsiX, psiY are respectively the angle of the pixel primary optical axis unit vector of the corresponding image in each reference mark and satellite body coordinate system X-axis, Y-axis; U is scale factor; M₀for the installation matrix of satellite body coordinate system with respect to camera, before satellite launch, by ground survey, obtained; M₁for row corresponding to each reference mark, inscribe constantly satellite to orbital coordinate system rotation matrix, by the attitude angle of measuring on star, formed; M₂for row corresponding to each reference mark inscribed track constantly to J2000.0 coordinate system rotation matrix; M₃during for this, inscribe J2000.0 to WGS84 coordinate system rotation matrix;

3) according to step 2) the strict imaging model that obtains, calculate the positioning error size and the direction that obtain each place, reference mark, preliminary excluding gross error reference mark, obtains initial control point pair;

4) set up how much calibration mathematical models of camera elements of exterior orientation, the initial control point pair obtaining according to step 3), calculates and obtains each reference mark processing opinion pointing vector and actual pointing vector;

41) in strict imaging model, introduce error matrix M, set up how much calibration mathematical models of camera elements of exterior orientation

42) to arbitrary group of reference mark, utilize place, elements of interior orientation (psiX, psiY) the calculating reference mark optical axis at place, reference mark to point to vector V in camera coordinates system, be normalized the rear actual pointing vector V as this place, reference mark_actual;

43) to arbitrary group of reference mark, utilize reference mark ground point location coordinate (lat, lon, height) to calculate place, reference mark optical axis and point to vector V in camera coordinates system, be normalized the rear theoretical pointing vector V as this place, reference mark_theoretical;

Whereina is semimajor axis of ellipsoid, and b is semiminor axis of ellipsoid;

5) the three axle rotation error equations that theorize between pointing vector and actual pointing vector, carry out linearization process to this three axles rotation error equation, obtain system of linear equations;

51) establishing error matrix M is three axle rotating orthogonal matrixes, the matrix conversion model theorizing between pointing vector and actual pointing vector:

Wherein

for theoretical pointing vector,for actual pointing vector;

52) model is changed, obtained each nonlinearity in parameters equation of each parameter of theoretical pointing vector and actual pointing vector;

53) nonlinear equation is located to carry out Taylor expansion at initial point (0,0,0), obtain the linear-apporximation equation of nonlinear equation;

54) for each, organize theoretical pointing vector and actual pointing vector, calculate linear equation parameter, obtain system of linear equations;

6) system of linear equations that adopts least square method Iterative step 5) to obtain, obtains and demarcates outer parameter.

2. a kind of remote sensing satellite linear array push according to claim 1 is swept optical camera outer orientation parameter calibration method in-orbit, it is characterized in that: the concrete grammar that gathers reference mark information in step 1) is:

11) based on SIFT algorithm, treat calibration original image I₁carry out feature point extraction, obtain M unique point PixI_1i(i=1,2,3 ..., M), M is positive integer; Record the SIFT proper vector of each unique point;

12) based on SIFT algorithm, high precision is controlled to image I₂carry out feature point extraction, obtain N unique point PixI_2i(i=1,2,3 ..., N), N is positive integer; Record the SIFT proper vector at each unique point place;

13) adopt Euclidean distance as similarity measurement criterion, the unique point of two width images to be mated, obtain T registering control points to { PixI_1i, PixI_2i(i=1,2,3 ..., T).

3. a kind of remote sensing satellite linear array push according to claim 1 is swept optical camera outer orientation parameter calibration method in-orbit, it is characterized in that: step 2) in set up the strict imaging model at place, each reference mark concrete grammar be:

21) according to dominating pair of vertices { PixI_1i, PixI_2iimage coordinate (i, j), while obtaining its corresponding row;

22) utilize Lagrange's interpolation algorithm interpolation calculation dominating pair of vertices { PixI_1i, PixI_2iorbital position (PX, PY, PZ, VX, VY, VZ) during corresponding row;

23) attitude quaternion satellite being passed down carries out coordinate system conversion process, utilizes Lagrangian Arithmetic interpolation calculation dominating pair of vertices { PixI_1i, PixI_2iconstantly locate camera with respect to the three-axis attitude angle (Roll, Pitch, Yaw) of orbital coordinate system during corresponding row;

24) according to camera laboratory measurement data or how much calibration results of elements of interior orientation, calculate dominating pair of vertices { PixI_1i, PixI_2icorresponding CCD visits first optical axis and points to angle (psiX, psiY);

25) according to step 21)-step 24) result that obtains, set up the strict imaging model at place, reference mark, certain on the image obtaining for a certain moment of linear array push-scanning image camera a bit, builds the strict imaging model of remote sensing image.

4. a kind of remote sensing satellite linear array push according to claim 1 is swept optical camera outer orientation parameter calibration method in-orbit, it is characterized in that: the concrete acquisition methods that in step 3), initial control point is right is:

31) according to step 2) the strict imaging model that obtains, calculates each dominating pair of vertices { PixI_1i, PixI_2iin treat calibration imaging point PixI_1iinitially millet cake position coordinates { lon ', lat ' };

32) according to each dominating pair of vertices { PixI_1i, PixI_2ithe middle imaging point PixI that controls_2ilocation coordinate information, calculates reference mark PixI_1iinitial alignment error { Δ xi, Δ yi};

33) calculate dominating pair of vertices { PixI_1i, PixI_2ipositioning error mean value { Δ x, Δ y};

34) to each dominating pair of vertices, judge its initial alignment error Δ xi, and Δ yi} whether Δ x ± 20% Δ x, within the scope of Δ y ± 20% Δ y}, will be not Δ x ± 20% Δ x, the dominating pair of vertices within the scope of Δ y ± 20% Δ y} is rejected;

35) reject after error dominating pair of vertices, obtain K group initial control point to { PixI_1i, PixI_2i(i=1,2,3 ..., K).

Images