Terrain matching aided navigation method without satellite navigation

Abstract

The application discloses a terrain matching aided navigation method without satellite navigation, and mainly solves the problem of the accumulated error of the inertial navigation in the existing technology when the satellite positioning cannot be used. The method is as follows: according to the current time inertial navigation estimated position longitude and latitude, the inertial navigation estimated position error and the corresponding elevation terrain map resolution, the range needing terrain matching is calculated; translation, convolution, height measurement and multiplication calculation are sequentially carried out in the range to obtain a posterior probability density array; the mean and mean square deviation of the array are calculated, and the left upper corner longitude and latitude of the matching range is combined to obtain the current time terrain matching estimated position longitude and latitude, the posterior standard deviation, the map height standard deviation and the confidence degree. The application can effectively reduce the accumulated error of the inertial navigation system drift under the condition that the aircraft flies for a long time or the satellite positioning fails, the positioning precision reaches within 100 meters, and can be used for the safe flight of the aircraft, the missile or the unmanned aerial vehicle without satellite navigation.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of navigation technology, and particularly relates to a terrain matching assisted navigation method, which can be used for the flight of an airplane, a missile or a UAV in the absence of satellite navigation. BACKGROUND

[0002] Navigation is widely used in modern military technology, and it provides key speed, position and time information for various military operations and weapon equipment, and occupies a very important position. The navigation technology commonly used on a UAV or a cruise missile is a global positioning system (GPS), an inertial navigation system (INS) and a combined navigation system (INS / GPS).

[0003] The inertial navigation system (INS) measures the current heading and speed of a flying vehicle according to its internal gyroscopes and accelerometers, and thus recursively estimates the position of the flying vehicle. However, the inertial navigation has a drift characteristic, and the voyage of a flying vehicle can be hundreds of kilometers to thousands of kilometers long. If an uncorrected inertial navigation system is used, the position error caused by the drift characteristic will gradually increase with the increase of the voyage, and accurate positioning or striking cannot be achieved.

[0004] The combined navigation system (INS / GPS) actually uses the global positioning system (GPS) to assist in correcting the INS, and it can achieve long-range accurate navigation in the case of normal GPS signals. However, it relies on GPS signals, and in the case of interference or forgery of GPS signals, it will appear false position information.

[0005] The global positioning system (GPS) is a high-precision radio navigation positioning system based on artificial earth satellites, and it can provide accurate geographic position, speed and precise time information anywhere in the world and near space. It is called satellite navigation, and it is also the currently commonly used accurate navigation system. However, in the case of human interference, tampering and forgery of GPS signals, satellite navigation GPS cannot be used, and the flying vehicle will lose its direction and even crash. SUMMARY

[0006] The present application aims at the deficiencies of the prior art, and provides a terrain matching assisted navigation method in the absence of satellite navigation, so as to obtain accurate position information, avoid the loss of direction of a flying vehicle and ensure its safe flight.

[0007] The technical scheme of the present application is as follows: when satellite positioning GPS cannot be used or a flying vehicle flies for a long time, the flying data of the flying vehicle and a flight area map are obtained, the flying data is continuously calculated and corrected, and thus more accurate position information is obtained for safe flight. The implementation steps include the following:

[0008] 1) The aircraft collects flight data and the resolution of the height map of its flight area, the flight data including the inertial navigation estimated position longitude J, the inertial navigation estimated position latitude W, the barometric height Q, the radar height R, the inertial navigation estimated position error E, and the resolution of the height map including the map longitudinal resolution FA and the map latitudinal resolution FB;

[0009] 2) According to the resolution of the height map, the inertial navigation estimated position error E m at the current time t m and the inertial navigation estimated position longitude J m and latitude W m , the initial range for terrain matching is calculated, m is 0-n, n is an integer and ≥1;

[0010] 3) The terrain matching calculation at the first time t0 is performed:

[0011] 3a) The area height map data within the terrain matching range is taken out from the flight control computer memory of the aircraft, and combined with the inertial navigation estimated position longitude J0 and latitude W0 at the time t0, the matching range upper left corner longitude PJ0 and latitude PW0 are calculated;

[0012] 3b) The numerical difference between the barometric height Q0 and the radar height R0 at the time t0 is taken as the measurement value m0, combined with the height error and the height value of the area height map, the height measurement probability density array H0 is obtained;

[0013] 3c) The height measurement probability density array H0 is multiplied by the prior probability density array P0 of the same size but with elements of 1 to obtain the posterior probability density array PA0 of the first terrain matching at the time t0;

[0014] 4) The terrain matching calculation at the subsequent time t n is performed:

[0015] 4a) According to the difference between the inertial navigation estimated position longitude and latitude at the current time t n and the previous time t n-1 , the terrain matching range at the previous time t n-1 is translated, and the posterior probability density array PA n-1 obtained by multiplying the previous time t n-1 is taken as the prior probability density array P n of the matching time t n ;

[0016] 4b) According to the known set navigation error and sharp coefficient, the convolution probability density array C n is obtained by convolution calculation on the prior probability density array P n ;

[0017] 4c) Repeat step 2) to obtain the current tn The terrain matching range at the moment is intersected with the terrain matching range at the previous moment t n-1 ;

[0018] 4d) The area elevation map of the intersection part of the matching range is called to carry out the height measurement calculation at the current moment t n , and the height measurement probability density array H n at the current moment t n is obtained.

[0019] 4e) The height measurement probability density array H n at the current moment t n is multiplied with the convolution probability density array C n correspondingly, and the posterior probability density array PA n at the current moment t n is obtained.

[0020] 5) The standard deviation of the map height in the matching range in steps 3a) and 4d) is calculated respectively, and the terrain height standard deviation HS m at the moment t m is obtained.

[0021] 6) The mean value and mean square deviation of the posterior probability density array PA m obtained in steps 3c) and 4e) are calculated, and the longitude EJ m and latitude EW m of the estimated position of the terrain matching at the moment t m , the posterior standard deviation BS m and the confidence Z m are calculated by combining the longitude and latitude of the upper left corner of the corresponding matching range.

[0022] 7) The calculation results of steps 5) and 6) are output, and the terrain matching at the moment t m is completed.

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

[0024] I. High positioning accuracy.

[0025] In the prior art, the integrated navigation system INS / GPS relies on GPS signals, and in the case of human interference, tampering and falsification of GPS signals, incorrect position information will be obtained, and the same is true for the global positioning system GPS; and the position error caused by the inherent drift characteristics of the inertial navigation system INS will gradually accumulate and become larger with the increase of the flight time of the aircraft, which can reach thousands or even tens of thousands of meters.

[0026] The application calculates the position information at the current time by comparing the height data on the high-precision height topographic map with the height of the terrain under the aircraft according to the measured height of the terrain under the aircraft, and can reduce the accumulated error of the drift of the inertial navigation system, obtain more accurate position information, and ensure the safe flight of the aircraft.

[0027] II. Fast calculation speed

[0028] The flight data is input once at different times t m , only the mathematical calculations of translation, convolution, multiplication and probability density are performed, and the estimated position of the terrain matching of this time can be obtained within 0.2s or even shorter time. BRIEF DESCRIPTION OF DRAWINGS

[0029] Figure 1 The flowchart for the realization of the application;

[0030] Figure 2 The comparison chart of the longitude and latitude of the predicted position of the terrain matching under a certain flight route, the longitude and latitude of the inertial navigation estimated position and the real flight longitude and latitude;

[0031] Figure 3 The comparison chart of the predicted position error and the inertial navigation estimated position error under a certain flight route.

[0032] Figure 4 The posterior standard deviation chart of the predicted position of a certain flight route of the application;

[0033] Figure 5 The terrain height standard deviation chart under a certain flight route calculated by the application. DETAILED DESCRIPTION

[0034] With reference Figure 1 to the implementation steps of the application:

[0035] Step 1: Obtain the flight data of the aircraft and the height topographic map data.

[0036] 1.1) Obtain the flight data from the flight control computer memory of the aircraft, and the data includes the current barometric height Q, the radar height R, the inertial navigation estimated position, the inertial navigation estimated position error E, the inertial navigation estimated position includes the inertial navigation estimated position longitude J and the inertial navigation estimated position latitude W;

[0037] 1.2) Obtain the elevation topographic map data of the flight area where the aircraft is located from the high-speed hard disk of the flight control computer. The map format is ".dxt" and includes the map starting longitude jb, map starting latitude wb, map longitude sampling number jn, map latitude sampling number wn, the reciprocal of the longitude difference between two adjacent sampling points in the longitude direction jd, the reciprocal of the latitude difference between two adjacent sampling points in the latitude direction wd, and map altitude data h of different sampling points, where jn and wn are positive integers.

[0038] Step 2: Calculate the range of terrain matching.

[0039] 2.1) Calculate the longitude resolution FA and latitudinal resolution FB of the map based on the reciprocal jd of the longitude difference between two adjacent sampling points in the longitude direction and the reciprocal wd of the latitude difference between two adjacent sampling points in the latitude direction of the topographic map data.

[0040]

[0041] 2.2) Calculate the latitude W of the position estimated by the inertial navigation system. m Longitude JM represented by 1 meter on the same latitude line m And the inertial navigation system predicts the longitude J m Latitude WM represented by 1 meter on the same meridian m :

[0042]

[0043]

[0044] In the formula, r is the Earth's radius, approximately 6,371,000 meters;

[0045] 2.3) Compare t m Inertial navigation system position prediction error E at time 1 m The position error ER is obtained by taking the larger of the set error SE. m By combining the map's meridional resolution FA and latitudinal resolution FB, a region with a radius of t is drawn on the elevation topographic map. m Inertial navigation system predicts position longitude J m Inertial navigation system predicts latitude and position W m Centered on, with twice the position error ER m The initial matching range is defined by the side length, specifically the longitude range from 0 to ro. m The latitudinal range is 0 to ra. m ,in:

[0046] ro m =ER m ×JM m ÷FA×2

[0047] ra m = ER m × WM m ÷ FB × 2

[0048] wherein, ro m , ra m are positive integers, m is 0~n, n is an integer and ≥1.

[0049] Step 3, the first time t0 topographic matching calculation is carried out.

[0050] 3.1) according to the inertial navigation system predicted position longitude J0, latitude W0, map longitude resolution FA, map latitude resolution FB and matching range at t0 time, the left upper corner longitude PJ0 and latitude PW0 of the matching range are calculated:

[0051] PJ0 = J0 - ro0 / 2 × FA

[0052] PW0 = W0 - ra0 / 2 × FB

[0053] In the formula, ro0, ra0 are respectively the maximum value of the longitude direction matching range and the maximum value of the latitude direction matching range at t0 time;

[0054] 3.2) a two-dimensional array with ra0 rows and ro0 columns and all elements being 1 is defined, which is the prior probability density array P0 at t0 time;

[0055] 3.3) the map height in the matching range is stored in the height two-dimensional array mh0[i][j] with ra0 rows and ro0 columns, and the conditional probability density q(m0|i,j) of the measured value m0 at t0 time under the map height mh0[i][j] of the current sampling point is calculated:

[0056]

[0057] In the formula, 0≤i<ro0, 0≤j<ra0, i, j are all integers, ro0, ra0 are respectively the maximum value of the longitude direction matching range and the maximum value of the latitude direction matching range at t0 time; m0 = Q0 - R0, Q0 is the barometric height at t0 time, R0 is the radar height at t0 time; mh0[i][j] represents the map height of the jth row and the ith column sampling point in the current matching range, and σ is the standard deviation of the current measured height value and σ = 30;

[0058] 3.4) a two-dimensional array with ra0 rows and ro0 columns is defined, the conditional probability density q(m0|i,j) of each sampling point on the matching range is stored in the array, and the height measurement probability density array H0 at t0 time is obtained.

[0059] 3.5) Multiply the height probability density array H0 with the prior probability density array P0 to get the posterior probability density array PA0 of the first terrain match at time t0.

[0060] Step 4, perform terrain match calculation at time t n .

[0061] 4.1) According to the flight direction of the aircraft, i.e. the current time t n and the previous time t n-1 , the difference between the latitude and longitude of the inertial navigation estimated position, the terrain matching range of the previous time t n-1 is translated, and the posterior probability density array PA n-1 multiplied by the previous time t n-1 is obtained as the prior probability density array P n of the matching time t n , n is an integer and ≥ 0;

[0062] 4.2) According to the known set navigation error DE and sharp coefficient JR, the prior probability density array P n is convoluted to obtain the convolution probability density array C n :

[0063] 4.2.1) Define a two-dimensional convolution kernel array Cor with 3 rows and 3 columns and all elements being 1;

[0064] 4.2.2) Set the navigation error DE = 1, respectively, add 1 row of elements with value 0 on the top, bottom, left and right four sides of the posterior probability density array PA n-1 of the previous time t n-1 , n is an integer and ≥ 1;

[0065] 4.2.3) Weighted convolution of the posterior probability density array PA n-1 with the convolution kernel array Cor to obtain the prior probability density array P n in the matching range of the current time t n :

[0066] P n [i][j] = PA n-1 [i-1][j-1] + PA n-1 [i-1][j] + PA n-1 [i-1][j+1] + PA n-1 [i][j-1] + (1 + JR) PA n-1 [i][j] + PA n-1 [i][j+1] + PA n-1 [i+1][j-1] + PA n-1 [i+1][j] + PAn-1 [i+1][j+1]

[0067] Among them, P n [i][j] represents the prior probability density of the sampling point in the j-th row and i-th column within the matching range, JR = 100, 0 ≤ i < r n +2*DE, 0≤j<ra n +2*DE, where i and j are integers, ro n ra n Representing the current t respectively n The maximum value of the initial matching range in the longitude direction and the maximum value of the initial matching range in the latitude direction;

[0068] 4.2.4) For the prior probability density array P n Perform normalization and remove the outer array sv. n line su n The elements with a column value of 0 are used to obtain the convolution probability density array C. n The array C n There is a v n line u n Column, where:

[0069] u n =ro n +2*DE-su n ,

[0070] v n =ra n +2*DE-sv n ;

[0071] 4.3) According to t n Inertial navigation system predicts longitude J at time. m Inertial navigation system predicts latitude and position W m and the inertial navigation system's predicted position error E m Repeat step 2) to obtain the current t. n The terrain matching range at each time point is then compared with the previous time point t. n-1 Intersect the terrain matching ranges;

[0072] 4.4) Retrieve the regional elevation map of the intersection of the matching ranges for the current t n The height measurement at time t is calculated to obtain the current t. n The probability density array H of the height measurement at time t is n :

[0073] 4.4.1) Calculate t n The range of time-based height measurement matching:

[0074] The current t nMaximum value of the initial matching range in the longitude direction at any given time (ro) n With convolution probability density array C n The number of columns u n Compare the values ​​and take the smaller value as the current t. n Maximum value of the longitude-direction altimetry matching range co n The longitude range is 0 to 0°. n ;

[0075] The current t n The initial maximum value of the matching range in the latitude direction at any given time. n With convolution probability density array C n The number of rows v n Compare the values ​​and take the smaller value as the current t. n Maximum value of altimetry matching range ca in latitude direction at any given time n The latitude range is 0~ca n ;

[0076] 4.4.2) t n The elevation of the area on the map at that moment is entered into the CA. n co n Column height two-dimensional array mh n [i][j] in;

[0077] 4.4.3) For the two-dimensional height array mh n Calculate [i][j] to obtain t n The measured value m at time n The conditional probability density q(m) at the current map height of the sampling point n |i,j):

[0078]

[0079] In the formula, 0 ≤ i < co n , 0≤j<ca n where i and j are integers, and m n =Q n -R n Q n For t n At what moment is the air pressure altitude, R n For t n radar altitude at any given time; mh n [i][j] represents the map height of the sampling point in the j-th row and i-th column within the current matching range, and σ is the standard deviation of the current measured height value and σ = 30;

[0080] 4.4.4) Define a ca n co na two-dimensional array of columns, the conditional probability density q(m n | i,j) of each sampling point in the matching range is stored in the array, obtaining t n the height measurement probability density array H n at the current time t

[0081] 4.5) Corresponding multiplication is performed between the height measurement probability density array H n at the current time t n and the convolution probability density array C n , obtaining the posterior probability density array PA n at the current time t n :

[0082]

[0083]

[0084] wherein H n [i][j] is the height measurement probability density of the jth row ith column sampling point in the height measurement matching range at the time t n ,

[0085] C n [i][j] is the convolution probability density of the jth row ith column sampling point in the convolution matching range at the time t n ,

[0086] PA n [i][j] is the posterior probability density of the jth row ith column sampling point at the time t n ;

[0087] a is the longitude direction coordinate corresponding to the longitude of the upper left corner sampling point of the intersection part of the height measurement matching range and the convolution matching range at the current time t n ;

[0088] b is the latitude direction coordinate corresponding to the latitude of the upper left corner sampling point of the intersection part of the height measurement matching range and the convolution matching range at the current time t n ;

[0089] A = a + xo n , xo n is the smaller value between the maximum value co n of the longitude direction of the height measurement matching range at the current time t n and the column number u n of the convolution probability density array C n ;

[0090] B = b + xa n , xa n is the smaller value between the maximum value co n of the latitude direction of the height measurement matching range at the current time tMaximum value of altimetry matching range ca in latitude direction at any given time n With convolution probability density array C n The number of rows v n The smaller value in the range.

[0091] Step 5, calculate t m The estimated latitude and longitude, posterior standard deviation, confidence level, and standard deviation of terrain height for each time-based terrain matching are calculated and output to complete the terrain matching.

[0092] 5.1) Based on the two-dimensional height array mh from steps 3.3) and 4.4) m [i][j]Calculate t m Standard deviation of terrain elevation at time HS m :

[0093]

[0094] in, The height is a two-dimensional array mh m The height mean of [i][j], mhE m ca m ,co m The two-dimensional arrays mh and h represent the height respectively. m The rows and columns of [i][j], where m is 0 to n, and n is an integer ≥ 1;

[0095] 5.2) According to t m The posterior probability density array PA at time step m Calculate the estimated location longitude EJ for terrain matching. m and latitude EW m :

[0096] 5.2.1) Calculate t m The posterior probability density array PA at time step m The expected value of the longitudinal direction EX m Dimensional mathematical expectation EY m :

[0097]

[0098] 5.2.2) Calculate t m The posterior probability density array PA at time step m longitudinal variance EX 2 m Latitudinal variance EY 2 m :

[0099]

[0100]

[0101] where PA m is the posterior probability density of the sample point at the i-th column and j-th row at time t m ; g is the posterior probability density array PA m at time t m ; p is the column number of the posterior probability density array PA m at time t m ; and q is the row number of the posterior probability density array PA

[0102] 5.2.3) Add all the numbers in the posterior probability density array PA m at time t m together as a normalization factor s m , and combine the longitudinal mathematical expectation EX m and the latitudinal mathematical expectation EY m to calculate the longitudinal mean evx m and the latitudinal mean evy m :

[0103]

[0104] 5.2.4) According to the longitudinal mean evx m and the latitudinal mean evy m , calculate the longitudinal standard deviation dex m and the latitudinal standard deviation dey m :

[0105]

[0106] 5.2.5) According to the longitudinal mean evx m and the latitudinal mean evy m , calculate the estimated position longitude EJ m and latitude EW m at time t m :

[0107] EJ m = PJ m + evx m · FA

[0108] EW m = PW m + evy m · FB

[0109] where PJ m is the longitude of the upper-left sample point of the intersection part of the current t m time's height matching range and the convolution matching range, and PW m is the latitude of the upper-left sample point of the intersection part of the current t mThe latitude of the upper left corner sample point of the intersection part of the time instant height measurement matching range and the convolution matching range, FA and FB represent the longitude resolution of the height topographic map and the latitude resolution of the height topographic map respectively;

[0110] 5.3) Calculate t m The terrain matching posterior standard deviation BS at time instant m :

[0111] 5.3.1) Calculate t m The latitude W of the inertial navigation estimated position at time instant m The number of meters MJ represented by 1 degree of longitude on the same latitude m :

[0112]

[0113] In the formula, r is the radius of the earth, which is approximately 6371000 meters;

[0114] 5.3.2) Calculate t m The longitude J of the inertial navigation estimated position at time instant m The number of meters MW represented by 1 degree of latitude on the same longitude m :

[0115]

[0116] 5.3.3) Calculate t m According to the longitude standard deviation dex m， and the latitude standard deviation dey m The posterior standard deviation BS of the terrain matching predicted position at time instant m :

[0117]

[0118] In the formula, FA and FB represent the longitude resolution of the height topographic map and the latitude resolution of the height topographic map respectively;

[0119] 5.4) Calculate t m The confidence Z of terrain matching at time instant m :

[0120] 5.4.1) Calculate the probability density sum fm m of the posterior probability density array PA m centered on the longitude mean evx m and the latitude mean evy m within the longitude standard deviation dex m and the latitude standard deviation dey m :

[0121]

[0122] In the formula, PA m [k][l] represents the posterior probability density array PA. m The posterior probability density of row l and column k in the middle;

[0123] 5.4.2) Calculate the mean value evx along the longitudinal direction. m Latitudinal mean evy m Centered on the posterior probability density array PA within the range of two sampling points above and below it. m probability density and fz m :

[0124]

[0125] 5.4.3) According to fm m and fz m Calculate confidence level Z m :

[0126]

[0127] 5.5) Output t m Standard deviation of terrain elevation at time HS m Estimated location longitude EJ for terrain matching m Latitude EW m Posterior standard deviation (BS) m and confidence level Z m , complete t m Terrain matching at any given time.

[0128] The effects of this invention can be further illustrated by the following simulation results:

[0129] I. Simulation Content

[0130] Based on the steps outlined above in this example, a program is written to input flight data along a specific flight path of the aircraft. Terrain matching is performed every 1 second, and the predicted latitude and longitude of the terrain matching position are compared with the inertial navigation system's estimated latitude and longitude and the actual flight latitude and longitude. Figure 2 As shown; the posterior standard deviation of terrain matching prediction is as follows. Figure 4 As shown, the standard deviation of terrain height within the terrain matching range is as follows: Figure 5 As shown.

[0131] Will Figure 2 The difference in latitude and longitude is converted into a real distance difference, resulting in Figure 3 ;

[0132] from Figure 2 As can be seen, compared with the latitude and longitude estimated by the inertial navigation system, the latitude and longitude of the position predicted by the present invention are closer to the actual flight position.

[0133] fromFigure 3 As can be seen from the table, when the flight time of the aircraft reaches 3000s, the inertial navigation drift cumulative error, i.e. the distance between the inertial navigation estimated longitude and latitude and the true flight longitude and latitude, is as high as 3300m, while the predicted position error, i.e. the distance between the predicted position longitude and latitude and the true flight longitude and latitude, has been kept at about 100m after 200s, which more directly shows that the present application can reduce the inertial navigation drift cumulative error.

[0134] As can be seen from the table, when the flight time of the aircraft reaches 3000s, the inertial navigation drift cumulative error, i.e. the distance between the inertial navigation estimated longitude and latitude and the true flight longitude and latitude, is as high as 3300m, while the predicted position error, i.e. the distance between the predicted position longitude and latitude and the true flight longitude and latitude, has been kept at about 100m after 200s, which more directly shows that the present application can reduce the inertial navigation drift cumulative error. Figure 4 As can be seen from the table, after the aircraft has flown for 200s, the predicted posterior standard deviation is basically kept within 200m.

[0135] As can be seen from the table, after the aircraft has flown for 200s, the predicted posterior standard deviation is basically kept within 200m. Figure 5 As can be seen from the table, the terrain height standard deviation is mostly greater than 50m during the whole flight, which shows that the terrain of the region flown by the aircraft this time is relatively large.

[0136] As can be seen from the table, the terrain height standard deviation is mostly greater than 50m during the whole flight, which shows that the terrain of the region flown by the aircraft this time is relatively large. Figure 4 As can be seen from the table, the terrain height standard deviation is mostly greater than 50m during the whole flight, which shows that the terrain of the region flown by the aircraft this time is relatively large. Figure 5 As can be seen from the table, the terrain height standard deviation is mostly greater than 50m during the whole flight, which shows that the terrain of the region flown by the aircraft this time is relatively large.

[0137] The above description is only one specific example of the present application and does not constitute any limitation on the present application. Obviously, for those skilled in the art, after understanding the content and principles of the present application, various modifications and changes in form and details can be made without departing from the principles and structures of the present application, but these modifications and changes based on the idea of the present application are still within the protection scope of the claims of the present application.

Claims

1. A terrain matching aided navigation method without satellite navigation, characterized by, The method comprises the following steps: 1) the aircraft collects flight data and an elevation map resolution of its flight area, the flight data comprising an inertial navigation estimated longitude J, an inertial navigation estimated latitude W, an air pressure height Q, a radar height R, and an inertial navigation estimated position error E, and the elevation map resolution comprising a map longitude resolution FA and a map latitude resolution FB; 2) according to the resolution of the height topographic map, the current t m moment inertial navigation estimation position error E m and the inertial navigation estimation position longitude J m , latitude W m , the initial range of terrain matching is calculated, m is 0~n, n is an integer and ≥1; 3) performing terrain matching calculation at a first time t0: 3a) obtaining the area elevation map data within the terrain matching range from the memory of the flight control computer of the aircraft, and combining the inertial navigation estimated longitude J0 and latitude W0 at the time t0 to calculate the longitude PJ0 and latitude PW0 at the top left corner of the matching range; 3b) taking the numerical difference between the air pressure height Q0 and the radar height R0 at the time t0 as a measurement value m0, and combining the height value of the area elevation map and the height measurement error to perform height measurement calculation to obtain a height probability density array H0; 3c) multiplying the height probability density array H0 and a prior probability density array P0 of the same size but with elements of 1 to obtain a posterior probability density array PA0 of the first terrain matching at the time t0; 4) Perform subsequent t n Topographic matching computation at the moment: 4a) According to the current t n moment and the previous moment t n-1 , the difference between the inertial navigation estimated position longitude and latitude, the terrain matching range of the previous moment t n-1 is translated, and the product of the previous moment t n-1 is multiplied to obtain the posterior probability density array PA n-1 as the prior probability density array P n of the matching moment t n ; 4b) Convolve the a priori probability density array P with the known set of navigation error and sharpness coefficients to obtain a convolved probability density array C n n ;​ 4c) repeating step 2) to obtain a current t n terrain matching range, and intersecting it with the terrain matching range of the previous time instant t n-1 . 4d) retrieve the regional elevation map of the intersection of the matching range parts for the current t n time, resulting in the altimetric probability density array H n for the current t n time; 4e) multiply the current t n time height probability density array H n with the convolution probability density array C n to obtain the current t n time posterior probability density array PA n ; 5) Calculate the standard deviation of the map height in the matching range for step 3a) and step 4d) respectively, to obtain t m the standard deviation of the terrain height at the time instant HS m ; 6) Calculate the mean, mean square deviation of the arrays of posterior probability densities PA resulting from step 3c) and step 4e) respectively, and combine with the longitude and latitude of the upper left corner of the corresponding matching range to calculate t m m the estimated position longitude EJ m , latitude EW m , posterior standard deviation BS m and confidence Z m at the time of terrain matching;​ 7) output the results of steps 5) and 6) to complete t m topographic matching at the instant.

2. The method according to claim 1, wherein the initial range for which terrain matching is performed in step 2) is calculated as follows: 2a) Calculate the latitude of the position estimated by the inertial navigation W m Longitude expressed in meters on the same parallel JM m : 2b) calculating the latitude WM represented in meters on the same meridian as the inertial predicted position longitude J m latitude WM represented in meters on the same meridian as the inertial predicted position longitude J m : 2c) the inertial navigation estimated position error E m The position error ER is obtained by comparing the set error SE with the inertial navigation estimated position error E and taking the greater value m ; 2d) according to the position error ER m determining an initial range of the terrain match, i.e. a range in the longitude direction of 0 to ro m and a range in the latitude direction of 0 to ra m wherein: ro m = ER m × JM m ÷ FA x 2 ra m = ER m x WM m ÷ FB x 2 where r is the earth radius, approximately 6371000 meters, ro m , ra m are positive integers, and FA, FB represent the longitudinal resolution of the elevation topographic map and the latitudinal resolution of the elevation topographic map, respectively.

3. The method according to claim 1, wherein the longitude PJ0 and latitude PW0 at the top left corner of the matching range in step 3a) are calculated according to the following formulas: PJ0 = J0 - ro02 × FA PW0 = W0 - ra02 × FB wherein J0 and W0 are the inertial navigation estimated longitude and latitude at the time t0, respectively, and ro0 and ra0 are the maximum matching range values in the longitude direction and the latitude direction at the time t0, respectively; and FA and FB represent the longitude resolution and the latitude resolution of the elevation map, respectively.

4. The method according to claim 1, wherein the height measurement calculation in step 3b) is performed as follows: 3b1) placing the height of the area elevation map at the time t0 into a height two-dimensional array mh0[i][j] with ra0 rows and ro0 columns; 3b2) calculating the height two-dimensional array mh0[i][j] to obtain the conditional probability density q(m0|i,j) of the measurement value m0 at the map height mh0[i][j] of the current sampling point at the time t0: wherein 0≤i<ro0 and 0≤j<ra0, i and j are integers, ro0 and ra0 are the maximum matching range values in the longitude direction and the latitude direction at the time t0, respectively; m0 = Q0 - R0, Q0 is the air pressure height at the time t0, and R0 is the radar height at the time t0; mh0[i][j] represents the map height of the jth row and ith column sampling point within the current matching range, and σ is the standard deviation of the current measurement height value and σ = 30; 3b3) defining a two-dimensional array with ra0 rows and ro0 columns, and storing the conditional probability density q(m0|i,j) of each sampling point in the matching range into the array to obtain the height probability density array H0 at the time t0.

5. The method according to claim 1, wherein the convolution calculation in step 4b) is performed as follows: 4b1) define a two-dimensional convolution kernel array Cor of 3 rows and 3 columns and with elements all being 1; 4b2) set navigation error DE = 1 at the previous time t n-1 posterior probability density array PA n-1 each of the four sides of the upper, lower, left and right is increased by one row of elements with a value of 0, n is an integer and ≥ 1; 4b3) an array PA of posterior probabilities of the increasing elements n-1 a weighted convolution with the array of convolution kernels Cor, resulting in t n an array P of prior probabilities of the increasing elements in the matching range at the current time instant n : P n [i][j] = PA n-1 [i-1][j-1] + PA n-1 [i-1][j] + PA n-1 [i-1][j+1] + PA n-1 [i][j-1] + (1 + JR)PA n-1 [i][j] + PA n-1 [i][j+1] + PA n-1 [i+1][j-1] + PA n-1 [i+1][j] + PA n-1 [i+1][j+1] wherein P n [i][j] represents the prior probability density of the jth row and ith column sampling point in the matching range, JR is the sharpness coefficient and has a value of 100, 0≤i<ro n +2*DE, 0≤j<ra n +2*DE, i, j are integers, ro n , ra n respectively represent the maximum value of the initial matching range in the longitude direction and the maximum value of the initial matching range in the latitude direction at the current t n time. 4b4) the a priori probability density array P n is normalized and the outer array of sv is deleted n row su n elements with column value 0, resulting in the convolution probability density array C n , which array C n has v n rows and u n columns, u n = ro n + 2*DE-su n , v n = ra n + 2*DE-sv n .

6. The method of claim 1, wherein the altimeter calculation at the current t n is performed in step 4d) to obtain an altimeter probability density array H n at the current t n , which is implemented as follows: 4d1) Calculate t n Range of time-height matches: The current t n Maximum value of the initial matching range in the longitude direction at any given time (ro) n With convolution probability density array C n The number of columns u n Compare the values ​​and take the smaller value as the current t. n Maximum value of the longitude-direction altimetry matching range co n The longitude range is 0 to 0°. n ; The current t n The initial maximum value of the matching range in the latitude direction at any given time. n With convolution probability density array C n The number of rows v n Compare the values ​​and take the smaller value as the current t. n Maximum value of altimetry matching range ca in latitude direction at any given time n The latitude range is 0~ca n ; 4d2) putting the height of the area at the time instant into ca n height of the area at the time instant into ca n row co n column height two-dimensional array mh n [i][j] 4d3) calculating for the highly two-dimensional array mh n [i][j] the t n measurement value m n at the current sampling point map height mh n [i][j] the conditional probability density q(m n |i,j) : wherein 0≤i n , 0≤j n , i, j are integers, m n = Q n -R n , Q n is the air pressure height at time t n , R n is the radar height at time t n ; mh n [i][j] represents the map height of the jth row and ith column sampling point in the current matching range, and σ is the standard deviation of the current measured height value and σ = 30. 4d4) define a ca n row co n a two-dimensional array of columns, storing the conditional probability density of each sampling point on the matching range into the array, obtaining the H n height probability density array H n at time t.

