A method for real-time positioning of an airborne moving target by dual-satellite cooperative observation

Abstract

The application discloses a kind of double star cooperative observation air movement target real-time positioning method, steps are as follows: 1, using target detection and tracking algorithm, respectively from the air movement target sequence image that two optical remote sensing satellites current time simultaneously shoots in detection target's particle coordinates;2, according to optical remote sensing satellite imaging mechanism and air target movement characteristics, constructs air movement target movement track fusion model;3, using movement track sliding solution method: using the particle coordinates of air movement target on two satellite sequence images, constructs standard error equation and iteratively solves air target movement track parameters;4, according to movement track parameters, solves the spatial position and speed of air movement target at current time.The application can eliminate the influence of double star time synchronization error on air movement target positioning accuracy from model itself, realizes the high-precision real-time positioning of double star cooperative observation air movement target.

Description

TECHNICAL FIELD

[0001] The present application belongs to the technical field of optical remote sensing satellite data geometric processing, and particularly relates to a method for real-time positioning of an aerial moving target based on cooperative observation of two satellites. BACKGROUND

[0002] From the perspective of high-precision positioning of moving targets, moving targets can be mainly divided into three categories: sea surface moving targets, ground moving targets and aerial moving targets. For aerial moving targets, the target flight height is constantly changing and unknown, and it is difficult to know the height value in advance, and it is difficult to achieve high-precision positioning of the target by single satellite single observation. Compared with static targets, aerial moving targets are in a moving state, and it is also difficult to achieve high-precision positioning by single satellite multiple observations. In view of the characteristics that the height of the aerial moving target is unknown and in a moving state, two or more remote sensing satellites need to be used for cooperative observation, and time synchronization and observation from different angles need to be achieved, so as to achieve high-precision positioning of the aerial moving target.

[0003] Due to the influence of factors such as imaging frame frequency and time control of optical remote sensing satellites, when two satellites are used for cooperative observation of aerial moving targets, it is difficult for the two satellites to strictly achieve time synchronization, that is, there is a time synchronization error. The conventional model and method for target positioning of optical remote sensing satellite images are mainly aimed at static targets, and do not need to consider the time synchronization error of two satellites. However, for aerial moving targets, the positioning model and method for conventional static targets are used for positioning processing, the time synchronization error of two satellites is ignored, and the positioning result will inevitably deviate from the true position. Therefore, the existing model and method have obvious deficiencies in high-precision real-time positioning of aerial moving targets, and it is difficult to achieve high-precision positioning of aerial moving targets. SUMMARY

[0004] The technical problem to be solved by the present application is to overcome the deficiencies of the prior art, and to provide a method for real-time positioning of an aerial moving target by constructing a target motion trajectory fusion model based on cooperative observation of two satellites, so as to achieve high-precision positioning of the aerial moving target.

[0005] The technical solution of the present application is a method for real-time positioning of an aerial moving target based on cooperative observation of two satellites, and the steps are as follows:

[0006] I. A target detection and tracking algorithm is used to detect the particle coordinates of the target in the image coordinate system from the aerial moving target sequence images taken by the first and second remote sensing satellites at the current time, respectively;

[0007] II. According to the imaging mechanism of optical remote sensing satellites and the position, velocity and acceleration information of the aerial target, an aerial moving target motion trajectory fusion model is constructed on the first and second remote sensing satellites, respectively;

[0008] III. Take the motion trajectory sliding solution method: for the current time first remote sensing satellite image on the detected air target particle coordinates, using the first remote sensing satellite before the time of continuous N frame and the second remote sensing satellite before the time of continuous M frame image on the air target particle coordinates, construct the standard error equation, iterative solution of the current time air target motion trajectory fusion model parameters, M>1, N>1;

[0009] IV. The air target motion trajectory fusion model parameters obtained in step III are substituted into the position equation to solve the spatial position (X, Y, Z) of the air moving target at the current time, and substituted into the velocity equation to solve the velocity (V X , V Y , V Z ) of the air moving target at the current time.

[0010] Further, the target detection and tracking algorithm in step I is a time domain target detection and tracking algorithm or a spatial domain target detection and tracking algorithm.

[0011] Further, the air moving target motion trajectory fusion model constructed for each remote sensing satellite in step II is as follows:

[0012]

[0013]

[0014] In the above formula, t is the imaging time; X0, δ1, δ2 are the initial value of the X direction position of the air target, the 1st order term of the X direction position about t and the 2nd order term of the X direction position about t parameters; Y0, β1, β2 are the initial value parameter of the Y direction position of the air target, the 1st order term parameter of the Y direction position about t and the 2nd order term parameter of the Y direction position about t; Z0, θ1, θ2 are the initial value of the Z direction position of the air target, the 1st order term of the Z direction position about t and the 2nd order term of the Z direction position about t parameters; (X S，1 , Y S，1 , Z S，1 ) is the space coordinate of the first remote sensing satellite in the WGS84 coordinate system, which is measured by the on-board GNSS measurement device, (X S，2 , Y S，2 , Z S，2 ) is the space coordinate of the second remote sensing satellite in the WGS84 coordinate system; λ1 is the first scale factor, and λ2 is the second scale factor;

[0015] is the rotation matrix from J2000 coordinate system to WGS84 coordinate system; is the rotation matrix from the first remote sensing satellite attitude measurement coordinate system to the J2000 coordinate system, is a rotation matrix of the second remote sensing satellite attitude measurement coordinate system to the J2000 coordinate system; a camera on the first remote sensing satellite is referred to as a first camera, and a camera on the second remote sensing satellite is referred to as a second camera, is a placement matrix of the first camera in the first remote sensing satellite attitude measurement coordinate system, is a placement matrix of the second camera in the second remote sensing satellite attitude measurement coordinate system;

[0016] are pointing angles of a first camera imaging element in the first remote sensing satellite attitude measurement coordinate system along the track direction and the cross-track direction, respectively, are pointing angles of a second camera imaging element in the second remote sensing satellite attitude measurement coordinate system along the track direction and the cross-track direction, respectively;

[0017] The origin of the first remote sensing satellite attitude measurement coordinate system is located at the center of mass of an attitude measurement device on the first remote sensing satellite, the X-axis is the flight direction of the first remote sensing satellite, the Z-axis points to the center of the Earth, and the Y-axis is perpendicular to the Z-axis and the X-axis to form a right-hand system; the origin of the second remote sensing satellite attitude measurement coordinate system is located at the center of mass of an attitude measurement device on the second remote sensing satellite, the X-axis is the flight direction of the second remote sensing satellite, the Z-axis points to the center of the Earth, and the Y-axis is perpendicular to the Z-axis and the X-axis to form a right-hand system.

[0018] Further, the pointing angle of the first camera imaging element in the first remote sensing satellite attitude measurement coordinate system along the track direction is and the pointing angle along the cross-track direction is The pointing angle of the second camera imaging element in the second remote sensing satellite attitude measurement coordinate system along the track direction is and the pointing angle along the cross-track direction is The pointing angles are calculated according to the following method:

[0019]

[0020]

[0021] In the above formula, (s1, l1) is the coordinate of a target mass point detected by the first remote sensing satellite in the image coordinate system, (s2, l2) is the coordinate of the target mass point detected by the second remote sensing satellite in the image coordinate system, s1 is the column number of the target mass point coordinate on the first remote sensing satellite, s2 is the column number of the target mass point coordinate on the second remote sensing satellite, l1 is the row number of the target mass point coordinate on the first remote sensing satellite, and l2 is the row number of the target mass point coordinate on the second remote sensing satellite;

[0022] h 0，1 is a constant term of the first remote sensing satellite pointing angle along the track direction h 1，1 , h 2，1 , h 3，1 , h4，1 h 5，1 h 6，1 h 7，1 h 8，1 h 9，1 These correspond to the pointing angles along the orbit of the first remote sensing satellite. In the direction s1, l1, s1l1, The coefficient, k 0，1 The pointing angle of the vertical orbit of the first remote sensing satellite The constant term in the direction, k 1，1 k 2，1 k 3，1 k 4，1 k 5，1 k 6，1 k 7，1 k 8，1 k 9，1 Corresponding to the pointing angles of the vertical orbit of the first remote sensing satellite In the direction s1, l1, s1l1, The coefficient;

[0023] h 0，2 The pointing angle along the orbit of the second remote sensing satellite The constant term in the direction; h 1，2 h 2，2 h 3，2 h 4，2 h 5，2 h 6，2 h 7，2 h 8，2 h 9，2 These correspond to the pointing angles along the orbit of the second remote sensing satellite. In the direction s2, l2, s2l2, The coefficient, k 0，2 The pointing angle of the vertical orbit of the second remote sensing satellite The constant term in the direction, k 1，2 k 2，2 k 3，2 k 4，2 k 5，2 k 6，2 k 7，2 k 8，2 k 9，2 Corresponding to the pointing angles of the vertical orbit of the second remote sensing satellite In the direction s2, l2, s2l2, The coefficient.

[0024] Further, definition

[0025]

