Space grid method for interconversion of geodetic height and atmospheric pressure height

Abstract

The invention discloses a space grid method for interconversion of geodetic height and barometric height, and belongs to the technical field of aviation navigation and geographic information systems. According to the method, the mutual conversion of geodetic height and barometric height is realized through establishing a global space-time grid reference system, establishing a space-time grid data base, converting barometric height to geodetic height, inversely converting the geodetic height to the barometric height, and carrying out reference conversion, data fusion and grid compression between the barometric heights. According to the method, the vertical availability of the height estimation value obtained through air pressure measurement within the height range within the stratosphere can be improved, the method is suitable for most weather conditions, the height estimation value with the precision reaching the ten-meter level can be achieved, and compared with an existing ISA-based method, the precision is higher, the application range is wider, and the application range is wider. The method can effectively improve the reliability of flight control systems of various aircrafts, and has wide application prospects in urban air traffic operation.

Description

Technical Field

[0001] The present invention belongs to the technical field of aviation navigation and geographic information system, and in particular relates to a spatial grid method for mutual conversion between geodetic height and atmospheric pressure height. Background Art

[0002] Currently, the majority of aircraft operating in urban air traffic (UAM) include general aviation helicopters, electric vertical takeoff and landing (EVTL) aircraft, and drones. These aircraft can achieve reliable, stable, and continuous navigation and positioning services horizontally via satellite, airborne, or ground-based navigation systems. However, in the vertical range, the uncertainty and accuracy limitations of common altitude measurement methods require a large separation between aircraft during low-altitude airspace management, severely impacting airspace operational efficiency and flight safety. Currently, in air traffic management, aircraft use two widely used methods for measuring altitude, each with a distinct reference base: the barometric method and the geometric method. The barometric method uses a barometer to measure altitude (i.e., pressure altitude). Based on the reference base provided by the ICAO Standard Atmosphere (ISA), it converts pressure readings taken at a specific location into altitude above mean sea level using a globally accepted conversion relationship between pressure and temperature. This provides a more intuitive representation for aircraft flight and airspace management. Its main drawback is that the altitude measured by pressure can vary with weather. Using ISA to measure altitude can result in errors exceeding hundreds of meters between warmer and colder climates. Existing temperature compensation schemes only correct for temperature differences between the ISA reference level and the aircraft's flight environment, failing to account for altitude variations caused by real-time pressure and temperature. They also ignore changes in pressure and temperature profiles, as well as the effects of adverse weather conditions such as weather fronts and severe convection on pressure altitude. Furthermore, current temperature compensation methods are only applicable below transition altitude and near airports. Therefore, above transition altitude, the deviation between ISA and the true atmospheric state results in a significant difference between atmospheric altitude and actual flight altitude. Furthermore, pilots must manually switch the atmospheric altitude reference from ISA to the airport during approach or takeoff (below transition altitude) to ensure the aircraft adheres to an isobaric surface. In contrast, geodetic height is primarily derived from Global Navigation Satellite System (GNSS) positioning, using the WGS84 or CGCS2000 ellipsoid as a reference, or from radio altimeter measurements, using sea level or the local terrain surface as a reference. GNSS-based altitude measurements are unbiased and can currently achieve an accuracy of several meters. Their advantage is that they are unaffected by atmospheric variations. However, their representation of altitude in flight is less intuitive, making them less convenient for airspace management.

[0003] Therefore, it is necessary to design a universal method for converting between geodetic altitude and pressure altitude, so as to provide sufficient performance to support all stages of aircraft flight and digital airspace management, especially in take-off and landing, approach and regional flight, as well as meet the flight performance needs of manned and unmanned aircraft in urban air traffic. This will greatly improve the efficiency of airspace use and management, and ensure the safety and order of high-density, heterogeneous aircraft integrated flight. Summary of the Invention

[0004] In view of this, the present invention discloses a spatial grid method for converting geoid height and air pressure height.

[0005] The present invention adopts the following technical solutions:

[0006] A spatial grid method for converting geodetic height and air pressure height, for use in urban air traffic, comprising the following steps:

[0007] S1. Establishing a global space-time grid reference system: Based on the Geodetic 2000 coordinate system, using the gridding of geographic space and time, a global space-time grid reference system is established in the plane, altitude, and time dimensions.

[0008] S2 establish spatiotemporal grid data base: based on the global spatiotemporal grid reference system, associated static and dynamic parameters for each grid;

[0009] S3. Pressure height to geoid conversion: Based on the global space-time grid reference system and the space-time grid data base, pressure height to geoid conversion is performed;

[0010] S4. Inverse conversion of geodetic height to pressure height: inverse calculation of pressure height and pressure value based on the geodetic height;

[0011] S5. Pressure-altitude base conversion: Based on the inverse conversion from geodetic height to pressure altitude, pressure-altitude base conversion is performed. The conversion includes pressure-altitude base conversion under standard atmospheric conditions and pressure-altitude base conversion under actual atmospheric conditions.

