JP2018151300A - Position estimation method, position estimation device and position estimation program - Google Patents

Abstract

PROBLEM TO BE SOLVED: To provide a position estimation method capable of reducing the time required until the result of a calculation on the position of an aircraft converges.SOLUTION: The position estimation method includes a series of processing of: acquiring a piece of observation information of an aircraft; and calculating an estimation position of the aircraft by applying a nonlinear Bayesian filter based on a state equation which includes a term of the system error in which an initial value of the covariance is set to a first predetermined value and an observation equation which includes a term of an observation residual in which an initial value of the covariance is set to a second predetermined value.SELECTED DRAWING: Figure 5

Description

  The present invention relates to a technique for estimating the position of a flying object.

  There are techniques for estimating the position of a flying object (such as a ballistic bullet, rocket, or artificial satellite). For example, a document discloses a technique for estimating a target position based on a batch filter based on a least square method using observation information of a flying target.

  In the above-described technology, when the calculation result converges or when a time that can be considered that the calculation result has converged has elapsed, the flying object's drop position is determined using the estimated position and estimated speed of the flying object. Calculation to predict is started. Therefore, if the time until the calculation result converges can be shortened, the predicted position of the flying object can be known from an earlier stage. As a result, it is possible to take measures in preparation for the falling of the flying object at an earlier stage.

  However, when the calculation based on the least square method is executed, the calculation amount increases when the number of samples is increased by shortening the sampling interval of the observation values. As a result, the time until the calculation result converges becomes longer.

JP 2002-202367 A

  In one aspect, an object of the present invention is to provide a technique for shortening a time required for a calculation result related to a position of a flying object to converge.

  A position estimation method according to an aspect includes a state equation including a system error term in which observation information of a flying object is acquired and an initial value of covariance is set to a first predetermined value, and an initial value of covariance is 2 including a process of calculating an estimated position of the flying object by applying a nonlinear Bayes filter based on an observation equation including an observation residual term set to a predetermined value of 2 to the observation information.

  In one aspect, the time required for the calculation result regarding the position of the flying object to converge can be shortened.

FIG. 1 is a diagram for explaining the calculation of the position of a ballistic bullet that is a flying object. FIG. 2 is a diagram for explaining the process of predicting the position of the ballistic bullet based on the equation of motion. FIG. 3 is a diagram for explaining the prediction of the landing position based on the orientation process. FIG. 4 is a configuration diagram of a system according to the present embodiment. FIG. 5 is a functional block diagram of the information processing apparatus. FIG. 6 is a diagram illustrating a processing flow of processing executed by the sensor signal processing unit. FIG. 7 is a diagram illustrating a processing flow of processing executed by the orientation processing unit and the landing prediction unit. FIG. 8 is a diagram showing the origin in the local coordinate system. FIG. 9 is a diagram showing the relationship between the origin and the position of the ballistic bullet in the local coordinate system. FIG. 10 is a diagram showing the relationship between the ECEF coordinate system and the NED coordinate system. FIG. 11 is a diagram illustrating a relationship among the NED coordinate system, the body coordinate system, and the visual axis coordinate system. FIG. 12 is a diagram illustrating a processing flow of the orientation processing. FIG. 13 is a diagram illustrating the relationship between the position of the ballistic bullet and the position of the observation device. FIG. 14 is a diagram illustrating a processing flow of landing prediction processing. FIG. 15 is a diagram for explaining calculation of the predicted landing range. FIG. 16 is a diagram showing ballistic bullets in the simulation. FIG. 17 is a diagram illustrating a flight scenario. FIG. 18 is a diagram showing simulation conditions. FIG. 19 is a diagram showing simulation conditions. FIG. 20 is a diagram illustrating a processing flow of simulation. FIG. 21 is a diagram illustrating the time until the calculation result converges. FIG. 22 is a diagram illustrating the orientation error for each sampling interval. FIG. 23 is a diagram illustrating the relationship between the orientation error and the observation time. FIG. 24 is a diagram illustrating the orientation error when the setting value of the covariance of the system error is changed. FIG. 25 is a diagram illustrating the orientation error when the second set value is used. FIG. 26 is a diagram illustrating the prediction error of the landing position. FIG. 27 is a diagram illustrating the relationship between the prediction error of the landing position and the observation time.

FIG. 1 is a diagram for explaining the calculation of the position of the ballistic bullet 3 that is a flying object. The observation device 1 that executes the main processing of the present embodiment acquires the position information of the ballistic bullet 3 as the two-dimensional camera coordinates (H, V) (that is, the pixel position on the image). Position (a, b, c) of the observation device 1 in the ECEF (Earth Center Earth Fixed) coordinate system, visual axis direction (θ _e , Φ _e , ε _e ) and NED (North-East-Down) in the body coordinate system The attitude angles (α, β, γ) in the coordinate system can be measured by the observation device 1.

  The target azimuth (θ, Φ, ε) in the local coordinate system of the ballistic bullet 3 can be calculated using the above two-dimensional camera coordinates, posture angle, and visual axis azimuth. However, the position (x, y, z) of the ballistic bullet 3 cannot be calculated unless the distance from the observation device 1 to the ballistic bullet 3 is known.

  Therefore, in the present embodiment, the position (x, y, z) of the ballistic bullet 3 is calculated on the assumption that the distance from the observer 1 to the ballistic bullet 3 is a predetermined distance.

FIG. 2 is a diagram for explaining processing for predicting the position of the ballistic bullet 3 based on the equation of motion. When the target is the ballistic bullet 3 as in the present embodiment, the ballistic bullet 3 falls free when the boost phase ends after being fired from the launch pad. Therefore, the movement of the ballistic bullet 3 can be expressed by one equation of motion. Specifically, if the calculation is performed based on the equation of motion starting from the position (x _t1 , y _t1 , z _t1 ) calculated on the assumption that the distance from the observation device 1 to the ballistic bullet 3 is a predetermined distance, The position (x _t2 , y _t2 , z _t2 ) after the unit time can be predicted.

  FIG. 3 is a diagram for explaining the prediction of the landing position based on the orientation process. In the present embodiment, the orientation process is started when the boost phase of the ballistic bullet 3 ends. In the orientation process, the estimated position and estimated speed of the ballistic ammunition 3 are corrected so that the error between the target direction calculated at each sampling interval and its predicted value becomes small. In this embodiment, when a time that can be considered that the calculation result has converged has elapsed (for example, several tens of seconds), the predicted landing according to the equation of motion is performed using the estimated position and the estimated speed that are finally calculated. A position (more specifically, a predicted landing range) is calculated.

  FIG. 4 is a configuration diagram of a system according to the present embodiment. A GPS (Global Positioning System) antenna 1000, a search sensor 1001, a tracking sensor 1002, a gimbal 1003, and a position / orientation reference device 1004 are installed outside the observation device 1 in this embodiment. An information processing apparatus 10 is installed in the observation machine 1.

  Information on the position of the observation device 1 in the ECEF coordinate system is acquired by the GPS antenna 1000. Search sensor 1001 and tracking sensor 1002 are infrared sensors that are passive two-dimensional imaging devices. The search sensor 1001 detects the ballistic bullet 3 and acquires the two-dimensional camera coordinates of the ballistic bullet 3. The tracking sensor 1002 acquires the two-dimensional camera coordinates of the ballistic bullet 3 in accordance with the flight of the ballistic bullet 3 detected by the search sensor 1001. The gimbal 1003 controls the direction of the search sensor 1001 and the tracking sensor 1002. The position / orientation reference unit 1004 acquires body information (for example, information such as an attitude angle and a visual axis direction).

  The information processing apparatus 10 includes one or more CPUs (Central Processing Units) 151, one or more memories 153 (for example, DRAM (Dynamic Random Access Memory)), and one or more HDDs (Hard Disk Drive) 155. Have.

  The information processing apparatus 10 executes processing based on the acquired camera coordinates, position information, posture angle, visual axis direction, and the like. Further, the information processing apparatus 10 controls the gimbal 1003 based on the processing result.

  Thus, the system according to the present embodiment has the function of an IRST (Infra-Red Search and Track) system.

  FIG. 5 is a functional block diagram of the information processing apparatus 10. The information processing apparatus 10 includes a gimbal control unit 101, a sensor signal processing unit 103, an orientation processing unit 105, an impact prediction unit 107, an orientation result storage unit 111, and an output data storage unit 113.

  The gimbal control unit 101, the sensor signal processing unit 103, the orientation processing unit 105, and the landing prediction unit 107 are realized, for example, by causing the CPU 151 to execute a program loaded in the memory 153 illustrated in FIG. The orientation result storage unit 111 and the output data storage unit 113 are provided in the memory 153 or the HDD 155, for example.

  The gimbal control unit 101 controls the operation of the gimbal 1003 based on the result of the orientation process described later. The sensor signal processing unit 103 calculates the azimuth of the ballistic bullet 3 (hereinafter referred to as a target azimuth) and outputs the calculated target azimuth information to the orientation processing unit 105. The orientation processing unit 105 executes orientation processing described later based on the processing result of the sensor signal processing unit 103 and other information (for example, a set value), and stores the result of orientation processing in the orientation result storage unit 111. The landing prediction unit 107 executes processing based on the data stored in the orientation result storage unit 111 and stores the processing result in the output data storage unit 113.

  Next, the operation of the information processing apparatus 10 will be described with reference to FIGS.

  FIG. 6 is a diagram illustrating a processing flow of processing executed by the sensor signal processing unit 103. This process is executed every sampling interval (for example, 500 milliseconds).

  The sensor signal processing unit 103 acquires sensor data (specifically, the two-dimensional camera coordinates of the ballistic bullet 3) from the search sensor 1001 and the tracking sensor 1002. Further, the sensor signal processing unit 103 acquires body information (for example, information such as an attitude angle and a visual axis direction) from the position / direction reference unit 1004 via the gimbal control unit 101. Further, the sensor signal processing unit 103 acquires the position information of the observation device 1 from the GPS antenna 1000 (FIG. 6: Step S1).

