Three-dimensional wind speed and direction detection method and system based on GNSS tomography inversion

Abstract

The invention relates to the technical field of meteorological remote sensing, in particular to a three-dimensional wind speed and wind direction detection method and system based on GNSS tomography inversion, and the method comprises the steps: building a three-dimensional water vapor density tomography inversion model through obtaining observation data and meteorological observation data of a plurality of GNSS observation stations, obtaining water vapor density distribution of different height layers, and obtaining a three-dimensional wind speed and wind direction detection model; further obtaining the PWV (atmospheric precipitable water) of each height layer; a three-dimensional wind field inversion algorithm is established by utilizing a coupling relation between water vapor transportation and a wind field and combining a spatial position of a GNSS observation station and spatial and temporal variation characteristics of PWV of different height layers, so that quantitative estimation of wind speeds and wind directions of the different height layers is realized, and three-dimensional wind speed and wind direction quantitative detection based on GNSS tomography inversion can be realized. And important observation data support is provided for meteorological monitoring and numerical weather forecast.

Description

Technical Field

[0001] The present invention relates to the field of meteorological remote sensing technology, and in particular to a three-dimensional wind speed and direction detection method and system based on GNSS tomographic inversion. Background Art

[0002] Atmospheric water vapor is a crucial component of Earth's atmospheric system, and its spatiotemporal distribution is closely related to atmospheric motion. In atmospheric dynamics, water vapor transport is primarily controlled by wind fields, and there is a clear correlation between the two. In particular, the spatiotemporal variation of atmospheric precipitable water (PWV), an integral of water vapor content, is closely related to wind field variations, providing a theoretical basis for inferring wind fields based on water vapor.

[0003] Currently, wind field detection relies primarily on sounding balloons, wind profiler radars, meteorological satellites, and numerical weather prediction models. However, these methods all have varying degrees of limitations in terms of temporal resolution, spatial coverage, observation continuity, and adaptability to meteorological conditions. Sounding balloon observations are infrequent and offer insufficient spatial coverage; wind profiler radars are expensive to deploy and have limited spatial representation; meteorological satellites have poor vertical resolution and are significantly affected by cloud interference; and while numerical weather prediction models offer comprehensive coverage, their forecast accuracy is limited by the initial field and parameterization schemes, making it difficult to accurately reflect local wind field characteristics.

[0004] The application of global navigation satellite systems in atmospheric remote sensing provides a new approach for wind field detection. GNSS signals are all-weather, unaffected by weather conditions such as clouds and precipitation, enabling continuous and stable observations. Furthermore, utilizing existing GNSS infrastructure for atmospheric detection avoids the duplication of specialized equipment, offering significant cost advantages. Three-dimensional water vapor tomography technology enables high-temporal and spatial resolution atmospheric water vapor inversion, providing a reliable data foundation for wind field estimation.

[0005] However, the existing GNSS-based three-dimensional wind field detection technology is still in the exploratory stage. There is an urgent need to develop a high-precision three-dimensional wind field detection method based on GNSS three-dimensional water vapor density inversion, establish a more accurate parameterized tomography inversion model, and fully consider the temporal and spatial distribution characteristics of water vapor. This has important scientific significance and application value for improving the atmospheric wind field observation capability, improving the accuracy of numerical weather forecasts, and enhancing the extreme weather monitoring and early warning capabilities. Summary of the Invention

[0006] The main purpose of the present invention is to provide a three-dimensional wind speed and direction detection method and system based on GNSS tomographic inversion to solve the above technical problems.

[0007] In order to achieve the above object, the present invention provides a three-dimensional wind speed and direction detection method based on GNSS tomographic inversion, comprising the steps of: S1, obtaining GNSS observation data of multiple GNSS observation stations in the area to be measured and meteorological observation data of the corresponding GNSS observation stations; S2, divide the area to be measured into three-dimensional grids and establish the three-dimensional water vapor density tomography inversion observation equation: ,in, is the slant path water vapor content vector of the GNSS signal for all GNSS signals, is a coefficient matrix, representing the length of the signal path in the three-dimensional grid, is the three-dimensional water vapor density parameter vector to be determined; S3, establishing constraints, including horizontal constraints on water vapor density, vertical constraints on water vapor density, and prior constraints based on GNSS observation data and meteorological observation data; S4, build a complete 3D tomographic inversion model: ,in, is the prior constraint coefficient matrix, is the horizontal constraint coefficient matrix, is the vertical constraint coefficient matrix, is the prior constraint vector; S5, solving the three-dimensional tomographic inversion model in step S4 to obtain a three-dimensional water vapor density distribution, and then calculating the atmospheric precipitable water volume PWV at each altitude layer; S6, calculate three-dimensional wind speed and direction: For the k The wind speed and direction estimation model is established based on the temporal and spatial variation differences of atmospheric precipitable water volume PWV between adjacent GNSS observation stations. i and j The relationship between the PWV changes is expressed as: ,in, is the adjacent GNSS observation station at the height of the kth layer i and j The time difference between the PWV changes, is the distance between adjacent GNSS observation stations i and j, is the azimuth of the line connecting adjacent GNSS observation stations, is the wind direction at the height of the kth layer, is the wind speed at the height of the kth floor; The observation equations are established by combining multiple GNSS observation stations, and the least squares method is used to solve the wind speed at each altitude layer. and wind direction .

[0008] Preferably, the meteorological observation data includes ground temperature T, ground air pressure P and relative humidity RH, and the GNSS observation data includes pseudorange observation values, carrier phase observation values ​​and satellite orbit data; In the step S2: The expression for the slant path water vapor content of the GNSS signal is: , where SWV represents the water vapor content of the oblique path corresponding to a valid signal, represents the conversion factor, represents the wet mapping function, represents the tropospheric zenith wet delay, represents the atmospheric horizontal gradient mapping function, 、 Represent the north-south and east-west gradient components, respectively. 、 Represent the satellite azimuth and satellite altitude angle respectively, Represents the post-test residual correction of the observation value; In the step S3: The expression of the water vapor density level constraint is: , where represents the water vapor density of the i-th grid in a certain altitude layer, i = 1, 2, 3, ..., q, q represents the number of grids on the same altitude plane, Represents the distance weight between the other (q-1) grids and the i-th grid; , is the distance between grid i and grid j, and ii is the other grids in the same layer as j; is the smoothing factor; The expression of the vertical constraint condition of water vapor density is: , where and Horizontal position The water vapor density at the kth layer and the k+1th layer, and is the height of the corresponding k-th layer and k+1-th layer, is the atmospheric elevation; The priori constraints include: based on the ground water vapor density , using distance weighted interpolation or Kriging interpolation method, the ground water vapor density at each station Interpolate to each grid point, and then use the exponential decrease model to extrapolate to each altitude layer as the initial prior value of the water vapor density in each grid at different altitude layers, where the ground water vapor density is , , e is the surface water vapor pressure, ground temperature T, is the relative humidity, is the saturated water vapor pressure.

[0009] Preferably, in step S5, the multiplication algebraic reconstruction algorithm MART is used to iteratively solve the three-dimensional tomography inversion model, and the iterative formula is: Where, Indicates the k During the iteration j The water vapor density parameter value of each grid, The coefficient matrix representing the water vapor tomography model A No. i Row, No. j The element value of the column, Indicates the first i The observed values ​​of the rows, represents the relaxation factor of the iterative algorithm, n is the total number of grids.

[0010] Preferably, in step S5, the calculation formula of the atmospheric precipitation water volume PWV at each altitude layer is: ,in, is the 2D grid point (x,y) in the k The amount of precipitable water in the atmosphere, is the density of liquid water, For the The water vapor density of the layer, For the The height of the layer.

[0011] Preferably, step S6 includes: S61, extracting extreme feature points of the PWV time series at each altitude layer; S62, obtaining the time difference between adjacent GNSS observation stations when the same extreme value feature point appears , establish a space-time propagation model, for the k Layer, adjacent GNSS observation stations i and j The relationship between the PWV changes is: ; S63, establishes an observation equation group by combining multiple GNSS observation stations, and uses the least squares method to solve the wind speed at each altitude layer and wind direction .

[0012] Preferably, in step S61, the calculation formula for the time rate of change of PWV is: , where represents the time rate of change of PWV, represents the PWV value at a certain time t, Indicates a time interval; The method for determining extreme feature points is: ,in, is the threshold parameter; if it is greater than the threshold parameter, it is determined to be an extreme feature point.

[0013] Preferably, the step S63 includes: using the least square method to simultaneously solve the wind field parameters of multiple height layers: .

[0014] Preferably, in step S63, , , for the formula The observation equation obtained after linearization is: ; Solving for wind field components by iterative least squares method and ; Calculate wind speed ,wind direction , where wind direction The unit is degrees, with north wind being 0° and clockwise being positive.

[0015] Preferably, the distance between adjacent GNSS observation stations i and j is , calculated using the following formula: , where represents the radius of the Earth, 、 Respectively represent measuring stations i and j Latitude, 、 Respectively represent measuring stations i and j longitude; Azimuth of the line connecting adjacent GNSS observation stations , calculated using the following formula: .

[0016] Preferably, the wind speed at each altitude layer calculated in step S6 is also included. and wind direction The steps for verifying the results include: Determine whether the wind speed range meets the requirements m / s, whether the wind direction range complies with .

[0017] The present invention also provides a three-dimensional wind speed and direction detection system based on GNSS tomographic inversion, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the system implements the steps of any one of the above-mentioned three-dimensional wind speed and direction detection methods based on GNSS tomographic inversion.

[0018] This application proposes a three-dimensional wind speed and direction detection method and system based on GNSS tomographic inversion. By acquiring observation data and meteorological observation data from multiple GNSS observation stations, a three-dimensional water vapor density tomographic inversion model is established to obtain the water vapor density distribution at different altitudes, and further obtain the atmospheric precipitable water volume (PWV) at each altitude. By utilizing the coupling relationship between water vapor transport and wind field, combined with the spatial position of the GNSS observation station and the spatiotemporal variation characteristics of PWV at different altitudes, a three-dimensional wind field inversion algorithm is established to achieve quantitative estimation of wind speed and direction at different altitudes. This method can achieve quantitative detection of three-dimensional wind speed and direction based on GNSS tomographic inversion, providing important observation data support for meteorological monitoring and numerical weather forecasting. The present invention has at least the following advantages: (1) This method is not affected by weather conditions such as clouds and rainfall, and can achieve 24-hour continuous monitoring; (2) Utilize existing GNSS infrastructure to avoid duplication of specialized equipment and significantly reduce observation costs; (3) Obtain three-dimensional wind speed and direction information at different altitudes to compensate for the insufficient resolution of traditional observation methods; (4) Based on mature GNSS technology and tomographic inversion theory, it has the conditions for practical engineering application. BRIEF DESCRIPTION OF THE DRAWINGS

[0019] Figure 1 A schematic flow chart of a three-dimensional wind speed and direction detection method based on GNSS tomographic inversion is provided in this embodiment.

[0020] Figure 2 The spatial distribution diagrams of PWV at different altitudes at two adjacent moments are provided for this example.

[0021] Figure 3 A schematic diagram of the hardware structure of a three-dimensional wind speed and direction detection system based on GNSS tomographic inversion is provided in this embodiment. DETAILED DESCRIPTION

[0022] In order to better explain the present invention and facilitate understanding, the present invention is described in detail below through specific implementation methods in conjunction with the accompanying drawings.

[0023] like Figure 1 As shown, this embodiment provides a three-dimensional wind speed and direction detection method based on GNSS tomographic inversion, which includes steps: S1 to S6.

[0024] S1, obtain GNSS observation data of multiple GNSS observation stations in the test area and meteorological observation data of the corresponding GNSS observation stations. Specifically, the meteorological observation data includes ground temperature T, ground pressure P and relative humidity RH, and the GNSS observation data includes pseudorange observation values, carrier phase observation values ​​and satellite orbit data: Preferably, in this embodiment, observation data from multiple GNSS continuously operating reference stations within the measured area are obtained to provide basic data support for the subsequent water vapor density inversion. This includes pseudorange observations, carrier phase observations, satellite orbit data, etc. The tropospheric delay and corresponding satellite position in the signal path of each satellite above the station are calculated using a precise single point positioning algorithm to provide basic data support for the subsequent water vapor density inversion. At the same time, meteorological observation data from the corresponding GNSS observation station are obtained, including ground temperature T, ground air pressure P, and relative humidity RH; these meteorological data will be used to calculate the ground water vapor density and provide boundary conditions and prior constraint information for three-dimensional tomography inversion.

[0025] Specifically, the basic GNSS observation data and meteorological observation data are used to calculate known parameters, which will not be described in detail in this embodiment.

[0026] S2, divide the area to be measured into three-dimensional grids and establish the three-dimensional water vapor density tomography inversion observation equation: ,in, is the slant path water vapor content vector of all GNSS signals, is a coefficient matrix, representing the length of the signal path in the three-dimensional grid, is the three-dimensional water vapor density parameter vector to be determined.

[0027] In this embodiment, the measured area is divided into a three-dimensional grid, ensuring that each grid cell has multiple signal rays passing through it. Each GNSS signal penetrating the model from the top of the grid layer is segmented into small segments distributed within each unit cell within the grid model. Specifically, because water vapor changes rapidly in the lower atmosphere and slowly in the upper atmosphere, a vertical stratification scheme is designed accordingly: a resolution of 200-400m is set within the 0-2km range, 400-600m is used within the 2-5km range, and 600-1000m is used above 5km. Horizontal grid division uses a grid size of 0.1°×0.1° resolution.

[0028] Each GNSS signal penetrating the model from the top of the grid model layer will be divided into small segments distributed within each unit grid within the grid model. Based on this, the basic observation equation of the GNSS signal oblique path precipitable water and three-dimensional water vapor density distribution is established: ,in, is the oblique path water vapor content vector, is a coefficient matrix, which represents the length of the signal path in each grid. is the three-dimensional water vapor density parameter vector to be determined.

[0029] The expression for the slant path water vapor content of the GNSS signal is: , where SWV represents the water vapor content of the oblique path corresponding to a valid signal, represents the conversion factor, represents the wet mapping function, represents the tropospheric zenith wet delay, represents the atmospheric horizontal gradient mapping function, 、 Represent the north-south and east-west gradient components, respectively. 、 Represent the satellite azimuth and satellite altitude angle respectively, Represents the post-test residual correction of the observation value; S3, establish constraints, including: horizontal constraints on water vapor density, vertical constraints on water vapor density, and prior constraints based on GNSS observation data and meteorological observation data.

[0030] Preferably, the horizontal constraint condition means that the atmospheric water vapor density in the same horizontal direction is related, that is, the grid i The water vapor density of the grid can be calculated from the water vapor density of other grids on the same horizontal plane.

[0031] The expression of the water vapor density level constraint is: , where represents the water vapor density of the i-th grid in a certain altitude layer, j = 1, 2, 3, ..., q, q represents the number of grids on the same altitude plane, Represents the distance weight between the other (q-1) grids j and the i-th grid; , is the distance between grid i and grid j. If i and j are not in the same layer, .

[0032] Preferably, water vapor density has the physical property of decreasing exponentially with altitude, so the water vapor densities of adjacent altitude layers along the same vertical direction satisfy an exponential relationship. The expression of the vertical constraint condition of water vapor density is: , where and Horizontal position The water vapor density at the kth layer and the k+1th layer, and is the height of the corresponding k-th layer and k+1-th layer, is the atmospheric scale height, usually 1-2 km.

[0033] Preferably, the initial water vapor density value of each grid point is obtained as a priori constraint by spatial interpolation method based on the ground meteorological data provided by the GNSS observation station. Specifically, the priori constraint conditions include: , using distance weighted interpolation or Kriging interpolation method, the ground water vapor density at each station Interpolate to each grid point, and then use the exponential decrease model to extrapolate to each altitude layer as the initial prior value of the water vapor density in each grid at different altitude layers, where the ground water vapor density is , , e is the surface water vapor pressure, ground temperature T, is the relative humidity, is the saturated water vapor pressure.

[0034] S4, build a complete 3D tomographic inversion model: ,in, is the prior constraint coefficient matrix, is the horizontal constraint coefficient matrix, is the vertical constraint coefficient matrix, is the prior constraint vector.

[0035] S5, solving the three-dimensional tomographic inversion model in step S4 to obtain a three-dimensional water vapor density distribution, and then calculating the atmospheric precipitable water volume PWV at each altitude layer.

[0036] Preferably, in step S5, the multiplication algebraic reconstruction algorithm MART is used to iteratively solve the three-dimensional tomography inversion model, and the iterative formula is: Where, Indicates the k During the iteration j The water vapor density parameter value of each grid, The coefficient matrix representing the water vapor tomography model A No. i Row, No. j The element value of the column, Indicates the first i The observed values ​​of the rows, represents the relaxation factor of the iterative algorithm, n is the total number of grids.

[0037] Preferably, in step S5, based on the three-dimensional water vapor density obtained above , the atmospheric precipitable water volume (PWV) at different altitudes of each horizontal grid can be calculated by vertical integration. The calculation formula for the atmospheric precipitable water volume (PWV) at each altitude is: ,in, is the 2D grid point (x,y) in the k The amount of precipitable water in the atmosphere, is the density of liquid water, For the The water vapor density of the layer, For the The height of the layer.

[0038] S6, calculate three-dimensional wind speed and direction: For the k The wind speed and direction estimation model is established based on the temporal and spatial variation differences of atmospheric precipitable water volume PWV between adjacent GNSS observation stations. i and j The relationship between the PWV changes is expressed as: ,in, is the adjacent GNSS observation station at the height of the kth layer i and j The time difference between the PWV changes, is the distance between adjacent GNSS observation stations i and j, is the azimuth of the line connecting adjacent GNSS observation stations, is the wind direction at the height of the kth layer, is the wind speed at the height of the kth floor; The observation equations are established by combining multiple GNSS observation stations, and the least squares method is used to solve the wind speed at each altitude layer. and wind direction .

[0039] Preferably, step S6 includes: S61, extracting extreme feature points of the PWV time series at each altitude layer; S62, obtaining the time difference between adjacent GNSS observation stations when the same extreme value feature point appears , establish a space-time propagation model, for the k Layer, adjacent GNSS observation stations i and j The relationship between the PWV changes is: ; S63, establishes an observation equation group by combining multiple GNSS observation stations, and uses the least squares method to solve the wind speed at each altitude layer and wind direction .

[0040] Preferably, in step S61, the calculation formula for the time rate of change of PWV is: , where represents the time rate of change of PWV, represents the PWV value at a certain time t, Indicates a time interval; The method for determining extreme feature points is: ,in, is the threshold parameter. If it is greater than the threshold parameter, it can be determined as an extreme feature point.

[0041] Preferably, the step S63 includes: using the least square method to simultaneously solve the wind field parameters of multiple height layers: .

[0042] Preferably, in step S63, in order to linearize the solution formula ,set up , , for the formula The observation equation obtained after linearization is: ; Solved by iterative least squares method and ; Calculate wind speed ,wind direction , where wind direction The unit is degrees, with north wind being 0° and clockwise being positive.

[0043] Preferably, the distance between adjacent GNSS observation stations i and j is , calculated using the following formula: , where represents the radius of the Earth, 、 Respectively represent measuring stations i and j Latitude, 、 Respectively represent measuring stations i and j longitude; Azimuth of the line connecting adjacent GNSS observation stations , calculated using the following formula: .

[0044] Preferably, the wind speed at each altitude layer calculated in step S6 is also included. and wind direction The steps for verifying the results include: Determine whether the wind speed range meets the requirements m / s, whether the wind direction range complies with .

[0045] Preferably, the verification method can also be calculated by using a root mean square error test method, specifically including: using ERA5 to reanalyze data at different altitudes as a reference, and calculating the root mean square error of wind speed and wind direction: , ,in, and are the estimated values ​​of wind speed and direction, and are the reference values ​​of wind speed and direction obtained from ERA5.

[0046] Based on the observation data of CORS stations in a certain area, the three-dimensional wind speed and direction detection method based on GNSS tomographic inversion proposed in this paper can achieve quantitative estimation of wind speed and direction at different altitudes. The spatial distribution of CORS stations in the area, the tomographic area is 113.86°E – 114.36°E, 22.18° - 22.58°, and the horizontal spatial resolution is 0.1° The observation stations are evenly distributed, covering the entire study area and providing continuous and stable atmospheric water vapor data. Figure 2 The spatial distribution of atmospheric precipitable water content (PWV) at different altitudes at two adjacent time points is calculated using the method of the present invention. By calculating the PWV change rate at two adjacent time points, three-dimensional wind speed and direction detection can be achieved. The following are the statistical results of the average wind speed and average wind direction at different altitudes:

[0047] Comparing the average wind speed and direction data at different altitudes with the ERA5 reanalysis data, the RMSE for average wind speed across all altitudes was 3.51 km / h, and the RMSE for average wind direction was 6.96°. This demonstrates that the proposed method can quantitatively estimate wind speed and direction at different altitudes.

[0048] like Figure 3 As shown, the present application also provides a three-dimensional wind speed and direction detection system based on GNSS tomographic inversion, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the steps of a three-dimensional wind speed and direction detection method based on GNSS tomographic inversion as described in any one of the above items are implemented.

[0049] Figure 3 FIG. 1 is a schematic diagram of the hardware structure for running a three-dimensional wind speed and direction detection method based on GNSS tomographic inversion, provided by an embodiment of the present invention. Figure 3As shown, this embodiment / computer 6 includes: a processor 60, a memory 61, and a computer program 62 stored in the memory 61 and executable on the processor 60, such as a program for executing a three-dimensional wind speed and direction detection method based on GNSS tomographic inversion. When the processor 60 executes the computer program 62, it implements the steps described above in the embodiments for executing a three-dimensional wind speed and direction detection method based on GNSS tomographic inversion. Alternatively, when the processor 60 executes the computer program 62, it implements the functions of the modules / units described above in the various device embodiments.

[0050] For example, the computer program 62 may be divided into one or more modules / units, which are stored in the memory 61 and executed by the processor 60 to implement the present invention. The one or more modules / units may be a series of computer program instruction segments capable of implementing specific functions, and the instruction segments are used to describe the execution process of the computer program 62 in the computer 6.

[0051] The computer 6 can be a computing device such as a desktop computer, a notebook, a PDA, or a cloud server. The computer 6 device can include, but is not limited to, a processor 60 and a memory 61. Those skilled in the art will understand that Figure 3 This is only an example of computer 6 and does not constitute a limitation on computer 6. It may include more or fewer components than shown in the figure, or a combination of certain components, or different components. For example, the computer 6 may also include input and output devices, network access devices, buses, etc.

[0052] The processor 60 may be a central processing unit (CPU), other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA), other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The general-purpose processor may be a microprocessor or any conventional processor.

[0053] The memory 61 can be an internal storage unit of the computer 6, such as the computer's hard drive or memory. Alternatively, the memory 61 can be an external storage device of the computer 6, such as a plug-in hard drive, a Smart Media Card (SMC), a Secure Digital (SD) card, a flash memory card, etc., equipped on the terminal device. Furthermore, the memory 61 can include both the internal storage unit of the computer 6 and an external storage device. The memory 61 is used to store the computer program and other programs and data required by the terminal device. The memory 61 can also be used to temporarily store data that has been output or is about to be output.

[0054] Those skilled in the art can clearly understand that, for the convenience and brevity of description, only the division of the above-mentioned functional units and modules is used as an example for illustration. In actual applications, the above-mentioned functions can be distributed and completed by different functional units and modules as needed, that is, the internal structure of the device can be divided into different functional units or modules to complete all or part of the functions described above. The functional units and modules in the embodiment can be integrated into one processing unit, or each unit can exist physically alone, or two or more units can be integrated into one unit. The above-mentioned integrated unit can be implemented in the form of hardware or in the form of software functional units. In addition, the specific names of the functional units and modules are only for the convenience of distinguishing each other, and are not used to limit the scope of protection of this application. The specific working process of the units and modules in the above-mentioned system can refer to the corresponding process in the aforementioned method embodiment, and will not be repeated here.

[0055] In the above embodiments, the description of each embodiment has its own focus. For parts that are not described or recorded in detail in a certain embodiment, reference can be made to the relevant description of other embodiments.

[0056] Those skilled in the art will appreciate that the units and algorithm steps of each example described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are performed in hardware or software depends on the specific application and design constraints of the technical solution. Professionals and technicians can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of the present invention.

[0057] In the embodiments provided by the present invention, it should be understood that the disclosed devices / terminal equipment and methods can be implemented in other ways. For example, the device / terminal equipment embodiments described above are merely illustrative. For example, the division of the modules or units is merely a logical function division. In actual implementation, there may be other division methods, such as multiple units or components can be combined or integrated into another system, or some features can be ignored or not executed. Another point is that the mutual coupling or direct coupling or communication connection shown or discussed can be through some interfaces, indirect coupling or communication connection of devices or units, which can be electrical, mechanical or other forms.

[0058] The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units, that is, they may be located in one place or distributed across multiple network units. Some or all of these units may be selected to achieve the purpose of this embodiment according to actual needs.

[0059] In addition, the functional units in the various embodiments of the present invention may be integrated into a single processing unit, each unit may exist physically separately, or two or more units may be integrated into a single unit. The aforementioned integrated units may be implemented in the form of hardware or software functional units.

[0060] If the integrated module / unit is implemented as a software functional unit and sold or used as a standalone product, it can be stored in a computer-readable storage medium. Based on this understanding, the present invention can implement all or part of the process steps in the above-mentioned method embodiments by using a computer program to instruct the relevant hardware. The computer program can be stored in a computer-readable storage medium. When executed by a processor, the computer program can implement the steps of each of the above-mentioned method embodiments. The computer program includes computer program code, which can be in source code form, object code form, executable file, or some intermediate form. The computer-readable medium can include: any entity or device capable of carrying the computer program code, recording medium, USB flash drive, mobile hard drive, magnetic disk, optical disk, computer memory, read-only memory (ROM), random access memory (RAM), electric carrier signal, telecommunication signal, and software distribution medium. It should be noted that the content of the computer-readable medium can be appropriately increased or decreased based on the requirements of legislation and patent practice in a jurisdiction. For example, in some jurisdictions, based on legislation and patent practice, computer-readable media does not include electric carrier signals and telecommunication signals.

[0061] The above are merely specific application examples of the present invention and do not constitute any limitation on the scope of protection of the present invention. In addition to the above embodiments, the present invention may also have other implementations. Any technical solution formed by equivalent replacement or equivalent transformation falls within the scope of protection claimed by the present invention.

Claims

1. A three-dimensional wind speed and direction detection method based on GNSS tomographic inversion, characterized in that: The following steps are involved: S1, obtaining GNSS observation data of multiple GNSS observation stations in the area to be measured and meteorological observation data of the corresponding GNSS observation stations; S2, divide the area to be measured into three-dimensional grids and establish the three-dimensional water vapor density tomography inversion observation equation: ,in, is the slant path water vapor content vector of all GNSS signals, is a coefficient matrix, representing the length of the signal path in the three-dimensional grid, is the three-dimensional water vapor density parameter vector to be determined; S3, establishing constraints, including horizontal constraints on water vapor density, vertical constraints on water vapor density, and prior constraints based on GNSS observation data and meteorological observation data; S4, build a complete 3D tomographic inversion model: ,in, is the prior constraint coefficient matrix, is the horizontal constraint coefficient matrix, is the vertical constraint coefficient matrix, is the prior constraint vector; S5, solving the three-dimensional tomographic inversion model in step S4 to obtain a three-dimensional water vapor density distribution, and then calculating the atmospheric precipitable water volume PWV at each altitude layer; S6, calculate three-dimensional wind speed and direction: For the k The wind speed and direction estimation model is established based on the temporal and spatial variation differences of atmospheric precipitable water volume PWV between adjacent GNSS observation stations. i and j The relationship between the PWV changes is expressed as: ,in, is the adjacent GNSS observation station at the height of the kth layer i and j The time difference between the PWV changes, is the distance between adjacent GNSS observation stations i and j, is the azimuth of the line connecting adjacent GNSS observation stations, is the wind direction at the height of the kth layer, is the wind speed at the height of the kth floor; The observation equations are established by combining multiple GNSS observation stations, and the least squares method is used to solve the wind speed at each altitude layer. and wind direction .

2. The three-dimensional wind speed and direction detection method based on GNSS tomographic inversion according to claim 1 is characterized in that: The meteorological observation data includes ground temperature T, ground pressure P and relative humidity RH, and the GNSS observation data includes pseudorange observation values, carrier phase observation values ​​and satellite orbit data; In the step S2: The expression for the slant path water vapor content of the GNSS signal is: , where SWV represents the water vapor content of the oblique path corresponding to a valid signal, represents the conversion factor, represents the wet mapping function, represents the tropospheric zenith wet delay, represents the atmospheric horizontal gradient mapping function, 、 Represent the north-south and east-west gradient components, respectively. 、 Represent the satellite azimuth and satellite altitude angle respectively, Represents the post-test residual correction of the observation value; In the step S3: The expression of the water vapor density level constraint is: , where represents the water vapor density of the i-th grid in a certain altitude layer, j = 1, 2, 3, ..., q, q represents the number of grids on the same altitude plane, Represents the distance weight between the other (q-1) grids j and the i-th grid; , is the distance between grid i and grid j; The expression of the vertical constraint condition of water vapor density is: , where and Horizontal position The water vapor density at the kth layer and the k+1th layer, and is the height of the corresponding k-th layer and k+1-th layer, is the atmospheric elevation; The priori constraints include: based on the ground water vapor density , using distance weighted interpolation or Kriging interpolation method, the ground water vapor density at each station Interpolate to each grid point, and then use the exponential decrease model to extrapolate to each altitude layer as the initial prior value of the water vapor density in each grid at different altitude layers, where the ground water vapor density is , , e is the surface water vapor pressure, ground temperature T, is the relative humidity, is the saturated water vapor pressure.

3. The three-dimensional wind speed and direction detection method based on GNSS tomographic inversion according to claim 1 is characterized in that: In step S5, the multiplication algebraic reconstruction algorithm MART is used to iteratively solve the three-dimensional tomographic inversion model, and the iterative formula is: Where, Indicates the k During the iteration j The water vapor density parameter value of each grid, The coefficient matrix representing the water vapor tomography model A No. i Row, No. j The element value of the column, Indicates the first i The observed values ​​of the rows, represents the relaxation factor of the iterative algorithm, n is the total number of grids.

4. The three-dimensional wind speed and direction detection method based on GNSS tomographic inversion according to claim 3 is characterized in that: In step S5, the calculation formula of the atmospheric precipitable water volume PWV at each altitude layer is: ,in, is the 2D grid point (x,y) in the k The amount of precipitable water in the atmosphere, is the density of liquid water, For the The water vapor density of the layer, For the The height of the layer.

5. The three-dimensional wind speed and direction detection method based on GNSS tomographic inversion according to claim 1 is characterized in that: The step S6 comprises: S61, extracting extreme feature points of the PWV time series at each altitude layer; S62, obtaining the time difference between adjacent GNSS observation stations when the same extreme value feature point appears , establish a space-time propagation model, for the k Layer, adjacent GNSS observation stations i and j The relationship between the PWV changes is: ; S63, establishes an observation equation group by combining multiple GNSS observation stations, and uses the least squares method to solve the wind speed at each altitude layer and wind direction .

6. The three-dimensional wind speed and direction detection method based on GNSS tomographic inversion according to claim 5 is characterized in that: In step S61, the calculation formula for the time rate of change of PWV is: , where represents the time rate of change of PWV, represents the PWV value at a certain time t, Indicates a time interval; The method for determining extreme feature points is: ,in, is the threshold parameter; if it is greater than the threshold parameter, it is determined to be an extreme feature point.

7. The three-dimensional wind speed and direction detection method based on GNSS tomographic inversion according to claim 6 is characterized in that: The step S63 includes: The least squares method is used to solve the wind field parameters of multiple height layers at the same time: ; set up , , for the formula The observation equation obtained after linearization is: ; Solving for wind field components by iterative least squares method and ; Calculate wind speed ,wind direction , where wind direction The unit is degrees, with north wind being 0° and clockwise being positive.

8. The three-dimensional wind speed and direction detection method based on GNSS tomographic inversion according to claim 5 is characterized in that: The distance between adjacent GNSS observation stations i and j , calculated using the following formula: , where represents the radius of the Earth, 、 Respectively represent measuring stations i and j Latitude, 、 Respectively represent measuring stations i and j longitude; Azimuth of the line connecting adjacent GNSS observation stations , calculated using the following formula: .

9. A three-dimensional wind speed and direction detection method based on GNSS tomographic inversion according to any one of claims 1 to 8, characterized in that: Also includes the wind speed of each altitude layer calculated in step S6 and wind direction The steps for verifying the results include: Determine whether the wind speed range meets the requirements m / s, whether the wind direction range complies with .

10. A three-dimensional wind speed and direction detection system based on GNSS tomographic inversion, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the computer program, the steps of the three-dimensional wind speed and direction detection method based on GNSS tomographic inversion according to any one of claims 1 to 9 are implemented.

Images