[0012] S6. Data fusion and grid compression: Based on the benchmark conversion, the conversion results are subjected to data fusion and grid compression. The model output geodetic height and the measured value of the global navigation satellite system are dynamically corrected through Kalman filtering. The conversion parameters are compressed into a binary data stream according to the grid coding rules and embedded into the avionics system database.

[0013] Furthermore, the plane dimension in step S1 includes: the plane range adopts a basic area grid size of 30′×30′, and the range is between 46-55 kilometers, while being compatible with the hierarchical division and mutual conversion of imaging and observation in the fields of meteorology and surveying and mapping; the altitude dimension includes: dividing the flight airspace into multiple layers of three-dimensional grids, and dynamically adjusting the grid resolution according to the altitude range: the low-altitude layer of 0-500 meters adopts a high-resolution grid, and the vertical interval of the high-resolution grid is ≤1 meter; the medium and high-altitude layer of 500 meters-20 kilometers adopts a medium and low-resolution grid, and the vertical interval of the medium and low-resolution grid is ≥10 meters; the time dimension includes: for static planning application scenarios, the hour level is adopted; for dynamic monitoring application scenarios, the minute level is adopted.

[0014] Furthermore, the generation and encoding method of the three-dimensional grid includes:

[0015] The earth space is divided into a quadtree structure, and each grid is encoded as:

[0016]

[0017] Code is the grid code, n is the level, n is a positive integer, b k ∈{0, 1, 2, 3} represents the partition index of level k.

[0018] Furthermore, the static parameters in step S2 include the geoid height of the grid center point, and the dynamic parameters include regional corrected sea level pressure, temperature gradient and humidity distribution.

[0019] Furthermore, the updating and fusion process of the dynamic parameters includes:

[0020] S2.1. Spatial interpolation of regional corrected sea level pressure. For non-grid center points between two adjacent grid cells, the geoid height is calculated using bilinear interpolation. For cross-grid flight scenarios, the regional corrected sea level pressure at the current point is calculated using inverse distance weighted interpolation.

[0021] S2.2. Dynamic parameter update uses sliding window filtering algorithm,

[0022] For meteorological data time series {QNH t , T t , RH t} Perform Kalman filtering, where t is the time, QNH t T is the corrected sea level pressure for the region at time t, t is the temperature at time t, RH t is the humidity at time t, and the prediction grid parameters are:

[0023]

[0024] Where, is the grid parameter predicted at time t+1, is the grid parameter at time t, QNH g Corrected pressure for the area at the center of the grid, is the temperature gradient, F is the state transfer matrix, B is the external influencing factor, u(t) is the external input, and w(t) is the noise disturbance.

[0025] Furthermore, step S3 includes:

[0026] S3.1 Inputting airborne sensor data and grid index code, wherein the airborne sensor data includes air pressure, temperature, and humidity;

[0027] S3.2 obtaining local parameters according to the grid index code: regional corrected air pressure at the grid center, geoid height at the grid center, and temperature gradient at the grid center;

[0028] S3.3 Calculate the virtual temperature of moist air:

[0029] Where e is the actual water vapor pressure, unit: Pa, RH is humidity, T is the virtual temperature,

[0030] Where q is the specific humidity, e is the actual water vapor pressure,

[0031] T v =T·(1+0.608·q), where T v is the virtual temperature of wet air, T is the virtual temperature, and q is the specific humidity; S3.4 introduces an equivalent gas constant in the virtual temperature correction to replace the standard gas constant to improve the humidity correction accuracy;

[0032] S3.5 obtains the regional corrected sea level pressure value and temperature gradient of the current grid, and calculates the corrected pressure altitude; queries the geoid height of the center point of the grid parameter grid, corrects the geoid, and calculates the geoid altitude.

[0033] Furthermore, step S4 includes:

[0034] S4.1 Input target geodetic height and grid parameters;

[0035] S4.2 Reverse thrust pressure altitude:

[0036] H QNH =H GNSS +N g , where H QNH To correct the pressure altitude, H GNSS is the earth's height, N g is the geoid height of the grid center point;

[0037] S4.3 Reverse calculation of air pressure P:

[0038] , where P is the measured air pressure, QNH g is the corrected air pressure for the area at the center of the grid, L is the temperature gradient, H QNH To correct the pressure altitude, T v is the virtual temperature, R is the gas constant of dry air, and g is the acceleration due to gravity.

[0039] Furthermore, the conversion between pressure altitudes under standard atmospheric conditions in step S5 includes: obtaining the difference between different altitudes based on the pressure difference between the surface pressure QFE, the corrected sea level pressure QNH and the standard atmospheric pressure QNE and the vertical pressure lapse rate under standard atmospheric conditions; according to the conversion relationship: surface pressure QFE = standard atmospheric pressure QNE - corrected sea level pressure QNH, combined with the airport elevation H S And the pressure lapse rate correction value β under the current temperature conditions, calculate the airport zero altitude: H0=H S -QFE×β, where H0 is the zero altitude of the airport, H S is the airport elevation, β is the pressure lapse rate correction value, and QFE is the surface air pressure.

[0040] Furthermore, the conversion between pressure-altitude references under actual atmospheric conditions in step S5 includes: using bilinear interpolation meteorological forecast data at the actual aircraft position, and comparing the aircraft's pressure measurement with the forecast data at the grid level for each altitude level, and using the surrounding grid data level as a reference surface for the corresponding level to determine the actual temperature decay rate: Where β is the temperature decay rate, m1 is the height grid of the vertical range, and T′ m1 is the temperature of the upper layer of the grid where the vertical range is located, T m1 is the temperature of the grid where the vertical range is located, h′ m1 The height of the previous grid layer where the vertical range is located, h m1 The height of the previous grid layer where the vertical range is located.

[0041] Furthermore, the Kalman filter dynamic correction in step S6 includes: In the formula is the output of the model after correction, K k is the Kalman gain, H model is the model output before correction, H GNSS is the geoid height; the conversion parameters include: the geoid height of the grid center point, the regional corrected air pressure and temperature gradient of the grid center point.

[0042] Beneficial effects of the present invention:

[0043] Through gridded local correction, the method of the present invention can reduce the meter-level error of traditional methods to the decimeter level (error ≤ 0.3 meters in low-altitude airspace), effectively improving accuracy. The method of the present invention uses grid indexing to reduce computational time by more than 50%, meeting the millisecond-level response requirements of urban air traffic and effectively improving the real-time performance of conversion. The method of the present invention can be seamlessly integrated into existing avionics systems, supporting data protocols such as ARINC 429 and AFDX, and has excellent compatibility.

[0044] This method has been verified by simulation. In a typical urban environment, the standard deviation of the height conversion error is 0.21 meters, which is significantly better than the traditional ISA model (1.5 meters).

[0045] The present invention provides a universal spatial grid application method for mutual conversion between geodetic altitude and pressure altitude, which can improve the vertical usability of altitude estimates derived from pressure measurements within the stratosphere. Through two-way comparative verification and analysis of historical observation data and actual flight data, this method is applicable to most weather conditions and can achieve altitude estimates with an accuracy of up to ten meters. Compared with existing ISA-based methods, it has higher accuracy, a wider range of adaptability, and reduces altitude errors by multiple orders of magnitude. This method can effectively improve the reliability of flight control systems for various types of aircraft, has high scalability for manned and unmanned aircraft, and can provide a more feasible solution for vertical RNP and urban air traffic operations based on pressure-measured altitude. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0047] Figure 1 It is a flow chart of a spatial grid method for converting geoid height and air pressure height according to the present invention. DETAILED DESCRIPTION

[0048] In order to better understand the technical solution of the present invention, the embodiments of the present invention are described in detail below with reference to the accompanying drawings.

[0049] It should be understood that the embodiments described are only a portion of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by persons of ordinary skill in the art without creative work are within the scope of protection of the present invention.

[0050] Example 1

[0051] A spatial grid method for converting geodetic height and pressure height, such as Figure 1 As shown, for urban air traffic, the method includes the following steps:

[0052] S1. Establish a global space-time grid reference system: Based on the Geodetic 2000 Coordinate System (CGCS-2000), a global space-time grid reference system is established in the plane, altitude, and time dimensions using the gridding of geographic space and time.

[0053] S2 establish spatiotemporal grid data base: based on the global spatiotemporal grid reference system, associated static and dynamic parameters for each grid;

[0054] S3. Pressure height to geoid conversion: Based on the global space-time grid reference system and the space-time grid data base, pressure height to geoid conversion is performed;

[0055] S4. Inverse conversion of geodetic height to pressure height: inverse calculation of pressure height and pressure value based on the geodetic height;

[0056] S5. Pressure-altitude base conversion: Based on the inverse conversion from geodetic height to pressure altitude, pressure-altitude base conversion is performed. The conversion includes pressure-altitude base conversion under standard atmospheric conditions and pressure-altitude base conversion under actual atmospheric conditions.

[0057] S6. Data fusion and grid compression: Based on the benchmark conversion, the conversion results are subjected to data fusion and grid compression. The model output geodetic height and the measured value of the global navigation satellite system are dynamically corrected through Kalman filtering. The conversion parameters are compressed into a binary data stream according to the grid coding rules and embedded into the avionics system database.