The sensor signal processing unit 103 calculates the target azimuth of the ballistic bullet 3 based on the information acquired in step S1 (step S3) and stores it in a storage device such as the memory 153. Specifically, the azimuth angle AZ _L (ε = 0, Φ) and the elevation angle EL _L (ε = 0, θ) are calculated. The coordinate system is a local coordinate system and may be the same coordinate system as the visual axis direction or a different coordinate system. Then, the process ends.

  If the target orientation of the ballistic bullet 3 calculated by the above processing is used, the position of the ballistic bullet 3 can be calculated in the subsequent processing.

  FIG. 7 is a diagram illustrating a processing flow of processing executed by the orientation processing unit 105 and the landing prediction unit 107. This process is executed after the ballistic bullet 3 is detected by the search sensor 1001 and tracking by the tracking sensor 1002 is started.

  First, the system model of this embodiment is shown. The system model of the present embodiment is based on the extended Kalman filter, and the state equation is expressed as follows.

  Here, v (t) is a system error, and follows a distribution in which the covariance matrix is Q (t) and 0 is an average.

  It is assumed that the observation equation in the present embodiment is the following observation equation.

  Here, w (t) is an observation residual, follows a distribution in which the covariance matrix is R (t) and 0 is an average.

  In the following description, the covariance matrix may be simply referred to as covariance.

  The time transition model F (t) is expressed as a matrix of 6 rows and 6 columns as follows.

  The hat symbol represents the predicted value, and the bar symbol represents the estimated value.

  The state variable x (t) is expressed as follows.

  The observation model H (t) is expressed as a 2 × 6 matrix as follows.

  Then, the orientation processing unit 105 determines whether or not the boost phase of the ballistic bullet 3 has ended (FIG. 7: Step S11). Whether or not the boost phase has ended is determined based on, for example, acceleration estimated from the transition of the two-dimensional camera coordinates.

  If the boost phase of the ballistic bullet 3 has not ended (step S11: No route), the process returns to step S11.

  On the other hand, when the boost phase of the ballistic bullet 3 is completed (step S11: Yes route), the orientation processing unit 105 calculates the initial position of the ballistic bullet 3 in the local coordinate system (step S13). In step S13, for example, the initial position is calculated according to the following equation.

  R represents a set value of the distance from the observation device 1 to the position where the ballistic bullet 3 is detected for the first time. In the present embodiment, for example, R is set as R = 1200 (km).

  FIG. 8 is a diagram showing the origin in the local coordinate system. Here, the local coordinate system is a three-dimensional orthogonal coordinate system, and the positions of the search sensor 1001 and the tracking sensor 1002 (that is, the position of the observation device 1) are the origin O. The positive direction of the X axis is the north direction, the positive direction of the Y axis is the east direction, and the positive direction of the Z axis is the zenith direction.

  FIG. 9 is a diagram illustrating the relationship between the origin and the position of the ballistic bullet 3 in the local coordinate system. O represents the origin and corresponds to the position of the observation device 1. A position 901 whose distance from the origin is R corresponds to the position of the ballistic bullet 3. The position 901 of the ballistic bullet 3 is specified by the azimuth angle Φ and the elevation angle θ with respect to the origin.

  The orientation processing unit 105 finally converts the coordinates of the initial position of the ballistic bullet 3 into coordinates in the ECEF coordinate system.

  Here, the transformation of the coordinate system will be described.

FIG. 10 is a diagram showing the relationship between the ECEF coordinate system and the NED coordinate system. The ECEF coordinate system is an orthogonal coordinate system, and the origin is at the center of the earth. On the other hand, the NED coordinate system is an orthogonal coordinate system whose origin is the position of the sensor (observer 1 in the present embodiment) (the position where the longitude is λ and the latitude is δ in FIG. 10). As described above, the positive direction of the X axis is the north direction, the positive direction of the Y axis is the east direction, and the positive direction of the Z axis is the zenith direction. If R _{ECEF_NED} is a determinant that _converts the ECEF coordinate system into the coordinates of the local coordinate system with the position of the observation device 1 as the origin, R _{NED_ECEF} = R _{ECEF_NED} ⁻¹ is satisfied.

FIG. 11 is a diagram illustrating a relationship among the NED coordinate system, the body coordinate system, and the visual axis coordinate system. In FIG. 11, R _{NED_AC} is a determinant that transforms the NED coordinate system into the airframe coordinate system using posture angles (roll, pitch, yaw), and R _{AC_VI} uses the visual axis direction to convert the airframe coordinate system into the visual axis coordinates It is a determinant to convert to a system. In this case, R _{AC_NED} = R _{NED_AC} ⁻¹ is satisfied, and R _{VI_AC} = R _{AC_VI} ⁻¹ is satisfied. In the present embodiment, when the rotation θ is 0 in the visual axis coordinate system, ε = AZ _L and Φ = EL _L.

  The relationship between the ECEF coordinate system and the geodetic coordinate system is as follows.

  Here, X, Y, and Z represent the coordinates of the ECEF coordinate system, respectively. λ, Φ, and h represent latitude, longitude, and altitude, respectively. a represents the equator radius and is 6,378,137 (m). f represents the flatness ratio and is 1 / 298.257223563.

  The ECI (Earth-Centered Inertial) coordinate system and the ECEF coordinate system have a difference of ΔΦ with respect to the rotation angle around the Z axis due to the earth rotation, and the relationship between the ECI coordinate system and the ECEF coordinate system is as follows.

  Here, x, y, and z each represent a coordinate in the ECI coordinate system. ΔΦ represents a rotation angle around the z-axis (Z-axis) axis of the ECEF coordinate system with respect to the ECI coordinate system.

  Conversion from the geodetic coordinate system to the ECI coordinate system or the ECEF coordinate system can be performed by substituting into the equation shown above.

  The coordinate conversion from the ECI coordinate system or the ECEF coordinate system to the geodetic coordinate system is as follows.

  Here, a represents the Earth's long radius, which is 6378137.0 (m). e represents the eccentricity and is 0.0818191910428158. f represents the aspect ratio, which is 0.003352811.

  Returning to the explanation of FIG. 7, the orientation processing unit 105 calculates the initial velocity of the ballistic bullet 3 (step S15). The initial velocity of the ballistic bullet 3 is calculated from, for example, the initial position of the ballistic bullet 3 and the motion model. However, it may be calculated by other methods, or the initial speed may be set in advance.

  The orientation processing unit 105 sets the initial covariance of the system error and the initial covariance of the observation residual (step S17). The initial covariance of the system error is based on, for example, the sum of the residues of all singular points when it is assumed that the ballistic bullet 3 moves at equal acceleration in a minute time. The initial covariance of the observation residual is based on, for example, the maximum error that can occur during observation. However, the initial covariance of the system error and the initial covariance of the observation residual may be set on the assumption that the ballistic bullet 3 moves at equal acceleration in a very short time.

  The orientation processing unit 105 determines whether the altitude of the ballistic bullet 3 is smaller than 0 (step S19).

  When the altitude of the ballistic bullet 3 is smaller than 0 (step S19: Yes route), the process ends. When the altitude of the ballistic bullet 3 is not smaller than 0 (step S19: No route), the orientation processing unit 105 executes orientation processing (step S21). The orientation process will be described with reference to FIGS. 12 and 13.

  First, the orientation processing unit 105 converts the coordinate system from an ECEF coordinate system to a rotating coordinate system (that is, a non-inertial coordinate system) (FIG. 12: step S31).

  The orientation processing unit 105 calculates the predicted position and predicted speed of the ballistic bullet 3 (step S33). In step S33, the predicted position is calculated based on the following equation.

  Here, x hat (t + Δt) is the predicted value of the x coordinate at time (t + Δt), y hat (t + Δt) is the predicted value of the y coordinate at time (t + Δt), and z hat (t + Δt) is the time (t + Δt). ) Is the predicted value of the z coordinate. x bar (t) is an estimated value of x coordinate at time t, y bar (t) is an estimated value of y coordinate at time t, and z bar (t) is an estimated value of z coordinate at time t. Is the initial position of the ballistic bullet 3 when t = 0. Vx bar (t) is an estimated value of velocity in the x direction at time t, Vy bar (t) is an estimated value of velocity in the y direction at time t, and Vz bar (t) is an estimated value of z direction at time t. The estimated value of the velocity is the initial velocity of the ballistic bullet 3 when t = 0. Δt is a sampling interval (second).

Note that the predicted speed of the ballistic bullet 3 is calculated based on, for example, an estimated value of the velocity of the ballistic bullet 3 at time t, gravity acceleration (m / s ² ), Coriolis force, and centrifugal force.

  The orientation processing unit 105 calculates a predicted value of system error covariance (step S35). The predicted value of the system error covariance is calculated as follows.

  Here, P hat (t + Δt | t) represents a predicted value of the covariance of the system error at time t. P bar (t | t) represents an estimated value of the covariance of the system error at time t.

  The orientation processing unit 105 converts the coordinate system from a rotating coordinate system to a local coordinate system (that is, an inertial coordinate system) (step S37). In step S37, conversion is performed based on the position of the observation device 1 at time (t + Δt).

  The orientation processing unit 105 calculates a predicted value of the target azimuth based on the position at time (t + Δt) (step S39). The predicted value of the target direction is calculated as follows.

