CN113962057A - Remote missile active section motion parameter correction method based on time sequence intersection - Google Patents

Remote missile active section motion parameter correction method based on time sequence intersection Download PDF

Info

Publication number
CN113962057A
CN113962057A CN202110727620.5A CN202110727620A CN113962057A CN 113962057 A CN113962057 A CN 113962057A CN 202110727620 A CN202110727620 A CN 202110727620A CN 113962057 A CN113962057 A CN 113962057A
Authority
CN
China
Prior art keywords
missile
ballistic
point
satellite
time
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110727620.5A
Other languages
Chinese (zh)
Other versions
CN113962057B (en
Inventor
盛庆红
任瑞
王博
吴凡
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University of Aeronautics and Astronautics
Original Assignee
Nanjing University of Aeronautics and Astronautics
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University of Aeronautics and Astronautics filed Critical Nanjing University of Aeronautics and Astronautics
Priority to CN202110727620.5A priority Critical patent/CN113962057B/en
Publication of CN113962057A publication Critical patent/CN113962057A/en
Application granted granted Critical
Publication of CN113962057B publication Critical patent/CN113962057B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/20—Analysis of motion
    • G06T7/277—Analysis of motion involving stochastic approaches, e.g. using Kalman filters
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/04—Constraint-based CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30241—Trajectory

Abstract

The invention provides a remote missile active section motion parameter correction method based on time sequence intersection. The method comprises the following steps: extracting a unit observation sight vector from the center of the satellite to the target ballistic missile in the direction unit under the earth center fixed coordinate system according to the conversion relation between the target point coordinate of the ballistic missile and the image point coordinate of the early warning satellite infrared sensor; analyzing the geometrical relationship existing among all ballistic planes, and cutting the extracted observation sight vector by using the ballistic planes to establish a ballistic plane cutting model; constructing an active section dynamic model based on momentum conservation, and screening alternative ballistic planes; designing a conversion matrix according to the image information, and converting all image planes to obtain a plurality of virtual image plane coordinates; performing Kalman prediction on image points on a plurality of virtual image planes; and correcting the constraint according to the observation sight line formed by the prediction result. The method can effectively realize the estimation of the active section parameters of the ballistic missile under the single-satellite observation background.

Description

Remote missile active section motion parameter correction method based on time sequence intersection
Technical Field
The invention relates to the technical field of photogrammetry, infrared image processing and trajectory prediction, in particular to a remote missile active section motion parameter correction method based on time sequence intersection.
Background
The space-based infrared early warning system mainly adopts a satellite-borne infrared sensor to realize the functions of monitoring, tracking, parameter estimation and the like of ballistic missile target launching. The infrared sensor obtains the angle observation quantity of the ballistic missile through passive detection, cannot measure the distance and the speed, and belongs to incomplete observation for missile positioning. In the development process of the high orbit space-based early warning system, due to the stage of early warning constellation deployment and verification and test in the early stage of system construction, early warning detection is in single-satellite detection in most cases.
The active section is the key of the estimation of the space motion parameters of the ballistic missile, and for the active parameter estimation model of the ballistic missile, the main research method is ballistic modeling based on a ballistic template library by using prior knowledge. Due to the limited types of Ballistic missiles, different types of target active segment Ballistic templates can be stored in a database in advance, called a standard Ballistic Profile (Nominal Ballistic Profile). The parameters of the alternative trajectory in the existing method are usually set according to experience or traversal, and the problems of difficulty in determining the alternative trajectory plane and low calculation efficiency exist.
Disclosure of Invention
The purpose of the invention is as follows: the invention provides a remote missile active section motion parameter correction method based on time sequence rendezvous, which solves the problem that the center of the prior art is difficult to determine an alternative trajectory plane by applying Kalman filtering and a time sequence rendezvous principle.
The technical scheme is as follows: the invention relates to a remote missile active section motion parameter correction method based on time sequence intersection, which comprises the following steps of:
(1) extracting a unit observation sight vector from the center of a satellite to a target ballistic missile in a direction unit under a geocentric fixed coordinate system according to the conversion relation between the target point coordinate of the ballistic missile and the early warning image point coordinate of an early warning satellite infrared sensor;
(2) analyzing the geometrical relationship existing among all ballistic planes, and cutting the extracted observation sight vector by using the ballistic planes to establish a ballistic plane cutting model;
(3) constructing an active section dynamic model based on momentum conservation according to the stress condition and the motion characteristic of the ballistic missile during active flight and the common characteristic parameters of the ballistic missile, and screening alternative ballistic planes;
(4) converting the image plane coordinate system according to the time sequence to obtain a plurality of virtual image plane coordinate systems;
(5) performing time sequence two-dimensional track prediction on target information of an image plane coordinate system by using Kalman filtering;
(6) and correcting the track according to the time sequence observation sight of the predicted two-dimensional track.
Further, the step (1) includes the steps of:
(11) based on a collinear condition equation, establishing a strict geometric imaging model for an original image of an early warning image point coordinate of an early warning satellite infrared sensor:
Figure BDA0003138093620000021
wherein (x y) is the image plane coordinate obtained by the imaging model, (x)0 y0) F is the focal length of the camera, (X) for early warning the coordinates of the principal point of the early warning image of the satellite infrared sensors Ys Zs) For the coordinates of the early warning satellite in the earth center fixed coordinate system, (X Y Z) is the coordinates of the ballistic missile in the earth center fixed coordinate system, ak,bk,ck(k is 1,2,3) is a function of the early warning satellite system parameters, and is specifically expressed as follows:
Figure BDA0003138093620000022
in the formula, RJ2000-WGS84A rotation matrix from a J2000 coordinate system to a WGS84 coordinate system; rorbit-J2000Rotation of an orbital measurement coordinate system to a J2000 coordinate system for early warning satellitesRotating the matrix; rbody-orbitA transformation matrix from the body coordinate system of the early warning satellite to the orbit measurement coordinate system; rcamera-bodyA transformation matrix from a sensor coordinate system to a body coordinate system;
(12) establishing a direction unit observation sight vector (uv w) from the center of an observation satellite to a target trajectory missileTMathematical functional relationship with the image plane coordinates (x, y) of the missile:
Figure BDA0003138093620000023
wherein m and lambda are proportionality coefficients; (x, y)TImage plane coordinates for ballistic missile imaging; (x)0,y0,-f)TThe image points are corresponding to the internal orientation elements during imaging; obtaining an observation sight line vector (uv w)TU, v, and w are observation visual line vectors in the X direction, the Y direction, and the Z direction, respectively.
Further, the step (2) comprises the steps of:
(21) introducing coordinates of a first starting point and an end point of the early warning satellite in a geocentric fixed coordinate system, wherein the coordinates are respectively (X)1 Y1 Z1)TAnd (X)N YN ZN)TThe first starting point is the position of the early warning satellite for detecting and finding the ballistic missile for the first time, and the end point is the position of the last detected and found ballistic missile, and the method comprises the following steps:
Figure BDA0003138093620000031
wherein (u)1 v1 w1)TAs the observation sight line vector corresponding to the initial point, (u)N vN wN)TIs the observation sight line vector corresponding to the terminal point, (X)S YS ZS)TFor early warning of the coordinates of the satellite in the earth-centered fixed coordinate system, reIs the average earth radius, h1And hNAssuming observation of the satellite, respectively initial and final elevationDetecting tail flame of the missile engine with the active section trajectory when the satellite-borne infrared sensor monitors, wherein N groups of observed values are provided and correspond to the observation time t respectively1,t2,……,tN；
(22) Solving for (X) according to the visibility constraint in the basic ballistic constraint1 Y1 Z1)TAnd (X)N YN ZN)TThe perspective constraint means that the observation sight line vector can not pass through the earth before the ballistic missile when the ballistic plane is cut, and then a ballistic cutting plane passing through the center of the earth is determined.
Further, the step (3) includes the steps of:
(31) neglecting the influence of air action, Coriolis inertia force and a traction inertia force, and constructing an acceleration model of the ballistic missile according to the impulse of the combined external force borne by the ballistic missile after the momentum change of the ballistic missile in the active section:
Figure BDA0003138093620000032
wherein p '(t) is the acceleration of the missile at the time t, M (t) is the missile mass at the time t, M' (t) is the change rate of the missile mass at the time t, and g0The gravity acceleration is adopted, and u is the relative effective injection speed of the high-temperature gas of the missile;
(32) introducing the second consumption a of the rocket engine on the basis of the step (31), and establishing a dynamic model of the active section of the ballistic missile by taking s as a/M:
Figure BDA0003138093620000033
where p (t) denotes the position of the missile at time t, p0,v0Initial conditions representing position and velocity, respectively; suppose the high-temperature fuel gas has a constant effective injection speed u relative to the missilec，gcThe subscript c represents an assumed constant value as the gravitational acceleration constant;
(33) and screening the alternative trajectory planes by using an active section dynamic model of the trajectory missile, and removing trajectories which do not accord with the dynamic model to obtain trajectory objects which accord with the characteristics of the missile.
Further, in the step (4), the coordinate system is converted by using the following conversion matrix:
Figure BDA0003138093620000041
wherein the parameters are respectively as follows:
Figure BDA0003138093620000042
Figure BDA0003138093620000043
Figure BDA0003138093620000044
b1＝cosωsinκ
b2＝cosωcosκ
b3＝-sinω
Figure BDA0003138093620000045
Figure BDA0003138093620000046
Figure BDA0003138093620000047
wherein
Figure BDA0003138093620000048
ω and κ are three attitude angles of the satellite: pitch angle, roll angle andand (4) yaw angle.
Further, the step (5) includes the steps of:
(51) performing Kalman prediction on the virtual image plane obtained in the step (4), wherein a Kalman filtering recursion formula is as follows:
equation of state prediction
Figure BDA0003138093620000049
Equation of state measurement
Figure BDA00031380936200000410
Mean square error prediction equation
Figure BDA00031380936200000411
Kalman gain equation
Figure BDA00031380936200000412
Equation of state modification
Figure BDA00031380936200000413
Mean square error modification equation
Figure BDA00031380936200000414
Where k represents a time point, matrix
Figure BDA0003138093620000051
State, matrix, representing k time points without considering estimation errors
Figure BDA0003138093620000052
Representing the state at k time points, matrix zkRepresenting data measured at k time points, matrix AkA state transition matrix representing states from a point in time k-1 to a point in time k
Figure BDA0003138093620000053
Mean square error matrix, matrix H, representing k time pointskRepresenting a Kalman gain matrix;
(52) defining Kalman filter tracker system states
Figure BDA0003138093620000054
Is a four-dimensional vector (P)x,Vx,Py,Vy)T，Px，Vx，Py，VyThe position and the speed of the target in the X-axis direction and the Y-axis direction under the image plane coordinate system respectively; defining an observation vector zk＝(Px,Py)T；
Defining a state transition matrix Ak：
Figure BDA0003138093620000055
Wherein Δ T is a satellite camera shooting time interval;
the observation matrix C is known from the relationship between the system state and the observation statekComprises the following steps:
Figure BDA0003138093620000056
let recurrence formula:
Figure BDA0003138093620000057
wherein delta is 0.01, r is 100;
kalman filtering is carried out according to the final X after multiple times of calculation of all image plane coordinates in a certain virtual image planekA predicted image plane coordinate point is obtained.
Further, the step (6) comprises the steps of:
(61) the prediction point calculated by the Kalman filtering is (x)0,y0,z0) The actual point observed is (x'0,y′0,z′0) Setting the observation vector as a parameter equation:
Figure BDA0003138093620000061
wherein i1、j1、k1、i2、j2、k2The vector parameters of the satellite and the missile are obtained by the coordinates of the satellite and the missile;
(62) according to the forward intersection principle, the intersection point of the predicted time sequence observation vector and the real time sequence observation vector is the missile position, and in order to accurately calculate the intersection point of the two spatial straight lines, the intersection point of a common perpendicular line of the two straight lines and the two straight lines is set as A (i)1m+x0,j1m+y0,k1m+z0)，B(i2n+x′0,j2n+y′0,k2n+z′0) Wherein m and n are parameters to be solved, and according to the properties of the plumb line, the following are obtained:
Figure BDA0003138093620000062
(63) will i1、j1、k1、i2、j2、k2Predicted point (x)0,y0,z0) And actual point (x'0,y′0,z′0) Substituting into the calculation formula in the step (62) to obtain the values of the parameters m and n; substituting m and n into the calculation formula in the step (61) to obtain A, B two-point three-dimensional information by the formula;
(64) calculating the midpoint C of the intersection point A, B, where C is a coordinate set C at all time pointsiConstraining the rear track with CiAnd modifying the parameters after the comparison to realize the track modification.
Has the advantages that:
1. according to the method, a missile movement physical model is established according to the geometric relation between the cutting trajectory plane and the observation sight and the momentum theorem, a strongly-constrained trajectory missile active section parameter estimation model is established, the geometric change rule of the early warning satellite in the process of observing a space target is analyzed, the dependence of the parameter estimation model on prior information can be eliminated, and a plurality of accurate alternative trajectory planes can be calculated only by depending on image plane coordinates and satellite information.
2. The method provided by the invention corrects the trajectory track by applying Kalman filtering and time sequence intersection principles and using a Kalman filtering spatial track correction model, and breaks through the problem of inaccurate parameter estimation caused by lack of prior information support.
Drawings
FIG. 1 is a flow chart of a remote missile active section motion parameter correction method based on time sequence intersection.
Detailed Description
The present invention is described in further detail below with reference to the attached drawing figures. The invention provides a remote missile active section motion parameter correction method based on time sequence intersection, which comprises the following steps as shown in figure 1:
the method comprises the following steps: and extracting a unit observation sight vector from the satellite center to the target ballistic missile in the direction under the earth center fixed coordinate system according to the conversion relation between the target point coordinate of the ballistic missile and the early warning image point coordinate of the early warning satellite infrared sensor. The method specifically comprises the following steps:
(11) based on a collinear condition equation, establishing a strict geometric imaging model for an original image of an early warning image point coordinate of an early warning satellite infrared sensor:
Figure BDA0003138093620000071
wherein (x y) is the image plane coordinate obtained by the imaging model, (x)0 y0) Is the coordinate of the principal point of the image, f is the focal length of the camera, (X)s Ys Zs) The method comprises the following steps of (1) pre-warning coordinates of a satellite center under a geocentric fixed coordinate system, wherein the geocentric coordinates are obtained by measurement of an orbit measurement device, and the orbit model corresponding to a platform relates to aspects such as orbit dynamics and satellite perturbation; (X Y Z) is the coordinate of the ballistic missile in the earth center fixed coordinate system, ak,bk,ck(k＝1,2,3)The function of the early warning satellite system parameters is specifically expressed as follows:
Figure BDA0003138093620000072
in the formula, RJ2000-WGS84A rotation matrix from a J2000 coordinate system to a WGS84 coordinate system; rorbit-J2000Acquiring a rotation matrix from a coordinate system to a J2000 coordinate system for the orbit measurement of the early warning satellite by attitude measurement equipment; rbody-orbitA transformation matrix from the body coordinate system of the early warning satellite to the orbit measurement coordinate system; rcamera-bodyA transformation matrix from a sensor coordinate system to a body coordinate system, namely a camera installation matrix;
(12) establishing a direction unit observation sight vector (uv w) from the center of an observation satellite to a target trajectory missileTMathematical functional relationship between (i.e. direction finding information) and image plane coordinates (x, y) of the missile:
Figure BDA0003138093620000081
wherein m and lambda are proportional coefficients which are equivalent to scaling coefficients; (x, y)TImage plane coordinates for ballistic missile imaging; (x)0,y0,-f)TThe image point corresponds to the internal orientation element when imaging. Obtaining the observation sight line vector (uv w)TWherein u, v and w are observation sight line vectors in the X direction, the Y direction and the Z direction respectively.
Step two: and analyzing the geometrical relationship existing among all the ballistic planes, and cutting the extracted observation sight line vector by using the ballistic planes to establish a ballistic plane cutting model. The method specifically comprises the following steps:
(21) introducing a GEO early warning satellite to detect the coordinates of the ballistic missile (namely the initial point) for the first time and the ballistic missile (namely the terminal point) for the last time under the geocentric fixed coordinate system, wherein the coordinates are respectively (X)1 Y1 Z1)TAnd (X)N YN ZN)TThen, there are:
Figure BDA0003138093620000082
wherein (u)1 v1 w1)TAs the observation sight line vector corresponding to the initial point, (u)N vN wN)TIs the observation sight line vector corresponding to the terminal point, (X)S YS ZS)TCoordinates r of the GEO early warning satellite in the geocentric fixed coordinate systemeIs the average earth radius, h1And hNRespectively as a starting point and an end point, detecting the tail flame of the active section trajectory missile engine when the satellite-borne infrared sensor of the GEO observation satellite is assumed to monitor, wherein N groups of observed values are provided and respectively correspond to the observation time t1,t2,……,tN；
(22) According to the visibility constraint in the basic trajectory constraint, namely that when the trajectory plane is cut, the observation sight vector cannot pass through the earth before the trajectory missile, the (X) can be calculated1 Y1 Z1)TAnd (X)N YN ZN)TAnd a through-the-earth ballistic cutting plane can be determined.
The rapid cutting model method based on the trajectory plane traversal cutting assumes that an observed trajectory missile is not maneuvered, and the active section of the trajectory missile always moves in a trajectory plane which is approximately over the center of the earth sphere, namely, the plane assumption is met. Cutting the viewing line of sight vector (u) with these ballistic planesi vi wi)TN, all the cut curves are likely to be ballistic curves, i.e. alternative ballistic curves.
(a) And determining the interval range of the first starting point elevation and the end point elevation. Elevation interval of first issue point [ h ]1min，h1max]The early warning image data can be obtained according to the data of the early warning image of the initial point and the local meteorological data; the lower limit boundary of the terminal elevation interval can be selected to be the same as the upper limit boundary of the initial elevation interval, namely hNmin＝h1maxWhose upper bound can be hNmax＝-g0·T2/2+T·umaxThus obtaining the product. Wherein T is TN-t1To warn of the total observed time length, g, of the satellite from the first origin to the end point0As acceleration of gravity, umaxThe maximum gas spraying speed value is obtained.
(b) Set up h1And hNAnd traversing all possible candidate ballistic cutting planes, and solving candidate ballistic trajectory data in each cutting plane.
Step three: and (4) analyzing the movement characteristics of the ballistic missiles to screen the trajectories of the ballistic missiles.
(31) Neglecting the influence of air action, Coriolis inertia force and a traction inertia force, and constructing an acceleration model of the ballistic missile according to the impulse of the combined external force borne by the ballistic missile after the momentum change of the ballistic missile in the active section:
Figure BDA0003138093620000091
wherein, p '(t) is the acceleration of the missile at the moment t, M (t) is the missile mass at the moment t, M' (t) is the change rate of the missile mass at the moment t, and the change rate refers to the ratio of the change of the missile mass to the time infinitesimal after the time infinitesimal passes at the moment t. g0And u is the relative effective jet velocity of the high-temperature gas of the missile.
(32) Introducing the second consumption a of the rocket engine on the basis of the step (31), and establishing a dynamic model of the active section of the ballistic missile by taking s as a/M:
Figure BDA0003138093620000092
where p (t) denotes the position of the missile at time t, i.e. the ballistic model, p0,v0Initial conditions representing position and velocity, respectively; suppose the high-temperature fuel gas has a constant effective injection speed u relative to the missilec，gcThe subscript c represents an assumed constant value as the gravitational acceleration constant; in the active section and flight time section of the ballistic missile, the change of the space position of the ballistic missile can be ignoredThe gravity acceleration can be ignored, and can be processed as constant, and recorded as gc。
(3) And strictly constraining the alternative trajectory plane by using a dynamical model of the ballistic missile, and eliminating trajectories which do not accord with the dynamical model of the ballistic missile to obtain trajectory objects which accord with the characteristics of the missile.
Step four: and designing a conversion matrix according to the image information, and converting all the image planes to obtain a plurality of virtual image plane coordinates.
And converting the multi-time-point image plane coordinate system to obtain a plurality of virtual image plane coordinate systems. The specific transformation matrix is:
Figure BDA0003138093620000101
wherein the parameters are respectively as follows:
Figure BDA0003138093620000102
Figure BDA0003138093620000103
Figure BDA0003138093620000104
b1＝cosωsinκ
b2＝cosωcosκ
b3＝-sinω
Figure BDA0003138093620000105
Figure BDA0003138093620000106
Figure BDA0003138093620000107
wherein
Figure BDA0003138093620000108
ω and κ are three attitude angles of the satellite: pitch angle, roll angle and yaw angle, three rotation angles can be modified according to actual need.
Modifying according to the actual virtual image plane requirement
Figure BDA0003138093620000109
The values of omega and kappa can make the image plane coordinate [ x ] in the step one0 y0 -f]TConverted into virtual image plane coordinates [ x 'y' -f]T：
Figure BDA00031380936200001010
Each virtual image plane has the same number of coordinate points, which is about 10-20, depending on the actual situation.
Step five: and performing two-dimensional track prediction on target information of the image plane coordinate system by using Kalman filtering. The method specifically comprises the following steps:
(51) and performing Kalman prediction on the virtual image plane obtained in the step four. Kalman filtering is an algorithm for performing optimal estimation on the system state by using a linear system state equation and inputting and outputting observation data through a system, is the most widely applied filtering method at present, and is better applied to multiple fields such as communication, navigation, guidance and control. The image on the focal plane may be approximated as a point, the position of which is the system input to the kalman system. By means of kalman filtering, system noise can be suppressed and an optimal position estimate is obtained.
The Kalman filtering does not store the past data, and after new data are obtained each time, the new parameter estimation value can be calculated according to the new data and the parameter estimation value at the previous moment and a set of deduced recursion formula by means of the state transition equation of the system. The storage and calculation amount of the Kalman filtering device are greatly reduced due to the characteristic of using algorithm, and the limit of a stable random process is broken through.
The n-dimensional state equation and the m-dimensional measurement equation of a discrete dynamic system are respectively:
xk＝Akxk-1+wk-1 (1.1)
yk＝Ckxk+vk (1.2)
in the formulae (1.1) and (1.2), the matrix AkA state transition matrix representing the state from the point of time k-1 to the point of time k, CkCalled observation matrix, ykIs measured data, wk-1And vkAre noise error matrices. A. thek，Ck，ykIs known, using xk-1Is estimated value of
Figure BDA0003138093620000111
Determining an estimated value for a k time point
Figure BDA0003138093620000112
If there is no wk-1And vkThen, x can be immediately obtained from the formula (1.1) and the formula (1.2)kI.e. there is no estimation problem. The estimation problem arises because the signal and noise are added together so that the true value of the signal is estimated. Temporarily disregarding wk-1And vxIn this case x is obtained according to formulae (1.1) and (1.2)kAnd ykAre used separately
Figure BDA0003138093620000113
And zkThis means that there are:
Figure BDA0003138093620000114
Figure BDA0003138093620000115
will zkAnd ykAre compared, and the difference is used
Figure BDA0003138093620000116
It shows, as follows:
Figure BDA0003138093620000117
zkdeviation ykTo produce
Figure BDA0003138093620000118
Apparently due to neglect of wk-1And vkCaused, that is to say
Figure BDA0003138093620000119
Implies wk-1And vkOf (2), in other words
Figure BDA00031380936200001110
Implying the current observed value ykThe information of (1).
If it will be
Figure BDA00031380936200001111
Multiplied by a gain HkTo correct the original
Figure BDA00031380936200001112
The values, a better estimate is obtained:
Figure BDA00031380936200001113
it can therefore be said that
Figure BDA00031380936200001114
And true value xkThe mean square error of (a) is an error square matrix. If H under the minimum condition of the error can be obtainedkThen this H is addedkSubstitution of formula (1.6) to give
Figure BDA00031380936200001115
Is to xkLinear optimal estimation of (2).
The kalman filter recursion formula is as follows:
equation of state prediction
Figure BDA0003138093620000121
Equation of state measurement
Figure BDA0003138093620000122
Mean square error prediction equation
Figure BDA0003138093620000123
Kalman gain equation
Figure BDA0003138093620000124
Equation of state modification
Figure BDA0003138093620000125
Mean square error modification equation
Figure BDA0003138093620000126
Where k represents a time point, the matrix
Figure BDA0003138093620000127
State, matrix, representing k time points without considering estimation errors
Figure BDA0003138093620000128
Representing the state at k time points, matrix zkRepresenting data measured at k time points, matrix AkTo representState transition matrix from k-1 time point to k time point, matrix
Figure BDA0003138093620000129
Mean square error matrix, matrix H, representing k time pointskRepresenting a kalman gain matrix.
(52) Defining Kalman filter tracker system states
Figure BDA00031380936200001210
Is a four-dimensional vector (P)x,Vx,Py,Vy)T。Px，Vx，Py，VyThe position and velocity of the target in the X-axis and Y-axis directions, respectively, in the image plane coordinate system. Since only the two-dimensional position of the target can be determined on the image, the observation vector z is defined herek＝(Px,Py)T。
Since the motion of the object of interest within a unit time interval is regarded as uniform motion, the state transition matrix AkIs defined as:
Figure BDA00031380936200001211
where Δ T is the satellite camera capture interval.
The observation matrix C is known from the relationship between the system state and the observation statekComprises the following steps:
Figure BDA00031380936200001212
further, assume that in the recurrence formula:
Figure BDA0003138093620000131
wherein delta is 0.01 and r is 100.
Kalman filtering is based on the coordinates of all image planes in a virtual image planeAfter multiple calculations, the final X can be obtainedkTo obtain a predicted image plane coordinate point.
Step six: and correcting the track according to the observation sight of the predicted two-dimensional track. The method specifically comprises the following steps:
(61) assume that the prediction point calculated from the time-series Kalman filtering is (x)0,y0,z0) The actual point observed is (x'0,y′0,z′0). The observation vector can be set to the parametric equation:
Figure BDA0003138093620000132
wherein i1、j1、k1、i2、j2、k2The vector parameter of the satellite-missile can be obtained by the coordinates of the satellite and the missile.
(62) According to the forward intersection principle, the intersection point of the predicted time sequence observation vector and the real time sequence observation vector is the missile position. In order to accurately calculate the intersection point of two spatial straight lines, the intersection point of the common perpendicular line of the two straight lines and the two straight lines is set as A (i)1m+x0,j1m+y0,k1m+z0)，B(i2n+x′0,j2n+y′0,k2n+z′0). Wherein m and n are parameters to be solved. According to the properties of the plumb line, the following can be obtained:
Figure BDA0003138093620000133
(63) will i1、j1、k1、i2、j2、k2Predicted point (x)0,y0,z0) And actual point (x'0,y′0,z′0) The values of the parameters m and n can be obtained by substituting the formula (11). Two-point three-dimensional information A, B can be obtained by substituting m and n into the formula (10).
(64) The midpoint C, of intersection A, B is calculated as (A + B)/2. All deficiencyAll the sets of C obtained by pseudo-image plane calculation are Ci. From the method of determining a plane at three points, from CiOptionally three coordinates in the set define a plane D, the set of all different planes D being Di。
Plane D is defined as:
CS＝OX+PY+LZ (12)
wherein O, P, L, CS is a plane parameter.
If step three alternative ballistic planes lack plane set DiThe missing plane is added to the alternative ballistic plane.
The method researches a trajectory plane cutting prediction model under a single star visual axis, analyzes the motion characteristics of a trajectory active section, constructs a constraint plane by using Kalman filtering, determines the trajectory space trajectory, and can solve the problem that the center of the prior art is difficult to determine an alternative trajectory plane.
The foregoing is only a preferred embodiment of this invention and it should be noted that modifications can be made by those skilled in the art without departing from the principle of the invention and these modifications should also be considered as the protection scope of the invention.

Claims (8)

1. A remote missile active section motion parameter correction method based on time sequence intersection is characterized by comprising the following steps:
(1) extracting a unit observation sight vector from the center of a satellite to a target ballistic missile in a direction unit under a geocentric fixed coordinate system according to the conversion relation between the target point coordinate of the ballistic missile and the early warning image point coordinate of an early warning satellite infrared sensor;
(2) analyzing the geometrical relationship existing among all ballistic planes, and cutting the extracted observation sight vector by using the ballistic planes to establish a ballistic plane cutting model;
(3) constructing an active section dynamic model based on momentum conservation according to the stress condition and the motion characteristic of the ballistic missile during active flight and the common characteristic parameters of the ballistic missile, and screening alternative ballistic planes;
(4) converting the image plane coordinate system according to the time sequence to obtain a plurality of virtual image plane coordinate systems;
(5) performing time sequence two-dimensional track prediction on target information of an image plane coordinate system by using Kalman filtering;
(6) and correcting the track according to the time sequence observation sight of the predicted two-dimensional track.
2. The remote missile active segment motion parameter modification method based on time-series rendezvous as claimed in claim 1, wherein the step (1) comprises the following steps:
(11) based on a collinear condition equation, establishing a strict geometric imaging model for an original image of an early warning image point coordinate of an early warning satellite infrared sensor:
Figure FDA0003138093610000011
wherein (x y) is the image plane coordinate obtained by the imaging model, (x)0 y0) F is the focal length of the camera, (X) for early warning the coordinates of the principal point of the early warning image of the satellite infrared sensors Ys Zs) For the coordinates of the early warning satellite in the earth center fixed coordinate system, (X Y Z) is the coordinates of the ballistic missile in the earth center fixed coordinate system, ak,bk,ck(k is 1,2,3) is a function of the early warning satellite system parameters, and is specifically expressed as follows:
Figure FDA0003138093610000012
in the formula, RJ2000-WGS84A rotation matrix from a J2000 coordinate system to a WGS84 coordinate system; rorbit-J2000Measuring a rotation matrix from a coordinate system to a J2000 coordinate system for the orbit of the early warning satellite; rbody-orbitA transformation matrix from the body coordinate system of the early warning satellite to the orbit measurement coordinate system; rcamera-bodyA transformation matrix from a sensor coordinate system to a body coordinate system;
(12) establishing an observation satellite center-to-target missileDirection unit observation sight vector of road missile (uv w)TMathematical functional relationship with the image plane coordinates (x, y) of the missile:
Figure FDA0003138093610000021
wherein m and lambda are proportionality coefficients; (x, y)TImage plane coordinates for ballistic missile imaging; (x)0,y0,-f)TThe image points are corresponding to the internal orientation elements during imaging; obtaining an observation sight line vector (uv w)TU, v, and w are observation visual line vectors in the X direction, the Y direction, and the Z direction, respectively.
3. The remote missile active segment motion parameter modification method based on time-series rendezvous as claimed in claim 1, wherein the step (2) comprises the following steps:
(21) introducing coordinates of a first starting point and an end point of the early warning satellite in a geocentric fixed coordinate system, wherein the coordinates are respectively (X)1 Y1 Z1)TAnd (X)NYN ZN)TThe first starting point is the position of the early warning satellite for detecting and finding the ballistic missile for the first time, and the end point is the position of the last detected and found ballistic missile, and the method comprises the following steps:
Figure FDA0003138093610000022
wherein (u)1 v1 w1)TAs the observation sight line vector corresponding to the initial point, (u)N vN wN)TIs the observation sight line vector corresponding to the terminal point, (X)S YS ZS)TFor early warning of the coordinates of the satellite in the earth-centered fixed coordinate system, reIs the average earth radius, h1And hNRespectively as a starting point and an end point elevation, detecting the tail flame of the missile engine at the active section trajectory on the assumption that a satellite-borne infrared sensor of an observation satellite monitors,a total of N groups of observed values respectively corresponding to the observation time t1,t2,……,tN；
(22) Solving for (X) according to the visibility constraint in the basic ballistic constraint1 Y1 Z1)TAnd (X)N YN ZN)TThe perspective constraint means that the observation sight line vector can not pass through the earth before the ballistic missile when the ballistic plane is cut, and then a ballistic cutting plane passing through the center of the earth is determined.
4. The remote missile active segment motion parameter modification method based on time-series rendezvous as claimed in claim 3, wherein the step (22) comprises the following steps:
(a) determining the interval range of the first-transmitting-point elevation and the terminal elevation, wherein the first-transmitting-point elevation interval [ h1min，h1max]Obtaining the early warning image data according to the first sending point and the local meteorological data; the lower limit boundary of the terminal elevation interval is selected to be the same as the upper limit boundary of the first-generation elevation interval, namely hNmin＝h1maxIts upper bound is hNmax＝-g0·T2/2+T·umaxObtained, wherein T is TN-t1To warn of the total observed time length, g, of the satellite from the first origin to the end point0As acceleration of gravity, umaxThe maximum gas spraying speed value is obtained;
(b) set up h1And hNAnd traversing all possible candidate ballistic cutting planes, and solving candidate ballistic trajectory data in each cutting plane.
5. The remote missile active segment motion parameter modification method based on time-series rendezvous as claimed in claim 1, wherein the step (3) comprises the following steps:
(31) neglecting the influence of air action, Coriolis inertia force and a traction inertia force, and constructing an acceleration model of the ballistic missile according to the impulse of the combined external force borne by the ballistic missile after the momentum change of the ballistic missile in the active section:
Figure FDA0003138093610000031
wherein p '(t) is the acceleration of the missile at the time t, M (t) is the missile mass at the time t, M' (t) is the change rate of the missile mass at the time t, and g0The gravity acceleration is adopted, and u is the relative effective injection speed of the high-temperature gas of the missile;
(32) introducing the second consumption a of the rocket engine on the basis of the step (31), and establishing a dynamic model of the active section of the ballistic missile by taking s as a/M:
Figure FDA0003138093610000032
where p (t) denotes the position of the missile at time t, p0,v0Initial conditions representing position and velocity, respectively; suppose the high-temperature fuel gas has a constant effective injection speed u relative to the missilec，gcThe subscript c represents an assumed constant value as the gravitational acceleration constant;
(33) and screening the alternative trajectory planes by using an active section dynamic model of the trajectory missile, and removing trajectories which do not accord with the dynamic model to obtain trajectory objects which accord with the characteristics of the missile.
6. The method for modifying motion parameters of an active segment of a remote missile based on time-series rendezvous as claimed in claim 1, wherein the step (4) is implemented by converting a coordinate system by using a conversion matrix as follows:
Figure FDA0003138093610000041
wherein the parameters are respectively as follows:
Figure FDA0003138093610000042
Figure FDA0003138093610000043
Figure FDA0003138093610000044
b1＝cosωsinκ
b2＝cosωcosκ
b3＝-sinω
Figure FDA0003138093610000045
Figure FDA0003138093610000046
Figure FDA0003138093610000047
wherein
Figure FDA0003138093610000048
ω and κ are three attitude angles of the satellite: pitch angle, roll angle, and yaw angle.
7. The remote missile active segment motion parameter modification method based on time-series rendezvous as claimed in claim 1, wherein the step (5) comprises the following steps:
(51) performing Kalman prediction on the virtual image plane obtained in the step (4), wherein a Kalman filtering recursion formula is as follows:
equation of state prediction
Figure FDA0003138093610000049
Equation of state measurement
Figure FDA00031380936100000410
Mean square error prediction equation
Figure FDA00031380936100000411
Kalman gain equation
Figure FDA00031380936100000412
Equation of state modification
Figure FDA00031380936100000413
Mean square error modification equation
Figure FDA00031380936100000414
Where k represents a time point, matrix
Figure FDA00031380936100000415
State, matrix, representing k time points without considering estimation errors
Figure FDA00031380936100000416
Representing the state at k time points, matrix zkRepresenting data measured at k time points, matrix AkA state transition matrix representing states from a point in time k-1 to a point in time k
Figure FDA0003138093610000051
Mean square error matrix, matrix H, representing k time pointskRepresenting a Kalman gain matrix;
(52) defining Kalman filter tracker system states
Figure FDA0003138093610000052
Is a four-dimensional vector (P)x,Vx,Py,Vy)T，Px，Vx，Py，VyThe position and the speed of the target in the X-axis direction and the Y-axis direction under the image plane coordinate system respectively; defining an observation vector zk＝(Px,Py)T；
Defining a state transition matrix Ak：
Figure FDA0003138093610000053
Wherein Δ T is a satellite camera shooting time interval;
the observation matrix C is known from the relationship between the system state and the observation statekComprises the following steps:
Figure FDA0003138093610000054
let recurrence formula:
Figure FDA0003138093610000055
wherein delta is 0.01, r is 100;
kalman filtering is carried out according to the final X after multiple times of calculation of all image plane coordinates in a certain virtual image planekA predicted image plane coordinate point is obtained.
8. The remote missile active segment motion parameter modification method based on time-series rendezvous as claimed in claim 1, wherein the step (6) comprises the following steps:
(61) the prediction point calculated by the Kalman filtering is (x)0,y0,z0) The actual point observed is (x'0,y′0,z′0) Setting the observation vector as a parameter equation:
Figure FDA0003138093610000061
wherein i1、j1、k1、i2、j2、k2The vector parameters of the satellite and the missile are obtained by the coordinates of the satellite and the missile;
(62) according to the forward intersection principle, the intersection point of the predicted time sequence observation vector and the real time sequence observation vector is the missile position, and in order to accurately calculate the intersection point of the two spatial straight lines, the intersection point of a common perpendicular line of the two straight lines and the two straight lines is set as A (i)1m+x0,j1m+y0,k1m+z0)，B(i2n+x′0,j2n+y′0,k2n+z′0) Wherein m and n are parameters to be solved, and according to the properties of the plumb line, the following are obtained:
Figure FDA0003138093610000062
(63) will i1、j1、k1、i2、j2、k2Predicted point (x)0,y0,z0) And actual point (x'0,y′0,z′0) Substituting into the calculation formula in the step (62) to obtain the values of the parameters m and n; substituting m and n into the calculation formula in the step (61) to obtain A, B two-point three-dimensional information by the formula;
(64) calculating the midpoint C of the intersection point A, B, where C is a coordinate set C at all time pointsiConstraining the rear track with CiAnd modifying the parameters after the comparison to realize the track modification.
CN202110727620.5A 2021-06-29 2021-06-29 Remote missile active section motion parameter correction method based on time sequence intersection Active CN113962057B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110727620.5A CN113962057B (en) 2021-06-29 2021-06-29 Remote missile active section motion parameter correction method based on time sequence intersection

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110727620.5A CN113962057B (en) 2021-06-29 2021-06-29 Remote missile active section motion parameter correction method based on time sequence intersection

Publications (2)

Publication Number Publication Date
CN113962057A true CN113962057A (en) 2022-01-21
CN113962057B CN113962057B (en) 2022-06-24

Family

ID=79460314

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110727620.5A Active CN113962057B (en) 2021-06-29 2021-06-29 Remote missile active section motion parameter correction method based on time sequence intersection

Country Status (1)

Country Link
CN (1) CN113962057B (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117052542A (en) * 2023-10-13 2023-11-14 太仓点石航空动力有限公司 Propulsion control optimizing method and system for aeroengine

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1596361A (en) * 2001-11-28 2005-03-16 福图尔特克股份公司 Projectile having a high penetrating action and lateral action and equipped with an integrated fracturing device
US20120316819A1 (en) * 2008-04-22 2012-12-13 United States Government, As Represented By The Secretary Of The Navy Process for estimation of ballistic missile boost state
CN105022035A (en) * 2015-07-31 2015-11-04 中国电子科技集团公司第三十八研究所 Trajectory target launch point estimate apparatus based on model updating and method
CN106643297A (en) * 2016-12-23 2017-05-10 南京长峰航天电子科技有限公司 Estimation and correction method for vector miss distance parameters of motion platform
CN107607947A (en) * 2017-08-22 2018-01-19 西安电子科技大学 Spaceborne radar imaging parameters On-line Estimation method based on Kalman filtering
CN109933847A (en) * 2019-01-30 2019-06-25 中国人民解放军战略支援部队信息工程大学 A kind of improved boost phase trajectory algorithm for estimating
CN111462182A (en) * 2020-03-31 2020-07-28 南京航空航天大学 Trajectory missile three-dimensional trajectory estimation method based on infrared early warning image

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1596361A (en) * 2001-11-28 2005-03-16 福图尔特克股份公司 Projectile having a high penetrating action and lateral action and equipped with an integrated fracturing device
US20120316819A1 (en) * 2008-04-22 2012-12-13 United States Government, As Represented By The Secretary Of The Navy Process for estimation of ballistic missile boost state
CN105022035A (en) * 2015-07-31 2015-11-04 中国电子科技集团公司第三十八研究所 Trajectory target launch point estimate apparatus based on model updating and method
CN106643297A (en) * 2016-12-23 2017-05-10 南京长峰航天电子科技有限公司 Estimation and correction method for vector miss distance parameters of motion platform
CN107607947A (en) * 2017-08-22 2018-01-19 西安电子科技大学 Spaceborne radar imaging parameters On-line Estimation method based on Kalman filtering
CN109933847A (en) * 2019-01-30 2019-06-25 中国人民解放军战略支援部队信息工程大学 A kind of improved boost phase trajectory algorithm for estimating
CN111462182A (en) * 2020-03-31 2020-07-28 南京航空航天大学 Trajectory missile three-dimensional trajectory estimation method based on infrared early warning image

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
YICONG LI ETC.: ""Trajectory and Launch Point Estimation for Ballistic Missiles from Boost Phase LOS Measurements"", 《1999 IEEE AEROSPACE CONFERENCE. PROCEEDINGS (CAT. NO.99TH8403)》 *
YICONG LI ETC.: ""Trajectory and Launch Point Estimation for Ballistic Missiles from Boost Phase LOS Measurements"", 《1999 IEEE AEROSPACE CONFERENCE. PROCEEDINGS (CAT. NO.99TH8403)》, 6 August 2002 (2002-08-06) *
申镇 等: ""单星无源探测弹道导弹射向估计新方法"", 《宇航学报》 *
申镇 等: ""单星无源探测弹道导弹射向估计新方法"", 《宇航学报》, 31 July 2011 (2011-07-31) *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117052542A (en) * 2023-10-13 2023-11-14 太仓点石航空动力有限公司 Propulsion control optimizing method and system for aeroengine
CN117052542B (en) * 2023-10-13 2023-12-08 太仓点石航空动力有限公司 Propulsion control optimizing method and system for aeroengine

Also Published As

Publication number Publication date
CN113962057B (en) 2022-06-24

Similar Documents

Publication Publication Date Title
US9073648B2 (en) Star tracker rate estimation with kalman filter enhancement
CN110081881B (en) Carrier landing guiding method based on unmanned aerial vehicle multi-sensor information fusion technology
CN107741229B (en) Photoelectric/radar/inertia combined carrier-based aircraft landing guiding method
US6639553B2 (en) Passive/ranging/tracking processing method for collision avoidance guidance
Mammarella et al. Machine vision/GPS integration using EKF for the UAV aerial refueling problem
CN109059907B (en) Trajectory data processing method and device, computer equipment and storage medium
US20080195316A1 (en) System and method for motion estimation using vision sensors
Zhao et al. Vision-aided estimation of attitude, velocity, and inertial measurement bias for UAV stabilization
CN114216454B (en) Unmanned aerial vehicle autonomous navigation positioning method based on heterogeneous image matching in GPS refusing environment
Gur fil et al. Partial aircraft state estimation from visual motion using the subspace constraints approach
Johnson et al. Combining stereo vision and inertial navigation for automated aerial refueling
CN109211231B (en) Cannonball attitude estimation method based on Newton iteration method
US9816786B2 (en) Method for automatically generating a three-dimensional reference model as terrain information for an imaging device
CN113962057B (en) Remote missile active section motion parameter correction method based on time sequence intersection
Helgesen et al. Tracking of ocean surface objects from unmanned aerial vehicles with a pan/tilt unit using a thermal camera
Campbell et al. Vision-based geolocation tracking system for uninhabited aerial vehicles
Moore et al. Vision-only estimation of wind field strength and direction from an aerial platform
Zsedrovits et al. Performance analysis of camera rotation estimation algorithms in multi-sensor fusion for unmanned aircraft attitude estimation
KR102050995B1 (en) Apparatus and method for reliability evaluation of spatial coordinates
Veth et al. Tightly-coupled ins, gps, and imaging sensors for precision geolocation
CN115388890A (en) Visual sense-based multi-unmanned aerial vehicle cooperative ground target positioning method
Yuan et al. Dynamic initial alignment of the MEMS-based low-cost SINS for AUV based on unscented Kalman filter
Yun et al. Range/optical flow-aided integrated navigation system in a strapdown sensor configuration
Bashir et al. Kalman Filter Based Sensor Fusion for Altitude Estimation of Aerial Vehicle
Hoshizaki et al. Performance of Integrated Electro‐Optical Navigation Systems

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant