Three-dimensional wind field retrieval method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing

Abstract

The application discloses a three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing, and the method comprises the following steps: obtaining radar site geographic information; using a triangulation method to perform radar networking according to the radar site geographic information; performing inversion grid division after the radar networking; obtaining radial velocity data through radar detection and performing interpolation; obtaining radial velocity two-dimensional filtering results; and performing three-dimensional wind field inversion through the inversion grid, the radar site geographic information and the radial velocity two-dimensional filtering results. Compared with the existing three-dimensional wind field inversion algorithm taking the variational method as the core, the application focuses on the optimization of the radial velocity of each radar before the wind field inversion, thereby reducing the calculation burden in the inversion process, optimizing the radar networking method and the interpolation algorithm in the overall three-dimensional wind field inversion, improving the inversion accuracy while reducing the calculation complexity, and being more suitable for the promotion in actual business.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of meteorological detection and analysis, and particularly relates to a three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing. BACKGROUND

[0002] Whether it is a phased array radar or a Doppler radar, the atmospheric wind field information obtained by them is relatively limited, usually being reflectivity, radial velocity and velocity spectrum width, and they cannot directly complete the detection of complete three-dimensional wind field. Therefore, using multiple radars to network for overall three-dimensional wind field inversion has always been a topic of interest to meteorologists.

[0003] When a three-dimensional wind field is constructed according to the positional relationship of the radial velocity measured by the radar in the radial direction and the three-dimensional components of the wind field vector in the atmosphere, there is a large vertical error in the low layer result in the wind field result. In order to improve this problem, most of the current three-dimensional wind field inversion methods use three-dimensional variational method as the core algorithm. The use of the cost function of this algorithm lacks standard, the optimal solution is difficult to judge, the vertical velocity error problem is limited to improve, and the calculation complexity is high, which leads to the fact that the three-dimensional variational method has not been widely promoted in actual meteorological business. SUMMARY

[0004] The purpose of the present application is to provide a three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing, so as to solve the problem that the three-dimensional variational method as the core algorithm in the existing three-dimensional wind field inversion method has not been widely promoted in actual meteorological business.

[0005] The embodiments of the present application are implemented in the following manner. A three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing is provided, which comprises the following steps.

[0006] Using a triangulation method to perform radar networking according to the radar site geographical information;

[0007] Performing inversion grid division after radar networking;

[0008] Obtaining radial velocity data through radar detection and performing interpolation;

[0009] Performing two-dimensional filtering on the interpolation result of the radial velocity data obtained through radar detection by using the height of the divided inversion grid and the result obtained by averaging the signal-to-noise ratio obtained through radar detection, to obtain a radial velocity two-dimensional filtering result;

[0010] Performing three-dimensional wind field inversion by using the inversion grid, the radar site geographical information and the radial velocity two-dimensional filtering result.

[0011] In some embodiments,

[0012] In some embodiments, the radar site geographic information comprises longitude, latitude and altitude; and / or

[0013] The radar networking obtains a radar network, in which the geographic information data of each radar site is converted from a polar coordinate system to a Cartesian coordinate system, and the conversion formula is as follows:

[0014] x = lcosεsinθ + x0

[0015] y = lcosεcosθ + y0

[0016] z = lsine + z0

[0017] In the formula, x is the x-axis coordinate in the converted Cartesian coordinate system, y is the y-axis coordinate in the converted Cartesian coordinate system, z is the z-axis coordinate in the converted Cartesian coordinate system, x0 is the x-axis coordinate of the radar site position in the Cartesian coordinate system, y0 is the y-axis coordinate of the radar site position in the Cartesian coordinate system, z0 is the altitude of the radar site position in the Cartesian coordinate system, l is the radial distance of the radar, θ is the elevation angle of the radar, and ε is the azimuth angle of the radar.

[0018] In some embodiments, after converting all the geographic position data of the radar sites from the polar coordinate system to the Cartesian coordinate system, an inversion grid is divided from the radar network; and / or

[0019] Using a triangulation method, the radial velocity data obtained by the radar is triangulated according to the relative Cartesian coordinates to obtain a plurality of tetrahedrons whose vertices are composed of point coordinates of the radial velocity data obtained by the radar and are closely connected, and the radial velocity data obtained by the radar used is determined by judging the distribution of the inversion grid center P in the tetrahedron; and / or

[0020] According to the average results of the signal-to-noise ratio of each radar in the triangular sub-region in the radar network and the radial velocity data in the different inversion grid height layers in the three-dimensional space which have been interpolated, corresponding frequency domain two-dimensional filters are respectively set, and after all the filters are set, a filter bank is obtained.

[0021] In some embodiments, the method for judging the distribution of the inversion grid center P in the tetrahedron comprises:

[0022] The volume of the tetrahedron ABCD is calculated, and the formula is as follows:

[0023]

[0024] wherein x1 is the x-axis coordinate of the tetrahedron vertex A in space, y1 is the y-axis coordinate of the tetrahedron vertex A in space, z1 is the z-axis coordinate of the tetrahedron vertex A in space, x2 is the x-axis coordinate of the tetrahedron vertex B in space, y2 is the y-axis coordinate of the tetrahedron vertex B in space, z2 is the z-axis coordinate of the tetrahedron vertex B in space, x3 is the x-axis coordinate of the tetrahedron vertex C in space, y3 is the y-axis coordinate of the tetrahedron vertex C in space, z3 is the z-axis coordinate of the tetrahedron vertex C in space, x4 is the x-axis coordinate of the tetrahedron vertex D in space, y4 is the y-axis coordinate of the tetrahedron vertex D in space, and z4 is the z-axis coordinate of the tetrahedron vertex D in space; V ΔABCD is the directed volume of the tetrahedron ABCD.

[0025] The volumes of the tetrahedron formed by the inversion grid center P and each vertex of the tetrahedron ABCD are obtained respectively, and the formula is as follows:

[0026]

[0027]

[0028] wherein x is the x-axis coordinate of the inversion grid center P in space, y is the y-axis coordinate of the inversion grid center P in space, z is the z-axis coordinate of the inversion grid center P in space, x1 is the x-axis coordinate of the tetrahedron vertex A in space, y1 is the y-axis coordinate of the tetrahedron vertex A in space, z1 is the z-axis coordinate of the tetrahedron vertex A in space, x2 is the x-axis coordinate of the tetrahedron vertex B in space, y2 is the y-axis coordinate of the tetrahedron vertex B in space, z2 is the z-axis coordinate of the tetrahedron vertex B in space, x3 is the x-axis coordinate of the tetrahedron vertex C in space, y3 is the y-axis coordinate of the tetrahedron vertex C in space, z3 is the z-axis coordinate of the tetrahedron vertex C in space, x4 is the x-axis coordinate of the tetrahedron vertex D in space, y4 is the y-axis coordinate of the tetrahedron vertex D in space, and z4 is the z-axis coordinate of the tetrahedron vertex D in space; V ΔPBCD is the directed volume of the tetrahedron PBCD, V ΔAPCD is the directed volume of the tetrahedron APCD, V ΔABPD is the directed volume of the tetrahedron ABCD, V ΔABCP is the directed volume of the tetrahedron ABCP.

[0029] It is judged whether the inversion grid center P is in the tetrahedron ABCD, and when the condition V ΔABCD = V ΔPBCD + V ΔAPCD + V ΔABPD + V ΔABCPIf the condition is satisfied, it means that the tetrahedron ABCD contains the inversion grid center P; if the condition is not satisfied, it means that the tetrahedron ABCD does not contain the inversion grid center P, and another tetrahedron meeting the condition needs to be found for the inversion grid center P;

[0030] According to the volume relationship between the inversion grid center P and the sub-tetrahedron and the tetrahedron ABCD, the weight used in the interpolation process is obtained, and then the radial velocity V of the grid is obtained by interpolation Pr

[0031]

[0032] In the formula, V Pr is the radial velocity after interpolation, V ΔABCD is the directed volume of the tetrahedron ABCD, V ΔPBCD is the directed volume of the tetrahedron PBCD, V ΔAPCD is the directed volume of the tetrahedron APCD, V ΔABPD is the directed volume of the tetrahedron ABPD, V ΔABCP is the directed volume of the tetrahedron ABCP; V Ar , V Br , V Cr , V Dr are respectively the radial velocities of the radar detection at the vertices of the tetrahedron ABCD; and / or

[0033] The formula of the frequency response of the two-dimensional filter is as follows:

[0034]

[0035] f(h)=k*h+SNR*s

[0036] In the formula, |H h (ω1,ω2)| is the frequency response of the two-dimensional filter, D(ω1,ω2) is the distance from each point of the radial velocity of each radar in the frequency domain to the 0 frequency point, and is expressed as ω1 is the horizontal coordinate of the two-dimensional frequency axis, and ω2 is the vertical coordinate of the two-dimensional frequency axis; f(h) is a function related to the height of the inversion grid, h is the height of the wind field inversion; k is a height influence factor, ranging from 0 to 1; SNR is the average result of the signal-to-noise ratio of the triangular sub-region located at the vertex obtained by radar detection in the inversion process; s is a basic passband influence factor, ranging from 0 to 1.

[0037] ​In some embodiments, the interpolated radial velocity is subjected to a two-dimensional fast Fourier transform to convert to a frequency domain, multiplied by a two-dimensional filter frequency response corresponding to the height of the inversion grid in the filter bank to obtain a frequency domain filter result, and finally subjected to a two-dimensional inverse Fourier transform to obtain a two-dimensional filtered radial velocity result.

[0038] In some embodiments, the selection of the three-dimensional wind field inversion data includes: calculating the directed area of a triangular sub-region in the radar network, and then calculating the directed area of a triangle formed by the projection P'(x, y) of the inversion grid center P on the xy plane and the vertices A, B and C of the triangular sub-region, respectively, according to the following formula:

[0039]

[0040] In the formula, x is the x-axis coordinate of the inversion grid center P in space, y is the y-axis coordinate of the inversion grid center P in space, x1 is the x-axis coordinate of the vertex A of the triangle in space, y1 is the y-axis coordinate of the vertex A of the triangle in space, x2 is the x-axis coordinate of the vertex C of the triangle in space, y2 is the y-axis coordinate of the vertex B of the triangle in space, x3 is the x-axis coordinate of the vertex C of the triangle in space, and y3 is the y-axis coordinate of the vertex C of the triangle in space; and S ΔP BC is the directed area of the triangle P'BC, S ΔAP′C is the directed area of the triangle AP'C, and S ΔABP′ is the directed area of the triangle ABP'.

[0041] When the condition S ΔP +S BC +S ΔAP′C +S ΔABP′ =S ΔABC is satisfied, the radial velocity data detected by the radar located at the vertices of the triangular sub-region ABC can be selected as the three-dimensional wind field inversion data.

[0042] In some embodiments, the inversion grid center P is located in the triangular sub-region ABC, the coordinates of the inversion grid center P are (x, y, z), the radar coordinates of the vertices of the triangular sub-region ABC are (x1, y1, z1), (x2, y2, z2) and (x3, y3, z3) respectively, and a direction vector is constructed according to the following formula:

[0043]

[0044] In the formula, d1 is a direction vector composed of the first part of the radar coordinates in the triangular sub-region and the inversion grid center P, ​​the direction vector composed of the third part of radar coordinates in the triangular sub-region and the inversion grid center P, the direction vector composed of the third part of radar coordinates in the triangular sub-region and the inversion grid center P,

[0045] The relationship between the two-dimensional filtering result of radial velocity and the wind speed component is obtained by vector projection, as follows:

[0046]

[0047] In the formula, V r1 the radial velocity of the first part of radar in the triangular sub-region after two-dimensional filtering at the inversion grid center P, V r2 the radial velocity of the second part of radar in the triangular sub-region after two-dimensional filtering at the inversion grid center P, V r3 the radial velocity of the third part of radar in the triangular sub-region after two-dimensional filtering at the inversion grid center P; u is the to-be-solved north-south direction horizontal wind speed component; v is the to-be-solved east-west direction horizontal wind speed component, and w is the to-be-solved vertical wind speed component; the direction vector composed of the first part of radar coordinates in the triangular sub-region and the inversion grid center P, the direction vector composed of the second part of radar coordinates in the triangular sub-region and the inversion grid center P, the direction vector composed of the third part of radar coordinates in the triangular sub-region and the inversion grid center P, the 2-norm of the direction vector the 2-norm of the direction vector the 2-norm of the direction vector the 2-norm of the direction vector the 2-norm of the direction vector the 2-norm of the direction vector

[0048] The matrix of the unknown speed component uvw inversion equation calculation process is as follows:

[0049]

[0050] In the formula, V r1 the radial velocity of the first part of radar in the triangular sub-region after two-dimensional filtering at the inversion grid center P, V r2 the radial velocity of the second part of radar in the triangular sub-region after two-dimensional filtering at the inversion grid center P, V r3Let be the radial velocity of the third radar in the triangular subregion, obtained by two-dimensional filtering at the centroid P of the inversion grid; x is the x-axis coordinate of the centroid P in space; y is the y-axis coordinate of the centroid P in space; z is the z-axis coordinate of the centroid P in space; x1 is the x-axis coordinate of the first radar position in the triangular subregion in Cartesian coordinates; y1 is the y-axis coordinate of the first radar position in the triangular subregion in Cartesian coordinates; and z1 is the altitude of the first radar position in the triangular subregion in Cartesian coordinates. x2 represents the x-coordinate of the second radar location in the triangular sub-region on the Cartesian coordinate system; y2 represents the y-coordinate of the second radar location in the triangular sub-region on the Cartesian coordinate system; z2 represents the altitude of the second radar location in the triangular sub-region on the Cartesian coordinate system; x3 represents the x-coordinate of the third radar location in the triangular sub-region on the Cartesian coordinate system; y3 represents the y-coordinate of the third radar location in the triangular sub-region on the Cartesian coordinate system; z3 represents the altitude of the third radar location in the triangular sub-region on the Cartesian coordinate system. The direction vector formed by the first radar coordinates in the triangular sub-region and the centroid P of the inversion grid. The direction vector formed by the second radar coordinates in the triangular subregion and the centroid P of the inversion grid. The direction vector formed by the third radar coordinates in the triangular sub-region and the centroid P of the inversion grid. Represented as direction vector 2 Fan, Represented as direction vector 2 Fan, Represented as direction vector In the formula, u represents the horizontal wind speed component in the north-south direction, v represents the horizontal wind speed component in the east-west direction, and w represents the vertical wind speed component; A is a simplified explanation of the matrix in the formula. - 1 represents the inversion matrix.

[0051] In some embodiments, the intensity data obtained by radar detection is then interpolated, and the interpolated intensity data is fused, as shown in the following formula:

[0052]

[0053] In the formula, Z represents the fused intensity data, Z1 represents the intensity obtained by interpolating the first radar in the triangular sub-region at the centroid P of the inversion grid, Z2 represents the intensity obtained by interpolating the first radar in the triangular sub-region at the centroid P of the inversion grid, and Z3 represents the intensity obtained by interpolating the first radar in the triangular sub-region at the centroid P of the inversion grid. The direction vector formed by the first radar coordinates in the triangular sub-region and the centroid P of the inversion grid. The direction vector formed by the second radar coordinates in the triangular subregion and the centroid P of the inversion grid. The direction vector formed by the third radar coordinates in the triangular sub-region and the centroid P of the inversion grid. Represented as direction vector 2 Fan, Represented as direction vector 2 Fan, Represented as direction vector 2.

[0054] In some embodiments, the true vertical component of the wind vector is as follows:

[0055] w ′ =ww t

[0056] w t =2.65(ρ0 / ρ)Z 0.114

[0057] In the formula, w′ is the true vertical component of the wind vector; w t ρ represents the falling velocity of precipitation particles; w represents the vertical velocity component obtained from the 3D wind field inversion; Z represents the fused intensity data; ρ0 represents the air density at the inversion height plane; and ρ represents the average air density.

[0058] In some embodiments, an inversion grid on the boundary of an adjacent radar network triangular sub-region is selected. Simultaneously, wind field inversion is performed on this inversion grid using the radars belonging to the two adjacent triangular sub-regions of the radar network, yielding two results for the wind speed components belonging to the boundary grid, denoted as u1, v1, w′1 and u2, v2, w′2. These components are then synthesized by averaging, as shown in the following formula:

[0059]

[0060] In the formula, u1 is the north-south horizontal component of the inversion grid result of the first triangular sub-region in the adjacent radar network triangular sub-region, v1 is the east-west horizontal component of the inversion grid result of the first triangular sub-region in the adjacent radar network triangular sub-region, w′1 is the vertical component of the inversion grid result of the first triangular sub-region in the adjacent radar network triangular sub-region; u2 is the north-south horizontal component of the inversion grid result of the second triangular sub-region in the adjacent radar network triangular sub-region, v2 is the east-west horizontal component of the inversion grid result of the second triangular sub-region in the adjacent radar network triangular sub-region, w′2 is the vertical component of the inversion grid result of the second triangular sub-region in the adjacent radar network triangular sub-region; u is the synthesized result of the north-south horizontal component of the inversion grid on the boundary of the adjacent radar network triangular sub-region, v is the synthesized result of the east-west horizontal component of the inversion grid on the boundary of the adjacent radar network triangular sub-region, w′ is the synthesized result of the vertical component of the inversion grid on the boundary of the adjacent radar network triangular sub-region;

[0061] obtaining the wind speed size v′, the horizontal direction θ and the vertical direction The formula is as follows:

[0062]

[0063] In the formula, u is the synthesized result of the north-south horizontal component of the inversion grid, v is the synthesized result of the east-west horizontal component of the inversion grid, w is the synthesized result of the vertical component of the inversion grid, and w′ is the synthesized result of the vertical component of the inversion grid on the boundary of the adjacent radar network triangular sub-region.

[0064] Correspondingly, the embodiment of the application further provides a three-dimensional wind field inversion system based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing, comprising:

[0065] A radar site geographic information system is configured to obtain radar site geographic information.

[0066] A radar networking system is configured to use a triangular partitioning method to perform radar networking according to the radar site geographic information.

[0067] An inversion grid division module is configured to perform inversion grid division after radar networking.

[0068] An interpolation module is configured to obtain radial velocity data through radar detection and perform interpolation.

[0069] A radial velocity two-dimensional filtering result module is configured to perform two-dimensional filtering on the interpolation result of the radial velocity data obtained through radar detection by using the height of the divided inversion grid and the result obtained by averaging the signal-to-noise ratio obtained through radar detection, to obtain a radial velocity two-dimensional filtering result.

[0070] A three-dimensional wind field inversion module is configured to perform three-dimensional wind field inversion by inverting a grid, geographical information of radar sites and two-dimensional filtering results of radial velocities.

[0071] Correspondingly, the embodiments of the present application further provide a computer device, comprising a storage and a processor, wherein the storage stores a computer program, and the computer program is executed by the processor to make the processor execute the steps of the above method.

[0072] Correspondingly, the embodiments of the present application further provide a computer readable storage medium, which stores a computer program, and the computer program is executed by a processor to make the processor execute the steps of the above method.

[0073] In summary, due to the adoption of the above technical solutions, the present application has the following advantages:

[0074] Compared with the existing three-dimensional wind field inversion algorithm based on variational method, the present application focuses on the optimization of radial velocities of each radar before wind field inversion, thereby reducing the calculation burden in the inversion process, optimizing the radar networking method and interpolation algorithm in the whole three-dimensional wind field inversion, improving the inversion accuracy while reducing the calculation complexity, and being more suitable for practical business promotion. BRIEF DESCRIPTION OF DRAWINGS

[0075] Figure 1 A flowchart of the three-dimensional wind field inversion method based on two-dimensional filtering of radial velocities and smoothing of wind vectors on the edges of a radar network is provided for the embodiments of the present application;

[0076] Figure 2 A geometric relationship diagram between the inversion grid center P and the tetrahedron obtained by partitioning is provided for the embodiments of the present application;

[0077] Figure 3 A data diagram existing on the edges of the divided radar network is provided for the embodiments of the present application;

[0078] Figure 4 A radar network diagram obtained by partitioning is provided for the embodiments of the present application;

[0079] Figure 5 A partial two-dimensional filter display diagram is provided for the embodiments of the present application;

[0080] Figure 6 A grid area diagram of the edges of the divided radar network is provided for the embodiments of the present application;

[0081] Figure 7 A comparison diagram of two-dimensional filtering results before and after filtering is provided for the embodiments of the present application;

[0082] Figure 8 A comparison diagram of accuracy before and after two-dimensional filtering is provided for the embodiments of the present application;

[0083] Figure 9 A comparison chart of actual wind field inversion results before and after filtering is provided for the embodiments of the present application.

[0084] Figure 10 An additional radar verification result chart is provided for the embodiments of the present application. DETAILED DESCRIPTION

[0085] In order to make the purpose, technical solutions and advantages of the present application clearer, the present application is further described in detail below in combination with the drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and do not limit the present application.

[0086] The technical solutions of the present application are as follows:

[0087] In a first aspect, the embodiments of the present application provide a three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing, comprising:

[0088] S01, obtaining radar site geographic information;

[0089] S02, using a triangulation method to perform radar networking according to the radar site geographic information;

[0090] S03, performing inversion grid division after radar networking;

[0091] S04, obtaining radial velocity data by radar detection and performing interpolation;

[0092] S05, performing two-dimensional filtering on the interpolation results of the radial velocity data obtained by radar detection by the height of the divided inversion grid and the results obtained by averaging the signal-to-noise ratio of radar detection;

[0093] S06, performing three-dimensional wind field inversion by the inversion grid, radar site geographic information and radial velocity two-dimensional filtering results.

[0094] Compared with the existing three-dimensional wind field inversion algorithm with variational method as the core, the present application focuses on the optimization of the radial velocity of each radar before wind field inversion, thereby reducing the calculation burden in the inversion process, and optimizing the radar networking method and interpolation algorithm in the whole three-dimensional wind field inversion. In addition to further improving the inversion accuracy, it also reduces the calculation complexity and is more suitable for popularization in practical business. In addition, it also solves the vertical velocity error of the low-altitude wind field in the inversion result, and to some extent, it also optimizes the inversion error of the horizontal velocity component.

[0095] The application reconsiders the whole wind field inversion process, and proposes a three-dimensional wind field inversion construction method based on two-dimensional filtering of radial velocity at different heights according to a large amount of error analysis. The algorithm with two-dimensional filtering of radial velocity at different inversion heights is proposed to solve the problem of obvious error of low-altitude inversion result.

[0096] Finally, the overall wind field inversion process is applied to the three-dimensional wind field inversion experiment of actual detection data. The effectiveness of the overall three-dimensional wind field inversion process optimized by the application is verified through additional S-band Doppler radar detection data. Meanwhile, the necessity and reliability of the two-dimensional filtering of radial velocity are proved through the comparison of the effects of two-dimensional filtering of radial velocity at different heights.

[0097] In the S01, the radar site geographical information includes longitude, latitude and altitude.

[0098] In some embodiments, the radar site geographical information includes longitude, latitude and altitude.

[0099] In the S02, the radar site geographical information includes longitude, latitude and altitude.

[0100] It can be understood that, by regarding the radars as endpoints on a plane, the Bowyer-Watson algorithm can make the sum of the minimum internal angles of the divided triangular sub-regions greater than the sum of the minimum internal angles of the triangles formed by any non-Delaunay triangulation, without the singularity of the divided triangular regions. This property makes the sub-regions surrounded by the divided radar network closest to equilateral triangles. Moreover, it can ensure that the divided regions are unique, so that the divided sub-regions are completely determined after one network formation, without the need to perform triangulation every time the wind field is inverted. Moreover, when a radar site is added, deleted or moved, it will only affect the adjacent sub-regions, without affecting the network formation result of the original overall radar, which is more conducive to the use and update of the three-dimensional wind field inversion algorithm in actual business.

[0101] In the S03, the radar site geographical information includes longitude, latitude and altitude.

[0102] It can be understood that, after the radars are networked, the inversion grid division can be performed. First, the coordinate system needs to be converted according to the geographical positions of the radars. The geographical position of each radar in the radar network can be represented as (x0, y0, z0), which respectively represent the relative coordinates (x0, y0) of the radar about the distance converted from the absolute latitude and longitude coordinates to the relative Cartesian coordinate system and the altitude z0 thereof.

[0103] In some embodiments, the radar networking obtains a radar network, and in the radar network, the geographical information data of each radar site is converted from the polar coordinate system to the Cartesian coordinate system, and the conversion formula is as follows:

[0104] x = lcosεsinθ + x0

[0105] y = lcosεcosθ + y0

[0106] z = lεsin + z0

[0107] In the formula, x is the x-axis coordinate in the converted Cartesian coordinate system, y is the y-axis coordinate in the converted Cartesian coordinate system, z is the z-axis coordinate in the converted Cartesian coordinate system, x0 is the x-axis coordinate of the radar site position in the Cartesian coordinate system, y0 is the y-axis coordinate of the radar site position in the Cartesian coordinate system, z0 is the altitude of the radar site position in the Cartesian coordinate system, l is the radial distance of the radar, θ is the elevation angle of the radar, and ε is the azimuth angle of the radar.

[0108] Further, after converting all radar site geographical position data from the polar coordinate system to the Cartesian coordinate system, the inversion grid is divided from the radar network.

[0109] In S04, the following steps are performed:

[0110] It can be understood that the regular inversion grid and the detection data of each radar after conversion of the coordinate system cannot be aligned, and the distribution of the detection data is also uneven. The distribution of the radial velocity data will be more dense closer to the radar, and the farther away from the radar, the more sparse the data distribution. Therefore, the scattered radial velocity data in the coordinate system needs to be uniformly distributed on the regular inversion grid through interpolation.

[0111] In some embodiments, the triangular partitioning method is adopted, and the radial velocity data obtained by radar detection is triangularly partitioned according to the relative Cartesian coordinates to obtain a plurality of tetrahedrons whose vertexes are composed of point coordinates of the radial velocity data obtained by radar detection and are closely connected. The radial velocity data obtained by radar detection used is determined by judging the distribution of the inversion grid center P in the tetrahedron.

[0112] Further, the method for judging the distribution of the inversion grid center P in the tetrahedron includes:

[0113] S041, the volume of the tetrahedron ABCD is calculated, and the formula is as follows:

[0114]

[0115] wherein x1 is the x-axis coordinate of the tetrahedron vertex A in space, y1 is the y-axis coordinate of the tetrahedron vertex A in space, z1 is the z-axis coordinate of the tetrahedron vertex A in space, x2 is the x-axis coordinate of the tetrahedron vertex B in space, y2 is the y-axis coordinate of the tetrahedron vertex B in space, z2 is the z-axis coordinate of the tetrahedron vertex B in space, x3 is the x-axis coordinate of the tetrahedron vertex C in space, y3 is the y-axis coordinate of the tetrahedron vertex C in space, z3 is the z-axis coordinate of the tetrahedron vertex C in space, x4 is the x-axis coordinate of the tetrahedron vertex D in space, y4 is the y-axis coordinate of the tetrahedron vertex D in space, and z4 is the z-axis coordinate of the tetrahedron vertex D in space; V ΔABCD is the directed volume of the tetrahedron ABCD;

[0116] S042, the volume surrounded by the inversion grid center P and each vertex of the tetrahedron ABCD is obtained respectively, and the formula is as follows:

[0117]

[0118] wherein x is the x-axis coordinate of the inversion grid center P in space, y is the y-axis coordinate of the inversion grid center P in space, z is the z-axis coordinate of the inversion grid center P in space, x1 is the x-axis coordinate of the tetrahedron vertex A in space, y1 is the y-axis coordinate of the tetrahedron vertex A in space, z1 is the z-axis coordinate of the tetrahedron vertex A in space, x2 is the x-axis coordinate of the tetrahedron vertex B in space, y2 is the y-axis coordinate of the tetrahedron vertex B in space, z2 is the z-axis coordinate of the tetrahedron vertex B in space, x3 is the x-axis coordinate of the tetrahedron vertex C in space, y3 is the y-axis coordinate of the tetrahedron vertex C in space, z3 is the z-axis coordinate of the tetrahedron vertex C in space, x4 is the x-axis coordinate of the tetrahedron vertex D in space, y4 is the y-axis coordinate of the tetrahedron vertex D in space, and z4 is the z-axis coordinate of the tetrahedron vertex D in space; V ΔPBCD is the directed volume of the tetrahedron PBCD, V ΔAPCD is the directed volume of the tetrahedron APCD, V ΔABPD is the directed volume of the tetrahedron ABPD, V ΔABCP is the directed volume of the tetrahedron ABCP;

[0119] S043, whether the inversion grid center P is in the tetrahedron ABCD is judged, and when the condition V ΔABCD = V ΔPBCD + V ΔAPCD + V ΔABPD + V ΔABCPIf the condition is met, it means that the tetrahedron ABCD contains the inversion grid center P; if the condition is not met, it means that the tetrahedron ABCD does not contain the inversion grid center P, and another tetrahedron meeting the condition needs to be found for the inversion grid center P;

[0120] S044、According to the volume relationship between the sub-tetrahedron surrounded by the inversion grid center P and the vertex and the tetrahedron ABCD, the weight needed in the interpolation process is obtained, and then the radial velocity V of the grid is obtained by interpolation Pr

[0121]

[0122] In the formula, V Pr is the radial velocity after interpolation, V ΔABCD is the directed volume of the tetrahedron ABCD, V ΔPBCD is the directed volume of the tetrahedron PBCD, V ΔAPCD is the directed volume of the tetrahedron APCD, V ΔABPD is the directed volume of the tetrahedron ABPD, V ΔABCP is the directed volume of the tetrahedron ABCP; V Ar , V Br , V Cr , V Dr are respectively the radial velocities of the radar detection at the vertices of the tetrahedron ABCD.

[0123] In some embodiments, the radar detects intensity data, which is then interpolated by the method of step S04.

[0124] It can be understood that the intensity data detected by the radar and the method of interpolating the radial velocity detected by the radar are the same, and only need to be replaced during calculation.

[0125] It can be understood that in the present application, the weight is calculated only once during the initial operation of the algorithm, and the subsequent interpolation work only needs to be directly looked up according to the weight that has been determined, which saves a lot of calculation cost, avoids additional errors caused by manual parameter setting, and is more conducive to business promotion.

[0126] Exemplarily, please refer to Figure 2 , the case where the inversion grid center P is surrounded by the tetrahedron ABCD.

[0127] In the S05:

[0128] It can be understood that through error analysis, the wind speed error in the inversion result is related to the height of the inversion grid, the lower the height of the inversion grid, the greater the error introduced in the inversion result; and when the height of the inversion grid gradually increases, the error introduced in the inversion gradually decreases. ​

[0129] In some embodiments, the averaged signal-to-noise ratio obtained from each radar in the triangular sub-region in the radar network and the radial velocity data interpolated in the three-dimensional space are respectively set with corresponding frequency-domain two-dimensional filters at different inversion grid height layers, and after all the filters are set, a filter bank is obtained.

[0130] It can be understood that the two-dimensional filter of the present application should be analytical in the design process, and the designed two-dimensional filter corresponds to the inversion grid height layer one by one.

[0131] It can also be understood that all two-dimensional filters in the filter bank should increase the bandwidth with the increase of height.

[0132] Further, the formula of the frequency response of the two-dimensional filter is as follows:

[0133]

[0134] f(h)=k*h+SNR*s

[0135] In the formula, |H h (ω1,ω2)| is the frequency response of the two-dimensional filter, D(ω1,ω2) is the distance from each point of the radial velocity of each radar in the frequency domain to the 0 frequency point, and is expressed as ω1 is the horizontal coordinate of the two-dimensional frequency axis, ω2 is the vertical coordinate of the two-dimensional frequency axis; f(h) is a function related to the height of the inversion grid, h is the height of the wind field inversion; k is a height influencing factor, ranging from 0 to 1; SNR is the averaged signal-to-noise ratio obtained from the radar detection at the top of the triangular sub-region in the inversion process; s is a basic passband influencing factor, ranging from 0 to 1.

[0136] It can be understood that k is a height influencing factor, ranging from 0 to 1, which is used to adjust the sensitivity of the radial velocity two-dimensional filter to height.

[0137] It can be understood that the SNR*s term in the formula is used to provide a basic passband, which is often related to the averaged signal-to-noise ratio obtained from the radar detection. Because the error influence of SNR on the wind field inversion result is the largest, if the signal-to-noise ratio is smaller, the two-dimensional filter should be designed more strictly. And k*h will gradually reduce the strictness of the two-dimensional filter with the increase of height, because the amplification ability of the inversion process to noise decreases at higher inversion height, at this time, the error suppression is not the primary task, but to preserve weather features as much as possible to deal with extreme weather conditions.

[0138] Further, the interpolated radial velocity is subjected to two-dimensional fast Fourier transform to convert to frequency domain, and the frequency domain result of the interpolated radial velocity is multiplied by the two-dimensional filter frequency response in the filter set corresponding to the height of the inversion grid to obtain a frequency domain filtering result, and finally the filtering result is subjected to two-dimensional inverse Fourier transform to obtain a radial velocity two-dimensional filtering result.

[0139] It can be understood that the two-dimensional filter set can effectively correct the vertical velocity component in the final wind field inversion result, and at the same time can further optimize the horizontal component.

[0140] In S06, the following steps are performed:

[0141] It can be understood that, in order to better cooperate with the radar network, before performing three-dimensional wind field inversion, the radial velocity two-dimensional filtering result of the radar belonging to the actual participation in the three-dimensional wind field inversion of the triangular sub-region in the radar network obtained after the radar networking is selected. The selection method is similar to the method of judging whether the grid center is in the tetrahedron, and the difference is that the height information is omitted, and only the projection of each inversion grid center on the xy plane is concerned.

[0142] In some embodiments, the selection of the three-dimensional wind field inversion data includes: calculating the directed area of the triangular sub-region in the radar network, and then calculating the directed area of the triangle surrounded by the projection P'(x, y) of the inversion grid center P on the xy plane and the vertices ABC of the triangular sub-region, respectively, and the formula is as follows:

[0143]

[0144] In the formula, x is the x-axis coordinate of the inversion grid center P in space, y is the y-axis coordinate of the inversion grid center P in space, x1 is the x-axis coordinate of the vertex A of the triangle in space, y1 is the y-axis coordinate of the vertex A of the triangle in space, x2 is the x-axis coordinate of the vertex C of the triangle in space, y2 is the y-axis coordinate of the vertex B of the triangle in space, x3 is the x-axis coordinate of the vertex C of the triangle in space, and y3 is the y-axis coordinate of the vertex C of the triangle in space; and S ΔP BC ΔAP′C ΔABP′

[0145] When the condition S ΔP′BC +S ΔAP′C +S ΔABP′ =S ΔABC is met, the radial velocity data detected by the radar located at the vertices of the triangular sub-region ABC is selected as the three-dimensional wind field inversion data. ​​​​

[0146] It can be understood that when the condition S ΔP′BC +S ΔAP′C +S ΔABP′ =S ΔABC , it indicates that P' is located in the triangular sub-region ABC, then the radial velocity data of the radar probe located at the vertex of the triangular sub-region ABC can be selected for subsequent inversion work corresponding to the inversion grid center P.

[0147] It can be understood that when each grid finds its corresponding radar network, the three-dimensional wind field inversion work can be carried out.

[0148] Further, the inversion grid center P is located in the triangular sub-region ABC, the coordinates of the inversion grid center P are (x, y, z), and the radar coordinates of the vertices of the triangular sub-region ABC are (x1, y1, z1), (x2, y2, z2), and (x3, y3, z3). The direction vectors are constructed, and the formula is as follows:

[0149]

[0150] In the formula, is the direction vector composed of the first part of the radar coordinates in the triangular sub-region and the inversion grid center P, is the direction vector composed of the second part of the radar coordinates in the triangular sub-region and the inversion grid center P, is the direction vector composed of the third part of the radar coordinates in the triangular sub-region and the inversion grid center P;

[0151] The relationship between the two-dimensional filtering results of the radial velocity and the wind speed components is obtained through the vector projection relationship, as follows:

[0152]

[0153] In the formula, V r1 is the radial velocity of the first part of the radar in the triangular sub-region obtained by two-dimensional filtering at the inversion grid center P, V r2 is the radial velocity of the second part of the radar in the triangular sub-region obtained by two-dimensional filtering at the inversion grid center P, V r3 is the radial velocity of the third part of the radar in the triangular sub-region obtained by two-dimensional filtering at the inversion grid center P; u is the to-be-solved north-south direction horizontal wind speed component; v is the to-be-solved east-west direction horizontal wind speed component, and w is the to-be-solved vertical wind speed component; is the direction vector composed of the first part of the radar coordinates in the triangular sub-region and the inversion grid center P, is the direction vector composed of the second part of the radar coordinates in the triangular sub-region and the inversion grid center P, a direction vector composed of the third radar coordinate in the triangular sub-region and the inversion grid center P; 2-norm of the direction vector 2-norm of the direction vector 2-norm of the direction vector 2-norm of the direction vector 2-norm of the direction vector 2-norm of the direction vector

[0154] It can be understood that the inversion equation can be constructed by the above relationship for solving the unknown velocity component uvw.

[0155] The matrix of the unknown velocity component uvw inversion equation calculation process is as follows:

[0156]

[0157] It can be understood that the wind speed component uvw can be obtained by solving the inverse of the matrix A.

[0158] In the formula, V r1 is the radial velocity of the first radar in the triangular sub-region obtained by two-dimensional filtering at the inversion grid center P, V r2 is the radial velocity of the second radar in the triangular sub-region obtained by two-dimensional filtering at the inversion grid center P, V r3 is the radial velocity of the third radar in the triangular sub-region obtained by two-dimensional filtering at the inversion grid center P; x is the x-axis coordinate of the inversion grid center P in space, y is the y-axis coordinate of the inversion grid center P in space, z is the z-axis coordinate of the inversion grid center P in space, x1 is the x-axis coordinate of the first radar position in the triangular sub-region in the Cartesian coordinate system, y1 is the y-axis coordinate of the first radar position in the triangular sub-region in the Cartesian coordinate system, z1 is the altitude of the first radar position in the triangular sub-region in the Cartesian coordinate system, x2 is the x-axis coordinate of the second radar position in the triangular sub-region in the Cartesian coordinate system, y2 is the y-axis coordinate of the second radar position in the triangular sub-region in the Cartesian coordinate system, z2 is the altitude of the second radar position in the triangular sub-region in the Cartesian coordinate system, x3 is the x-axis coordinate of the third radar position in the triangular sub-region in the Cartesian coordinate system, y3 is the y-axis coordinate of the third radar position in the triangular sub-region in the Cartesian coordinate system, z3 is the altitude of the third radar position in the triangular sub-region in the Cartesian coordinate system, a direction vector composed of the first radar coordinate in the triangular sub-region and the inversion grid center P, a direction vector composed of the second radar coordinate in the triangular sub-region and the inversion grid center P, a direction vector composed of the third radar coordinate in the triangular sub-region and the inversion grid center P, 2-norm of the direction vector 2-norm of the direction vector 2-norm of the direction vector - 1 is the inversion matrix.

[0159] It can be understood that the Euclidean distance is used as the weight of the intensity data in the fusion of the three radars.

[0160] The interpolated intensity data is fused, and the formula is as follows:

[0161]

[0162] In the formula, Z is the fused intensity data, Z1 is the intensity of the first radar in the triangular sub-region interpolated at the inversion grid center P, Z2 is the intensity of the first radar in the triangular sub-region interpolated at the inversion grid center P, and Z3 is the intensity of the first radar in the triangular sub-region interpolated at the inversion grid center P. is a direction vector composed of the coordinates of the first radar in the triangular sub-region and the inversion grid center P, is a direction vector composed of the coordinates of the second radar in the triangular sub-region and the inversion grid center P, is a direction vector composed of the coordinates of the third radar in the triangular sub-region and the inversion grid center P, 2-norm of the direction vector 2-norm of the direction vector 2-norm of the direction vector

[0163] It can be understood that the existing way of calculating the weight is used to fuse the interpolated intensity data.

[0164] It can be understood that the weight composed of the intensity data of the radar detection of the inversion grid in the radar network and the Euclidean distance can complete the intensity fusion of the three-dimensional wind field inversion.

[0165] It can be understood that in the wind field inversion process, the true vertical component of the wind vector is also needed to be obtained by subtracting the falling speed of the precipitation particles.

[0166] Further, the true vertical component of the wind vector is as follows:

[0167] w​​​​​​′ =ww t

[0168] w t =2.65(ρ0 / ρ)Z 0 .114

[0169] In the formula, w′ is the true vertical component of the wind vector; w t ρ is the falling velocity of precipitation particles; w is the vertical velocity component obtained from the 3D wind field inversion; Z is the fused intensity data; ρ0 is the air density at the inversion height plane; and is the average air density.

[0170] It is understandable that, regardless of the networking method, data will exist at the edges of the divided radar network, such as... Figure 3 In the process of this invention, this situation is manifested as a directed area S ΔP′BC S ΔAP′C S ΔABP′ The sum of these values ​​is 0, indicating that the centroid P of the inversion grid is located at the edge of two triangular sub-regions. In this case, the selection of radars participating in the inversion involves subjective factors, resulting in a suboptimal selection and often leading to wind vector results that do not meet expectations. Therefore, this invention proposes a method for jointly inverting smoothed wind vectors from adjacent triangular sub-regions.

[0171] Furthermore, inversion grids on the boundaries of adjacent radar network triangular sub-regions are selected, and wind field inversions are performed on these grids using the radars belonging to the adjacent triangular sub-regions of the radar networks on both sides. This yields two results for the wind speed components belonging to the boundary grids, denoted as u1, v1, w′1 and u2, v2, w′2. These components are then combined by averaging, as shown in the following formula:

[0172]

[0173] In the formula, u1 is the north-south horizontal component of the inversion grid result of the first triangular sub-region in the adjacent radar network triangular sub-region, v1 is the east-west horizontal component of the inversion grid result of the first triangular sub-region in the adjacent radar network triangular sub-region, w′1 is the vertical horizontal component of the inversion grid result of the first triangular sub-region in the adjacent radar network triangular sub-region; u2 is the north-south horizontal component of the inversion grid result of the second triangular sub-region in the adjacent radar network triangular sub-region, v2 is the east-west horizontal component of the inversion grid result of the second triangular sub-region in the adjacent radar network triangular sub-region, w′2 is the vertical horizontal component of the inversion grid result of the second triangular sub-region in the adjacent radar network triangular sub-region; u is the composite result of the north-south horizontal components of the inversion grid on the boundary of the adjacent radar network triangular sub-region, v is the composite result of the east-west horizontal components of the inversion grid on the boundary of the adjacent radar network triangular sub-region, and w′ is the composite result of the vertical components of the inversion grid on the boundary of the adjacent radar network triangular sub-region.

[0174] Obtain wind speed magnitude v′, horizontal direction θ, and vertical direction. The formula is as follows:

[0175]

[0176] In the formula, u is the composite result of the north-south horizontal components of the inverted grid, v is the composite result of the east-west horizontal components of the inverted grid, w is the composite result of the vertical components of the inverted grid, and w′ is the composite result of the vertical components of the inverted grid on the boundary of the triangular sub-region of the adjacent radar network.

[0177] It is understandable that, ultimately, by obtaining the wind speed components from each inversion grid, one can directly obtain the wind speed magnitude v′, the horizontal direction θ, and the vertical direction θ.

[0178] Secondly, embodiments of this application provide a three-dimensional wind field inversion system based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing, including:

[0179] A radar site geographic information system is used to obtain geographic information about radar sites.

[0180] A radar networking system is used to network radars based on the geographical information of radar sites using a triangulation method.

[0181] The inversion grid division module is used for inversion grid division after radar networking;

[0182] The interpolation module is used to interpolate the radial velocity data obtained from radar detection.

[0183] a radial velocity two-dimensional filtering result module configured to perform two-dimensional filtering on the interpolation result of the radial velocity data obtained by the radar probe by using the height of the divided inversion grid and the signal-to-noise ratio averaged result of the radar probe to obtain a radial velocity two-dimensional filtering result;

[0184] a three-dimensional wind field inversion module configured to perform three-dimensional wind field inversion by using the inversion grid, the geographic information of the radar site and the radial velocity two-dimensional filtering result.

[0185] In a third aspect, the present application provides a computer device, comprising a storage and a processor, wherein the storage stores a computer program, and the computer program, when executed by the processor, causes the processor to perform the steps of the three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing.

[0186] The computer device can be a desktop computer, a notebook computer, a palm computer, a cloud server or the like. The computer device can interact with a user through a keyboard, a mouse, a remote controller, a touchpad, a voice control device or the like.

[0187] The storage comprises at least one type of readable storage medium, such as a flash memory, a hard disk, a multimedia card, a card-type memory (e.g., an SD or D interface display memory), a random access memory (RAM), a static random access memory (SRAM), a read-only memory (ROM), an electrically erasable programmable read-only memory (EEPROM), a programmable read-only memory (PROM), a magnetic memory, a magnetic disk, an optical disk or the like. In some embodiments, the storage can be an internal storage unit of the computer device, such as a hard disk or a memory of the computer device. In other embodiments, the storage can also be an external storage device of the computer device, such as a plug-in hard disk, a smart media card (SMC), a secure digital (SD) card, a flash card or the like. Of course, the storage can also comprise both the internal storage unit and the external storage device of the computer device. In the present embodiment, the storage is usually used to store an operating system and various application software installed in the computer device, such as program codes of the three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing. In addition, the storage can also be used to temporarily store various data that have been output or will be output.

[0188] In some embodiments, the processor may be a central processing unit (CPU), a controller, a microcontroller, a microprocessor, or other data processing chip. The processor is typically used to control the overall operation of the computer device. In this embodiment, the processor is used to run program code stored in the memory or process data, for example, to run the program code for the three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing.

[0189] Fourthly, this application provides a computer-readable storage medium storing a computer program, which, when executed by a processor, causes the processor to perform the steps of the three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing as described above.

[0190] The computer-readable storage medium stores an interface display program that can be executed by at least one processor to enable the at least one processor to perform the steps of the three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing as described above.

[0191] Application examples

[0192] This application example uses six radars in Guangzhou. The radar network result after triangulation is as follows: Figure 4 As shown.

[0193] After interpolating the radial velocity of each radar, a two-dimensional radial velocity filter was applied to each altitude level. Figure 5 The frequency responses of the two-dimensional filter in this application example at inversion altitudes of 1km, 5km, 10km, 15km, and 20km are respectively... Figure 5 Examples (a), (b), (c), (d), and (e) are shown in the diagram. In this case, the height influence factor k of the radial velocity two-dimensional filter is 0.3, the basic passband influence factor s is 0.2, and the average signal-to-noise ratio of the phased array radar is assumed to be 20 dB. In this filter, the passband portion is sufficient to satisfy most of the actual wind speed variation frequencies.

[0194] During wind field inversion, the edge inversion grids located in the radar network sub-regions need to be averaged according to the region using the method of jointly inverting smooth wind vectors from adjacent triangular sub-regions of this invention. These grids are as follows: Figure 6As shown. After the wind field inversion is completed in each sub-region, the edge grid of region 1 needs to be vector-averaged using the wind field inversion results of the sub-regions enclosed by radars 235 and 231; the edge grid of region 2 needs to be vector-averaged using the wind field inversion results of the sub-regions enclosed by radars 231 and 134; the edge grid of region 3 needs to be vector-averaged using the wind field inversion results of the sub-regions enclosed by radars 134 and 436; the edge grid of region 4 needs to be vector-averaged using the wind field inversion results of the sub-regions enclosed by radars 436 and 635; and the edge grid of region 5 needs to be vector-averaged using the wind field inversion results of the sub-regions enclosed by radars 635 and 532.

[0195] The comparison of the effects before and after filtering at inversion heights of 1km, 5km, and 10km, with the vertical velocity component of the inversion results at 20dB, is shown in the figure. Figure 7 , Figure 7 In the dataset, (a), (b), and (c) represent results over 1 km; (d), (e), and (f) represent results over 5 km; and (g), (h), and (i) represent results over 10 km. Figure 7 The results show that the radial velocity two-dimensional filter has a very good error correction effect in the low-altitude wind field inversion results. Although there is some distortion in the results at an inversion height of 1 km, the overall results after using the radial velocity two-dimensional filter are closer to the ideal prior artificial wind field than the results without filtering. The effect of the radial velocity two-dimensional filter can also be demonstrated by comparing the RMSE at various altitude levels, such as... Figure 8 As shown in the figure, the effect of radial velocity two-dimensional filtering on the uvw components at various heights is compared when the signal-to-noise ratio of the probe data is 20dB and 10dB, respectively. Figure 8 Figures (a) and (b) show the probe signal-to-noise ratio (SNR) at 20 dB; figures (c) and (d) show the probe SNR at 10 dB. Furthermore, the optimized overall 3D wind field inversion process was applied to the actual probe data. The final inversion result of the overall 3D wind field is shown below. Figure 9 As shown, the overall wind field inversion range resembles a pentagon, encompassing the main areas of Guangzhou. Results with inversion heights of 1km, 5km, 9km, 13km, and 17km are displayed. Figure 9The paper demonstrates the improvement in overall wind field inversion results after using a two-dimensional radial velocity filter. It is evident from the low-altitude inversion results at 1km and 5km altitudes that, after filter optimization, the abruptly disappearing vertical component in the low-altitude wind field results is completed, and the discontinuous vertical components are made more continuous. Furthermore, numerous spikes in the descending airflow structure within the wind field are smoothed out after filter optimization, while maintaining the original descending airflow trend. Simultaneously, the wind vector direction, which initially showed obvious errors, is also significantly improved after optimization. Figure 9 It can be seen that although the design of the radial velocity two-dimensional filter requires a more stringent passband in the low-altitude wind field obtained by inversion than in the high-altitude wind field, the actual experimental comparison shows that the optimization results still retain the continuity trend of the original wind field, indicating that the method of the present invention is relatively reliable.

[0196] The three-dimensional wind field inversion results at various altitudes at four different time points were compared with the verification radar detection data for validation. Figure 10 As shown, the effect before and after using radial velocity two-dimensional filtering to correct errors is also demonstrated. After actual radar data acquisition, the inverted three-dimensional wind field before and after radial velocity two-dimensional filtering was verified and compared at four time points, which can be described as time points 1-4, (a) is time point 1, (b) is time point 2, (c) is time point 3, and (d) is time point 4. This can be seen from... Figure 10 The results show that the error in the inversion results initially decreases gradually and then increases as the inversion height increases. This is because the signal-to-noise ratio of the radar detection data decreases with increasing inversion height, and the distribution of detection points becomes sparser, leading to increased errors. Therefore, when comparing the effects of filtering, the analysis of the graphs should focus on the verification results at inversion heights below 15km. After two-dimensional filtering, the errors at each height in the results at these four time points are improved to varying degrees, especially the accuracy of the low-level wind field at inversion heights below 5km, which is significantly improved.

[0197] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for 3D wind field retrieval based on radial velocity 2D filtering and radar network edge wind vector smoothing, characterized in that, The method comprises the following steps: obtaining radar station geographic information; using a triangulation method to perform radar networking according to the radar station geographic information; performing inversion grid division after the radar networking; obtaining radial velocity data through radar detection and performing interpolation; performing two-dimensional filtering on the interpolation result of the radial velocity data obtained through radar detection by using the height of the divided inversion grid and the average result of the signal-to-noise ratio obtained through radar detection, to obtain a radial velocity two-dimensional filtering result; performing three-dimensional wind field inversion by using the inversion grid, the radar station geographic information and the radial velocity two-dimensional filtering result; wherein, according to the average result of the signal-to-noise ratio obtained through each radar in the triangular sub-region in the radar network and the radial velocity data which has been interpolated in the three-dimensional space, a corresponding frequency domain two-dimensional filter is set at the height layer of the inversion grid, and after all the filters are set, a filter bank is obtained; the formula of the frequency response of the two-dimensional filter is as follows: f(h)=k*h+SNR*s In the formula, |H h (ω1,ω2) is a two-dimensional filter frequency response, D(ω1,ω2) is the distance from each point of the radial velocity of each part of the radar in the frequency domain to the 0 point, and is expressed as ω1 is the horizontal coordinate of the two-dimensional frequency axis, ω2 is the vertical coordinate of the two-dimensional frequency axis; f(h) is a function related to the height of the inversion grid, h is the height of the wind field inversion; k is a height influence factor, ranging from 0 to 1; SNR is the average result of the signal-to-noise ratio obtained by the radar detection at the top of the triangular sub-region during the inversion process; s is a basic passband influence factor, and the range is between 0 and 1.

2. The three-dimensional wind field retrieval method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing according to claim 1, characterized in that, The radar station geographic information comprises longitude, latitude and altitude; and / or radar networking obtains a radar network, and in the radar network, the geographic information data of each radar station is converted from a polar coordinate system to a Cartesian coordinate system, and the conversion formula is as follows: x=lcosεsinθ+x0 y=lcosεcosθ+y0 z=lsinε+z0 In the formula, x is the x-axis coordinate in the converted Cartesian coordinate system, y is the y-axis coordinate in the converted Cartesian coordinate system, z is the z-axis coordinate in the converted Cartesian coordinate system, x0 is the x-axis coordinate of the position of the radar station in the Cartesian coordinate system, y0 is the y-axis coordinate of the position of the radar station in the Cartesian coordinate system, z0 is the altitude of the position of the radar station in the Cartesian coordinate system, l is the radial distance of the radar, θ is the elevation angle of the radar, and ε is the azimuth angle of the radar.

3. The three-dimensional wind field retrieval method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing according to claim 2, characterized in that, After converting all the geographic position data of the radar stations from the polar coordinate system to the Cartesian coordinate system, the inversion grid is divided from the radar network; and / or the triangulation method is adopted to triangulate the radial velocity data obtained through radar detection according to the relative Cartesian coordinates, to obtain a plurality of tetrahedrons whose vertices are composed of the point coordinates of the radial velocity data obtained through radar detection and which are closely connected, and the radial velocity data obtained through radar detection used is determined by judging the distribution of the inversion grid barycenter P in the tetrahedrons.

4. The three-dimensional wind field retrieval method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing according to claim 3, characterized in that, The method for judging the distribution of the inversion grid barycenter P in the tetrahedrons comprises the following steps: calculating the volume of the tetrahedron ABCD, and the formula is as follows: In the formula, x1 is the x-axis coordinate of tetrahedron vertex A in space, y1 is the y-axis coordinate of tetrahedron vertex A in space, z1 is the z-axis coordinate of tetrahedron vertex A in space, x2 is the x-axis coordinate of tetrahedron vertex B in space, y2 is the y-axis coordinate of tetrahedron vertex B in space, z2 is the z-axis coordinate of tetrahedron vertex B in space, x3 is the x-axis coordinate of tetrahedron vertex C in space, y3 is the y-axis coordinate of tetrahedron vertex C in space, z3 is the z-axis coordinate of tetrahedron vertex C in space, x4 is the x-axis coordinate of tetrahedron vertex D in space, y4 is the y-axis coordinate of tetrahedron vertex D in space, and z4 is the z-axis coordinate of tetrahedron vertex D in space. ΔABCD is the directed volume of the tetrahedron ABCD. respectively obtaining the volume surrounded by the inversion grid barycenter P and each vertex of the tetrahedron ABCD, and the formula is as follows: In the formula, x is the x-axis coordinate of the inversion grid center P in space, y is the y-axis coordinate of the inversion grid center P in space, z is the z-axis coordinate of the inversion grid center P in space, x1 is the x-axis coordinate of the tetrahedron vertex A in space, y1 is the y-axis coordinate of the tetrahedron vertex A in space, z1 is the z-axis coordinate of the tetrahedron vertex A in space, x2 is the x-axis coordinate of the tetrahedron vertex B in space, y2 is the y-axis coordinate of the tetrahedron vertex B in space, z2 is the z-axis coordinate of the tetrahedron vertex B in space, x3 is the x-axis coordinate of the tetrahedron vertex C in space, y3 is the y-axis coordinate of the tetrahedron vertex C in space, z3 is the z-axis coordinate of the tetrahedron vertex C in space, x4 is the x-axis coordinate of the tetrahedron vertex D in space, y4 is the y-axis coordinate of the tetrahedron vertex D in space, and z4 is the z-axis coordinate of the tetrahedron vertex D in space; V ΔPBCD is the directed volume of the tetrahedron PBCD, V ΔAPCD is the directed volume of the tetrahedron APCD, V ΔABPD is the directed volume of the tetrahedron ABPD, V ΔABCP is the directed volume of the tetrahedron ABCP. Determine whether the inversion grid center P is in the tetrahedron ABCD. When the condition V ΔABCD =V ΔPBCD +V ΔAPCD +V ΔABPD +V ΔABCP is met, it means that the tetrahedron ABCD contains the inversion grid center P; if the condition is not met, it means that the tetrahedron ABCD does not contain the inversion grid center P, and another tetrahedron meeting the condition needs to be found for the inversion grid center P. The weight used in the interpolation process is obtained according to the volume relationship between the sub-tetrahedron surrounded by the inversion grid center P and the vertices and the tetrahedron ABCD, and then the radial velocity V of the grid is obtained by interpolation Pr is: where V Pr is the interpolated radial velocity, V ΔABCD is the directed volume of tetrahedron ABCD, V ΔPBCD is the directed volume of tetrahedron PBCD, V ΔAPCD is the directed volume of tetrahedron APCD, V ΔABPD is the directed volume of tetrahedron ABPD, V ΔABCP is the directed volume of tetrahedron ABCP; V Ar , V Br , V Cr , V Dr are the radial velocities at the vertices of tetrahedron ABCD obtained from radar.