Here, EL _L (t | t) and AZ _L (t | t) are as follows in the second quadrant and the fourth quadrant.

  FIG. 13 is a diagram illustrating the relationship between the position of the ballistic bullet 3 and the position of the observation device 1. O represents the origin and corresponds to the position of the observation device 1. The position 1301 of the ballistic bullet 3 is specified by the azimuth angle Φ and the elevation angle θ with respect to the origin. Vx, Vy, and Vz each represent a velocity vector. The direction of the combined vector of Vx, Vy, and Vz corresponds to the movement direction of the ballistic bullet 3.

  The orientation processing unit 105 calculates a predicted value of covariance of the observation residual (step S41). The predicted value of the covariance of the observation residual is calculated as follows.

  Here, H (t + Δt) represents an observation matrix at time (t + Δt).

  The orientation processing unit 105 calculates the correction amount (that is, Kalman gain) of the predicted value (step S43). The Kalman gain is calculated as follows.

  Here, K (t + Δt) represents the Kalman gain at time (t + Δt). S (t + Δt) represents the covariance of the observation residual at time (t + Δt).

  The orientation processing unit 105 calculates the estimated position and the estimated speed (step S45), and stores the estimated position and the estimated speed in the orientation result storage unit 111. The estimated position and the estimated speed are calculated as follows.

  Here, z (t + Δt) represents an observed value (including information on the target direction) at time (t + Δt), and is acquired from the sensor signal processing unit 103.

  The orientation processing unit 105 calculates an estimated value of system error covariance (step S47). Processing then returns to the caller. The estimated value of the system error covariance is calculated as follows.

  Here, P-bar (t + Δt | t + Δt) represents an estimated value of the system error covariance at time (t + Δt).

  In this way, using the extended Kalman filter based on the state equation and observation equation as described above (the extended Kalman filter is a type of Bayes filter) reduces the amount of calculation, thereby shortening the sampling interval and increasing the number of samples. Even so, it is possible to suppress an increase in the amount of calculation. Note that the extended Kalman filter is a technique that can shorten the time until convergence by applying it to the motion in accordance with the laws of physics.

  Returning to the description of FIG. 7, the orientation processing unit 105 determines whether a predetermined time (for example, 200 seconds) has elapsed since the boost phase ended (step S <b> 23).

  If the predetermined time has not elapsed (step S23: No route), the process returns to step S19. When the predetermined time has elapsed (step S23: Yes route), the orientation processing unit 105 notifies the landing prediction unit 107 that the orientation processing has been completed. In response, the landing prediction unit 107 executes a landing prediction process (step S25). The landing prediction process will be described with reference to FIGS. 14 and 15.

  First, the landing prediction unit 107 reads the result of the orientation process (that is, the estimated position and estimated speed calculated last in the orientation process) from the orientation result storage unit 111 (FIG. 14: step S51).

The landing prediction unit 107 calculates the predicted position and predicted speed of the ballistic bullet 3 using the estimated position and estimated speed read in step S51 (step S53). In step S53, the predicted position and the predicted speed are calculated based on Expression (14). However, t is the time corresponding to the orientation results storage section 111 the estimated position and estimated position is read (hereinafter, the time t _s). The initial value of Δt is a predetermined unit time, and is updated as 2Δt, 3Δt, 4Δt,... Each time the process of step S53 is executed.

  The landing prediction unit 107 determines whether the altitude of the predicted position of the ballistic bullet 3 is smaller than 0 (step S55). When the altitude of the predicted position of the ballistic bullet 3 is not smaller than 0 (step S55: No route), the landing prediction unit 107 updates Δt. Then, the process returns to step S53.

  On the other hand, when the altitude of the predicted position of the ballistic bullet 3 is smaller than 0 (step S55: Yes route), the landing prediction unit 107 calculates a predicted landing range (step S57).

  The landing prediction unit 107 stores the data of the predicted landing range calculated in step S57 in the output data storage unit 113 (step S59). Then, the process returns to the caller and ends.

  Note that the predicted landing range data is transmitted to, for example, a computer on the ground and displayed on the display device or the like, and is confirmed by the administrator.

FIG. 15 is a diagram for explaining calculation of the predicted landing range. The position of the ballistic bullet 3 when the altitude of the predicted position becomes 0 is the predicted landing position. The predicted landing range corresponds to an ellipse when an error ellipsoid centered on the predicted landing position is projected onto the ground surface in the velocity vector direction. In addition, the time t _f is the predicted impact time.

  Next, simulation based on the method of the present embodiment will be described with reference to FIGS.

