阅读说明：本技术 利用rpc参数建立高分辨率卫星影像等效几何成像模型的方法 (The method for establishing the equivalent geometry imaging model of high-resolution satellite image using RPC parameter ) 是由 曹辉 李海鸿 陶鹏杰 于 2019-07-12 设计创作，主要内容包括：本发明涉及一种利用RPC参数建立高分辨率卫星影像等效几何成像模型的方法。首先选择高分辨率卫星影像中任意扫描行列出投影中心、像点及其对应物方点的共线条件方程,利用高分辨率卫星运行特点和几何特征,抽象出等效主距、主点坐标,偏置矩阵中的三个旋转角以及六个轨道根数,构成了等效几何成像模型的11个独立参数；然后通过利用RPC参数,直接计算等效几何成像模型的参数。本发明解决了RPC存在过度参数化、参数间强相关以及拟合严密成像模型具有残余误差等内在问题,其参数具有明显的几何和物理意义,方便对各种成像误差源进行建模分析,本发明解决了RPC参数相互转换,有利于提高高分辨率卫星影像的定位精度。(The present invention relates to a kind of methods for establishing the equivalent geometry imaging model of high-resolution satellite image using RPC parameter.Arbitrary scan ranks in high-resolution satellite image are selected to go out projection centre, picture point and its collinearity condition equation for corresponding to object space point first, utilize high-resolution satellite operation characteristic and geometrical characteristic, equivalent master is taken out away from, principal point coordinate, three rotation angles and six orbital trackings in bias matrix, constitute 11 independent parameters of equivalent geometry imaging model；Then by utilizing RPC parameter, the parameter of equivalent geometry imaging model is directly calculated.The present invention solves RPC, and there are strong correlations between overparameterization, parameter and the tight imaging model of fitting to have the built in problems such as residual error, its parameter has apparent geometry and physical significance, it is convenient that modeling analysis is carried out to various image error sources, the present invention solves RPC parameter and mutually converts, and is conducive to the positioning accuracy for improving high-resolution satellite image.)

1. the method for establishing the equivalent geometry imaging model of high-resolution satellite image using RPC parameter, which is characterized in that including such as Lower step:

Step 1, linear array push sweeps the building of the equivalent geometry imaging model of satellite image, and specific implementation is as follows,

Select the earth's core rectangular coordinate system as object coordinates system, projection centre is as satellite body coordinate origin, in object space Position vector under coordinate system in the plane that x-axis is located at satellite flight motion vector and z-axis is constituted and is directed toward movement as z-axis Direction vector, y-axis are determined perpendicular to xz plane, direction by right-hand rule；The imaging of any scan line correspondence in image Geometry can be expressed as the collinearity condition equation of instantaneous projection centre, image picture point and its corresponding object space point three, equivalent geometry Imaging model general type is expressed as following equation:

Wherein [l, s] is picpointed coordinate, and l is row, and s is column, y₀It is principal point pixel coordinate, f is equivalent master away from Xs (l), Ys (l), Zs (l) is the corresponding projection centre instantaneous position of scan line l；R_sat(l) be satellite platform transient posture matrix；R_sensor It (l) is bias matrix of the sensor relative to satellite platform；Three rotations of the equivalent master away from, principal point coordinate, in bias matrix Angle and six orbital trackings, constitute the independent parameter of equivalent geometry imaging model, and by RPC parameter calculation, these parameters are built Found corresponding equivalent geometry imaging model；

Step 2, the equivalent geometry imaging model parameter calculation based on RPC parameter, specifically includes following sub-step:

Step 2.1, the corresponding object coordinates of picture point are calculated based on RPC parameter；

Step 2.2, imaging sensor geometrical characteristic is swept using high-resolution satellite linear array push, calculates equivalent master away from sitting with principal point Mark；

Step 2.3, projection centre is intersected at using the corresponding imaging ray of picture points all on same one scan line, calculates arbitrary scan The corresponding instantaneous projection centre of row；

Step 2.4, using chebyshev approximating polynomial satellite motion track；

Step 2.5, satellite platform coordinate system is calculated relative to object space according to the corresponding unit vector of three reference axis of satellite platform The transient posture matrix of coordinate system；

Step 2.6, using four element Iterative sensor bias matrix of unit.

2. the method for establishing the equivalent geometry imaging model of high-resolution satellite image using RPC parameter as described in claim 1, It is characterized by: the specific implementation of step 2.1 is as follows,

RFM general expression is as follows,

Wherein, (e, n, h) is the geodetic coordinates of ground point, and Num (e, n, h) and Den (e, n, h) are cubic polynomial, and (L, S) is Normalized picpointed coordinate；

Take up an official post the pixel coordinate (l, s) of one picture point of meaning to fixing, normalized coordinate is calculated as follows:

Wherein line_offset, line_scale, sample_offset, sample_scale are rational polynominal model normalizing Change parameter；

The corresponding ground elevation h of the picture point is further given, then the picture point is obtained by equation (2) iterative calculation and corresponds to object space point Longitude and latitude (e, n), using WGS84 reference ellipsoid parameter, major semiaxis a and eccentric ratio e s, as the following formula by geodetic coordinates (e, n, h) Be converted to geocentric rectangular coordinate system coordinate:

3. the method for establishing the equivalent geometry imaging model of high-resolution satellite image using RPC parameter as claimed in claim 2, It is characterized by: the specific implementation of step 2.2 is as follows,

According to equivalent geometrical model it is assumed that ideal line array CCD Pixel size is consistent, it is flat that camera coke is placed in an ideal line On face, conformation corresponds to a scan line of satellite image, and first pixel of the scan line is selected to make on any scan line For picture point c (l, 0), any pixel is as picture point a (l, s_a)；According to geometrical relationship, to picture point c (l, 0) and projection centre line, Picture point a (l, s_a) and projection centre line between the angle α that constitutes, be listed below equation,