[0058] Furthermore, the plane dimension in step S1 includes: the plane range adopts a basic area grid size of 30′×30′, the range is between 46-55 kilometers, and is compatible with the hierarchical division and mutual conversion of imaging and observation in the fields of meteorology and surveying and mapping, 1°=60′; the altitude dimension includes: dividing the flight airspace into multiple layers of three-dimensional grids, and dynamically adjusting the grid resolution according to the altitude range: the low-altitude layer of 0-500 meters adopts a high-resolution grid, and the vertical spacing of the high-resolution grid is ≤1 meter; the medium and high-altitude layer of 500 meters-20 kilometers adopts a medium and low-resolution grid, and the vertical spacing of the medium and low-resolution grid is ≥10 meters; the time dimension includes: for static planning application scenarios, the hour level is adopted; for dynamic monitoring application scenarios, the minute level is adopted.

[0059] Furthermore, the generation and encoding method of the three-dimensional grid includes:

[0060] The earth space is divided into a quadtree structure, and each grid is encoded as:

[0061]

[0062] Code is the grid code, n is the level, n is a positive integer, b k ∈{0,1,2,3} represents the partition index of level k.

[0063] Furthermore, the static parameters in step S2 include the geoid height of the grid center point, and the dynamic parameters include regional corrected sea level pressure, temperature gradient and humidity distribution.

[0064] Furthermore, the updating and fusion process of the dynamic parameters includes:

[0065] S2.1. Spatial interpolation of regional corrected sea level pressure,

[0066] For non-grid center points between two adjacent grids, the geoid height is calculated by bilinear interpolation:

[0067]

[0068] Where N(x,y) is the geoid height of the location, x is the horizontal distance of the location relative to the center of the grid, y is the vertical distance of the location relative to the center of the grid, and x is the vertical distance of the location relative to the center of the grid. i is the horizontal distance between the center point of the grid and the boundary of the grid, y j is the vertical distance between the center point of the grid and the boundary of the grid, N ij is the geoid height of the grid center point;

[0069] For cross-grid flight scenarios, the corrected sea level pressure in the current point area is calculated by inverse distance weighted interpolation:

[0070]

[0071] Where QNH c Corrected sea level pressure for the current point area, QNH i Corrected sea level pressure for the area of ​​the i-th neighboring grid, d i is the distance from the current point to the center of the adjacent grid, k is the number of adjacent grids, and k is usually 4;

[0072] S2.2. Dynamic parameter update uses a sliding window filtering algorithm, where i is the horizontal coordinate of the grid or adjacent grid, and j is the vertical coordinate of the grid or adjacent grid.

[0073] For meteorological data time series {QNH t , T t, RH t} Perform Kalman filtering, where t is the time, QNH t T is the corrected sea level pressure for the region at time t, t is the temperature at time t, RH t is the humidity at time t, and the prediction grid parameters are:

[0074]

[0075] Where, is the grid parameter predicted at time t+1, is the grid parameter at time t, QNH g Corrected pressure for the area at the center of the grid, is the temperature gradient, F is the state transfer matrix, B is the external influencing factor, u(t) is the external input, and w(t) is the noise disturbance.

[0076] Furthermore, step S3 includes:

[0077] S3.1 Inputting airborne sensor data and grid index code, wherein the airborne sensor data includes air pressure, temperature, and humidity;

[0078] S3.2 obtaining local parameters according to the grid index code: regional corrected air pressure at the grid center, geoid height at the grid center, and temperature gradient at the grid center;

[0079] S3.3 Calculate the virtual temperature of moist air:

[0080] Where e is the actual water vapor pressure, unit: Pa, RH is humidity, T is the virtual temperature,

[0081] Where q is the specific humidity, e is the actual water vapor pressure,

[0082] T v =T·(1+0.608·q), where T v is the virtual temperature of wet air, T is the virtual temperature, and q is the specific humidity; S3.4 introduces an equivalent gas constant in the virtual temperature correction to replace the standard gas constant to improve the humidity correction accuracy;

[0083] S3.5 obtains the corrected sea level pressure and temperature gradient for the current grid area and calculates the corrected pressure altitude:

[0084] Where L is the temperature gradient, p is the measured air pressure (Pa), QNH gis the corrected air pressure for the area at the center of the grid, R is the gas constant for dry air (287.05 J / (kg·K)), g is the acceleration due to gravity (9.80665 m / s 2 ), is the temperature gradient at the center of the grid, where R = 287.05;

[0085] Query the geoid height of the center point of the grid parameter grid, correct the geoid, and calculate the geoid height.

[0086] Furthermore, step S4 includes:

[0087] S4.1 Input target geodetic height and grid parameters;

[0088] S4.2 Reverse thrust pressure altitude:

[0089] H QNH =H GNSS +N g , where H QNH To correct the pressure altitude, H GBSS is the earth's height, N g is the geoid height of the grid center point;

[0090] S4.3 Reverse calculation of air pressure P:

[0091] , where P is the measured air pressure, QNH g is the corrected air pressure for the area at the center of the grid, L is the temperature gradient, H QNH To correct the pressure altitude, T v is the virtual temperature, R is the gas constant of dry air, and g is the acceleration due to gravity.

[0092] Furthermore, the conversion between pressure altitudes under standard atmospheric conditions in step S5 includes: obtaining the difference between different altitudes based on the pressure difference between the surface pressure QFE, the corrected sea level pressure QNH and the standard atmospheric pressure QNE and the vertical pressure lapse rate under standard atmospheric conditions; according to the conversion relationship: surface pressure QFE = standard atmospheric pressure QNE - corrected sea level pressure QNH, combined with the airport elevation H S And the pressure lapse rate correction value β under the current temperature conditions, calculate the airport zero altitude: H0=H S -QFE×β, where H0 is the zero altitude of the airport, H S is the airport elevation, β is the pressure lapse rate correction value, and QFE is the surface air pressure.

[0093] Furthermore, the conversion between pressure-altitude references under actual atmospheric conditions in step S5 includes: using bilinear interpolation meteorological forecast data at the actual aircraft position, and comparing the aircraft's pressure measurement with the forecast data at the grid level for each altitude level, and using the surrounding grid data level as a reference surface for the corresponding level to determine the actual temperature decay rate: Where β is the temperature decay rate, m1 is the height grid of the vertical range, and T′ m1 is the temperature of the upper layer of the grid where the vertical range is located, T m1 is the temperature of the grid where the vertical range is located, h′ m1 The height of the previous grid layer where the vertical range is located, h m1 The height of the previous grid layer where the vertical range is located.

[0094] Furthermore, the Kalman filter dynamic correction in step S6 includes: In the formula is the output of the model after correction, K k is the Kalman gain, H model is the model output before correction, H GNSS is the geoid height; the conversion parameters include: the geoid height of the grid center point, the regional corrected air pressure and temperature gradient of the grid center point.

[0095] Example 2

[0096] Construct a dynamic space-time grid model.

[0097] The urban airspace is divided into multiple layers of three-dimensional grids (based on the global space-time grid framework), and the following parameters are pre-stored or calculated in real time within each grid:

[0098] Static parameters: geoid height Ng, benchmark terrain elevation;

[0099] Dynamic parameters: regional pressure reference (QNH), temperature gradient Humidity distribution;

[0100] Grid coding: Compresses and stores data using spatiotemporal grid coding rules, supporting fast indexing and transmission.

[0101] By integrating the local physical model within the grid with real-time data, bidirectional conversion between geodetic height and atmospheric pressure height is achieved. The specific process includes:

[0102] Forward conversion: air pressure height → virtual temperature correction → QNH adjustment → geoid correction → geoid height;

[0103] Reverse conversion: geoid height → geoid correction → QNH inverse → pressure height.

[0104] Example 3

[0105] A spatial grid method for converting geodetic height and pressure height, such as Figure 1 As shown, for urban air traffic, the method includes the following steps:

[0106] Step S1 establishes a global space-time grid reference system.

[0107] Based on the Geodetic 2000 Coordinate System (CGCS-2000), the grid division principles of geographic space and time are used to establish a global space-time grid reference system from the plane, altitude and time dimensions. The plane range adopts a basic area grid size of 30′×30′, which is between 46-55 kilometers. At the same time, it is compatible with the hierarchical division and mutual conversion of imaging and observation in professional fields such as meteorology and surveying and mapping. In terms of altitude range, the flight airspace is divided into multiple layers of three-dimensional grids, and the levels are dynamically adjusted according to the altitude range: the low-altitude layer (0-500 meters) adopts a high-resolution grid (vertical interval ≤1 meter), and the medium- and high-altitude layer (500 meters-20 kilometers) adopts a medium- and low-resolution grid (vertical interval ≥100 meters). In terms of time dimension, for application scenarios such as static planning, the hour level is adopted; for application scenarios such as dynamic monitoring, the minute level is adopted. The three-dimensional grid generation and encoding method is as follows:

[0108] The earth space is divided into a quadtree structure, and each grid is encoded as:

[0109]

[0110] Where n is the level, n is a positive integer, b k ∈{0, 1, 2, 3} represents the partition index of level k.

[0111] Step S2 establishes a spatiotemporal grid data base.

[0112] The following static and dynamic parameters are associated with each grid. Pre-stored static parameters include the geoid height N of the grid center point. g , which can be pre-calculated using the EGM-2008 model spherical harmonic coefficients and stored as floating-point arrays; real-time update of dynamic parameters including regional corrected sea level pressure (QNH), temperature gradient Humidity distribution can be obtained through the meteorological data link and stored as a time series queue. The dynamic parameter update and fusion process is as follows:

[0113] Step S2.1: QNH space interpolation.

[0114] For non-grid center points between two adjacent grids, the geoid height is calculated by bilinear interpolation:

[0115]

[0116] Where N(x,y) is the geoid height of the location, x is the horizontal distance of the location relative to the center of the grid, y is the vertical distance of the location relative to the center of the grid, and x is the vertical distance of the location relative to the center of the grid. i is the horizontal distance between the center point of the grid and the boundary of the grid, y j is the vertical distance between the center point of the grid and the boundary of the grid, N ij is the geoid height of the grid center point.

[0117] For cross-grid flight scenarios, the inverse distance weighted interpolation (IDW) is used to calculate the QNH of the current point:

[0118] Where QNH c Corrected sea level pressure for the current point area, QNH i Corrected sea level pressure for the area of ​​the i-th neighboring grid, d i is the distance from the current point to the center of the adjacent grid, k is the number of adjacent grids (usually 4), i is the horizontal coordinate of the grid or the adjacent grid, and j is the vertical coordinate of the grid or the adjacent grid.

[0119] Step S2.2: Dynamic parameter update adopts sliding window filtering algorithm.

[0120] For meteorological data time series {QNH t , T t , RH t} Perform Kalman filtering, where t is the time, QNH t T is the corrected sea level pressure for the region at time t, t is the temperature at time t, RH t is the humidity at time t, and the prediction grid parameters are:

[0121]

[0122] Where, is the grid parameter predicted at time t+1, is the grid parameter at time t, QNH g Corrected pressure for the area at the center of the grid, is the temperature gradient, F is the state transfer matrix, B is the external influencing factor, u(t) is the external input (such as weather forecast, etc.), and w(t) is the noise disturbance.

[0123] Step S3: Forward conversion from pressure altitude to geodetic altitude.

[0124] Step S3.1: Input airborne sensor data (air pressure P, temperature T, humidity RH) and grid index code;

[0125] Step S3.2 obtains local parameters based on the grid index: the regional corrected pressure QNH at the grid center point g , geoid height N of the grid center point g and the temperature gradient at the center of the grid

[0126] Step S3.3 calculates the wet air virtual temperature T according to the following formula: v :

[0127] Where e is the actual water vapor pressure, unit: Pa, RH is the humidity, T is the virtual temperature;

[0128] Where q is the specific humidity and e is the actual water vapor pressure;

[0129] T v =T·(1+0.608·q), where T v is the virtual temperature of moist air, T is the virtual temperature, and q is the specific humidity;

[0130] Step S3.4 introduces the equivalent gas constant R in the virtual temperature correction eff , used to replace the standard R dry To improve humidity correction accuracy.

[0131] Where R eff is the equivalent gas constant, R dry is the standard gas constant, q is the specific humidity;

[0132] Step S3.5 obtains the QNH value QNH of the current grid g and temperature gradient L, calculate the corrected pressure altitude H QNH :

[0133]

[0134] Where L is the temperature gradient, p is the measured air pressure (Pa), QNH g is the corrected air pressure for the area at the center of the grid, R is the gas constant for dry air (287.05 J / (kg·K)), g is the acceleration due to gravity (9.80665 m / s 2 ), is the temperature gradient at the center of the grid, where R = 287.05;

[0135] Query grid parameter N g , correct the geoid and calculate the geoid height H GNSS :

[0136] H GNSS =H QNH -Ng ,

[0137] Where H GNSS is the earth's height, H QNN To correct the pressure altitude, N g is the geoid height of the grid center point.

[0138] Step S4: reverse conversion from geoid altitude to pressure altitude.

[0139] Step S4.1 Input target earth height H GNSS and grid parameters;

[0140] Step S4.2 Reverse pressure height H QNH :

[0141] H QNH =H GNSS +N g , where H QNH To correct the pressure altitude, H GNSS is the earth's height, N g is the geoid height of the grid center point;

[0142] Step S4.3: Back-calculate the air pressure value P:

[0143] , where P is the measured air pressure (Pa), QNH g is the corrected air pressure of the area at the center of the grid, L is the temperature gradient (0.0065K / m), H QNH To correct the pressure altitude, T v is the virtual temperature, R is the gas constant of dry air (287.05 J / (kg·K)), g is the acceleration due to gravity (9.80665 m / s 2 ).

[0144] Step S5: switching between high and low pressure references.

[0145] Step S5.1: conversion between pressure and altitude under standard atmospheric conditions;

[0146] Based on the pressure difference between the surface pressure QFE, the corrected sea level pressure QNH and the standard atmospheric pressure QNE and the vertical pressure lapse rate under standard atmospheric conditions -8.25m / hPa, the difference between different heights is calculated:

[0147] Δh=Δp×8.25m / hPa, where Δp is the pressure difference and Δh is the height difference;

[0148] Based on the conversion relationship QFE=QNE-QNH, combined with the airport elevation H S And the pressure lapse rate correction value β under the current temperature conditions, the airport zero altitude can be calculated:

[0149] H0=H S -QFE×β, where H0 is the zero altitude of the airport, H S is the airport elevation, β is the pressure lapse rate correction value, and QFE is the surface air pressure.

[0150] Step S5.2: conversion between pressure and altitude under actual atmospheric conditions;