FIG. 16 is a diagram showing the ballistic bullet 3 in this simulation. The ballistic bullet 3 is a 1300 km class medium-range ballistic missile. The total length of the ballistic bullet 3 is 16.0 (m) and the width is 1.32 (m). The cross-sectional area is 13700 (cm ² ) and the resistance coefficient is 0.1. The mass of the payload is 1000 (kg), and the total mass at the time of launch is 16250 (kg). A portion represented by a one-dot chain line represents a plume.

  FIG. 17 is a diagram illustrating a flight scenario. In this simulation, three scenarios are tried. The scenario for the observation device 1a is set as scenario A, the scenario for the observation device 1b is set as scenario B, and the scenario for the observation device 1c is set as scenario C. The flying altitude of each observation aircraft is 25000 feet. The observation devices 1a to 1c stay at a position where the distance from the position where the ballistic bullet 3 is first detected is 800 (km). However, the observation device 1a stays at a position immediately below the trajectory of the ballistic bullet 3, the observation device 1b stays at a position 200 (km) away from just below the orbit, and the observation device 1c is at a position 400 (km) away from just below the orbit. Stay in the air.

  The simulation conditions are as shown in FIGS. The distance from the observation devices 1a to 1c to the position where the ballistic bullet 3 is first detected is 800 (km) as described above, but a relative distance error of 400 (km) is superimposed as shown in FIG. Therefore, the distance R is set to 1200 (km).

  FIG. 20 is a diagram showing a processing flow of this simulation. Here, the operation of the observation device 1a will be described as an example, but the observation device 1b and the observation device 1c operate similarly.

  The observation device 1a generates a flight path (step S1901), and generates a trajectory of the target ballistic bullet 3 (hereinafter referred to as a target trajectory) (step S1902).

  The observation device 1a calculates an observation angle (the observation angle corresponds to the target azimuth described above) based on the generated flight path and target trajectory (step S1903). The observation angle calculated in step S1903 is treated as a true value.

  The observation device 1a generates a position error of the observation device 1a from the generated flight path (step S1904). Further, the observation device 1a generates an angle measurement error from the observation angle that is a true value (step S1905).

  The observation device 1a executes the orientation process based on the observation device position (including the error) and the observation angle (including the error) of the observation device 1a (step S1906).

  An orientation error (in this embodiment, a positional deviation) is calculated from the orientation result generated by the orientation processing in step S1906 and the true value of the target position based on the target trajectory. The above is the processing flow of this simulation.

  FIG. 21 is a diagram showing time until calculation results converge for each scenario. The sampling interval is 500 (milliseconds). The upper part of the table shows the time when the simulation is executed under the same simulation conditions based on the algorithm of the least square method described in Patent Document 1, and the lower part shows the time when the simulation based on the method of the present embodiment is executed. Show. As shown in FIG. 21, in this embodiment, the time is shortened in all scenarios, and in particular, scenario B and scenario C are shortened by about 50%.

  As for the calculation amount of one cycle (here, the calculation amount when the orientation process is executed once), the calculation amount when the method of the present embodiment is executed is the least square method described in Patent Document 1. It has been confirmed that the amount of calculation is considerably smaller than the amount of calculation when the algorithm is used.

  FIG. 22 is a diagram illustrating the time until convergence for each sampling interval. The upper row shows the time when the sampling interval is 500 (milliseconds), and the lower row shows the time when the sampling interval is 100 (milliseconds). Since the calculation amount is reduced by the method of the present embodiment, the sampling interval can be shortened. It was confirmed that the time until convergence is shortened by about 10% by shortening the sampling interval.

  FIG. 23 is a diagram illustrating the relationship between the orientation error and the observation time. FIG. 23 shows a graph for scenario A, where the sampling interval is 100 (milliseconds). The vertical axis represents the orientation error, and the horizontal axis represents the observation time. As shown in FIG. 23, it was confirmed that when the observation time is longer than 100 seconds, the orientation error is rapidly reduced, and when the observation time is about 170 seconds, the orientation error is 10 (km) or less.

  FIG. 24 is a diagram illustrating time to convergence when the setting value of the initial covariance of the system error is changed. The sampling interval is 100 (milliseconds). The upper part of the table of FIG. 24 shows the time when the simulation is executed with the first set value, and the lower part of the table of FIG. 24 is the time when the simulation is executed with the second set value larger than the first set value. Indicates. By increasing the set value, the difference from the actual observed value variance increases, so the amount of correction of the position and velocity increases. Therefore, for scenario A, it was confirmed that the time until convergence was shortened by using the second set value. However, the results for Scenario B and Scenario C were almost the same.

  FIG. 25 is a diagram illustrating the orientation error when the second set value is used. FIG. 25 shows a graph for scenario A, where the sampling interval is 100 (milliseconds). The vertical axis represents the orientation error, and the horizontal axis represents the observation time. When the second set value was used, it was confirmed that the orientation error rapidly decreased before the observation time reached 100 (seconds).

  FIG. 26 is a diagram illustrating the prediction error of the landing position. The sampling interval is 100 (milliseconds). The upper part of the table in FIG. 26 shows the prediction error when the simulation is executed under the same simulation condition based on the algorithm of the least square method described in Patent Document 1, and the lower part shows the simulation based on the method of the present embodiment. The prediction error in the case is shown. There was no significant improvement in accuracy in each scenario.

  FIG. 27 is a diagram illustrating the relationship between the prediction error of the landing position and the observation time. FIG. 27 shows a graph for scenario A, where the sampling interval is 100 (milliseconds). The vertical axis represents the landing position prediction error, and the horizontal axis represents the observation time. As shown in FIG. 27, it was confirmed that the landing position prediction error was 10 (kmCEP) or less when the observation time was about 130 seconds.

  In Bayesian estimation, the true value is treated as a random variable, and the true value distribution is updated based on the sample. Therefore, if the distribution of the sample can be assumed in advance, it can be expected that convergence will be accelerated by incorporating it in the design. In the present embodiment, the amount of calculation can be reduced by appropriately setting the initial covariance of the system error and the initial covariance of the observation residual in advance. In addition, when an algorithm based on the least square method is used, the amount of calculation in one cycle may increase in proportion to the number of samples. However, in the method according to the present embodiment, the amount of calculation does not matter regardless of the number of samples. Can be made constant.

  As described above, since the sampling interval can be shortened, the number of samples that can be acquired within the same observation period can be increased. As a result, the landing position of the ballistic bullet 3 can be predicted at an earlier stage.

  Although one embodiment of the present invention has been described above, the present invention is not limited to this. For example, the functional block configuration of the information processing apparatus 10 described above may not match the actual program module configuration.

  Also in the processing flow, if the processing result does not change, the processing order can be changed. Further, it may be executed in parallel.

  It should be noted that this embodiment can be applied to flying objects other than the ballistic bullet 3. Further, the calculation may be performed using non-inertia system coordinates instead of inertia system coordinates.

  The embodiment of the present invention described above is summarized as follows.

  In the position estimation method according to the first aspect of the present embodiment, (A) the observation information of the flying object is acquired, and (B) the system error term in which the initial value of the covariance is set to the first predetermined value. Is applied to the observation information by applying a non-linear Bayesian filter to the observation information, and an observation equation including an observation residual term whose initial covariance value is set to the second predetermined value. Including the process of calculating.

  Since the non-linear Bayes filter as described above is used, it is possible to shorten the time required until the calculation result related to the position of the flying object converges.

  In addition, the first predetermined value may be a value based on the sum of the residues of all singular points, and the second predetermined value may be a value based on the maximum error that can occur in the observation information. .

  Given a reasonable covariance, the time to convergence is shortened.

  The nonlinear Bayes filter may be an extended Kalman filter.

  Further, in the process of calculating the estimated position, (b1) the predicted position of the flying object is calculated according to the equation of motion in the non-inertial system based on the observation information or the previously calculated estimated value, and (b2) the predicted position and the observation information (B3) the estimated position may be calculated based on the predicted position and the Kalman gain.

  For example, Coriolis force and centrifugal force can be reflected in the estimated position.

  In the process of calculating the predicted position, (b11) the coordinate system of the predicted position is converted from the rotating coordinate system to the local coordinate system, and in the process of calculating the estimated position, (b31) the predicted position and the Kalman in the local coordinate system. Based on the gain, an estimated position in the local coordinate system may be calculated.

  The value of the position coordinate is calculated appropriately.

  In addition, the observation information may include information on a two-dimensional image including a flying object acquired by an infrared detection device. In the process of calculating the estimated position, (b311) the direction of the flying object may be calculated based on at least the information of the two-dimensional image, and the observed position of the flying object may be calculated based on the direction.

  Even when the flying object is located far away, the flying object may be captured by the infrared device. Therefore, the position of a distant projectile can be estimated by using the information of the two-dimensional image acquired by the infrared detection device (for example, information on the pixel position of the projectile).

  Further, in the process of calculating the estimated position, (b4) the estimated speed of the flying object may be further calculated from the observation information. The position estimation method may further include (C) a process of calculating the range of the flying object's fall position from the estimated position and the estimated speed, and (D) outputting information on the range of the fall position.

  It will be possible to take measures against flying objects falling.

  Further, the covariance of the system error may be represented by a 6 × 6 matrix including elements for each component of the three-dimensional position and each component of the three-dimensional velocity.

  The state equation may include a term for position and velocity in an orthogonal coordinate system with the center of the earth as the origin, and the observation equation may include a term for an angle in the local coordinate system.

  The position estimation apparatus according to the second aspect of the present embodiment includes (E) an acquisition unit that acquires observation information of the flying object (search sensor 1001 and tracking sensor 1002 of the embodiment are examples of an acquisition unit). And (F) a state equation including a system error term in which the initial value of the covariance is set to the first predetermined value, and an observation residual term in which the initial value of the covariance is set to the second predetermined value. Is applied to the observation information acquired by the acquisition unit by applying a nonlinear Bayes filter based on the observation equation including the calculation unit (the orientation processing unit 105 of the embodiment is an example of a calculation unit) It is).

  A program for causing a computer to perform the processing according to the above method can be created. The program can be a computer-readable storage medium such as a flexible disk, a CD-ROM, a magneto-optical disk, a semiconductor memory, a hard disk, or the like. It is stored in a storage device. The intermediate processing result is temporarily stored in a storage device such as a main memory.

  The following supplementary notes are further disclosed with respect to the embodiments including the above examples.