Tan (α)=tan (α₁+α₂)=(tan (α₁)+tan(α₂))/(1-tan(α₁)tan(α₂)) (5)

Wherein α₁For picture point c (l, 0) and projection centre line, principal point o (l, y₀) and projection centre line between the folder that constitutes Angle, α₂For picture point a (l, s_a) and projection centre line, principal point o (l, y₀) and projection centre line between the angle that constitutes, side After journey further arranges, obtain:

Tan (α)=f*s_a/[f²-y₀*(s_a-y₀)] (6)

The highest elevation h of RPC parameter expression is given respectively₁With minimum elevation h₂, according to the method for step 2.1, picture point is calculated The object coordinates of the corresponding object space point A1 and A2 of a；

To picture point c (l, 0) and picture point b (l, s_b) same equation is listed, quadratic term is eliminated using two equations and calculates equivalent master Away from f and principal point coordinate y₀Initial value；And by selecting multi-strip scanning line on a scape image, selected in every scan line multiple Picture point by formula (6) column error equation and linearizes, and thus forms normal equation and iterates to calculate f and y₀Optimal solution.

4. the method for establishing the equivalent geometry imaging model of high-resolution satellite image using RPC parameter as claimed in claim 3, It is characterized by: the specific implementation of step 2.3 is as follows,

In equivalent geometry imaging model, projected centered on the corresponding imaging geometry of the every one scan line of satellite image, so same The corresponding imaging ray of all picture points intersects at projection centre on scan line, and picture point a (l, s are selected on any scan line_a), point The highest elevation h of the point RPC parameter expression is not given₁With minimum elevation h₂, according to the method for step 2.1, picture point a is calculated The object coordinates of corresponding object space point A1 and A2, projection centre P (X_s, Y_s, Z_s) it is located at straight line A₁_A₂On, it thus lists with lower section Journey (7), similarly, to picture point b (l, s_b) also calculate object space point B₁And B₂And corresponding equation (8) are listed,

According to the principle of least square, A is calculated by equation group (7) and (8)₁_A₂And B₁_B₂That is, the intersection point of two space lines should The corresponding projection centre of scan line；Stablize to calculate, multiple points is selected to calculate the scanning with least square adjustment on scan line The projection centre coordinate Xs (l), Ys (l), Zs (l) at row corresponding moment.

5. the method for establishing the equivalent geometry imaging model of high-resolution satellite image using RPC parameter as claimed in claim 4, It is characterized by: the specific implementation of step 2.4 is as follows,

The instantaneous position of satellite flight movement is given a definition in geocentric rectangular coordinate system, and orbit equation passes through Chebyshev polynomials Fitting obtains, it is assumed that needs to be fitted image scan row [l₀, l_m] in satellite orbit, using formula (9) first by variable l ∈ [l₀, l_m] It is transformed into τ ∈ [- 1,1],

τ=2 (l-l₀)/(l_m-l₀) -1, l ∈ [l₀, l_m] (9)

Then satellite orbit coordinate, i.e. projection centre coordinate Xs (l), the Chebyshev polynomials expression of Ys (l), Zs (l) are as follows:

In formula, n is the order of Chebyshev polynomials,Respectively X, Y, the Chebyshev of Z coordinate component are multinomial Formula coefficient, T_i(τ) is determined according to recurrence formula (11),

T₀(τ)=1, T₁(τ)=τ, T_n+1(τ)=2 τ T_n(τ)-T_n-1(τ) (11)

Then the projection centre coordinate being calculated using step 2.3 is calculated satellite orbit by least square adjustment and is being imaged Chebyshev polynomials coefficient in period

6. the method for establishing the equivalent geometry imaging model of high-resolution satellite image using RPC parameter as claimed in claim 5, It is characterized by: the specific implementation of step 2.5 is as follows,

The given corresponding instantaneous projection centre position vector of arbitrary scan lineIt is calculated pair using formula (10) The velocity vector answeredThe corresponding unit vector of three reference axis of satellite platform is calculated as follows:

Thus the transient posture matrix for obtaining satellite platform coordinate system relative to object coordinates system is as follows:

7. the method for establishing the equivalent geometry imaging model of high-resolution satellite image using RPC parameter as claimed in claim 6, It is characterized by: the specific implementation of step 2.6 is as follows,

When known to projection centre instantaneous position and satellite platform transient posture matrix, equivalent geometry imaging model simplifies are as follows:

Wherein, the virtual object space point coordinate and projection centre instantaneous position that [u (l), v (l), w (l)] is calculated according to RPC parameter With corresponding transient posture matrix R_sat(l) it is calculated, using picpointed coordinate as observation, above formula is directed to spin matrix R_sensorParameter linearisation, and by the principle of least square form normal equation, Iterative sensor bias matrix element；

Using unit quaternion [q₀, q₁, q₂, q₃] expression spin matrix R_sensor, unit quaternion satisfactionIts corresponding spin matrix is calculated by following formula,

According to quaternion differential equation, attitudes vibration expressed by four elements is expressed as coordinate system and turns around three of them reference axis Dynamic, the update of quaternary number can be calculated with following formula；

Selection Chebyshev's orthogonal polynomial first expresses attitudes vibration item w_x, w_y, w_z, by spin matrix (15) to unit quaternary Number derivation is simultaneously substituted into (16) formula, is obtained linearized expression of the spin matrix relative to three independent parameters, is then substituted into (14) formula and proportionality coefficient λ is eliminated, obtains the error equation for solving sensor offset angle；Then every scape image multi-strip scanning is selected Line, the multiple picture points of every scan line arrange vertical error equation by formula (14) and linearize, thus form normal equation and iterate to calculate w_x, w_y, w_zCorresponding Chebyshev's orthogonal polynomial coefficient optimal solution, and the corresponding rotation of arbitrary scan line is calculated according to above formula Torque battle array R_sensor(l)。