[0026] u1, v1, w1 are respectively the first row, the second row, the third row of the column vector obtained by the matrix multiplication of the above formula; the first scale factor λ1 is calculated according to the following ellipsoid equation:

[0027]

[0028] In the above formula, A=a e +h, B=b e +h, a e =6378137.0m, b e =6356752.3m, h is the ellipsoid height.

[0029] Further, define

[0030]

[0031] u2, v2, w2 are respectively the first row, the second row, the third row of the column vector obtained by the matrix multiplication of the above formula; the second scale factor λ2 is calculated according to the following ellipsoid equation:

[0032]

[0033] In the above formula, A=a e +h, B=b e +h, a e =6378137.0m, b e =6356752.3m, h is the ellipsoid height.

[0034] Further, the sum of the number of images before the first remote sensing satellite and the second remote sensing satellite required in step three satisfies: M+N≥9.

[0035] Further, the motion trajectory sliding solving method of step three comprises the following steps:

[0036] S1, define:

[0037] R1, R2 are orthogonal matrices,

[0038] The motion trajectory fusion model of the first and second remote sensing satellites is deformed as follows:

[0039]

[0040]

[0041] The above two equations are left multiplied by R -1 , and the following formula can be obtained by four arithmetic operations:

[0042]

[0043]

[0044] In the above formula: F x，1 , F x，2 are the pointing angle residual function models of the first remote sensing satellite and the second remote sensing satellite along the orbit direction respectively, F y，1 , F y，2 are the pointing angle residual function models of the first remote sensing satellite and the second remote sensing satellite along the orbit direction respectively;

[0045] S2, for each target particle coordinate of the aerial moving target on the sequence images of the first remote sensing satellite and the second remote sensing satellite, a standard error equation is constructed, which is in the following form:

[0046] V = AX - L P

[0047] In the above formula: V is the residual vector of X, which is in the following form:

[0048]

[0049] is the calculation result of the right side of the standard error equation;

[0050] A is a design matrix composed of unknown partial derivatives, which is in the following form:

[0051]

[0052] Wherein, the subscripts 1 and 2 respectively represent the first remote sensing satellite and the second remote sensing satellite, the subscript i is the frame number of the first remote sensing satellite sequence image where the aerial moving target is located, i = 1, 2,..., N, and the subscript j is the frame number of the second remote sensing satellite sequence image where the aerial moving target is located, j = 1, 2,..., M;

[0053] X is an unknown matrix, which is iteratively calculated according to the least square adjustment principle, and the calculation formula is:

[0054] X = (A T PA) -1 (A T PL)

[0055] X is in the following form:

[0056] X = [X0, δ1, δ2, Y0, β1, β2, Z0, θ1, θ2]

[0057] When calculating the standard error equation for the first time, the initial value of X is assigned;

[0058] L is the current calculation value of the unknown matrix X substituted into F x，1 , Fx，2 , F y，1 , F y，2 The constant term matrix obtained by the expression of F

[0059] L=[F x，1，i , F y，1，i , F x，2，j , F y，2，j ] T

[0060] P is the observation weight of the target point coordinate, and is a unit matrix;

[0061] S3, according to the least square adjustment principle, the motion trajectory parameters X0, δ1, δ2, Y0, β1, β2, Z0, θ1, θ2 of the aerial target at the iteration stop are solved.

[0062] Further, according to the least square adjustment principle, the condition for X iteration convergence is that the difference between X0, δ1, δ2, Y0, β1, β2, Z0, θ1, θ2 obtained by the current and the last two times is within 1x10 -6 , it is considered that the iteration converges, and the iteration is stopped; and the value of the unknown matrix X obtained by the last iteration is taken as the obtained aerial target motion trajectory fusion model parameters.

[0063] Further, the method for solving the spatial position (X, Y, Z) and the velocity (V X , V Y , V Z ) of the aerial motion target at the current time in the fourth step is: according to the obtained aerial target motion trajectory fusion model parameters X0, δ1, δ2, Y0, β1, β2, Z0, θ1, θ2, the following position equation and velocity equation are substituted to obtain:

[0064]

[0065] The present application expands and constructs the aerial motion target motion trajectory fusion model on the basis of the conventional optical remote sensing satellite image static target positioning model, eliminates the influence of the double-star time synchronization error on the target positioning accuracy from the model itself, and adopts a motion trajectory sliding solving strategy to obtain the spatial position and the velocity of the aerial motion target at the current time in real time, and realizes the high-precision real-time positioning of the aerial motion target observed by the double satellites.

[0066] Compared with the prior art, the present application has the following advantages:

[0067] (1) The application expands and constructs an air moving target motion trajectory fusion model on the basis of a conventional static target positioning model, can eliminate the influence of double-star time synchronization error on the positioning precision of the air moving target from the model itself, and realizes high-precision positioning of the air moving target in the optical remote sensing satellite image.

[0068] (2) The application adopts a motion trajectory sliding solution method, uses the particle information of the air moving target at a plurality of previous time points, increases the redundancy observation of the current time motion trajectory solution, can improve the positioning precision of the air moving target, realizes real-time positioning of the air moving target, and meets the high timeliness demand of the air moving target positioning processing. BRIEF DESCRIPTION OF DRAWINGS

[0069] Figure 1 is a specific flowchart of the embodiment of the application.

[0070] Figure 2 is a sliding solution schematic diagram of the air moving target motion trajectory of the application. DETAILED DESCRIPTION

[0071] In order to more clearly illustrate the technical solutions in the embodiments of the application and / or the prior art, the specific embodiments of the application will be described below with reference to the drawings.

[0072] The air moving target real-time positioning method provided by the embodiment of the application is as shown in Figure 1 The specific steps are as follows:

[0073] I. The target detection and tracking algorithm is used to detect the particle coordinates of the target in the image coordinate system from the air moving target sequence images taken by the first and second remote sensing satellites at the current time respectively; in the embodiment, the time domain target detection and tracking algorithm is used to obtain the particle coordinates of the target in the image coordinate system.

[0074] II. According to the imaging mechanism of the optical remote sensing satellite, that is, the collinear condition equation in photogrammetry and the optical physical imaging principle, a rigorous optical imaging geometric model of the optical remote sensing satellite and the position, velocity and acceleration motion characteristics of the air target are constructed, and the air moving target motion trajectory fusion model of the first and second remote sensing satellites is constructed as follows:

[0075]

[0076]

[0077] In the above formula, t is the imaging time; X0, δ1, δ2 are the initial value, the first-order term of t and the second-order term of t parameters of the target X direction in the air, respectively; Y0, β1, β2 are the initial value, the first-order term of t and the second-order term of t parameters of the target Y direction in the air, respectively; Z0, θ1, θ2 are the initial value, the first-order term of t and the second-order term of t parameters of the target Z direction in the air, respectively; (X S，1 , Y S，1 , Z S，1 ) are the space coordinates of the first remote sensing satellite in the WGS84 coordinate system, which are measured by the on-board GNSS measurement device, (X S，2 , Y S，2 , Z S，2 ) are the space coordinates of the second remote sensing satellite in the WGS84 coordinate system; λ1 is the first scale factor, and λ2 is the second scale factor; is the rotation matrix from the J2000 coordinate system to the WGS84 coordinate system; the current satellite attitude measurement coordinate system is defined with the origin at the center of mass of the attitude measurement device, the X axis in the direction of flight of the satellite, the Z axis pointing to the center of the earth, and the Y axis perpendicular to the Z axis and the X axis to form a right-handed system, is the rotation matrix from the first remote sensing satellite attitude measurement coordinate system to the J2000 coordinate system, is the rotation matrix from the second remote sensing satellite attitude measurement coordinate system to the J2000 coordinate system; is the installation matrix of the first camera in the first remote sensing satellite attitude measurement coordinate system, is the installation matrix of the second camera in the second remote sensing satellite attitude measurement coordinate system; are the pointing angles of the first camera imaging probe element in the first remote sensing satellite attitude measurement coordinate system along the track direction and the cross-track direction, respectively, are the pointing angles of the second camera imaging probe element in the second remote sensing satellite attitude measurement coordinate system along the track direction and the cross-track direction, respectively.

[0078] the pointing angle along the track direction and the pointing angle along the cross-track direction are calculated according to the following method:

[0079]

[0080]

[0081] In the above formula, (s1, l1) is the coordinate of the target particle detected by the first remote sensing satellite in the image coordinate system, (s2, l2) is the coordinate of the target particle detected by the second remote sensing satellite in the image coordinate system, s1 is the column number of the target particle coordinate on the first remote sensing satellite, s2 is the column number of the target particle coordinate on the second remote sensing satellite, l1 is the row number of the target particle coordinate on the first remote sensing satellite, and l2 is the row number of the target particle coordinate on the second remote sensing satellite.