(Appendix 1)
On the computer,
Obtain observation information of the flying object,
A state equation including a system error term with an initial covariance value set to a first predetermined value, and an observation equation including an observation residual term with an initial covariance value set to a second predetermined value; Applying a non-linear Bayesian filter based on the observation information to calculate an estimated position of the flying object;
A position estimation program that executes processing.

(Appendix 2)
The first predetermined value is a value based on the sum of the residues of all singular points;
The second predetermined value is a value based on a maximum error that can occur in the observation information.
The position estimation program according to attachment 1.

(Appendix 3)
The nonlinear Bayes filter is an extended Kalman filter.
The position estimation program according to appendix 1 or 2.

(Appendix 4)
In the process of calculating the estimated position,
Based on the observation information or the previously calculated estimated value, calculate the predicted position of the flying object according to the equation of motion in a non-inertial system,
Based on the predicted position and the observation information, Kalman gain is calculated,
Calculating the estimated position based on the predicted position and the Kalman gain;
The position estimation program according to attachment 3.

(Appendix 5)
In the process of calculating the predicted position,
Converting the coordinate system of the predicted position from a rotating coordinate system to a local coordinate system;
In the process of calculating the estimated position,
Calculating the estimated position in the local coordinate system based on the predicted position in the local coordinate system and the Kalman gain;
The position estimation program according to attachment 4.