8. the method for carrying out block adjustment using equivalent geometry imaging model as claimed in claim 7, it is characterised in that: specific Implementation is as follows,

Instantaneous position and satellite platform transient posture are moved by disturbing influence in view of satellite orbits the earth, and with additional parameter Form expresses picpointed coordinate systematic error, based on the general of equivalent geometry imaging model satellite image block adjustment function model Form is as follows:

Wherein, Δ x, Δ y are the systematic error indicated in the form of additional parameter and model error, [dXs (l), dYs (l), dZs It (l)] is disturbance part of the satellite position with sporting flying, D_oriIt (l) is influence of the various disturbances to the attitude of satellite, this three parts Constitute the unknown number of self-calibration block adjustment system；

After eliminating λ, the corresponding collinearity equation of each picture point (17) formula is expressed as two error equations of following form:

Wherein (v_x, v_y) it is that picture point corrects residual error, (x, y, f) is picture point observation and equivalent master away from (u, v, w) is according to object space The auxiliary coordinates coordinate that coordinate initial value, projection centre coordinate and spin matrix calculate, p are power；

U in above formula, v, w are linearized, the error equation of the abstract geometry imaging model of matrix form is obtained after arrangement:

V₁=B₁₁X₁+B₁₂X₂+l₁, P₁ (19)

Wherein V₁For the expression matrix of picture point residual error, B₁₁, B₁₂Error equation factor arrays, l₁For constant term matrix, P₁For weight matrix, X₁And X₂For unknown matrix number.

Technical field

The present invention is under the jurisdiction of remote sensing science and technology field, and the geometry location for being adapted to high-resolution satellite image calculates, main Syllabus is that there are reductions of geometric positioning accuracy caused by the factors such as systematic error etc. to ask in solution rational polynominal model (RFM) Topic.The method for proposing to establish the equivalent geometry imaging model of high-resolution satellite image using rational polynominal model parameter (RPC), And by equivalent geometry imaging model, the parameters such as Compensation for Model Errors are added, eliminate the geometry due to caused by model error in RFM The problems such as positioning accuracy reduces, to improve the geometric positioning accuracy of high-resolution satellite image.

Background technique

The geometry imaging model of high-resolution satellite image is generally divided into two classes: tight imaging model and universal imaging mould Type^[1][2].All parameters of tight imaging model have clear geometric meaning, therefore are highest based on its positioning accuracy^[3]。 But establish tight imaging model and need using satellite transit orbital position, posture observational record and sensor structure parameter, and Shortage versatility, therefore the less data publication form as satellite image.As the mathematical simulation of tight imaging model, it is based on three RFM (Rational Function Model) the universal imaging model of secondary rational polynominal has sensor parameters confidentiality, leads to With the advantages that property is strong, form is simple and convenience of calculation is quick^[4], it is widely studied and applied^[5-8].Currently, RFM model RPC (Rational Polynomial Coefficient) parameter have become the standard configuration of high-resolution satellite image.

Satellite image with RPC parameter is usually pre-processed by geometry, including in-orbit geometry calibration, multi-disc CCD Linear array splicing, sensing system error correction etc. have certain direct georeferencing precision, can be used for geometry location, just penetrating The subsequent processings such as image rectification, stereoplotting.However, RFM still has certain defect in practical applications, mainly it is reflected in Two aspect of RPC built in problem and external action.Firstly, RPC there are strong correlation between overparameterization, parameter and fitting tightly at As model has the built in problems such as residual error；Secondly, the stability of many satellite platforms is not by technical restrictions such as sensors Height, and transient posture measurement accuracy is limited, and these external factor lead to the presence of apparent system in satellite image and its RPC parameter Error caused by the error, such as satellite platform high frequency flutter of uniting.These problems will result directly in satellite image in stereoplotting and Vertical parallax in the application such as Image Matching is difficult to eliminate, and stereopsis is difficult, positioning accuracy declines.

Summary of the invention

The purpose of the invention is to effectively eliminate obvious systematic error present in satellite image and its RPC parameter, It establishes a set of independent of sensor structure information and appearance rail measurement data, is suitable for a variety of high-resolution satellite images, it can be right Various image error sources carry out modeling analysis, and the imaging model that can mutually convert with RPC parameter, pass through self calibration region Net adjusted data compensation model error, and RPC parameter is advanced optimized, it corrects in image and corresponding RPC parameter and there is apparent system System error.

In order to achieve the above object, present invention provide the technical scheme that establishing high-resolution satellite shadow using RPC parameter As the method for equivalent geometry imaging model, include the following steps:

Step 1, linear array push sweeps the building of the equivalent geometry imaging model of satellite image, and specific implementation is as follows,

Select the earth's core rectangular coordinate system as object coordinates system, projection centre as satellite body coordinate origin, As z-axis, x-axis is located in the plane that satellite flight motion vector is constituted with z-axis and is directed toward position vector under object coordinates system Motion vector direction, y-axis are determined perpendicular to xz plane, direction by right-hand rule；Any scan line corresponds to it in image Imaging geometry can be expressed as the collinearity condition equation of instantaneous projection centre, image picture point and its corresponding object space point three, equivalent Geometry imaging model general type is expressed as following equation:

Wherein [l, s] is picpointed coordinate, and l is row, and s is column, y₀It is principal point pixel coordinate, f is equivalent master away from, Xs (l), Ys (l), Zs (l) are the corresponding projection centre instantaneous positions of scan line l；R_sat(l) be satellite platform transient posture matrix； R_sensorIt (l) is bias matrix of the sensor relative to satellite platform；Equivalent master is away from, principal point coordinate, three in bias matrix Rotation angle and six orbital trackings, constitute the independent parameter of equivalent geometry imaging model, pass through these ginsengs of RPC parameter calculation Number establishes corresponding equivalent geometry imaging model；

Step 2, the equivalent geometry imaging model parameter calculation based on RPC parameter, specifically includes following sub-step:

Step 2.1, the corresponding object coordinates of picture point are calculated based on RPC parameter；

Step 2.2, using high-resolution linear array push-scanning image sensor geometrical characteristic, equivalent master is calculated away from sitting with principal point Mark；

Step 2.3, projection centre is intersected at using the corresponding imaging ray of picture points all on same one scan line, calculated any The corresponding instantaneous projection centre of scan line；

Step 2.4, using chebyshev approximating polynomial satellite motion track；

Step 2.5, according to the corresponding unit vector of three reference axis of satellite platform calculate satellite platform coordinate system relative to The transient posture matrix of object coordinates system；

Step 2.6, using four element Iterative sensor bias matrix of unit.

Further, the specific implementation of step 2.1 is as follows,

RFM general expression is as follows,

Wherein, (e, n, h) is the geodetic coordinates of ground point, and Num (e, n, h) and Den (e, n, h) are cubic polynomial, (L, It S) is normalized picpointed coordinate；

Take up an official post the pixel coordinate (l, s) of one picture point of meaning to fixing, normalized coordinate is calculated as follows:

Wherein line_offset, line_scale, sample_offset, sample_scale are rational polynominal model Normalized parameter；

The corresponding ground elevation h of the picture point is further given, then the picture point counterpart is obtained by equation (2) iterative calculation The longitude and latitude (e, n) just put, using WGS84 reference ellipsoid parameter, major semiaxis a and eccentric ratio e s press geodetic coordinates (e, n, h) Following formula is converted to geocentric rectangular coordinate system coordinate:

Further, the specific implementation of step 2.2 is as follows,

According to equivalent geometrical model it is assumed that ideal line array CCD Pixel size is consistent, camera is placed in an ideal line On focal plane, conformation corresponds to a scan line of satellite image, and first picture of the scan line is selected on any scan line Element is used as picture point c (l, 0), and any pixel is as picture point a (l, s_a)；According to geometrical relationship, to picture point c (l, 0) and projection centre Line, picture point a (l, s_a) and projection centre line between the angle α that constitutes, be listed below equation,

Tan (α)=tan (α₁+α₂)=(tan (α₁)+tan(α₂))/(1-tan(α₁)tan(α₂)) (5)

Wherein α₁For picture point c (l, 0) and projection centre line, principal point o (l, y₀) and projection centre line between constitute Angle, α₂For picture point a (l, s_a) and projection centre line, principal point o (l, y₀) and projection centre line between the angle that constitutes, It is available after equation further arranges:

Tan (α)=f*s_a/[f²-y₀*(s_a-y₀)] (6)

The highest elevation h of RPC parameter expression is given respectively₁With minimum elevation h₂, according to the method for step 2.1, it is calculated The object coordinates of the corresponding object space point A1 and A2 of picture point a；

To picture point c (l, 0) and picture point b (l, s_b) same equation is listed, quadratic term is eliminated using two equations and is calculated Equivalent master is away from f and principal point coordinate y₀Initial value；And by selecting multi-strip scanning line on a scape image, selected in every scan line Multiple picture points are selected, by formula (6) column error equation and are linearized, normal equation is thus formed and iterates to calculate f and y₀Optimal solution.

Further, the specific implementation of step 2.3 is as follows,

In equivalent geometry imaging model, projected centered on the corresponding imaging geometry of the every one scan line of satellite image, so With picture points all on one scan line, corresponding imaging ray intersects at projection centre, on any scan line select picture point a (l, s_a), the highest elevation h of the point RPC parameter expression is given respectively₁With minimum elevation h₂, according to the method for step 2.1, it is calculated The object coordinates of the corresponding object space point A1 and A2 of picture point a, projection centre P (X_s, Y_s, Z_s) it is located at straight line A₁_A₂On, thus list with Lower equation (7), similarly, to picture point b (l, s_b) also calculate object space point B₁And B₂And corresponding equation (8) are listed,

According to the principle of least square, A is calculated by equation group (7) and (8)₁_A₂And B₁_B₂The intersection point of two space lines, That is the corresponding projection centre of the scan line；Stablize to calculate, selecting multiple points to calculate on scan line with least square adjustment should Scan line corresponds to the projection centre coordinate Xs (l) at moment, Ys (l), Zs (l).

Further, the specific implementation of step 2.4 is as follows,

The instantaneous position of satellite flight movement is given a definition in geocentric rectangular coordinate system, and orbit equation is more by Chebyshev Item formula is fitted to obtain, it is assumed that needs to be fitted image scan row [l₀, l_m] in satellite orbit, using formula (9) first by variable l ∈ [l₀, l_m] it is transformed into τ ∈ [- 1,1],

τ=2 (l-l₀)/(l_m-l₀) -1, l ∈ [l₀, l_m] (9)

Then satellite orbit coordinate, i.e. projection centre coordinate Xs (l), the Chebyshev polynomials expression of Ys (l), Zs (l) are as follows:

In formula, n is the order of Chebyshev polynomials,The Qie Bixue of respectively X, Y, Z coordinate component Husband's multinomial coefficient, T_i(τ) is determined according to recurrence formula (11),

T₀(τ)=1, T₁(τ)=τ, T_n+1(τ)=2 τ T_n(τ)-T_n-1(τ) (11)

Then the projection centre coordinate being calculated using step 2.3 is calculated satellite orbit by least square adjustment and existed Chebyshev polynomials coefficient in imaging time section

Further, the specific implementation of step 2.5 is as follows,

The given corresponding instantaneous projection centre position vector of arbitrary scan lineIt is calculated using formula (10) To corresponding velocity vectorThe corresponding unit of three reference axis of satellite platform is calculated as follows Vector:

Thus the transient posture matrix for obtaining satellite platform coordinate system relative to object coordinates system is as follows:

Further, the specific implementation of step 2.6 is as follows,

When known to projection centre instantaneous position and satellite platform transient posture matrix, equivalent geometry imaging model simplifies Are as follows:

Wherein, the virtual object space point coordinate and projection centre that [u (l), v (l), w (l)] is calculated according to RPC parameter are instantaneous Position and corresponding transient posture matrix R_sat(l) it is calculated, using picpointed coordinate as observation, above formula is directed to spin moment Battle array R_sensorParameter linearisation, and by the principle of least square form normal equation, Iterative sensor bias matrix element；

Using unit quaternion [q₀, q₁, q₂, q₃] expression spin matrix R_sensor, unit quaternion satisfactionIts corresponding spin matrix is calculated by following formula,

According to quaternion differential equation, attitudes vibration expressed by four elements is expressed as coordinate system around three of them reference axis Rotation, the update of quaternary number can be calculated with following formula；

Selection Chebyshev's orthogonal polynomial first expresses attitudes vibration item w_x, w_y, w_z, by spin matrix (15) to unit The derivation of quaternary number is simultaneously substituted into (16) formula, obtains linearized expression of the spin matrix relative to three independent parameters, then generation Enter (14) formula and eliminate proportionality coefficient λ, obtains the error equation for solving sensor offset angle；Then select every scape image is a plurality of to sweep Line is retouched, the multiple picture points of every scan line arrange vertical error equation by formula (14) and linearize, thus form normal equation iterative calculation w_x, w_y, w_zCorresponding Chebyshev's orthogonal polynomial coefficient optimal solution, and arbitrary scan line is calculated according to above formula and corresponds to Spin matrix R_sensor(l)。

Further, block adjustment is carried out using the equivalent geometry imaging model of the method for the present invention building, it is specific real Existing mode is as follows,

Instantaneous position and satellite platform transient posture are moved by disturbing influence in view of satellite orbits the earth, and with additional Parametric form expresses picpointed coordinate systematic error, based on equivalent geometry imaging model satellite image block adjustment function model General type is as follows:

Wherein, Δ x, Δ y are the systematic error indicated in the form of additional parameter and model error, [dXs (l), dYs (l), DZs (l)] it is disturbance part of the satellite position with sporting flying, D_oriIt (l) is various influences of the disturbance to the attitude of satellite, this three Divide and constitutes the unknown number of self-calibration block adjustment system；

After eliminating λ, the corresponding collinearity equation of each picture point (17) formula is expressed as two error equations of following form:

Wherein (v_x, v_y) it is that picture point corrects residual error, (x, y, f) is picture point observation and equivalent master away from (u, v, w) is basis The auxiliary coordinates coordinate that object coordinates initial value, projection centre coordinate and spin matrix calculate, p are power；

U in above formula, v, w are linearized, the error equation of the abstract geometry imaging model of matrix form is obtained after arrangement:

V₁=B₁₁X₁+B₁₂X₂+l₁, P₁ (19)

Wherein V₁For the expression matrix of picture point residual error, B₁₁, B₁₂Error equation factor arrays, l₁For constant term matrix, P₁For power Matrix, X₁And X₂For unknown matrix number.

Compared with prior art, the advantages of the present invention:

The present invention has independent of sensor structure information and appearance rail measurement data, fits compared with rigorous geometry model For a variety of high-resolution satellite images, modeling analysis can be carried out to various image error sources, and can be mutual with RPC parameter Conversion advanced optimizes RPC parameter by self-calibration block adjustment compensation model error, corrects image and corresponding RPC ginseng There are apparent systematic errors in number.And compared with rational polynominal model, the present invention solve RPC there are overparameterization, Strong correlation and the tight imaging model of fitting have a built in problems such as residual error between parameter, parameter have apparent geometry and Physical significance, it is convenient that modeling analysis is carried out to various image error sources, and can directly optimize during block adjustment.This Invention solves RPC parameter and mutually converts, and is conducive to the positioning accuracy for improving high-resolution satellite image.

Detailed description of the invention

Fig. 1 is that scan line corresponds to projection centre resolving schematic diagram.

Based on Fig. 2 away from principal point resolve relation schematic diagram.

Specific embodiment

Technical solution of the present invention is described further with reference to the accompanying drawings and examples.

The present invention is to calculate equivalent geometry imaging model based on RPC parameter, and process is as follows:

1, linear array push sweeps the equivalent geometry imaging model of satellite image

Equivalent geometry imaging model (Equivalent Geometric Sensor Model, EGSM) is high-resolution light The mathematical abstractions that satellite image linear array push sweeps imaging geometry are learned, the corresponding target point object coordinates of image picpointed coordinate are described Between geometrical relationship.The foundation of equivalent geometry imaging model is based on following brass tacks (assumptions):

The image of every one scan line correspondence is obtained in a manner of central projection in satellite image；

Ideal line array CCD Pixel size is consistent, is placed in camera focal plane with an ideal line, and conformation is corresponding A scan line of satellite image；

The earth meets around the axis of rotation at the uniform velocity rotation, the sporting flying of the satellite orbit about earth in shorter imaging time Kepler's law, position and posture can be calculated by orbit equation or orbital tracking.

Select the earth's core rectangular coordinate system as object coordinates system, projection centre as satellite body coordinate origin, As z-axis, x-axis is located in the plane that satellite flight motion vector is constituted with z-axis and is directed toward position vector under object coordinates system Motion vector direction, y-axis are determined perpendicular to xz plane, direction by right-hand rule.Any scan line corresponds to it in image Imaging geometry can be expressed as the collinearity condition equation of instantaneous projection centre, image picture point and its corresponding object space point three.It is equivalent Geometry imaging model general type is expressed as following equation:

Wherein [l, s] is picpointed coordinate (l is row, and s is column), y₀Principal point pixel coordinate, f be equivalent master away from.Xs (l), Ys (l), Zs (l) are the corresponding projection centre instantaneous positions of scan line l, can be calculated according to orbit equation or orbital tracking It obtains.R_sat(l) be satellite platform transient posture matrix, express satellite platform coordinate system relative to geocentric rectangular coordinate system Transient posture satellite instantaneous position vector and its corresponding velocity vector (track can be passed through according to the definition of coordinate system The first derivative of equation) it is calculated.R_sensor(l) it is bias matrix of the sensor relative to satellite platform, has reacted image space Coordinate system is relative to the rotation between satellite platform coordinate system, for example, the pitch angle of forward sight, rear view camera, the side of side view camera The sensor of pivot angle, asynchronous imaging (quick satellite) is directed toward angle etc..Three rotations of the equivalent master away from, principal point coordinate, in bias matrix Corner and six orbital trackings, constitute the independent parameter of equivalent geometry imaging model, pass through these parameters of RPC parameter calculation Corresponding equivalent geometry imaging model can be established.

It is generally defined under inertial coodinate system with tight geometry imaging model and is either directly defined under object coordinates system Sensor attitude angle is different, and equivalent geometry imaging model is given a definition projection centre instantaneous position in geocentric rectangular coordinate system, respectively Satellite platform posture and sensor attitude are defined, and is further established conllinear between picpointed coordinate, projection centre and object space point Equation, it is fully equivalent in tight geometry imaging model in mathematical meaning.The foundation of EGSM does not need satellite platform and sensor Any structural parameters, whole model parameters are calculated by RPC coefficient, maintain the versatility of RFM.

2, the equivalent geometry imaging model parameter calculation based on RPC parameter

From formula (1) as can be seen that equivalent geometry imaging model foundation it needs to be determined that equivalent master away from f, principal point offset y₀, partially Set matrix R_sensor(l) three rotation angles and projection centre instantaneous position Xs (l), Ys (l), Zs (l) or satellite flight fortune in Dynamic equation.Internal and external orientation is integrally resolved with using space resection or bundle adjustment^[9]Difference, the present invention is according to pumping As the geometrical property of imaging model, corresponding parameter is directly calculated by RPC parameter, avoids the strong correlation between parameter, simultaneously Guarantee the stability calculated.It is worth noting that, in addition to RPC parameter, the foundation of equivalent geometry imaging model do not need it is any its Its information ensure that the versatility of model.

2.1, the corresponding object coordinates of picture point are calculated based on RPC parameter

As the mathematical simulation of tight geometry imaging model, rational function model (Rational Function Model, RFM picpointed coordinate and its corresponding object space point ground coordinate (generally indicating using WGS84 longitude, latitude and geodetic height)) are established Between functional relation.On ordinary meaning, first order expresses deformation of image caused by central projection, Atmosphere Refraction, optics Systematical distortion, earth curvature equal error are described by quadratic term, and various other systematic errors are simulated by cubic term.RFM is as high score Resolution satellite image general mathematical imaging model has become actual professional standard and is used by satellite image provider, and 80 Rational polynominal coefficient (constant term of two of them denominator is fixed as 1.0) and 10 normalized parameters are referred to as RPC (Rational Polynomial Coefficients) parameter, together with the satellite image Jing Guo sensing system error correction It mentions for users to use.RFM general expression is as follows,

Wherein, (e, n, h) is the geodetic coordinates of ground point, and Num (e, n, h) and Den (e, n, h) are cubic polynomial (packet The operation containing normalization).(L, S) is normalized picpointed coordinate.

Take up an official post the pixel coordinate (l, s) of one picture point of meaning to fixing, normalized coordinate is calculated as follows:

Wherein line_offset, line_scale, sample_offset, sample_scale are rational polynominal model Normalized parameter；

The corresponding ground elevation h of the picture point is further given, then can be iterated to calculate by equation (2) and obtain the picture point pair Answer the longitude and latitude (e, n) of object space point.Using WGS84 reference ellipsoid parameter, major semiaxis a and eccentric ratio e s can be by geodetic coordinates (e, n, h) is converted to geocentric rectangular coordinate system coordinate as the following formula:

2.2, equivalent master away from principal point coordinate calculate

According to equivalent geometrical model it is assumed that ideal line array CCD Pixel size is consistent, camera is placed in an ideal line On focal plane, conformation corresponds to a scan line of satellite image.As shown in Fig. 2, selecting the scanning on any scan line First pixel of row is as picture point c (l, 0), and any pixel is as picture point a (l, s_a).Geometrical relationship according to Fig.2, to picture point C (l, 0), picture point a (l, s_a) and projection centre between the angle that constitutes, equation can be listed below,

Tan (α)=tan (α₁+α₂)=(tan (α₁)+tan(α₂))/(1-tan(α₁)tan(α₂)) (5)

Wherein, α is picture point c (l, 0) and projection centre line, picture point a (l, s_a) and projection centre line between constitute Angle, α₁For picture point c (l, 0) and projection centre line, principal point o (l, y₀) and projection centre line between the angle that constitutes, α₂ For picture point a (l, s_a) and projection centre line, principal point o (l, y₀) and projection centre line between the angle that constitutes, equation is into one It is available after step arranges:

Tan (α)=f*s_a/[f²-y₀*(s_a-y₀)] (6)

As shown in Figure 1, giving the highest elevation h of RPC parameter expression respectively₁With minimum elevation h₂, according to the side of 2.2 introductions The corresponding object coordinates of picture point can be calculated in method.Geometrical relationship as shown in Figure 1 is appreciated that α is space line A₁_A₂ And C₁_C₂Between angle, tan (α) can be calculated according to their corresponding space coordinates by vector.To picture point c (l, And picture point b (l, s 0)_b) same equation can be listed, quadratic term can be eliminated using two equations and calculate equivalent master away from f and Principal point coordinate y₀Initial value.In order to guarantee the stability and precision that calculate, multi-strip scanning line, In can be selected on a scape image Select multiple picture points in every scan line, by formula (6) column error equation and linearize, thus form normal equation iterative calculation f and y₀Optimal solution.