[0082] h 0，1 The pointing angle along the orbit of the first remote sensing satellite The constant term in the direction, h 1，1 h 2，1 h 3，1 h 4，1 h 5，1 h 6，1 h 7，1 h 8，1 h 9，1 These correspond to the pointing angles along the orbit of the first remote sensing satellite. In the direction s1, l1, s1l1, The coefficient, k 0，1 The pointing angle of the vertical orbit of the first remote sensing satellite The constant term in the direction, k 1，1 k 2，1 k 3，1 k 4，1 k 5，1 k 6，1 k 7，1 k 8，1 k 9，1 Corresponding to the pointing angles of the vertical orbit of the first remote sensing satellite In the direction s1, l1, s1l1, The coefficient;

[0083] h 0，2 The pointing angle along the orbit of the second remote sensing satellite The constant term in the direction; h 1，2 h 2，2 h 3，2 h 4，2 h 5，2 h 6，2 h 7，2 h 8，2 h 9，2 These correspond to the pointing angles along the orbit of the second remote sensing satellite. In the direction s2, l2, s2l2, The coefficient, k 0，2 The pointing angle of the vertical orbit of the second remote sensing satellite The constant term in the direction, k 1，2 k 2，2 k 3，2 k 4，2 k 5，2 k 6，2 k 7，2 k 8，2 k 9，2 Corresponding to the pointing angles of the vertical orbit of the second remote sensing satellite s2, l 2， s2l2, coefficient of.

[0084] The calculation method of the first scale factor λ1 and the second scale factor λ2 is as follows:

[0085] Definition

[0086]

[0087] u1, v1, w1 are the first row, the second row and the third row of the column vector obtained by the matrix multiplication; the first scale factor λ1 is calculated according to the following ellipsoid equation:

[0088]

[0089] In the above formula, A=a e +h, B=b e +h, a e =6378137.0m, b e =6356752.3m, and h is the ellipsoid height.

[0090] Similarly, define

[0091]

[0092] u2, v2, w2 are the first row, the second row and the third row of the column vector obtained by the matrix multiplication; the second scale factor λ2 is calculated according to the following ellipsoid equation:

[0093]

[0094] In the above formula, A=a e +h, B=b e +h, a e =6378137.0m, b e =6356752.3m, and h is the ellipsoid height.

[0095] III. The motion trajectory sliding solving method is adopted: for the coordinates of the aerial target particles detected on the image of a satellite at the current time, the coordinates of the aerial target particles on the continuous N frames of images before the current time of the satellite and the continuous M frames of images before the current time of another satellite are used to construct a standard error equation, and the motion trajectory fusion model parameters X0, δ1, δ2, Y0, β1, β2, Z0, θ1, θ2 are iteratively solved; M+N≥9.

[0096] The motion trajectory sliding solving method comprises the following steps:

[0097] S1, according to the established motion trajectory fusion model and the pointing angle The calculation result defines R1, R2 are orthogonal matrices;

[0098] The following deformation is made to the motion trajectory fusion model of the first and second remote sensing satellites, to obtain

[0099]

[0100]

[0101] The left side of the above two equations is multiplied by R -1 The following equation can be obtained:

[0102]

[0103]

[0104] The above two formulas are converted to the following two formulas by four arithmetic operations:

[0105]

[0106]

[0107] In the above formula, F x，1 , F x，2 are the pointing angle residual function models of the first remote sensing satellite and the second remote sensing satellite along the orbit direction, respectively, and F y，1 , F y，2 are the pointing angle residual function models of the first remote sensing satellite and the second remote sensing satellite in the direction perpendicular to the orbit, respectively.

[0108] S2, for each mass point coordinate of the aerial moving target on the sequence images of the first remote sensing satellite and the second remote sensing satellite, a standard error equation is constructed, which has the following form:

[0109] V = AX - L P

[0110] In the above formula, V is the residual vector of X, which has the following form:

[0111]

[0112] is the calculation result of the right side of the standard error equation;

[0113] A is a design matrix composed of unknown partial derivatives, which has the following form:

[0114]

[0115] Wherein: subscript 1 and 2 respectively represent the first remote sensing satellite and the second remote sensing satellite, subscript i is the frame number of the first remote sensing satellite sequence image where the aerial moving target is located, i = 1, 2,..., N, and subscript j is the frame number of the second remote sensing satellite sequence image where the aerial moving target is located, j = 1, 2,..., M; i Wherein: subscript 1 and 2 respectively represent the first remote sensing satellite and the second remote sensing satellite, subscript i is the frame number of the first remote sensing satellite sequence image where the aerial moving target is located, i = 1, 2,..., N, and subscript j is the frame number of the second remote sensing satellite sequence image where the aerial moving target is located, j = 1, 2,..., M;

