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

Abstract

The invention discloses a three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing. The method comprises the following steps: obtaining geographic information of a radar station; radar networking is carried out by using a triangulation method according to the geographic information of the radar stations; performing inversion grid division after radar networking; radar detection is carried out to obtain radial speed data, and interpolation is carried out; obtaining a radial speed two-dimensional filtering result; and performing three-dimensional wind field inversion through the inversion grid, the radar station geographic information and the radial speed two-dimensional filtering result. Compared with an existing three-dimensional wind field inversion algorithm taking a variational method as a core, the center of gravity is put on radial speed optimization of each radar before wind field inversion, so that the calculation burden in the inversion process is reduced, a radar networking method and an interpolation algorithm in the overall three-dimensional wind field inversion are optimized, and the calculation efficiency is improved. The calculation complexity is reduced while the inversion precision is further improved, and the method is more suitable for popularization in actual business.

Description

Technical Field

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

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

[0003] When constructing a three-dimensional wind field based on the positional relationship between the radial velocity measured by the radar in the radial direction and the three-dimensional component of the wind field vector in the atmosphere, the wind field results obtained have large vertical errors in the low-level results. In order to improve this problem, most of the current three-dimensional wind field inversion methods use the three-dimensional variational method as the core algorithm. The three-dimensional variational method has not been widely promoted in actual meteorological services due to the lack of standardization in the use of the cost function of the algorithm, the difficulty in determining the optimal solution, the limited improvement of the vertical velocity error problem, and the high computational complexity. Summary of the invention

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

[0005] The embodiment of the present application is implemented by providing a three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing, including: obtaining radar site geographic information;

[0006] Radar networking is carried out using triangulation method based on the geographic information of radar sites;

[0007] After radar networking, inversion grid division is performed;

[0008] Radar detection obtains radial velocity data and performs interpolation;

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

[0010] The three-dimensional wind field inversion is performed by combining the inversion grid with the geographic information of the radar sites and the two-dimensional filtering results of the radial velocity.

[0011] In some embodiments,.

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

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

[0014] x=lcosεsinθ+x 0

[0015] y=lcosεcosθ+y 0

[0016] z=lsinε+z 0

[0017] Where x is the x-axis coordinate in the converted Cartesian coordinate system, y is the y-axis coordinate in the converted Cartesian coordinate system, and z is the z-axis coordinate in the converted Cartesian coordinate system. 0 is the x-axis coordinate of the radar site location in the Cartesian coordinate system, and y 0 is the y-axis coordinate of the radar site in the Cartesian coordinate system, z 0 is the altitude of the radar site 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 the geographic location data of all radar sites from a polar coordinate system to a 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 radar detection is triangulated according to relative Cartesian coordinates to obtain a tetrahedron whose vertices are formed by the coordinates of multiple points of the radial velocity data obtained by radar detection and are closely connected, and the radial velocity data obtained by radar detection to be used is determined by judging the distribution of the centroid P of the inversion grid in the tetrahedron; and / or

[0020] According to the average result of the signal-to-noise ratio obtained by the detection of each radar in the triangular sub-area of ​​the radar network and the radial velocity data that has been interpolated in the three-dimensional space, corresponding frequency domain two-dimensional filters are set at different height layers of the inversion grid. After all the settings are completed, a filter group is obtained.

[0021] In some embodiments, the method for determining whether the inversion grid centroid P is distributed in a tetrahedron includes:

[0022] Calculate the volume of tetrahedron ABCD using the following formula:

[0023]

[0024] In the formula, x1 is the x-axis coordinate of the tetrahedron vertex A in space, y 1 is the y-axis coordinate of the tetrahedron vertex A in space, z 1 is the z-axis coordinate of the tetrahedron vertex A in space, x 2 is the x-axis coordinate of the tetrahedron vertex B in space, y 2 is the y-axis coordinate of the tetrahedron vertex B in space, z 2 is the z-axis coordinate of the tetrahedron vertex B in space, x 3 is the x-axis coordinate of the tetrahedron vertex C in space, y 3 is the y-axis coordinate of the tetrahedron vertex C in space, z 3 is the z-axis coordinate of the tetrahedron vertex C in space, x 4 is the x-axis coordinate of the tetrahedron vertex D in space; y 4 is the y-axis coordinate of the tetrahedron vertex D in space, z 4 is the z-axis coordinate of the tetrahedron vertex D in space; V ΔABCD Then it is the directed volume of tetrahedron ABCD;

[0025] The volume enclosed by the centroid P of the inversion grid and each vertex of the tetrahedron ABCD is obtained respectively, and the formula is as follows:

[0026]

[0027]

[0028] Where x is the x-axis coordinate of the inversion grid gravity center P in space, y is the y-axis coordinate of the inversion grid gravity center P in space, z is the z-axis coordinate of the inversion grid gravity center P in space, x 1 is the x-axis coordinate of the tetrahedron vertex A in space, y 1 is the y-axis coordinate of the tetrahedron vertex A in space, z 1 is the z-axis coordinate of the tetrahedron vertex A in space, x 2 is the x-axis coordinate of the tetrahedron vertex B in space, y 2 is the y-axis coordinate of the tetrahedron vertex B in space, z 2 is the z-axis coordinate of the tetrahedron vertex B in space, x 3 is the x-axis coordinate of the tetrahedron vertex C in space, y 3 is the y-axis coordinate of the tetrahedron vertex C in space, z 3 is the z-axis coordinate of the tetrahedron vertex C in space, x 4 is the x-axis coordinate of the tetrahedron vertex D in space; y 4 is the y-axis coordinate of the tetrahedron vertex D in space, z 4 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 tetrahedral APCD, V ΔABPD is the directed volume of tetrahedron ABCD, V ΔABCP is the directed volume of the tetrahedron ABCP;

[0029] Determine whether the inversion grid centroid P is within the tetrahedron ABCD. If the condition V is met, ΔABCD =V ΔPBCD +V ΔAPCD +V ΔABPD +V ΔABCP When , it means that the tetrahedron ABCD contains the inversion grid gravity center P; if this condition is not met, it means that there is no inclusion relationship between the tetrahedron ABCD and the inversion grid gravity center P, and it is necessary to find another tetrahedron that meets the conditions for the inversion grid gravity center P;

[0030] According to the volume relationship between the inversion grid center P and the sub-tetrahedron and tetrahedron ABCD enclosed by the vertices, the weights to be used in the interpolation process are obtained. Then the radial velocity V of the grid is obtained by interpolation. Pr for:

[0031]

[0032] Where V Pr is the interpolated radial velocity, V ΔABCD is the directed volume of tetrahedron ABCD, V ΔPBCD is the directed volume of the tetrahedron PBCD, V ΔAPCD is the directed volume of the tetrahedral APCD, V ΔABPD is the directed volume of the tetrahedral ABPD, V ΔABCP is the directed volume of the tetrahedron ABCP; V Ar , V Br , V Cr , V Dr are respectively the radial velocities detected by radar at the vertices of tetrahedron ABCD; and / or