[0151] Using bilinearly interpolated weather forecast data at the actual aircraft location, the aircraft's pressure measurements are compared to the forecasted pressure data at the relevant grid level for each altitude level, using the surrounding grid data levels as reference surfaces for the corresponding levels to determine the actual temperature decay rate:

[0152] Where β is the temperature decay rate, m1 is the height grid of the vertical range, and T′ m1 is the temperature of the upper layer of the grid where the vertical range is located, T m1 is the temperature of the grid where the vertical range is located, h′ m1 The height of the previous grid layer where the vertical range is located, h m1 The height of the previous grid layer where the vertical range is located.

[0153] Step S6: grid coding and data fusion compression.

[0154] Output the model H GNSS Dynamic correction with GNSS measured values ​​through Kalman filtering:

[0155] In the formula is the output of the model after correction, K k is the Kalman gain, H model is the model output before correction, H GNSS is the ground height. Kalman gain K k Updated according to the model error covariance matrix.

[0156] Convert the parameters It is compressed into a binary data stream according to grid coding rules and embedded into the avionics system database.

[0157] Example 4

[0158] The present invention is related to mathematical analysis and experimental verification.

[0159] 1. Error transfer analysis of virtual temperature correction.

[0160] Assuming the temperature sensor error is ΔT, the virtual temperature error propagation is:

[0161]

[0162] Experiments show that when ΔT = 1K, the pressure height error ΔH QNH ≤0.3 meters.

[0163] 2. The impact of grid resolution on accuracy.

[0164] At low altitudes (≤1 meter resolution), the geoid height interpolation error is:

[0165] ∈ N =max(|N g -N 真实 |)≤0.1m,

[0166] The error of traditional spherical harmonic function calculation (order 2190) in the same area is 0.3 meters.

[0167] 3. Error analysis of height conversion model.

[0168] QNH correction error: mainly comes from the estimated deviation of the temperature gradient L. Let ΔL be the error, then the height error is:

[0169]

[0170] Humidity sensitivity: to virtual temperature T v The partial derivative of is:

[0171]

[0172] 4. Real-time optimization.

[0173] Grid preloading: EGM-2008 spherical harmonic coefficients are pre-stored by grid partitions. Only local coefficients are loaded during calculation. The time complexity can be reduced from O(n 2 ) is reduced to O(m 2 )(m<<n).

[0174] Parallel computing: Using GPU to calculate the spherical harmonic function ∑ n，m C nm P nm Perform parallel summation.

[0175] The embodiments of the present invention are described in detail above. Specific examples are used herein to illustrate the principles and implementation methods of the present invention. The description of the above embodiments is only used to help understand the method of the present invention and its core ideas. At the same time, for those skilled in the art, according to the ideas of the present invention, there may be changes in the specific implementation methods and application scopes. In summary, the contents of this specification should not be understood as limiting the present invention.

Claims

1. A spatial grid method for converting geodetic height and atmospheric pressure height for urban air traffic, characterized in that: The method comprises the following steps: S1. Establishing a global space-time grid reference system: Based on the Geodetic 2000 coordinate system, using the gridding of geographic space and time, a global space-time grid reference system is established in the plane, altitude, and time dimensions. S2 establish spatiotemporal grid data base: based on the global spatiotemporal grid reference system, associated static and dynamic parameters for each grid; S3. Pressure height to geoid conversion: Based on the global space-time grid reference system and the space-time grid data base, pressure height to geoid conversion is performed; S4. Inverse conversion of geodetic height to pressure height: inverse calculation of pressure height and pressure value based on the geodetic height; S5. Pressure-altitude base conversion: Based on the inverse conversion from geodetic height to pressure altitude, pressure-altitude base conversion is performed. The conversion includes pressure-altitude base conversion under standard atmospheric conditions and pressure-altitude base conversion under actual atmospheric conditions. S6. Data fusion and grid compression: Based on the benchmark conversion, the conversion results are subjected to data fusion and grid compression. The model output geodetic height and the measured value of the global navigation satellite system are dynamically corrected through Kalman filtering. The conversion parameters are compressed into a binary data stream according to the grid coding rules and embedded into the avionics system database.

2. The method according to claim 1, characterized in that The plane dimension in step S1 includes: the plane range adopts a basic area grid size of 30′×30′, and the range is between 46-55 kilometers, while being compatible with the hierarchical division and mutual conversion of imaging and observation in the fields of meteorology and surveying and mapping; the altitude dimension includes: dividing the flight airspace into multiple layers of three-dimensional grids, and dynamically adjusting the grid resolution according to the altitude range: the low-altitude layer of 0-500 meters adopts a high-resolution grid, and the vertical interval of the high-resolution grid is ≤1 meter; the medium and high-altitude layer of 500 meters-20 kilometers adopts a medium- and low-resolution grid, and the vertical interval of the medium- and low-resolution grid is ≥10 meters; the time dimension includes: for static planning application scenarios, the hour level is adopted; for dynamic monitoring application scenarios, the minute level is adopted.