7. The method of claim 1, wherein the a posteriori probability density array PA n at the current time instant t is obtained in step 4e) n : In the formula, H n [i][j] is the height measurement probability density of the i-th column and j-th row sampling point in the height measurement matching range at time t n n is an integer and ≥ 1. C n [i][j] is the convolution probability density of the i-th column and j-th row sample point in the matching range at time t n j-th row and i-th column sample point in the matching range at time t PA n [i][j] is the posterior probability density of the sample point at time t in row j and column i n [i][j] is the posterior probability density of the sample point at time t in row j and column i a is the current t n the longitude of the upper left corner sampling point of the intersection part of the height-moment matching range and the convolution matching range corresponds to the longitude direction coordinate on the entire height terrain map; b is the current t n the latitude direction coordinate corresponding to the latitude direction coordinate of the upper left corner sampling point of the intersection part of the current time height matching range and the convolution matching range on the whole height topographic map; A = a + xo n , xo n is the maximum value of the height matching range in the longitude direction at the current time t n co n is the smaller value of the column number u n of the convolution probability density array C n ​ B = b + xa n ,xa n For the current t n Maximum value of altimetry matching range ca in latitude direction at any given time n With convolution probability density array C n The number of rows v n The smaller value in the range.

8. The method of claim 1, wherein step 6) calculates t m the estimated position longitude EJ m , latitude EW m at the time of the terrain match, is accomplished as follows: 6a) compute t m the longitudinal mathematical expectation EX m of the array PA m of posterior probability densities at time t m and the latitudinal mathematical expectation EY 6b) compute t m the longitudinal variance EX m of the posterior probability density array PA 2 m the longitudinal variance EX 2 m : In the formula, PA m [i][j] is the posterior probability density of the sampling point at the i-th column and the j-th row at the t m time; g is the posterior probability density array PA m at the t m time; p is the column number of the posterior probability density array PA m at the t m time; m is 0~n, n is an integer and ≥1. 6c) add up all the numbers in the array PA m of posterior probability densities at time t m as a normalization factor s m , in combination with the longitudinal mathematical expectation EX m , the latitudinal mathematical expectation EY m , to calculate the longitudinal mean evx m , the latitudinal mean evy m : evx m = EX m / s m evy m = EY m / s m 6d) calculating the meridional standard deviation dex m , the weftwise standard deviation dey m m , the weftwise standard deviation dey m :​ dex m = sqrt(EX 2 m / s m -evx m 2 ) dey m = sqrt(EY 2 m / s m -evy m 2 ) 6e) Calculate the estimated position longitude EJ m , latitude EW m at time t m , based on the meridional mean evx m , the parallel mean evy m ​ EJ m = PJ m + evx m · FA EW m = PW m + evy m · FB In the formula, PJ m is the longitude of the upper left sampling point of the intersection part of the current t m is the latitude of the upper left sampling point of the intersection part of the current t m is the latitude of the upper left sampling point of the intersection part of the current t m FA and FB represent the longitude resolution and the latitude resolution of the height topographic map respectively.

9. The method of claim 1, wherein step 6) computes t m the terrain matching posterior standard deviation BS m as follows: 6f) Calculate t m Time and Inertial Navigation Estimated Latitude W m Meters on 1 degree of longitude MJ m : Wherein, r is the radius of the earth, approximately 6371000 meters, m is 0~n, n is an integer and ≥1; 6g) Calculate t m Time and Inertial Navigation Estimated Position Latitude J m Meters on the same meridian represented by 1 degree of latitude MW m : 6h) Calculate the meridional standard deviation dex m and the parallel standard deviation dey m Calculate t m the posterior standard deviation BS of the terrain-matching prediction position at time instant t m : Wherein, FA, FB represent the longitudinal resolution of the elevation topographic map and the latitudinal resolution of the elevation topographic map respectively.

10. The method of claim 1, wherein step 6) computes t m the confidence of the terrain match at time Z m is implemented as follows: 6i) calculating the equatorial mean evx m , the meridional standard deviation dex m , the standard deviation dex m , the standard deviation dex m the posterior probability density array PA m the probability density and fm m : where PA m [k][l] represents the posterior probability density array PA m where PA m is 0 to n, n is an integer and ≥ 1; 6j) Calculate the meridional mean evx m , the latitudinal mean evy m The posterior probability density array PA m with the probability density sum fz m : 6k) according to fm m and fz m calculate the confidence Z m :

Images