[0033] The formula for the frequency response of a 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 radar radial velocity in the frequency domain to the 0 frequency point, expressed as ω 1 is the horizontal coordinate of the two-dimensional frequency axis, ω 2 is the ordinate of the two-dimensional frequency axis; f(h) is a function related to the height of the inversion grid, h is the wind field inversion height; k is the height influence factor, ranging from 0 to 1; SNR is the average signal-to-noise ratio obtained by radar detection at the vertex of the triangular sub-area during the inversion process; s is the 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 and converted to the frequency domain. The frequency domain result of the interpolated radial velocity is multiplied by the two-dimensional filter frequency response in the filter group corresponding to the inversion grid height according to the inversion grid height to obtain a filtering result in the frequency domain. Finally, the filtering result is subjected to a two-dimensional inverse Fourier transform to obtain a two-dimensional filtering result of the radial velocity.

[0038] In some embodiments, the selection of three-dimensional wind field inversion data includes: calculating the directed area of ​​the triangular sub-area in the radar network, and then calculating the directed area of ​​the triangle enclosed by the projection P′(x, y) of the inversion grid centroid P on the xy axis plane and the vertices ABC of the triangular sub-area, respectively, according to the following formula:

[0039]

[0040] Where x is the x-axis coordinate of the inversion grid gravity center P in space, y is the y-axis coordinate of the inversion grid gravity center P in space, 1 is the x-axis coordinate of the triangle vertex A in space, y 1 is the y-axis coordinate of the triangle vertex A in space, x 2 is the x-axis coordinate of the triangle vertex C in space, y 2 is the y-axis coordinate of the triangle vertex B in space, x 3 is the x-axis coordinate of the triangle vertex C in space, y 3 is the y-axis coordinate of the triangle vertex C in space; and S ΔP ' BC is the directed area of ​​triangle P′BC, S ΔAP′C is the directed area of ​​triangle AP′C, S ΔABP′ is the directed area of ​​triangle ABP′;

[0041] When the condition S is met ΔP ' BC +S ΔAP′C +S ΔABP′ =S ΔABCWhen , the radar detection at the vertex of the triangular sub-area ABC can be selected to obtain the radial velocity data as the three-dimensional wind field inversion data.

[0042] In some embodiments, the inversion grid centroid P is located in the triangular sub-region ABC, the coordinates of the inversion grid centroid P are (x, y, z), and the radar coordinates of the vertices of the triangular sub-region ABC are (x, y, z). 1 ,y 1 ,z 1 ), (x 2 ,y 2 ,z 2 ), (x 3 ,y 3 ,z 3 ), construct the direction vector, the formula is as follows:

[0043]

[0044] In the formula, is the direction vector composed of the coordinates of the first radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the coordinates of the second radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the coordinates of the third radar in the triangular sub-area and the centroid P of the inversion grid;

[0045] The relationship between the radial velocity two-dimensional filtering result and the wind speed component is obtained through the relationship of vector projection, as shown in the following formula:

[0046]

[0047] Where V r1 is the radial velocity of the first radar in the triangular sub-area obtained by two-dimensional filtering on the inversion grid centroid P, V r2 is the radial velocity of the second radar in the triangular sub-area obtained by two-dimensional filtering on the inversion grid centroid P, V r3 is the radial velocity obtained by two-dimensional filtering of the third radar in the triangular sub-area on the inversion grid centroid P; 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; is the direction vector composed of the coordinates of the first radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the coordinates of the second radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the coordinates of the third radar in the triangular sub-area and the centroid P of the inversion grid; Represented as a direction vector 2 Fan, Represented as a direction vector 2 Fan, Represented as a direction vector 2Fan;

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

[0049]

[0050] Where V r1 is the radial velocity of the first radar in the triangular sub-area obtained by two-dimensional filtering on the inversion grid centroid P, V r2 is the radial velocity of the second radar in the triangular sub-area obtained by two-dimensional filtering on the inversion grid centroid P, V r3 is the radial velocity of the third radar in the triangular sub-area obtained by two-dimensional filtering on the inversion grid gravity center P; x is the x-axis coordinate of the inversion grid gravity center P in space, y is the y-axis coordinate of the inversion grid gravity center P in space, z is the z-axis coordinate of the inversion grid gravity center P in space, x 1 is the x-axis coordinate of the first radar position in the triangular sub-area in the Cartesian coordinate system, and y 1 is the y-axis coordinate of the first radar position in the triangular sub-area in the Cartesian coordinate system, z 1 is the altitude of the first radar position in the triangular sub-area in the Cartesian coordinate system, x 2 is the x-axis coordinate of the second radar position in the triangular sub-area in the Cartesian coordinate system, and y 2 is the y-axis coordinate of the second radar position in the triangular sub-area in the Cartesian coordinate system, z 2 is the altitude of the second radar position in the triangular sub-area in the Cartesian coordinate system, x 3 is the x-axis coordinate of the third radar position in the triangular sub-area in the Cartesian coordinate system, and y 3 is the y-axis coordinate of the third radar position in the triangular sub-area in the Cartesian coordinate system, z 3 is the altitude of the third radar position in the triangular sub-area in the Cartesian coordinate system, is the direction vector composed of the coordinates of the first radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the coordinates of the second radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the third radar coordinate in the triangular sub-area and the inversion grid centroid P, Represented as a direction vector 2 Fan, Represented as a direction vector 2 Fan, Represented as a direction vector The 2-norm of the formula, u is the horizontal wind speed component in the north-south direction, v is the horizontal wind speed component in the east-west direction, and w is the vertical wind speed component. A is a simplified description of the matrix in the formula. - 1 is the inversion matrix.

[0051] In some embodiments, radar detection obtains intensity data, which is then interpolated and fused using the following formula:

[0052]

[0053] Where Z is the fused intensity data, Z 1 is the intensity of the first radar in the triangular sub-area interpolated at the inversion grid centroid P, Z 2 is the intensity of the first radar in the triangular sub-area interpolated at the inversion grid centroid P, Z 3 is the intensity of the first radar in the triangular sub-area interpolated at the center of the inversion grid P; is the direction vector composed of the coordinates of the first radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the coordinates of the second radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the third radar coordinate in the triangular sub-area and the inversion grid centroid P, Represented as a direction vector 2 Fan, Represented as a direction vector 2 Fan, Represented as a direction vector 2Fan.

[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] Where 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 by inversion of the three-dimensional wind field; Z is the fused intensity data; ρ 0 is the air density at the inversion height plane; ρ is the average air density.

[0058] In some embodiments, an inversion grid on the boundary of an adjacent radar network triangular sub-area is selected, and the radars belonging to the adjacent radar network triangular sub-areas on both sides are used to perform wind field inversion on the inversion grid respectively, and two results of wind speed components belonging to the boundary grid are obtained, which are expressed as u 1 , v 1 , w′ 1 and u 2 , v 2 , w′ 2 , and then synthesize the components in an average manner, the formula is as follows:

[0059]