[0116] X is an unknown matrix, which is iteratively calculated according to the least square adjustment principle, and the representation of X is as follows:

[0117] X = [X0, δ1, δ2, Y0, β1, β2, Z0, θ1, θ2]

[0118] At the beginning of the calculation, the initial value of X needs to be assigned, and in the embodiment, the initial value of X is set to [0 0 0 0 0 0 0 0 0] T . Each time the iteration calculation of the standard error equation is performed, the value of the unknown matrix X changes;

[0119] L is a constant matrix obtained by substituting the calculated value of the unknown matrix X into the expression of F x，1 , F x，2 , F y，1 , F y，2 , and the representation is as follows:

[0120] L = [F x，1，i , F y，1，i , F x，2，j , F y，2，j ] T

[0121] At the beginning of the iteration, the initial value of X is substituted into F x，1 , F x，2 , F y，1 , F y，2 to obtain the initial value of L, and the next round of iteration calculation is started.

[0122] P is the observation weight of the target particle coordinates, which is a unit matrix.

[0123] In the iteration process, the calculation method of the unknown matrix X is iteratively calculated according to the least square adjustment principle:

[0124] X = (A T PA) -1 (A T PL)

[0125] S3, according to the least square adjustment principle, the aerial target motion trajectory parameters X0, δ1, δ2, Y0, β1, β2, Z0, θ1, θ2 at the time of stopping iteration are solved.

[0126] The condition for determining the convergence of the X iteration is: the difference between the X0, δ1, δ2, Y0, β1, β2, Z0, θ1, and θ2 obtained from the previous two calculations is all within 1×10. -6 Within a certain range, the iteration is considered to have converged, and the iteration stops; the value of the unknown matrix X obtained from the last iteration is used as the parameter of the obtained aerial target motion trajectory fusion model.

[0127] IV. Using the fusion model parameters of the aerial target's motion trajectory obtained in Step III, solve for the spatial position (X, Y, Z) and velocity (V) of the aerial target at the current moment. X V Y V Z The calculation method is as follows:

[0128]

[0129] By substituting the imaging time parameter t at each moment, the motion trajectory of the aerial target can be obtained.

[0130] like Figure 2 As shown, the principle of using the sliding motion trajectory solution method to achieve real-time positioning of aerial moving targets is as follows: For the aerial moving target particle p6 and p5 detected on the first remote sensing satellite image at the current time t6, the parameters X0, δ1, δ2, Y0, β1, β2, Z0, θ1, θ2 of the aerial target's motion trajectory L6 at time t6 can be accurately solved by combining the target particles p2 and p4 on several consecutive frames (only two frames are listed here) of the first remote sensing satellite image before the current time, and the target particles p1 and p3 on several consecutive frames (only two frames are listed here) of the second remote sensing satellite image before the current time. Similarly, for the aerial target particle p7 at time t7, the trajectory parameters of the motion trajectory L7 at time t7 can be jointly solved by using the target particles p3, p5, and p7, as well as p2, p4, and p6.

[0131] Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make possible changes and modifications to the technical solutions of the present invention by utilizing the methods and techniques disclosed above without departing from the spirit and scope of the present invention. Therefore, any simple modifications, equivalent changes and alterations made to the above embodiments based on the technical essence of the present invention without departing from the content of the technical solutions of the present invention shall fall within the protection scope of the technical solutions of the present invention.

Claims

1. A method for real-time positioning of an airborne moving target by dual-satellite cooperative observation, characterized in that, The method comprises the following steps: I. using a target detection and tracking algorithm to detect target particle coordinates in an image coordinate system from a sequence of images of the air target taken by the first and second remote sensing satellites at the current time; II. constructing an air target motion trajectory fusion model on the first and second remote sensing satellites according to the imaging mechanism of the optical remote sensing satellite and the position, velocity and acceleration information of the air target; III. using a motion trajectory sliding solution method: for the particle coordinates of the air target detected on the first remote sensing satellite image at the current time, the particle coordinates of the air target on the continuous N frames of images of the first remote sensing satellite at previous time and the continuous M frames of images of the second remote sensing satellite at previous time are used to construct a standard error equation, and the motion trajectory fusion model parameters of the air target at the current time are iteratively solved, M>1, N>1. Four, the position equation is substituted with the fusion model parameter of the motion trajectory of the aerial target obtained in step three, the spatial position (X, Y, Z) of the aerial moving target at the current time is solved, the velocity equation is substituted, and the velocity (V X , V Y , V Z ) of the aerial moving target at the current time is solved.

2. The method according to claim 1, wherein: The target detection and tracking algorithm in step I is a time domain target detection and tracking algorithm or a space domain target detection and tracking algorithm.