3. The method according to claim 2, characterized in that The generation and encoding method of the three-dimensional grid includes: The earth space is divided into a quadtree structure, and each grid is encoded as: Code is the grid code, n is the level, n is a positive integer, b k ∈{0, 1, 2, 3} represents the partition index of level k.

4. The method according to claim 3, characterized in that The static parameters in step S2 include the geoid height of the grid center point, and the dynamic parameters include the regional corrected sea level pressure, temperature gradient, and humidity distribution.

5. The method according to claim 4, characterized in that: The updating and fusion process of the dynamic parameters includes: S2.

1. Spatial interpolation of regional corrected sea level pressure. For non-grid center points between two adjacent grid cells, the geoid height is calculated using bilinear interpolation. For cross-grid flight scenarios, the regional corrected sea level pressure at the current point is calculated using inverse distance weighted interpolation. S2.

2. Dynamic parameter update uses sliding window filtering algorithm, For meteorological data time series {QNH t , T t , RH t } Perform Kalman filtering, where t is the time, QNH t T is the corrected sea level pressure for the region at time t, t is the temperature at time t, RH t is the humidity at time t, and the prediction grid parameters are: Where, is the grid parameter predicted at time t+1, is the grid parameter at time t, QNH g Corrected pressure for the area at the center of the grid, is the temperature gradient, F is the state transfer matrix, B is the external influencing factor, u(t) is the external input, and w(t) is the noise disturbance.

6. The method according to claim 5, characterized in that Step S3 includes: S3.1 Inputting airborne sensor data and grid index code, wherein the airborne sensor data includes air pressure, temperature, and humidity; S3.2 obtaining local parameters according to the grid index code: regional corrected air pressure at the grid center, geoid height at the grid center, and temperature gradient at the grid center; S3.3 Calculate the virtual temperature of moist air: , where e is the actual water vapor pressure, unit: Pa, RH is humidity, T is the virtual temperature, Where q is the specific humidity, T v =T·(1+0.608·q), where T v is the virtual temperature of moist air, q is the specific humidity; S3.4 introduces an equivalent gas constant in virtual temperature correction to replace the standard gas constant to improve humidity correction accuracy; S3.5 obtains the regional corrected sea level pressure value and temperature gradient of the current grid, and calculates the corrected pressure altitude; queries the geoid height of the center point of the grid parameter grid, corrects the geoid, and calculates the geoid altitude.

7. The method according to claim 6, characterized in that Step S4 includes: S4.1 Input target geodetic height and grid parameters; S4.2 Reverse thrust pressure altitude: H QNH =H GNSS +N g , where H QNH To correct the pressure altitude, H GNSS is the earth's height, N g is the geoid height of the grid center point; S4.3 Reverse calculation of air pressure P: , where P is the measured air pressure, QNH g is the corrected air pressure for the area at the center of the grid, L is the temperature gradient, H QNH To correct the pressure altitude, T v is the virtual temperature, R is the gas constant of dry air, and g is the acceleration due to gravity.

8. The method according to claim 7, characterized in that: The conversion between pressure altitudes under standard atmospheric conditions in step S5 includes: obtaining the difference between different altitudes based on the pressure difference between the surface pressure QFE, the corrected sea level pressure QNH and the standard atmospheric pressure QNE and the vertical pressure lapse rate under standard atmospheric conditions; according to the conversion relationship: surface pressure QFE = standard atmospheric pressure QNE - corrected sea level pressure QNH, combined with the airport elevation H S And the pressure lapse rate correction value β under the current temperature conditions, calculate the airport zero altitude: H0=H S -QFE×β, where H0 is the zero altitude of the airport, H S is the airport elevation, β is the pressure lapse rate correction value, and QFE is the surface air pressure.

9. The method according to claim 8, characterized in that The conversion between pressure and altitude references under actual atmospheric conditions in step S5 includes: using bilinear interpolation meteorological forecast data at the actual aircraft position, and comparing the aircraft's pressure measurement with the forecast data at the grid level for each altitude level, and using the surrounding grid data level as a reference surface for the corresponding level to determine the actual temperature decay rate: Where β is the temperature decay rate, m1 is the height grid of the vertical range, and T′ m1 is the temperature of the upper layer of the grid where the vertical range is located, T m1 is the temperature of the grid where the vertical range is located, h′ m1 The height of the previous grid layer where the vertical range is located, h m1 The height of the previous grid layer where the vertical range is located.

10. The method according to claim 9, characterized in that: The Kalman filter dynamic correction in step S6 includes: In the formula is the output of the model after correction, K k is the Kalman gain, H model is the model output before correction, H GNSS is the geoid height; the conversion parameters include: the geoid height of the grid center point, the regional corrected air pressure and temperature gradient of the grid center point.

Images