2.3, arbitrary scan row projection centre calculates

In equivalent geometry imaging model, projected centered on the corresponding imaging geometry of the every one scan line of satellite image, so With picture points all on one scan line, corresponding imaging ray intersects at projection centre.Fig. 1 is high-resolution satellite image scan line The resolving signal of corresponding projection centre, selects picture point a (l, s on any scan line_a), the point RPC parameter expression is given respectively Highest elevation h₁With minimum elevation h₂, according to the method for 2.1 introductions, the corresponding object space point A1 and A2 of picture point a can be calculated Object coordinates, as shown in Figure 1, projection centre P (X_s, Y_s, Z_s) it is located at straight line A₁_A₂On, it thus can list following equation (7). Similarly, to picture point b (l, s_b) object space point B can also be calculated₁And B₂And corresponding equation (8) are listed,

According to the principle of least square, A is calculated by equation group (7) and (8)₁_A₂And B₁_B₂The intersection point of two space lines, That is the corresponding projection centre of the scan line (with 4 equation solvers, 3 unknown numbers).It, can be on scan line in order to calculate stabilization Multiple points are selected to calculate the projection centre coordinate Xs (l) that the scan line corresponds to the moment, Ys (l), Zs (l) with least square adjustment, The projection centre of i.e. every scan line is exactly the motion profile of satellite.

2.4, satellite motion track fitting

When not considering disturbing influence, satellite orbit can be by six Kepler orbit elements (parameter) under inertial coodinate system Description only needs two points or a point and its corresponding velocity vector on known track to can determine, the satellite of any time Position can be calculated by orbital tracking and orbit equation.Due to the transcendental equation that orbit equation is made of trigonometric function, It resolves to generally require and be realized by iteration, to reduce computational complexity to improve computational stability and efficiency, can use more Item formula is fitted track^[1,2,3].In the present invention, the instantaneous position of satellite flight movement consolidates rectangular coordinate system in ground heart and gives a definition, Its orbit equation can be obtained by chebyshev approximating polynomial.

Assuming that needing to be fitted image scan row [l₀, l_m] in satellite orbit, using formula (9) first by variable l ∈ [l₀, l_m] It is transformed into T ∈ [- 1,1],

τ=2 (l-l₀)/(l_m-l₀) -1, l ∈ [l₀, l_m] (9)

The Chebyshev polynomials of then satellite orbit coordinate Xs (l), Ys (l), Zs (l) are expressed are as follows:

In formula, n is the order of Chebyshev polynomials,The system of polynomials of respectively X, Y, Z coordinate component Number.T_i(τ) is determined according to recurrence formula (11).

T₀(τ)=1, T₁(τ)=τ, T_n+1(τ)=2 τ T_n(τ)-T_n-1(τ) (11)

Keplerian orbit is the ellipse at the uniform velocity rotated under geocentric rectangular coordinate system, experiments have shown that, in single scape image Imaging time in 3 ranks (2 times) Chebyshev polynomials can provide good fitting effect, even if longer time section is general Need not exceed 5 ranks.M item (m > n) scan line is selected in a scape image, calculates their corresponding throwings using 2.3 the methods Then shadow centre coordinate calculates Chebyshev polynomials coefficient of the satellite orbit in imaging time section by least square adjustment Using chebyshev approximating polynomial track, avoid in orbit computation largely using trigonometric function, instantaneous position Calculating stability and high efficiency is set, without any satellite platform and track knowledge is understood, simplifies the realization of equivalent geometry imaging model And application.

2.5, the calculating of the instantaneous spin matrix of satellite platform

It is defined according to equivalent geometry imaging model, geocentric rectangular coordinate system is used as and defends as object coordinates system, projection centre Star body coordinate system origin, for the position vector under object coordinates system as z-axis, x-axis is located at satellite flight motion vector and z In the plane that axis is constituted and it is directed toward motion vector direction, y-axis is determined perpendicular to xz plane, direction by right-hand rule.

Given vectorAnd corresponding velocity vector is calculated using formula (10) The corresponding unit vector of three reference axis of satellite platform can be calculated as follows:

It is hereby achieved that satellite platform coordinate system is as follows relative to the instantaneous spin matrix of object coordinates system:

Wherein velocity vector is not write by existing velocity vector calculation method, the present invention.

2.6, the resolving of sensor bias matrix

When known to projection centre instantaneous position and satellite platform transient posture matrix, equivalent geometry imaging model can letter It turns to:

Wherein, the virtual object space point coordinate and projection centre that [u (l), v (l), w (l)] is calculated according to RPC parameter are instantaneous Position and corresponding attitude matrix R_sat(l) it is calculated.Using picpointed coordinate as observation, above formula is directed to spin matrix R_sensorParameter linearisation, and by the principle of least square form normal equation, Iterative sensor bias matrix element.Consider It orbits the earth movement to satellite, satellite sensor image space coordinate different with traditional aerophotogrammetry exterior orientation angle element Posture relative to platform coordinate system is not necessarily low-angle, and it is special that high-resolution satellite image photographic light flux is narrow and field angle is small etc. Point leads between its image orientation parameter that there are very strong correlations, to influence the stability of elements of exterior orientation resolving, is simultaneously The problems such as avoiding Eulerian angles from indicating singularity existing for three-dimensional rotation and discontinuous angle interpolation, this algorithm use unit quaternary Number [q₀, q₁, q₂, q₃] expression posture spin matrix^[10,11].Unit quaternion meetsIts corresponding rotation Torque battle array is calculated by following formula.

According to quaternion differential equation^[11], attitudes vibration expressed by four elements can be expressed as coordinate system around three of them The update of the rotation of reference axis, quaternary number can be calculated with following formula.