(Appendix 6)
The observation information includes information of a two-dimensional image including the flying object acquired by an infrared detection device,
In the process of calculating the estimated position,
Based on at least the information of the two-dimensional image, calculate the orientation of the flying object,
Based on the azimuth, calculate the flying object observation position,
The position estimation program according to attachment 4.

(Appendix 7)
In the process of calculating the estimated position,
Further calculating the estimated speed of the flying object from the observation information,
In the computer,
From the estimated position and the estimated speed, calculate the range of the falling position of the flying object,
Outputting information of the range of the drop position;
The position estimation program according to supplementary note 1, wherein the program is further executed.

(Appendix 8)
The covariance of the system error is represented by a 6-by-6 matrix including elements for each component of the three-dimensional position and each component of the three-dimensional velocity.
The position estimation program according to attachment 2.

(Appendix 9)
The equation of state includes terms for position and velocity in an orthogonal coordinate system with the center of the earth as the origin,
The observation equation includes a term for the angle in the local coordinate system,
The position estimation program according to attachment 1.

(Appendix 10)
Computer
Obtain observation information of the flying object,
A state equation including a system error term with an initial covariance value set to a first predetermined value, and an observation equation including an observation residual term with an initial covariance value set to a second predetermined value; Applying a non-linear Bayesian filter based on the observation information to calculate an estimated position of the flying object;
A position estimation method for executing processing.

(Appendix 11)
An acquisition unit for acquiring observation information of the flying object;
A state equation including a system error term with an initial covariance value set to a first predetermined value, and an observation equation including an observation residual term with an initial covariance value set to a second predetermined value; Applying a non-linear Bayesian filter based on the observation information acquired by the acquisition unit to calculate an estimated position of the flying object;
A position estimation apparatus.

DESCRIPTION OF SYMBOLS 1 Observation machine 10 Information processing apparatus 151 CPU 153 Memory 155 HDD 1000 GPS antenna 1001 Search sensor 1002 Tracking sensor 1003 Gimbal 1004 Position and orientation reference device 101 Gimbal control unit 103 Sensor signal processing unit 105 Orientation processing unit 107 Landing prediction unit 111 Standardization result storage Part 113 Output data storage part

Claims (11)

On the computer,
Obtain observation information of the flying object,
A state equation including a system error term with an initial covariance value set to a first predetermined value, and an observation equation including an observation residual term with an initial covariance value set to a second predetermined value; Applying a non-linear Bayesian filter based on the observation information to calculate an estimated position of the flying object;
A position estimation program that executes processing. The first predetermined value is a value based on the sum of the residues of all singular points;
The second predetermined value is a value based on a maximum error that can occur in the observation information.
The position estimation program according to claim 1. The nonlinear Bayes filter is an extended Kalman filter.
The position estimation program according to claim 1 or 2. In the process of calculating the estimated position,
Based on the observation information or the previously calculated estimated value, calculate the predicted position of the flying object according to the equation of motion in a non-inertial system,
Based on the predicted position and the observation information, Kalman gain is calculated,
Calculating the estimated position based on the predicted position and the Kalman gain;
The position estimation program according to claim 3. In the process of calculating the predicted position,
Converting the coordinate system of the predicted position from a rotating coordinate system to a local coordinate system;
In the process of calculating the estimated position,
Calculating the estimated position in the local coordinate system based on the predicted position in the local coordinate system and the Kalman gain;
The position estimation program according to claim 4. The observation information includes information of a two-dimensional image including the flying object acquired by an infrared detection device,
In the process of calculating the estimated position,
Based on at least the information of the two-dimensional image, calculate the orientation of the flying object,
Based on the azimuth, calculate the flying object observation position,
The position estimation program according to claim 4. In the process of calculating the estimated position,
Further calculating the estimated speed of the flying object from the observation information,
In the computer,
From the estimated position and the estimated speed, calculate the range of the falling position of the flying object,
Outputting information of the range of the drop position;
The position estimation program according to claim 1, further causing the process to be executed. The covariance of the system error is represented by a 6-by-6 matrix including elements for each component of the three-dimensional position and each component of the three-dimensional velocity.
The position estimation program according to claim 2. The equation of state includes terms for position and velocity in an orthogonal coordinate system with the center of the earth as the origin,
The observation equation includes a term for the angle in the local coordinate system,
The position estimation program according to claim 1. Computer
Obtain observation information of the flying object,
A state equation including a system error term with an initial covariance value set to a first predetermined value, and an observation equation including an observation residual term with an initial covariance value set to a second predetermined value; Applying a non-linear Bayesian filter based on the observation information to calculate an estimated position of the flying object;
A position estimation method for executing processing. An acquisition unit for acquiring observation information of the flying object;
A state equation including a system error term with an initial covariance value set to a first predetermined value, and an observation equation including an observation residual term with an initial covariance value set to a second predetermined value; Applying a non-linear Bayesian filter based on the observation information acquired by the acquisition unit to calculate an estimated position of the flying object;
A position estimation apparatus.

Images