3. The real-time positioning method for aerial moving targets based on dual-satellite collaborative observation according to claim 1, characterized in that: The air target motion trajectory fusion model constructed by the first remote sensing satellite in step II is as follows: The air target motion trajectory fusion model constructed by the second remote sensing satellite is as follows: In the above formula, t is the imaging time; X0, δ1, δ2 are respectively the initial value of the X direction position of the aerial target, the first order term of the X direction position with respect to t and the second order term of the X direction position with respect to t parameters; γ0, β1, β2 are respectively the initial value parameter of the Y direction position of the aerial target, the first order term of the Y direction position with respect to t and the second order term of the Y direction position with respect to t parameters; Z0, θ1, θ2 are respectively the initial value of the Z direction position of the aerial target, the first order term of the Z direction position with respect to t and the second order term of the Z direction position with respect to t parameters; (X S，1 , Y S，1 , Z S，1 ) are the space coordinates of the first remote sensing satellite in the WGS84 coordinate system, which are measured by the on-board GNSS measurement device, (X S，2 , Y S，2 , Z S，2 ) are the space coordinates of the second remote sensing satellite in the WGS84 coordinate system; λ1 is the first scale factor, and λ2 is the second scale factor; is a rotation matrix from the J2000 coordinate system to the WGS84 coordinate system; is a rotation matrix from the first remote sensing satellite attitude measurement coordinate system to the J2000 coordinate system, is a rotation matrix from the second remote sensing satellite attitude measurement coordinate system to the J2000 coordinate system; the camera on the first remote sensing satellite is referred to as the first camera, and the camera on the second remote sensing satellite is referred to as the second camera, is a placement matrix of the first camera under the first remote sensing satellite attitude measurement coordinate system, is a placement matrix of the second camera under the second remote sensing satellite attitude measurement coordinate system; respectively, are the pointing angles of the first camera imaging probe element along the along-track direction and the cross-track direction in the first remote sensing satellite attitude measurement coordinate system, respectively, are the pointing angles of the second camera imaging probe element along the along-track direction and the cross-track direction in the second remote sensing satellite attitude measurement coordinate system; The origin of the attitude measurement coordinate system of the first remote sensing satellite is located at the center of mass of the attitude measurement device on the first remote sensing satellite, the X axis is the flight direction of the first remote sensing satellite, the Z axis points to the center of the earth, and the Y axis is perpendicular to the Z axis and the X axis to form a right-handed system; The origin of the attitude measurement coordinate system of the second remote sensing satellite is located at the center of mass of the attitude measurement device on the second remote sensing satellite, the X axis is the flight direction of the second remote sensing satellite, the Z axis points to the center of the earth, and the Y axis is perpendicular to the Z axis and the X axis to form a right-handed system.

4. The method according to claim 3, wherein: the pointing angle of the first camera imaging sensor along the track direction in the first remote sensing satellite attitude measurement coordinate system and the pointing angle along the cross-track direction the pointing angle of the second camera imaging sensor along the track direction in the second remote sensing satellite attitude measurement coordinate system and the pointing angle along the cross-track direction calculated according to the following method: In the above formula: (s1, l1) is the coordinates of the target particle detected by the first remote sensing satellite in the image coordinate system, (s2, l2) is the coordinates of the target particle detected by the second remote sensing satellite in the image coordinate system, s1 is the column number of the target particle coordinates on the first remote sensing satellite, s2 is the column number of the target particle coordinates on the second remote sensing satellite, l1 is the row number of the target particle coordinates on the first remote sensing satellite, and l2 is the row number of the target particle coordinates on the second remote sensing satellite; h 0，1 the constant term in the direction of the along-track pointing angle of the first remote sensing satellite h 1，1 , h 2，1 , h 3，1 , h 4，1 , h 5，1 , h 6，1 , h 7，1 , h 8，1 , h 9，1 corresponding to the along-track pointing angle of the first remote sensing satellite the coefficient of s1, l1, s1l1, in the direction of the along-track pointing angle of the first remote sensing satellite 0，1 the constant term in the direction of the cross-track pointing angle of the first remote sensing satellite k 1，1 , k 2，1 , k 3，1 , k 4，1 , k 5，1 , k 6，1 , k 7，1 , k 8，1 , k 9，1 corresponding to the cross-track pointing angle of the first remote sensing satellite the coefficient of s1, l1, s1l1, in the direction of the cross-track pointing angle of the first remote sensing satellite h 0，2 the pointing angle of the second remote sensing satellite in the along-track direction the constant term in the along-track direction; h 1，2 , h 2，2 , h 3，2 , h 4，2 , h 5，2 , h 6，2 , h 7，2 , h 8，2 , h 9，2 corresponding to the pointing angle of the second remote sensing satellite in the along-track direction the coefficient of s2, l2, s2l2, in the along-track direction; k 0，2 the pointing angle of the second remote sensing satellite in the cross-track direction the constant term in the cross-track direction; k 1，2 , k 2，2 , k 3，2 , k 4，2 , k 5，2 , k 6，2 , k 7，2 , k 8，2 , k 9，2 corresponding to the pointing angle of the second remote sensing satellite in the cross-track direction the coefficient of s2, l2, s2l2, in the cross-track direction.