[0060] In the formula, u 1 is the north-south horizontal component of the inversion grid result of the first triangular sub-area in the adjacent radar network triangular sub-area, v 1 The east-west horizontal component of the inversion grid result of the first triangular sub-region of the adjacent radar network triangular sub-region, w′ 1 u is the vertical square component of the inversion grid result of the first triangular sub-region of the adjacent radar network triangular sub-region; 2 is the north-south horizontal component of the inversion grid result of the second triangular sub-area in the adjacent radar network triangular sub-area, v 2 The east-west horizontal component of the inversion grid result of the second triangular sub-region of the adjacent radar network triangular sub-region, w′ 2 is the vertical component of the inversion grid result of the second triangular sub-region of the adjacent radar network triangular sub-region; u is the synthesis 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 synthesis 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 synthesis result of the vertical component of the inversion grid on the boundary of the adjacent radar network triangular sub-region;

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

[0062]

[0063] Where u is the synthesis result of the north-south horizontal component of the inversion grid, v is the synthesis result of the east-west horizontal component of the inversion grid, w is the synthesis result of the vertical component of the inversion grid, and w′ is the synthesis result of the vertical component of the inversion grid on the boundary of the triangular sub-area of ​​the adjacent radar network.

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

[0065] Radar site geographic information system, used to obtain radar site geographic information;

[0066] A radar networking system, used for radar networking using a triangulation method based on the geographic information of radar sites;

[0067] Inversion grid division module, used for inversion grid division after radar networking;

[0068] Interpolation module, used for obtaining radial velocity data through radar detection and performing interpolation;

[0069] The radial velocity two-dimensional filtering result module is used to perform two-dimensional filtering on the interpolation result of the radial velocity data obtained by radar detection through the result of averaging the height of the divided inversion grid and the signal-to-noise ratio obtained by radar detection, so as to obtain the radial velocity two-dimensional filtering result;

[0070] The three-dimensional wind field inversion module is used to perform three-dimensional wind field inversion through the inversion grid, radar site geographic information and radial velocity two-dimensional filtering results.

[0071] Accordingly, an embodiment of the present application also provides a computer device, including a storage device and a processor, wherein the storage device stores a computer program, and when the computer program is executed by the processor, the processor executes the steps of the above method.

[0072] Accordingly, an embodiment of the present application also provides a computer-readable storage medium storing a computer program, and when the computer program is executed by a processor, the processor executes the steps of the above method.

[0073] In summary, due to the adoption of the above technical solution, the beneficial effects of the present invention are:

[0074] Compared with the existing three-dimensional wind field inversion algorithm based on the variational method, the present invention 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 also optimizes the radar networking method and interpolation algorithm in the overall three-dimensional wind field inversion, further improving the inversion accuracy while reducing the calculation complexity, and is more suitable for promotion in actual business. BRIEF DESCRIPTION OF THE DRAWINGS

[0075] Figure 1 A flow chart of a three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing provided by an embodiment of the present invention;

[0076] Figure 2 A geometric relationship diagram between the centroid P of the inversion grid provided in an embodiment of the present invention and the tetrahedron obtained by segmentation;

[0077] Figure 3A data graph present at the edge of a divided radar network provided by an embodiment of the present invention;

[0078] Figure 4 A radar network diagram obtained by division provided by an embodiment of the present invention;

[0079] Figure 5 A diagram showing some two-dimensional filters provided by an embodiment of the present invention;

[0080] Figure 6 A radar network edge grid area diagram divided according to an embodiment of the present invention;

[0081] Figure 7 A comparison diagram of the results before and after two-dimensional filtering provided by an embodiment of the present invention;

[0082] Figure 8 A comparison chart of the accuracy before and after two-dimensional filtering provided by an embodiment of the present invention;

[0083] Fig. 9 A comparison diagram of actual wind field inversion results before and after filtering provided by an embodiment of the present invention;

[0084] Fig.10 This is an additional radar verification result diagram provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0085] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0086] The technical solution of this application is as follows:

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

[0088] S01. Obtaining the geographic information of radar sites;

[0089] S02. Use triangulation method to network radars according to the geographic information of radar sites;

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

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

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

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

[0094] Compared with the existing three-dimensional wind field inversion algorithm based on the variational method, the present invention focuses on the radial velocity optimization of each radar before the wind field inversion, thereby reducing the computational burden in the inversion process, and also optimizes the radar networking method and interpolation algorithm in the overall three-dimensional wind field inversion, further improving the inversion accuracy while reducing the computational complexity, and is more suitable for promotion in actual business. In addition, the vertical velocity error of the low-altitude wind field in the inversion result is also solved, and the inversion error of the horizontal velocity component is optimized to a certain extent.

[0095] The present invention re-examines the entire wind field inversion process and proposes a three-dimensional wind field inversion construction method based on two-dimensional filtering of radial velocities at different heights based on a large number of error analyses. It not only proposes an algorithm with two-dimensional filtering of radial velocities at different inversion heights as the core to solve the problem of obvious errors in low-altitude inversion results.

[0096] Finally, the overall wind field inversion process was 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 present invention was verified through additional S-band Doppler radar detection data. At the same time, the necessity and reliability of the proposed radial velocity two-dimensional filtering were proved by comparing the effects after two-dimensional filtering of radial velocities at different heights.

[0097] In the S01:

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

[0099] In the S02:

[0100] It can be understood that if the radar is compared to an endpoint on a plane, the Bowyer-Watson algorithm can ensure that the sum of the minimum internal angles of the divided triangular sub-areas is greater than the sum of the minimum internal angles of any triangle formed by non-Delaunay triangulation without the occurrence of singularities in the divided triangular area. This property will make the sub-areas surrounded by the divided radar network closest to an equilateral triangle. It can also ensure that the divided area is unique, so the divided sub-areas are completely determined after a network is established, and there is no need to triangulate each time the wind field is inverted. And whether it is adding, deleting or moving radar sites in the future, it will only affect the adjacent sub-areas, and will not affect the networking results 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 said S03:

[0102] It can be understood that after the radars are networked, the inversion grid can be divided. First, the coordinate system needs to be converted according to the geographical location of each radar. The geographical location of each radar in the radar network can be expressed as (x 0 ,y 0 ,z 0 ), respectively representing the conversion of the absolute latitude and longitude coordinates of the radar to the relative coordinates of the distance in the relative Cartesian coordinate system (x 0 ,y 0 ) and its altitude z 0 .

[0103] In some embodiments, the radar network is obtained by radar networking. In the radar network, the geographic information data of each radar site is converted from a polar coordinate system to a Cartesian coordinate system. The conversion formula is as follows:

[0104] x=lcosεsinθ+x 0

[0105] y=lcosεcosθ+y 0

[0106] z=lsinε+z 0

[0107] Where x is the x-axis coordinate in the converted Cartesian coordinate system, y is the y-axis coordinate in the converted Cartesian coordinate system, and z is the z-axis coordinate in the converted Cartesian coordinate system. 0 is the x-axis coordinate of the radar site location in the Cartesian coordinate system, and y 0 is the y-axis coordinate of the radar site in the Cartesian coordinate system, z 0 is the altitude of the radar site 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] Furthermore, after converting the geographic location data of all radar sites from the polar coordinate system to the Cartesian coordinate system, the inversion grid is divided from the radar network.

[0109] In said S04:

[0110] It can be understood that the detection data of each radar cannot be aligned with the regular inversion grid after the coordinate system is converted, and the distribution of the detection data is also uneven. The distribution of radial velocity data will show that the closer to the radar, the denser the data distribution, and the farther away from the radar, the sparser the data distribution. Therefore, it is necessary to unify the scattered radial velocity data in the coordinate system to the regular inversion grid by interpolation.

[0111] In some embodiments, a triangulation method is used to triangulate the radial velocity data obtained by radar detection according to relative Cartesian coordinates to obtain a tetrahedron whose vertices are formed by the point coordinates of multiple radial velocity data obtained by radar detection and are closely connected. The radial velocity data obtained by radar detection to be used is determined by judging the distribution of the inversion grid centroid P in the tetrahedron.

[0112] Further, the method for determining whether the inversion grid centroid P is distributed in the tetrahedron includes:

[0113] S041. Calculate the volume of tetrahedron ABCD using the following formula:

[0114]

[0115] In the formula, x 1 is the x-axis coordinate of the tetrahedron vertex A in space, y 1 is the y-axis coordinate of the tetrahedron vertex A in space, z 1 is the z-axis coordinate of the tetrahedron vertex A in space, x 2 is the x-axis coordinate of the tetrahedron vertex B in space, y 2 is the y-axis coordinate of the tetrahedron vertex B in space, z 2 is the z-axis coordinate of the tetrahedron vertex B in space, x 3 is the x-axis coordinate of the tetrahedron vertex c in space, y 3 is the y-axis coordinate of the tetrahedron vertex C in space, z 3 is the z-axis coordinate of the tetrahedron vertex C in space, x 4 is the x-axis coordinate of the tetrahedron vertex D in space; y 4 is the y-axis coordinate of the tetrahedron vertex D in space, z 4 is the z-axis coordinate of the tetrahedron vertex D in space; V ΔABCD Then it is the directed volume of tetrahedron ABCD;

[0116] S042. Obtain the volume enclosed by the centroid P of the inversion grid and each vertex of the tetrahedron ABCD respectively, the formula is as follows:

[0117]

[0118] Where x is the x-axis coordinate of the inversion grid gravity center P in space, y is the y-axis coordinate of the inversion grid gravity center P in space, z is the z-axis coordinate of the inversion grid gravity center P in space, x 1 is the x-axis coordinate of the tetrahedron vertex A in space, y 1 is the y-axis coordinate of the tetrahedron vertex A in space, z 1 is the z-axis coordinate of the tetrahedron vertex A in space, x 2is the x-axis coordinate of the tetrahedron vertex B in space, y 2 is the y-axis coordinate of the tetrahedron vertex B in space, z 2 is the z-axis coordinate of the tetrahedron vertex B in space, x 3 is the x-axis coordinate of the tetrahedron vertex C in space, y 3 is the y-axis coordinate of the tetrahedron vertex C in space, z 3 is the z-axis coordinate of the tetrahedron vertex C in space, x 4 is the x-axis coordinate of the tetrahedron vertex D in space; y 4 is the y-axis coordinate of the tetrahedron vertex D in space, z 4 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 tetrahedral APCD, V ΔABPD is the directed volume of the tetrahedral ABPD, V ΔABCP is the directed volume of the tetrahedron ABCP;

[0119] S043. Determine whether the inversion grid centroid P is within the tetrahedron ABCD. If the condition V is met, ΔABCD =V ΔPBCD +V ΔAPCD +V ΔABPD +V ΔABCP When , it means that the tetrahedron ABCD contains the inversion grid gravity center P; if this condition is not met, it means that there is no inclusion relationship between the tetrahedron ABCD and the inversion grid gravity center P, and it is necessary to find another tetrahedron that meets the conditions for the inversion grid gravity center P;

[0120] S044. The weights to be used in the interpolation process are obtained according to the volume relationship between the inversion grid centroid P and the sub-tetrahedron and tetrahedron ABCD enclosed by the vertices. Then the radial velocity V of the grid is obtained by interpolation. Pr for:

[0121]

[0122] Where V Pr is the interpolated radial velocity, V ΔABCD is the directed volume of tetrahedron ABCD, V ΔPBCD is the directed volume of the tetrahedron PBCD, V ΔAPCD is the directed volume of the tetrahedral APCD, V ΔABPD is the directed volume of the tetrahedral ABPD, V ΔABCP is the directed volume of the tetrahedron ABCP; V Ar , V Br , V Cr , V DrThey are the radial velocities detected by the radar at the vertices of tetrahedron ABCD.

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

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

[0125] It can be understood that in the present invention, when the algorithm is run for the first time, the weights only need to be calculated once during the segmentation, and the subsequent interpolation work only needs to directly look up the table according to the already determined weights. While saving a lot of computing costs, it avoids the additional errors caused by manually setting parameters, which is more conducive to business promotion.

[0126] For example, see Figure 2 , the inversion grid centroid P is exactly surrounded by the tetrahedron ABCD.

[0127] In said S05:

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

[0129] In some embodiments, corresponding frequency domain two-dimensional filters are set at different height layers of the inversion grid according to the average result of the signal-to-noise ratio obtained by detection of each radar in the triangular sub-area of ​​the radar network and the radial velocity data that has been interpolated in the three-dimensional space. After all the settings are completed, a filter group is obtained.

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

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

[0132] Furthermore, the formula for the frequency response of a 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 radar radial velocity in the frequency domain to the 0 frequency point, expressed as ω 1 is the horizontal coordinate of the two-dimensional frequency axis, ω 2 is the ordinate of the two-dimensional frequency axis; f(h) is a function related to the height of the inversion grid, h is the wind field inversion height; k is the height influence factor, ranging from 0 to 1; SNR is the average signal-to-noise ratio obtained by radar detection at the vertex of the triangular sub-area during the inversion process; s is the basic passband influence factor, ranging from 0 to 1.

[0136] It can be understood that k is a height influence factor ranging from 0 to 1, and 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 result of averaging the signal-to-noise ratio obtained by radar detection. Because SNR has the greatest impact on the error of wind field inversion results, the smaller the signal-to-noise ratio, the stricter the two-dimensional filter should be designed. And k*h will gradually reduce the strictness of the two-dimensional filter as the altitude increases, because at higher inversion altitudes, the inversion process's ability to amplify noise is reduced. At this time, suppressing errors is no longer the primary task, but to retain weather characteristics as much as possible to cope with extreme weather conditions.

[0138] Furthermore, the interpolated radial velocity is subjected to a two-dimensional fast Fourier transform and converted to the frequency domain. The frequency domain result of the interpolated radial velocity is multiplied by the two-dimensional filter frequency response in the filter group corresponding to the inversion grid height according to the inversion grid height to obtain the filtering result in the frequency domain. Finally, the filtering result is subjected to a two-dimensional inverse Fourier transform to obtain the two-dimensional filtering result of the radial velocity.

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

[0140] In said S06:

[0141] It can be understood that in order to better cooperate with the radar network, before performing the three-dimensional wind field inversion, it is necessary to select the radial velocity two-dimensional filtering results of the radar that actually participates in the three-dimensional wind field inversion of the area in combination with the triangular sub-areas in the radar network obtained after the radar networking. The selection method is similar to the aforementioned method of determining whether the centroid of the grid is within the tetrahedron. The difference is that the height information is omitted, and only the projection of the centroid of each inversion grid on the xy plane needs to be paid attention to.

[0142] In some embodiments, the selection of three-dimensional wind field inversion data includes: calculating the directed area of ​​the triangular sub-area in the radar network, and then calculating the directed area of ​​the triangle enclosed by the projection P′(x, y) of the inversion grid centroid P on the xy axis plane and the vertices ABC of the triangular sub-area, respectively, according to the following formula:

[0143]

[0144] Where x is the x-axis coordinate of the inversion grid gravity center P in space, y is the y-axis coordinate of the inversion grid gravity center P in space, 1 is the x-axis coordinate of the triangle vertex A in space, y 1 is the y-axis coordinate of the triangle vertex A in space, x 2 is the x-axis coordinate of the triangle vertex C in space, y 2 is the y-axis coordinate of the triangle vertex B in space, x 3 is the x-axis coordinate of the triangle vertex C in space, y 3 is the y-axis coordinate of the triangle vertex C in space; and S ΔP ' BC is the directed area of ​​triangle P′BC, S ΔAP′C is the directed area of ​​triangle AP′C, S ΔABP′ is the directed area of ​​triangle ABP′;

[0145] When the condition S is met ΔP′BC +S ΔAP′C +S ΔABP′ =S ΔABC When , the radar detection at the vertex of the triangular sub-area ABC can be selected to obtain the radial velocity data as the three-dimensional wind field inversion data.

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

[0147] It can be understood that when all grids have found their corresponding radar networks, the three-dimensional wind field inversion work can be carried out.

[0148] Furthermore, the inversion grid centroid P is located in the triangular sub-region ABC, and the coordinates of the inversion grid centroid P are (x, y, z). The radar coordinates of the vertices of the triangular sub-region ABC are (x 1 ,y 1 ,z 1 ), (x2 ,y 2 ,z 2 ), (x 3 ,y 3 ,z 3 ), construct the direction vector, the formula is as follows:

[0149]

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

[0151] The relationship between the radial velocity two-dimensional filtering result and the wind speed component is obtained through the relationship of vector projection, as shown in the following formula:

[0152]

[0153] Where V r1 is the radial velocity of the first radar in the triangular sub-area obtained by two-dimensional filtering on the inversion grid centroid P, V r2 is the radial velocity of the second radar in the triangular sub-area obtained by two-dimensional filtering on the inversion grid centroid P, V r3 is the radial velocity obtained by two-dimensional filtering of the third radar in the triangular sub-area on the inversion grid centroid P; 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; is the direction vector composed of the coordinates of the first radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the coordinates of the second radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the coordinates of the third radar in the triangular sub-area and the centroid P of the inversion grid; Represented as a direction vector 2 Fan, Represented as a direction vector 2 Fan, Represented as a direction vector 2Fan;

[0154] It can be understood that the inversion equation can be constructed through the above relationship to solve 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 obtaining the inverse of the matrix A.

[0158] Where V r1 is the radial velocity of the first radar in the triangular sub-area obtained by two-dimensional filtering on the inversion grid centroid P, V r2 is the radial velocity of the second radar in the triangular sub-area obtained by two-dimensional filtering on the inversion grid centroid P, V r3 is the radial velocity of the third radar in the triangular sub-area obtained by two-dimensional filtering on the inversion grid gravity center P; x is the x-axis coordinate of the inversion grid gravity center P in space, y is the y-axis coordinate of the inversion grid gravity center P in space, z is the z-axis coordinate of the inversion grid gravity center P in space, x 1 is the x-axis coordinate of the first radar position in the triangular sub-area in the Cartesian coordinate system, and y 1 is the y-axis coordinate of the first radar position in the triangular sub-area in the Cartesian coordinate system, z 1 is the altitude of the first radar position in the triangular sub-area in the Cartesian coordinate system, x 2 is the x-axis coordinate of the second radar position in the triangular sub-area in the Cartesian coordinate system, and y 2 is the y-axis coordinate of the second radar position in the triangular sub-area in the Cartesian coordinate system, z 2 is the altitude of the second radar position in the triangular sub-area in the Cartesian coordinate system, x 3 is the x-axis coordinate of the third radar position in the triangular sub-area in the Cartesian coordinate system, and y 3 is the y-axis coordinate of the third radar position in the triangular sub-area in the Cartesian coordinate system, z 3 is the altitude of the third radar position in the triangular sub-area in the Cartesian coordinate system, is the direction vector composed of the coordinates of the first radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the coordinates of the second radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the third radar coordinate in the triangular sub-area and the inversion grid centroid P, Represented as a direction vector 2 Fan, Represented as a direction vector 2 Fan, Represented as a direction vector The 2-norm of the formula, u is the horizontal wind speed component in the north-south direction, v is the horizontal wind speed component in the east-west direction, and w is the vertical wind speed component. A is a simplified description of the matrix in the formula.- 1 is the inversion matrix.

[0159] It can be understood that the Euclidean distance is used as the weight when the intensity data is fused between the three radars.

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

[0161]

[0162] Where Z is the fused intensity data, Z 1 is the intensity of the first radar in the triangular sub-area interpolated at the inversion grid centroid P, Z 2 is the intensity of the first radar in the triangular sub-area interpolated at the inversion grid centroid P, Z 3 is the intensity of the first radar in the triangular sub-area interpolated at the center of the inversion grid P; is the direction vector composed of the coordinates of the first radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the coordinates of the second radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the third radar coordinate in the triangular sub-area and the inversion grid centroid P, Represented as a direction vector 2 Fan, Represented as a direction vector 2 Fan, Represented as a direction vector 2Fan.

[0163] It can be understood that the interpolated intensity data is fused using the existing weight calculation method.

[0164] It can be understood that the three-dimensional wind field inversion intensity fusion can be completed by selecting the intensity data obtained by radar detection of the inversion grid in the radar network and combining it with the weight composed of the Euclidean distance.

[0165] It can be understood that in the process of wind field inversion, it is also necessary to subtract the falling velocity of precipitation particles from the inversion result to obtain the true vertical component of the wind vector.

[0166] Furthermore, the true vertical component of the wind vector is given by:

[0167] w ′ =ww t

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

[0169] Where 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 by inversion of the three-dimensional wind field; Z is the fused intensity data; ρ 0 is the air density at the inversion height plane; is the average air density.

[0170] It can be understood that no matter how the network is formed, there will be data at the edge of the divided radar network, such as Figure 3 In the process of the present invention, this situation is reflected as the directed area S ΔP′BC , S ΔAP′C , S ΔABP′ The sum is 0, which means that the inversion grid center P is located at the edge of the two triangular sub-areas. In this case, there are subjective factors in selecting the radars to participate in the inversion, which causes the selection result to be not the optimal solution, and usually leads to the wind vector result obtained by inversion not meeting expectations. Therefore, the present invention proposes a method for jointly inverting smooth wind vectors in adjacent triangular sub-areas.

[0171] Furthermore, the inversion grid on the boundary of the adjacent radar network triangular sub-area is selected, and the radars belonging to the adjacent radar network triangular sub-areas on both sides are used to perform wind field inversion on the inversion grid respectively, and two results of wind speed components belonging to the boundary grid are obtained, which are expressed as u 1 , v 1 , w′ 1 and u 2 , v 2 , w′ 2 , and then synthesize the components in an average manner, the formula is as follows:

[0172]

[0173] In the formula, u 1 is the north-south horizontal component of the inversion grid result of the first triangular sub-area in the adjacent radar network triangular sub-area, v 1 The east-west horizontal component of the inversion grid result of the first triangular sub-region of the adjacent radar network triangular sub-region, w′ 1 u is the vertical square component of the inversion grid result of the first triangular sub-region of the adjacent radar network triangular sub-region; 2 is the north-south horizontal component of the inversion grid result of the second triangular sub-area in the adjacent radar network triangular sub-area, v 2 The east-west horizontal component of the inversion grid result of the second triangular sub-region of the adjacent radar network triangular sub-region, w′ 2is the vertical component of the inversion grid result of the second triangular sub-region of the adjacent radar network triangular sub-region; u is the synthetic 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 synthetic 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 synthetic result of the vertical component of the inversion grid on the boundary of the adjacent radar network triangular sub-region.

[0174] Get the wind speed v′, horizontal direction θ and vertical direction The formula is as follows:

[0175]

[0176] Where u is the synthesis result of the north-south horizontal component of the inversion grid, v is the synthesis result of the east-west horizontal component of the inversion grid, w is the synthesis result of the vertical component of the inversion grid, and w′ is the synthesis result of the vertical component of the inversion grid on the boundary of the triangular sub-area of ​​the adjacent radar network.

[0177] It can be understood that the wind speed v′, horizontal direction θ and vertical direction θ can be directly obtained by obtaining the wind speed components of each inversion grid.

[0178] In a second aspect, the embodiment of the present application provides a three-dimensional wind field inversion system based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing, including:

[0179] Radar site geographic information system, used to obtain radar site geographic information;

[0180] A radar networking system, used for radar networking using a triangulation method based on the geographic information of radar sites;

[0181] Inversion grid division module, used for inversion grid division after radar networking;

[0182] Interpolation module, used for obtaining radial velocity data through radar detection and performing interpolation;

[0183] The radial velocity two-dimensional filtering result module is used to perform two-dimensional filtering on the interpolation result of the radial velocity data obtained by radar detection through the result of averaging the height of the divided inversion grid and the signal-to-noise ratio obtained by radar detection, so as to obtain the radial velocity two-dimensional filtering result;

[0184] The three-dimensional wind field inversion module is used to perform three-dimensional wind field inversion through the inversion grid, radar site geographic information and radial velocity two-dimensional filtering results.

[0185] In a third aspect, the present application provides a computer device comprising a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor executes 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.

[0186] The computer device may be a desktop computer, a notebook, a PDA, a cloud server, etc. The computer device may interact with the user through a keyboard, a mouse, a remote control, a touch pad, or a voice control device.

[0187] The memory includes at least one type of readable storage medium, and the readable storage medium includes flash memory, hard disk, multimedia card, card-type memory (for example, SD or D interface display memory, etc.), random access memory (RAM), static random access memory (SRAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), programmable read-only memory (PROM), magnetic memory, disk, optical disk, etc. In some embodiments, the memory can be an internal storage unit of the computer device, such as a hard disk or memory of the computer device. In other embodiments, the memory can also be an external storage device of the computer device, such as a plug-in hard disk equipped on the computer device, a smart memory card (Smart Media Card, SMC), a secure digital (Secure Digital, SD) card, a flash card (Flash Card), etc. Of course, the memory can also include both the internal storage unit of the computer device and its external storage device. In this embodiment, the memory is often used to store the operating system and various application software installed on the computer device, such as the program code 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 memory may also be used to temporarily store various types of data that have been output or are to be output.

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

[0189] In a fourth aspect, the present application implements a computer-readable storage medium storing a computer program. When the computer program is executed by a processor, the processor executes 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] Wherein, the computer-readable storage medium stores an interface display program, and the interface display program can be executed by at least one processor so that the at least one processor performs 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 mentioned above.

[0191] Application Examples

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

[0193] After interpolating the radial velocity for each radar, a two-dimensional filter of the radial velocity is applied to each altitude level. Figure 5 The frequency responses of the two-dimensional filter in this application example at inversion heights of 1km, 5km, 10km, 15km and 20km correspond to Figure 5 (a), (b), (c), (d) and (e) in the figure. 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 20dB. In this filter, the passband part is sufficient to meet the variation frequency of most real wind speeds.

[0194] When inverting the wind field, the inversion grids at the edge of the radar network sub-area need to be averaged according to the area using the method of jointly inverting the smooth wind vector of adjacent triangular sub-areas of the present invention. Figure 6 As shown. After the wind field inversion is completed in each sub-area, the edge grid of area 1 needs to use the wind field inversion results of the sub-area enclosed by radar 235 and radar 231 for vector averaging; the edge grid of area 2 needs to use the wind field inversion results of the sub-area enclosed by radar 231 and radar 134 for vector averaging; the edge grid of area 3 needs to use the wind field inversion results of the sub-area enclosed by radar 134 and radar 436 for vector averaging; the edge grid of area 4 needs to use the wind field inversion results of the sub-area enclosed by radar 436 and radar 635 for vector averaging; the edge grid of area 5 needs to use the wind field inversion results of the sub-area enclosed by radar 635 and radar 532 for vector averaging.

[0195] The vertical velocity component of the inversion result is compared before and after filtering at an inversion height of 1km, 5km and 10km under the condition of 20dB. Figure 7 , Figure 7 In the figure, (a), (b) and (c) are the results for 1 km; (d), (e) and (f) are the results for 5 km; (g), (h) and (i) are the results for 10 km. Figure 7 It can be seen that the radial velocity two-dimensional filter has a very good error correction effect in the low-altitude wind field inversion results. Although the results at the inversion height of 1km are partially distorted, the overall results after using the radial velocity two-dimensional filter are closer to the ideal results of the prior artificial wind field than those without filtering. At the same time, the role of the radial velocity two-dimensional filter can also be reflected by comparing the RMSE of each altitude layer, such as Figure 8 As shown, the figure compares the effects of the uvw component before and after radial velocity two-dimensional filtering at various heights when the signal-to-noise ratio of the detection data is 20dB and 10dB respectively. Figure 8 (a) and (b) are the images with a detection signal-to-noise ratio of 20 dB; (c) and (d) are the images with a detection signal-to-noise ratio of 10 dB. In addition, the optimized overall three-dimensional wind field inversion process was applied to the actual detection data. The final inversion result of the overall three-dimensional wind field is shown in Figure 1. Fig. 9 As shown, the inversion range of the overall wind field is similar to a pentagon, which includes the main areas of Guangzhou City, and the results of inversion heights of 1km, 5km, 9km, 13km and 17km are selected for display. Fig. 9 The improvement of the overall wind field inversion results after using radial velocity two-dimensional filtering is demonstrated in the figure. It can be clearly seen from the low-altitude inversion results at altitudes of 1km and 5km that after the filter optimization, the suddenly disappearing vertical component in the low-altitude wind field results is supplemented, and the discontinuous vertical component is made more continuous. The numerous burrs in the downdraft structure in the wind field are also smoothed out while maintaining the original downdraft change trend after filtering optimization. At the same time, the wind vector direction that is obviously wrong in the result is also greatly improved after optimization. Fig. 9 It can be seen that although the design of the radial velocity two-dimensional filter requires that the low-altitude wind field obtained by inversion has a more stringent passband than the high-altitude wind field, it can be seen through actual experimental comparison that the continuity trend of the original wind field is still retained in the optimized result, which shows that the method of the present invention is more reliable.