To reduce the correlation between polynomial parameters to the greatest extent, the present invention selects Chebyshev's orthogonal polynomial expression posture to become Change item w_x, w_y, w_z.By spin matrix (15) to unit quaternion derivation and with the substitution of (16) formula, then it is opposite that spin matrix can be obtained It in the linearized expression of three independent parameters, then substitutes into (14) formula and eliminates proportionality coefficient λ, obtain solving sensor biasing The error equation at angle.The present invention selects every scape image multi-strip scanning line (more than the order of Chebyshev polynomials), every scanning The multiple picture points of line (more than two picture point) arrange vertical error equation by formula (14) and linearize, thus form normal equation iterative calculation w_x, w_y, w_zCorresponding Chebyshev's orthogonal polynomial coefficient optimal solution, and arbitrary scan line is calculated according to above formula and corresponds to Spin matrix R_sensor(l)。

Based on algorithm is given above, the internal and external orientation in equivalent geometry imaging model may be by RPC parameter It being calculated, the line element and angle element in elements of exterior orientation can also solve respectively, by using unit quaternion, even if Inclination angle is larger and without that can also be restrained with fast and stable in the case where initial value.

3. the block adjustment based on equivalent geometry imaging model

Instantaneous position and satellite platform transient posture are moved by disturbing influence in view of satellite orbits the earth, and with additional Parametric form expresses picpointed coordinate systematic error, based on equivalent geometry imaging model satellite image block adjustment function model General type is as follows:

Wherein, Δ x, Δ y are the systematic error indicated in the form of additional parameter and model error, [dXs (l), dYs (l), DZs (l)] it is disturbance part of the satellite position with sporting flying, D_ori(l) be various influences of the disturbance to the attitude of satellite, they It is numerically that in a small amount, polynomial approximation can be used, this three parts constitutes the unknown number of self-calibration block adjustment system.

After eliminating λ, the corresponding collinearity equation of each picture point (17) formula can be expressed as two error sides of following form Journey:

Wherein (v_x, v_y) it is that picture point corrects residual error, (x, y, f) is picture point observation and equivalent master away from (u, v, w) is basis The auxiliary coordinates coordinate that object coordinates initial value, projection centre coordinate and spin matrix calculate, p are power；

U in above formula, v, w are linearized, the error of the abstract geometry imaging model of available matrix form after arrangement Equation:

V₁=B₁₁X₁+B₁₂X₂+l₁, P₁ (19)

Wherein V₁For the expression matrix of picture point residual error, B₁₁, B₁₂Error equation factor arrays, l₁For constant term matrix, P₁For power Matrix, X₁And X₂For unknown matrix number.

It is described in the present invention that specific embodiments are merely illustrative of the spirit of the present invention.Technology belonging to the present invention The technical staff in field can make various modifications or additions to the described embodiments or by a similar method Substitution, however, it does not deviate from the spirit of the invention or beyond the scope of the appended claims.

Bibliography

[1] Toutin T., Geometric Processing of Remote Sensing Images:Models, Algorithms and Methods [J] .International Journal of Remote Sensing, 2004,25 (10): 1893-1924

[2] POLI D., TOUTIN T., Review of Developments in Geometric Modelling for High Resolution Satellite Pushbroom Sensors[J].The Photogrammetric Record, 2012,27 (137): 58-73

[3] POLI D., A rigorous model for spaceborne linear array sensors [J] .Photogrammetric Engineering and Remote Sensing, 2007,73 (2): 187-196.

[4] Tao C.V., Hu Y., A Comprehensive Study of the Rational Function Model for Photogrammetric Processing[J].Photogrammetric Engineering and Remote Sensing, 2001,67 (12): 1347-1357.

[5] Yuan Xiuxiao, Lin Xianyong, the rational polynominal parametric solution method [J] based on ridge estimaion, Wuhan University Journal: Information science version, 2008,33 (11): 1130-1133.

YUAN Xiuxiao, LIN Xianyong, A Method for Solving Rational Polynomial Coefficients Based on Ridge Estimation[J].Geomatics and Information Science Of Wuhan University, 2008,33 (11): 1130-1133.

[6] Zhang Yongjun, Wang Lei, Lu Yihui satellite remote-sensing image rational function model optimization method [J] survey and draw journal, 2011,40 (6): 756-761.

ZHANG Yongjun, WANGLei, LUYihui.Optimization of the Rational Function Model of Satellite Imagery [J] .Acta Geodaetica et CartographicaSinica, 2011,40 (6): 756-761.

[7] FRASER CS, HANLEY H B, Bias Compensation in Rational Functions for IKONOS Satellite Imagery [J] .Photogrammetric Engineering and Remote Sensing, 2003,69 (1): 53-57

[8] GRODECKI J., DIAL G., Block Adjustment of High-resolution Satellite Images Described by Rational Polynomials[J].Photogrammetric Engineering and Remote Sensing, 2003,69 (1): 59-68.

[9] DI Kaichang, MA Ruijin, LI Rongxing, Rational Functions and Potential for Rigorous Sensor Model Recovery[J].Photogrammetric Engineering&Remote Sensing, 2003,60 (1): 33-41.

[10] Wuhan space resection [J] that Yan Li, Nie Qian, Zhao Zhan state linear CCD image using quaternary sketch is big Learn journal: information science version, 2010,35 (2): 201-204.

YAN Li, NIE Qian, ZHAO Zhan.Space Resection of Line Scanner CCD Image Based on the Description of Quaternions[J].Geomatics and Information Science Of Wuhan University, 2010,35 (2): 201-204.

[11] Gong Hui, Jiang Ting, Jiang Gangwu wait the high resolution ratio satellite remote-sensing image foreign side bit of the tetra- element differential equation of Element solves [J] and surveys and draws journal, 2012,41 (3): 409-416.

GONG Hui, JIANG Ting, JIANG Gangwu, et al.Solution of Exterior Orientation Parameters for High-resolution Satellite Imagery Based on Quaternion Differential Equation [J] .Acta Geodaetica et CartographicaSinica, 2012,41 (3): 409-416.