5. The method according to claim 3, wherein: Definition u1, v1, w1 are the first row, second row and third row of the column vector obtained by matrix multiplication, respectively; the first scale factor λ1 is calculated according to the following ellipsoid equation: In the above formula, A = a e + h, B = b e + h, a e = 6378137.0 m, b e = 6356752.3 m, h is the ellipsoid height.

6. The method according to claim 3, wherein: Definition u2, v2, w2 are the first row, second row and third row of the column vector obtained by matrix multiplication, respectively; the second scale factor λ2 is calculated according to the following ellipsoid equation: In the above formula, A = a e + h, B = b e + h, a e = 6378137.0 m, b e = 6356752.3 m, h is the ellipsoid height.

7. The method of claim 1, wherein the method further comprises: The total number of images of the first remote sensing satellite and the second remote sensing satellite at previous time required by step III satisfies: M+N≥9.

8. The method according to claim 3, wherein: The motion trajectory sliding solution method of step III comprises the following steps: S1, definition: R1, R2 are orthogonal matrices, The motion trajectory fusion model of the first and second remote sensing satellites is deformed as follows: R is multiplied on both sides of the above two equations -1 By the four arithmetic operations, the following equation is obtained: In the above formula: define F x，1 , F x，2 are the pointing angle residual function models of the first remote sensing satellite and the second remote sensing satellite along the orbit direction respectively, F y，1 , F y，2 are the pointing angle residual function models of the first remote sensing satellite and the second remote sensing satellite along the orbit direction respectively; S2, for each target particle coordinate of the air target on the sequence images of the first and second remote sensing satellites, a standard error equation is constructed, which is as follows: V=AX-LP In the above formula: V is the residual vector of X, which is expressed as follows: is the result of the calculation on the right side of the standard error equation; A is a design matrix composed of unknown partial derivatives, which is expressed as follows: Wherein: subscript 1 and 2 respectively represent the first remote sensing satellite and the second remote sensing satellite, subscript i is the frame number of the first remote sensing satellite sequence image where the airborne moving target is located, i=1, 2, …, N, subscript j is the frame number of the second remote sensing satellite sequence image where the airborne moving target is located, j=1, 2, …, M; X is an unknown matrix, which is iteratively calculated according to the least square adjustment principle, and the calculation formula is as follows: X = (A T PA) -1 (A T PL) The representation form of X is as follows: X=[X0, δ1, δ2, Y0, β1, β2, Z0, θ1, θ2] When calculating the standard error equation for the first time, the initial value is assigned to X; L is a constant matrix obtained by substituting the current calculated value of the unknown matrix X into F x，1 , F x，2 , F y，1 , F y，2 , and the expression of the constant term matrix is as follows: L = [F x，1，i , F y，1，i , F x，2，j , F y，2，j ] T P is the observation weight of the target particle coordinates, which is a unit matrix; S3, according to the least square adjustment principle, the motion trajectory parameters X0, δ1, δ2, Y0, β1, β2, Z0, θ1, θ2 of the airborne target at the iteration stop are solved.

9. The method according to claim 8, wherein: The condition for judging the convergence of X iteration according to the least square adjustment principle is that the difference between X0, δ1, δ2, Y0, β1, β2, Z0, θ1, θ2 obtained by the current and previous calculation is within 1×10 -6 and the iteration is stopped; and the value of the unknown matrix X obtained by the last iteration is taken as the obtained air target motion trajectory fusion model parameter.

10. The method of claim 3, wherein the method further comprises: The method for solving the spatial position (X, Y, Z) and the velocity (V X , V Y , V Z ) of the aerial moving target at the current time in the fourth step is: according to the solved aerial target motion trajectory fusion model parameters X0, δ1, δ2, Y0, β1, β2, Z0, θ1, θ2, the following position equation and velocity equation are substituted to obtain:

Images