5. The three-dimensional wind field retrieval method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing according to claim 4, characterized in that, performing two-dimensional fast Fourier transform on the interpolated radial velocity, converting to the frequency domain, multiplying the frequency domain result of the interpolated radial velocity and the frequency response of the two-dimensional filter corresponding to the inversion grid height in the filter bank according to the inversion grid height, obtaining the frequency domain filtering result, and finally performing two-dimensional inverse Fourier transform on the filtering result to obtain the radial velocity two-dimensional filtering result.

6. The method according to claim 1, wherein, The selection of three-dimensional wind field inversion data includes: calculating the directed area of a triangular sub-region in the radar network, and then calculating the directed area of the triangle formed by the projection P'(x, y) of the inversion grid center P in the xy axis plane and the vertices A, B and C of the triangular sub-region, respectively, according to the following formula: wherein x is the x-axis coordinate of the center of gravity P of the inversion grid in the space, y is the y-axis coordinate of the center of gravity P of the inversion grid in the space, x1 is the x-axis coordinate of the vertex A of the triangle in the space, y1 is the y-axis coordinate of the vertex A of the triangle in the space, x2 is the x-axis coordinate of the vertex C of the triangle in the space, y2 is the y-axis coordinate of the vertex B of the triangle in the space, x3 is the x-axis coordinate of the vertex C of the triangle in the space, and y3 is the y-axis coordinate of the vertex C of the triangle in the space; and S ΔP′BC is the directed area of the triangle P'BC, S ΔAP′C is the directed area of the triangle AP' C, S ΔABP′ is the directed area of the triangle ABP'. When the condition S ΔP′BC + S ΔAP′C + S ΔABP′ = S ΔABC is satisfied, the radial velocity data obtained by the radar detection at the vertex of the triangular sub-region ABC is selected as the three-dimensional wind field inversion data.

7. The three-dimensional wind field retrieval method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing according to claim 6, characterized in that, The inversion grid center P is located in the triangular sub-region ABC, the coordinates of the inversion grid center P are (x, y, z), and the radar coordinates of the vertices of the triangular sub-region ABC are (x1, y1, z1), (x2, y2, z2) and (x3, y3, z3) respectively. The direction vector is constructed according to the following formula: wherein is a direction vector of the first part of the radar coordinates in the triangular sub-region and the inversion grid center P, is a direction vector of the second part of the radar coordinates in the triangular sub-region and the inversion grid center P, is a direction vector of the third part of the radar coordinates in the triangular sub-region and the inversion grid center P. The relationship between the two-dimensional filtering result of radial velocity and the wind speed component is obtained through the vector projection relationship, as follows: In the formula, V r1 V is the radial velocity obtained by two-dimensional filtering of the first radar in the triangular sub-region at the centroid P of the inversion grid. r2 V is the radial velocity obtained by two-dimensional filtering of the second radar in the triangular subregion at the centroid P of the inversion grid. r3 denoted as radial velocity obtained by two-dimensional filtering of the third radar in the triangular sub-region at the centroid P of the inversion grid; u is the horizontal wind speed component to be determined in the north-south direction; v is the horizontal wind speed component to be determined in the east-west direction; and w is the vertical wind speed component to be determined. The direction vector formed by the first radar coordinates in the triangular sub-region and the centroid P of the inversion grid. The direction vector formed by the second radar coordinates in the triangular subregion and the centroid P of the inversion grid. The direction vector formed by the third radar coordinates in the triangular sub-region and the centroid P of the inversion grid; Represented as direction vector 2 Fan, Represented as direction vector 2 Fan, Represented as direction vector 2. The matrix of the unknown speed component uvw inversion equation calculation process is as follows: In the formula, V r1 is the radial velocity of the first radar in the triangular sub-region obtained by two-dimensional filtering at the inversion grid center P, r2 is the radial velocity of the second radar in the triangular sub-region obtained by two-dimensional filtering at the inversion grid center P, r3 is the radial velocity of the third radar in the triangular sub-region obtained by two-dimensional filtering at the inversion grid center P; x is the x-axis coordinate of the inversion grid center P in space, y is the y-axis coordinate of the inversion grid center P in space, z is the z-axis coordinate of the inversion grid center P in space, x1 is the x-axis coordinate of the position of the first radar in the triangular sub-region in the Cartesian coordinate system, y1 is the y-axis coordinate of the position of the first radar in the triangular sub-region in the Cartesian coordinate system, z1 is the altitude of the position of the first radar in the triangular sub-region in the Cartesian coordinate system, x2 is the x-axis coordinate of the position of the second radar in the triangular sub-region in the Cartesian coordinate system, y2 is the y-axis coordinate of the position of the second radar in the triangular sub-region in the Cartesian coordinate system, z2 is the altitude of the position of the second radar in the triangular sub-region in the Cartesian coordinate system, x3 is the x-axis coordinate of the position of the third radar in the triangular sub-region in the Cartesian coordinate system, y3 is the y-axis coordinate of the position of the third radar in the triangular sub-region in the Cartesian coordinate system, z3 is the altitude of the position of the third radar in the triangular sub-region in the Cartesian coordinate system, is the direction vector composed of the coordinates of the first radar in the triangular sub-region and the inversion grid center P, is the direction vector composed of the coordinates of the second radar in the triangular sub-region and the inversion grid center P, is the direction vector composed of the coordinates of the third radar in the triangular sub-region and the inversion grid center P, represents the 2-norm of the direction vector represents the 2-norm of the direction vector represents the 2-norm of the direction vector u is the to-be-solved north-south direction horizontal wind speed component, v is the to-be-solved east-west direction horizontal wind speed component, w is the to-be-solved vertical wind speed component; A is a simplified description of the matrix in the formula, A -1 is the inversion matrix.​​ 8. The three-dimensional wind field retrieval method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing according to claim 7, characterized in that, The intensity data is obtained by radar detection, and then interpolation is performed. The interpolated intensity data is fused according to the following formula: In the formula, Z represents the fused intensity data, Z1 represents the intensity obtained by interpolating the first radar in the triangular sub-region at the centroid P of the inversion grid, Z2 represents the intensity obtained by interpolating the first radar in the triangular sub-region at the centroid P of the inversion grid, and Z3 represents the intensity obtained by interpolating the first radar in the triangular sub-region at the centroid P of the inversion grid. The direction vector formed by the first radar coordinates in the triangular sub-region and the centroid P of the inversion grid. The direction vector formed by the second radar coordinates in the triangular subregion and the centroid P of the inversion grid. The direction vector formed by the third radar coordinates in the triangular sub-region and the centroid P of the inversion grid. Represented as direction vector 2 Fan, Represented as direction vector 2 Fan, Represented as direction vector 2.

9. The three-dimensional wind field retrieval method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing according to claim 8, characterized in that, The true vertical component of the wind vector is as follows: w' = w - w t w t = 2.65 (p0 / p) Z 0.114 where w' is the true vertical component of the wind vector; w t is the fall velocity of the precipitation particles; w is the vertical velocity component obtained from the three-dimensional wind field retrieval; Z is the fused intensity data; p0 is the air density at the height of the inversion; p is the average air density.

10. The method of Claim 1, wherein, The inversion grid on the boundary of the adjacent radar network triangular sub-region is selected, and the two adjacent radar network triangular sub-regions are used to perform wind field inversion on the inversion grid respectively, to obtain two results of the wind speed component belonging to the boundary grid, which are expressed as u1, v1, w'1 and u2, v2, w'2. Then the components are synthesized in an average manner, as follows: In the formula, u1 is the north-south horizontal component of the inversion grid result of the first triangular sub-region in the adjacent radar network triangular sub-region, v1 is the east-west horizontal component of the inversion grid result of the first triangular sub-region in the adjacent radar network triangular sub-region, w'1 is the vertical component of the inversion grid result of the first triangular sub-region in the adjacent radar network triangular sub-region; u2 is the north-south horizontal component of the inversion grid result of the second triangular sub-region in the adjacent radar network triangular sub-region, v2 is the east-west horizontal component of the inversion grid result of the second triangular sub-region in the adjacent radar network triangular sub-region, w'2 is the vertical component of the inversion grid result of the second triangular sub-region in the adjacent radar network triangular sub-region; u is the north-south horizontal component synthesis result of the inversion grid on the boundary of the adjacent radar network triangular sub-region, v is the east-west horizontal component synthesis result of the inversion grid on the boundary of the adjacent radar network triangular sub-region, w' is the vertical component synthesis result of the inversion grid on the boundary of the adjacent radar network triangular sub-region; The wind speed magnitude v', the horizontal direction Θ, and the vertical direction of the formula, as follows:

Images