[0196] The three-dimensional wind field inversion results at each height at four moments were selected for verification with the verification radar detection data, such as Fig.10As shown in the figure, it also shows the effect before and after using the radial velocity two-dimensional filter to correct the error. After the actual radar acquisition, the inverted three-dimensional wind field before and after the radial velocity two-dimensional filter was verified and compared at four time points, which can be described as time points 1 to 4, (a) is time point 1, (b) is time point 2, (c) is time point 3, and (d) is time point 4. It can be seen from Fig.10 It can be seen that the error of the inversion result decreases gradually and then increases with the increase of the inversion height. This is because the signal-to-noise ratio of the radar detection data is lower as the inversion height increases, and the detection points are more sparsely distributed, resulting in an increase in errors. Therefore, when comparing the effects of filtering, the analysis of the chart can be focused on the verification results with an inversion height of less than 15km. After two-dimensional filtering, the errors at each height in the results of these four moments have been improved to varying degrees, especially for the low-level wind field with an inversion height of less than 5km.

[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 in the protection scope of the present invention.

Claims

1. A three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing, characterized in that: include: Obtain radar site geographic information; Radar networking is carried out using triangulation method based on the geographic information of radar sites; After radar networking, inversion grid division is performed; Radar detection obtains radial velocity data and performs interpolation; Perform two-dimensional filtering on the interpolation result of radial velocity data obtained by radar detection by averaging the height of the divided inversion grid and the signal-to-noise ratio obtained by radar detection to obtain a two-dimensional filtering result of radial velocity; The three-dimensional wind field inversion is performed by combining the inversion grid with the geographic information of the radar sites and the two-dimensional filtering results of the radial velocity.

2. The three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing according to claim 1 is characterized in that: The radar site geographical information includes longitude, latitude and altitude; and / or Radar networking obtains a radar network. In the radar network, the geographic information data of each radar site is converted from the polar coordinate system to the Cartesian coordinate system. The conversion formula is as follows: x=lcosεsinθ+x0 y=lcosεcosθ+y0 z=lsinε+z0 Wherein, 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 in the Cartesian coordinate system, y0 is the y-axis coordinate of the radar site in the Cartesian coordinate system, z0 is the altitude of the radar site 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 inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing according to claim 2 is characterized in that: After converting all radar site geolocation data from polar to Cartesian coordinates, the inversion grid is generated from the radar network; and / or Using a triangulation method, the radial velocity data obtained by radar detection is triangulated according to relative Cartesian coordinates to obtain a tetrahedron whose vertices are formed by the coordinates of multiple points of the radial velocity data obtained by radar detection and are closely connected, and the radial velocity data obtained by radar detection to be used is determined by judging the distribution of the centroid P of the inversion grid in the tetrahedron; and / or According to the average result of the signal-to-noise ratio obtained by the detection of each radar in the triangular sub-area of ​​the radar network and the radial velocity data that has been interpolated in the three-dimensional space, corresponding frequency domain two-dimensional filters are set at different height layers of the inversion grid. After all the settings are completed, a filter group is obtained.

4. The three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing according to claim 3 is characterized in that: Methods for determining whether the inversion grid centroid P is distributed in the tetrahedron include: Calculate the volume of tetrahedron ABCD using the following formula: Wherein, 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; V ΔABCD Then it is the directed volume of tetrahedron ABCD; The volume enclosed by the centroid P of the inversion grid and each vertex of the tetrahedron ABCD is obtained respectively, and the formula is as follows: Wherein, x is the x-axis coordinate of the inversion grid centroid P in space, y is the y-axis coordinate of the inversion grid centroid P in space, z is the z-axis coordinate of the inversion grid centroid 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 tetrahedral APCD, V ΔABPD is the directed volume of the tetrahedral ABPD, V ΔABCP is the directed volume of the tetrahedron ABCP; Determine whether the inversion grid centroid P is within the tetrahedron ABCD. If the condition V is met, ΔABCD =V ΔPBCD +V ΔAPCD +V ΔABPD +V ΔABCP When , it means that the tetrahedron ABCD contains the inversion grid gravity center P; if this condition is not met, it means that there is no inclusion relationship between the tetrahedron ABCD and the inversion grid gravity center P, and it is necessary to find another tetrahedron that meets the conditions for the inversion grid gravity center P; According to the volume relationship between the inversion grid center P and the sub-tetrahedron and tetrahedron ABCD enclosed by the vertices, the weights to be used in the interpolation process are obtained. Then the radial velocity V of the grid is obtained by interpolation. Pr for: Where V Pr is the interpolated radial velocity, V ΔABCD is the directed volume of tetrahedron ABCD, V ΔPBCD is the directed volume of the tetrahedron PBCD, V ΔAPCD is the directed volume of the tetrahedral APCD, V ΔABPD is the directed volume of the tetrahedral ABPD, V ΔABCP is the directed volume of the tetrahedron ABCP; V Ar , V Br , V Cr , V Dr are respectively the radial velocities detected by radar at the vertices of tetrahedron ABCD; and / or The formula for the frequency response of a two-dimensional filter is as follows: f(h)=k*h+SNR*s In the formula, |H h (ω1,ω2)| is the frequency response of the two-dimensional filter, and D(ω1,ω2) is the distance from each point of the radial velocity of each radar to the 0 frequency point in the frequency domain, expressed as ω1 is the abscissa of the two-dimensional frequency axis, and ω2 is the ordinate of the two-dimensional frequency axis; f(h) is a function related to the height of the inversion grid, and h is the wind field inversion height; k is the height influence factor, ranging from 0 to 1; SNR is the average signal-to-noise ratio obtained by radar detection at the vertex of the triangular sub-area during the inversion process; s is the basic passband influence factor, ranging from 0 to 1.

5. The three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing according to claim 4 is characterized in that: The interpolated radial velocity is subjected to a two-dimensional fast Fourier transform and converted to the frequency domain. According to the inversion grid height, the frequency domain result of the interpolated radial velocity is multiplied by the two-dimensional filter frequency response in the filter group corresponding to the inversion grid height to obtain the filtering result in the frequency domain. Finally, the filtering result is subjected to a two-dimensional inverse Fourier transform to obtain the two-dimensional filtering result of the radial velocity.

6. The three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing according to claim 1 is characterized in that: The selection of three-dimensional wind field inversion data includes: calculating the directed area of ​​the triangular sub-area in the radar network, and then calculating the directed area of ​​the triangle enclosed by the projection P′(x, y) of the inversion grid centroid P on the xy axis plane and the vertices ABC of the triangular sub-area, respectively. The formula is as follows: Where x is the x-axis coordinate of the inversion grid centroid P in space, y is the y-axis coordinate of the inversion grid centroid P in space, x1 is the x-axis coordinate of the triangle vertex A in space, y1 is the y-axis coordinate of the triangle vertex A in space, x2 is the x-axis coordinate of the triangle vertex C in space, y2 is the y-axis coordinate of the triangle vertex B in space, x3 is the x-axis coordinate of the triangle vertex C in space, y3 is the y-axis coordinate of the triangle vertex C in space; and S ΔP′BC is the directed area of ​​triangle P′BC, S ΔAP′C is the directed area of ​​triangle AP′C, S ΔABP′ is the directed area of ​​triangle ABP′; When the condition S is met ΔP′BC +S ΔAP′C +S ΔABP′ =S ΔABC When , the radar detection at the vertex of the triangular sub-area ABC can be selected to obtain the radial velocity data as the three-dimensional wind field inversion data.

7. The three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing according to claim 6 is characterized in that: The inversion grid centroid P is located in the triangular sub-area ABC. The coordinates of the inversion grid centroid P are (x, y, z). The radar coordinates of the vertices of the triangular sub-area ABC are (x1, y1, z1), (x2, y2, z2), and (x3, y3, z3). The direction vector is constructed as follows: In the formula, is the direction vector composed of the coordinates of the first radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the coordinates of the second radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the coordinates of the third radar in the triangular sub-area and the centroid P of the inversion grid; The relationship between the radial velocity two-dimensional filtering result and the wind speed component is obtained through the relationship of vector projection, as shown in the following formula: Where V r1 is the radial velocity of the first radar in the triangular sub-area obtained by two-dimensional filtering on the inversion grid centroid P, V r2 is the radial velocity of the second radar in the triangular sub-area obtained by two-dimensional filtering on the inversion grid centroid P, V t3 is the radial velocity obtained by two-dimensional filtering of the third radar in the triangular sub-area on the inversion grid centroid P; 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; is the direction vector composed of the coordinates of the first radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the coordinates of the second radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the coordinates of the third radar in the triangular sub-area and the centroid P of the inversion grid; Represented as a direction vector 2 Fan, Represented as a direction vector 2 Fan, Represented as a direction vector 2Fan; The matrix of the unknown velocity component uvw inversion equation calculation process is as follows: Where V r1 is the radial velocity of the first radar in the triangular sub-area obtained by two-dimensional filtering on the inversion grid centroid P, V r2 is the radial velocity of the second radar in the triangular sub-area obtained by two-dimensional filtering on the inversion grid centroid P, V r3 is the radial velocity of the third radar in the triangular sub-area obtained by two-dimensional filtering on the inversion grid gravity center P; x is the x-axis coordinate of the inversion grid gravity center P in space, y is the y-axis coordinate of the inversion grid gravity center P in space, z is the z-axis coordinate of the inversion grid gravity center P in space, x1 is the x-axis coordinate of the first radar position in the triangular sub-area in the Cartesian coordinate system, y1 is the y-axis coordinate of the first radar position in the triangular sub-area in the Cartesian coordinate system, z1 is the altitude of the first radar position in the triangular sub-area in the Cartesian coordinate system, x2 is the x-axis coordinate of the second radar position in the triangular sub-area in the Cartesian coordinate system, y2 is the y-axis coordinate of the second radar position in the triangular sub-area in the Cartesian coordinate system, z2 is the altitude of the second radar position in the triangular sub-area in the Cartesian coordinate system, x3 is the x-axis coordinate of the third radar position in the triangular sub-area in the Cartesian coordinate system, y3 is the y-axis coordinate of the third radar position in the triangular sub-area in the Cartesian coordinate system, and z3 is the altitude of the third radar position in the triangular sub-area in the Cartesian coordinate system. is the direction vector composed of the coordinates of the first radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the coordinates of the second radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the third radar coordinate in the triangular sub-area and the inversion grid centroid P, Represented as a direction vector 2 Fan, Represented as a direction vector 2 Fan, Represented as a direction vector The 2-norm of the formula, u is the horizontal wind speed component in the north-south direction, v is the horizontal wind speed component in the east-west direction, and w is the vertical wind speed component. A is a simplified description of the matrix in the formula. -1 is the inversion matrix.

8. The three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing according to claim 7 is characterized in that: The radar detects the intensity data, then interpolates and fuses the interpolated intensity data. The formula is as follows: Wherein, Z is the fused intensity data, Z1 is the intensity of the first radar in the triangular sub-area interpolated on the inversion grid gravity center P, Z2 is the intensity of the first radar in the triangular sub-area interpolated on the inversion grid gravity center P, and Z3 is the intensity of the first radar in the triangular sub-area interpolated on the inversion grid gravity center P; is the direction vector composed of the coordinates of the first radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the coordinates of the second radar in the triangular sub-area and the centroid P of the inversion grid, is the direction vector composed of the third radar coordinate in the triangular sub-area and the inversion grid centroid P, Represented as a direction vector 2 Fan, Represented as a direction vector 2 Fan, Represented as a direction vector 2Fan.

9. The three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing according to claim 8 is characterized in that: The true vertical component of the wind vector is given by: w′=ww t w t =2.65(ρ0 / ρ)Z 0.114 Where w′ is the true vertical component of the wind vector; w t is the falling speed of precipitation particles; w is the vertical velocity component obtained by three-dimensional wind field inversion; Z is the fused intensity data; ρ0 is the air density at the inversion height plane; ρ is the average air density.

10. The three-dimensional wind field inversion method based on radial velocity two-dimensional filtering and radar network edge wind vector smoothing according to claim 1 is characterized in that: The inversion grid on the boundary of the adjacent radar network triangular sub-area is selected, and the radars belonging to the adjacent radar network triangular sub-areas on both sides are used to perform wind field inversion on the inversion grid respectively, and two results of wind speed components belonging to the boundary grid are obtained, which are expressed as u1, v1, w′1 and u2, v2, w′2. The components are then synthesized in an average manner, and the formula is as follows: Wherein, 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 plane 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 plane component of the inversion grid result of the second triangular sub-region in the adjacent radar network triangular sub-region; u is the synthesis 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 synthesis result of the east-west horizontal component of the inversion grid on the boundary of the adjacent radar network triangular sub-region, and w′ is the synthesis result of the vertical component of the inversion grid on the boundary of the adjacent radar network triangular sub-region; Get the wind speed v′, horizontal direction θ and vertical direction The formula is as follows: Where u is the synthesis result of the north-south horizontal component of the inversion grid, v is the synthesis result of the east-west horizontal component of the inversion grid, w is the synthesis result of the vertical component of the inversion grid, and w′ is the synthesis result of the vertical component of the inversion grid on the boundary of the triangular sub-area of ​​the adjacent radar network.

Images