Dynamic resolution 3D structure reconstruction method based on weather radar data combined with terrain

Abstract

The invention discloses a method for reconstructing a three-dimensional structure with dynamic resolution of weather radar base data combined with terrain. The method realizes three-dimensional structure reconstruction by combining a digital elevation model (DEM) with weather radar base data. The method comprises the following steps: cutting out DEM terrain data and dumping it into grid data in binary form, generating a bitmap as a mapping source, establishing an elevation index array based on weather radar base data and an azimuth index array based on a conical scanning surface to quickly map angles, determining a rectangular reconstruction area and establishing a three-dimensional grid, traversing the three-dimensional grid to obtain elevation and azimuth angles of grid points and a distance combination array to obtain grid point data, smoothing the grid point data, resampling terrain grid data and constructing a terrain grid surface, obtaining a mapping image of the terrain grid surface, and performing three-dimensional volume rendering. The method has the advantages of not relying on online services, being able to realize dynamic high-precision three-dimensional structure reconstruction of radar echoes in a designated area while ensuring interactive performance.

Description

Technical Field

[0001] The invention relates to a three-dimensional structure reconstruction technology, and in particular to a three-dimensional structure reconstruction method of weather radar base data combined with terrain and with dynamic resolution. Background Art

[0002] Weather radar base data is a series of spatial slice data obtained by performing conical scanning along multiple elevation angles with the radar station as the center. These slice data record the radar echo values ​​of various observation elements, such as reflectivity, velocity, Doppler shift, etc., which together constitute the three-dimensional spatial data set of radar observation. In order to more vividly reflect the cloud water distribution characteristics in the actual atmospheric environment, these slice data need to be converted into continuous three-dimensional grid data that can be used for three-dimensional volume rendering through interpolation technology. This process is called three-dimensional structure reconstruction of weather radar base data.

[0003] When reconstructing a 3D structure, you first need to establish a grid of a target coordinate system, which can be a Cartesian coordinate system or a spherical coordinate system based on longitude and latitude. Then, in the target coordinate system, each grid point is interpolated using the surrounding slice data to estimate the value of the grid point. After the interpolation is completed, the obtained 3D grid data will be sent to the display module, and visualized through 3D volume rendering technology to obtain a 3D image.

[0004] In actual interpolation calculations, for a fixed range, if the grid resolution is too low, it will lack details and cannot accurately reflect the cloud water distribution characteristics in the atmospheric environment; if the grid resolution is too high, the interpolation speed will be slow, and the generated data volume will become huge, which may cause difficulties in processing and storage. In order to solve these problems, some strategies are usually adopted in meteorological business work. For example, data can be pre-calculated or generated, or an isosurface tracking algorithm can be used to calculate a three-dimensional object composed of isosurfaces that have lost details. In addition, fixed-resolution interpolation can be used to create a high-resolution three-dimensional matrix, and then the user-selected area can be intercepted for display during interaction. These methods can balance the needs of data details and processing speed to a certain extent.

[0005] At the same time, in meteorological data analysis, terrain, as an important factor, has a significant impact on the determination of meteorological situations. Adding terrain data (generally provided in the form of a digital elevation model (DEM)) to the visualization process can enhance the understanding and prediction of meteorological phenomena. DEM data is provided in the form of grid points, and is usually triangulated to form a terrain grid surface, thereby reflecting the undulating changes of terrain in three-dimensional space. In terrain analysis of a fixed area, DEM data can be pre-calculated to construct a terrain grid surface. However, in interactive systems, such as real-time meteorological data analysis, real-time calculations are required for unspecified areas. Traditional triangulation methods may be inefficient in real-time calculations and are difficult to meet the needs of real-time interactive systems.

[0006] Currently, there are few weather radar application software that can provide real-time interactive analysis and combine terrain and three-dimensional reconstruction. Most applications use Web systems and are implemented by calling external three-dimensional geographic information services, which usually require online connections and may have limitations in flexibility. Summary of the invention

[0007] The technical problem to be solved by the present invention is to provide a method for dynamic resolution three-dimensional structure reconstruction of weather radar base data combined with terrain, which does not need to rely on online services and can achieve dynamic high-precision three-dimensional structure reconstruction of radar echoes in a specified area while ensuring interactive performance.

[0008] The technical solution adopted by the present invention to solve the above technical problems is: a method for reconstructing a three-dimensional structure with dynamic resolution of weather radar base data combined with terrain, characterized by comprising the following steps:

[0009] Step 1: Select a data source in the form of DEM text, and clip the DEM terrain data. The clipping range is limited to the coverage of the radar of interest. Then, the clipped DEM terrain data is dumped into binary grid data, and the number of grid points, resolution, and latitude and longitude of the starting and ending grid points of the binary grid data are obtained. Then, the binary grid data is used, and the mapping relationship between elevation and color is combined to generate the corresponding bitmap image.

[0010] Step 2: For the weather radar base data to be processed, it contains N conical scanning surfaces with different elevation angles, each conical scanning surface contains M radial scans, wherein each conical scanning surface records the elevation angle, the length and number of the range library on the conical scanning surface, and each radial scan records its azimuth in the conical scanning surface to which it belongs. A serial number is provided for each conical scanning surface in the ascending order of the elevation angles of all conical scanning surfaces, and the serial number is 0 to N-1. A serial number is provided for each radial scan contained in any conical scanning surface in the ascending order of the azimuth angles of these radial scans, and the serial number is 0 to M-1;

[0011] Step 3: Based on the weather radar base data, a one-dimensional elevation index array A of size 1800 is established. Array A quickly maps the serial number of the nearest conical scanning surface to any of the 1800 angle values ​​from -90 degrees to 90 degrees with a resolution increment of 0.1 degrees; based on the conical scanning surface with serial number n, a one-dimensional azimuth index array B of size 3600 is established. n , array B n Rapidly map the number of the radial scan closest to any of the 3600 angle values ​​from 0 to 360 degrees with a resolution of 0.1 degree, where n = 0, 1, ..., N-1;

[0012] Step 4: Within the coverage of the radar of interest on the electronic map, a rectangular reconstruction area is determined by drawing a frame, and the coordinate positions of the upper left corner and the lower right corner of the rectangular reconstruction area are recorded; then a three-dimensional grid is established on the rectangular reconstruction area, and the parameters of the three-dimensional grid are determined, including the number of grid points in the horizontal longitude and latitude directions, the longitude and latitude of the starting grid point of the horizontal plane, the longitude and latitude of the ending grid point of the horizontal plane, the total height value, and the number of height layers;

[0013] Step 5: Traverse all grid points in the three-dimensional grid established in step 4. For the currently traversed grid point, calculate the elevation angle φ, azimuth angle γ and distance D of the grid point relative to the radar position based on the latitude, longitude and altitude values ​​of the grid point and the latitude, longitude and altitude values ​​of the radar of interest, where the accuracy of the elevation angle and azimuth angle is taken to 0.1 degree; then, based on the elevation angle φ, azimuth angle γ and distance D, and combined with A and B n , n=0,1,…,N-1, obtain the data of the distance library closest to the grid point on the four radial scans closest to the grid point, and then perform bilinear interpolation or inverse square interpolation to obtain the data of the grid point;

[0014] Step 6: Based on step 5, a smoothing range is preset. For the data of each grid point in the three-dimensional grid, if the data of the grid point is within the smoothing range, the data of the grid point is subjected to median filtering to obtain a three-dimensional grid with smooth edges; if the data of the grid point is not within the smoothing range, the data of the grid point is not processed;

[0015] Step 7: Set the number of grids of the terrain grid data, and set the area formed by the latitude and longitude of the horizontal plane starting grid point and the horizontal plane ending grid point of the three-dimensional grid as the latitude and longitude area of ​​the terrain grid data, and calculate the latitude and longitude resolution of the terrain grid data; then, according to the latitude and longitude area and the latitude and longitude resolution of the terrain grid data, resample the terrain grid data of the latitude and longitude area from the binary grid data; then use the terrain grid data to construct triangles according to the diagonal principle, and then form a terrain grid surface;

[0016] Step 8: According to the longitude and latitude of the starting grid point of the horizontal plane of the three-dimensional grid and the longitude and latitude of the ending grid point of the horizontal plane, a region consisting of the longitude and latitude of the starting grid point of the horizontal plane of the three-dimensional grid and the longitude and latitude of the ending grid point of the horizontal plane is cut out from the bitmap image as a mapping image of the terrain grid surface;

[0017] Step 9: Send the three-dimensional grid obtained in step 6, the terrain grid surface obtained in step 7, and the texture image of the terrain grid surface obtained in step 8 to the display engine to realize three-dimensional volume rendering and terrain construction to obtain a three-dimensional image.

[0018] In step 3, the process of establishing the one-dimensional elevation index array A is as follows:

[0019] Step 3.1a: The subscripts of all elements in array A are 0 to 1799, which correspond to 1800 angle values ​​from -90 degrees to 90 degrees with a resolution of 0.1 degree. The values ​​of all elements in array A are initialized to -1. When the value of an element in array A is -1, it means that the angle value corresponding to the subscript of the element has no corresponding conical scanning surface.

[0020] Step 3.2a: Traverse all conical scanning surfaces contained in the weather radar base data; set the serial number of the currently traversed conical scanning surface to n, record the elevation angle of the currently traversed conical scanning surface as α, and record the elevation angle of the conical scanning surface with serial number n+1 as β, where n=0,1,…,N-2; then calculate the subscript i corresponding to each angle θ with a resolution of 0.1 degrees in the interval [α,α+(β-α) / 2), i=int(θ×10)+90 0, and then assign the element A[i] with subscript i in the array A to the serial number of the currently traversed conical scanning surface, that is, n, where int() is a rounding function; and calculate the subscript j corresponding to each angle η with a resolution of 0.1 degrees in the interval [α+(β-α) / 2,β], j=int(η×10)+900; and then assign the element A[j] with subscript j in the array A to the serial number of the next conical scanning surface of the currently traversed conical scanning surface, that is, n+1.

[0021] In step 3, the one-dimensional azimuth index array B n The establishment process is:

[0022] Step 3.1b: Array B n The subscripts of all elements in are 0 to 3599, corresponding to 3600 angle values ​​from 0 to 360 degrees with a resolution of 0.1 degree. n All elements in array B are initialized to -1. n When the value of an element in is -1, it means that the angle value corresponding to the subscript of the element has no corresponding radial scan;

[0023] Step 3.2b: traverse all radial scans contained in the conical scanning surface with serial number n; set the serial number of the currently traversed radial scan to m, record the azimuth angle of the currently traversed radial scan as a, and record the azimuth angle of the radial scan with serial number m+1 as b, where m=0, 1, ..., M-2;

[0024] When abs(ba) is greater than the missing threshold, keep array B n The value of the element with subscript int(a×10) to the element with subscript int(b×10) remains unchanged at -1, indicating missing measurement, where abs() is the absolute value function;

[0025] When abs(ba) is less than or equal to the missing threshold, calculate the subscript g corresponding to each angle δ in the interval [a, a+(ba) / 2) with a resolution of 0.1 degrees, g = int(δ×10), and then convert the array B n The element B with subscript g n [g] is assigned the serial number of the currently traversed radial scan, i.e. m; and the subscript h corresponding to each angle σ with a resolution of 0.1 degrees in the interval [a+(ba) / 2,b] is calculated, where h=int(σ×10), and then the array B is n The element B with subscript h n [h] is assigned the serial number of the next radial scan of the currently traversed radial scan, i.e. m+1;

[0026] Step 3.3b: Based on the azimuth wraparound characteristic, the azimuth of the radial scan numbered M-1 is recorded as c, and the azimuth of the radial scan numbered 0 is recorded as d; then calculate the subscript e corresponding to each angle ρ with a resolution of 0.1 degrees in the interval [c, c+(d+360-c) / 2). If int(ρ×10) is greater than or equal to 3600, then e=int(ρ×10)-3600; if int(ρ×10) is less than 3600, then e=int(ρ×10), and then replace array B with n The element B with subscript e n[e] is assigned a value of M-1; and calculate the subscript f corresponding to each angle ω in the interval [c+(d+360-c) / 2, d+360] with a resolution of 0.1 degrees. If int(ω×10) is greater than or equal to 3600, then f=int(ω×10)-3600; if int(ω×10) is less than 3600, then f=int(ω×10). Then replace array B n The element B with subscript f n [f] is assigned a value of 0.

[0027] In step 5, the following process is performed for the currently traversed grid point:

[0028] Step 5.1: Based on the latitude, longitude and altitude values ​​of the grid point, as well as the latitude, longitude and altitude values ​​of the radar concerned, calculate the elevation angle φ, azimuth angle γ and distance D of the grid point relative to the radar position, where the accuracy of the elevation angle and azimuth angle is taken to 0.1 degree;

[0029] Step 5.2: According to the elevation angle φ, the azimuth angle γ and the distance D, obtain the data V1, V2, V3, V4 of the distance library closest to the grid point on the four radial scans closest to the grid point. The specific acquisition process is as follows:

[0030] Step 5.2.1: Calculate the subscript p corresponding to the elevation angle φ in array A, p = int(φ×10)+900, and then define the conical scanning surface numbered A[p] as the first conical scanning surface;

[0031] Step 5.2.2: When the elevation angle φ is lower than or equal to the elevation angle of the first conical scanning surface, if A[p]-1 is greater than or equal to 0, the conical scanning surface numbered A[p]-1 is defined as the second conical scanning surface, and then step 5.2.3 is executed; if A[p]-1 is less than 0, it means that the grid point is in the blind area, and the data of the grid point is set to an invalid value;

[0032] When the elevation angle φ is higher than the elevation angle of the first conical scanning surface, if A[p]+1 is less than or equal to N-1, the conical scanning surface numbered A[p]+1 is defined as the second conical scanning surface, and step 5.2.3 is executed again; if A[p]+1 is greater than N-1, it means that the grid point is in the blind area, and the data of the grid point is set to an invalid value;

[0033] Step 5.2.3: Calculate the one-dimensional azimuth index array B based on the first conical scanning surface A[p] The subscript q corresponding to the azimuth γ is q = int(γ×10). If B A[p] [q] is -1, which means that the radial scan in the first conical scan surface is no longer searched, V1 and V2 are set to invalid values, and step 5.2.6 is executed again; if B A[p] [q] is not -1, then the sequence number is B A[p]The radial scan of [q] is defined as the first radial scan, and then step 5.2.4 is performed;

[0034] Step 5.2.4: When the azimuth angle γ is less than or equal to the azimuth angle of the first radial scan, if B A[p] [q]-1 is greater than or equal to 0, then the sequence number is B A[p] The radial scan of [q]-1 is defined as the second radial scan, and then step 5.2.5 is executed; if B A[p] If [q]-1 is less than 0, the radial scan numbered M-1 is defined as the second radial scan, and step 5.2.5 is executed again;

[0035] When the azimuth angle γ is greater than the azimuth angle of the first radial scan, if B A[p] [q]+1 is less than or equal to M-1, then the sequence number is B A[p] The radial scan of [q]+1 is defined as the second radial scan, and then step 5.2.5 is executed; if B A[p] If [q]+1 is greater than M-1, the radial scan with sequence number 0 is defined as the second radial scan, and step 5.2.5 is executed again;

[0036] Step 5.2.5: Determine whether the absolute value of the difference between the azimuth angle of the first radial scan and the azimuth angle of the second radial scan is greater than the missing measurement threshold. If so, set V1 and V2 to invalid values, and then execute step 5.2.6; otherwise, calculate the index t of the distance library, t is equal to the distance D divided by the length of the distance library on the first conical scanning surface. When t is greater than the number of distance libraries on the first conical scanning surface, set V1 and V2 to invalid values, and then execute step 5.2.6; when t is less than or equal to the number of distance libraries on the first conical scanning surface, V1 is the data on the distance library with index t on the first radial scan, and V2 is the data on the distance library with index t on the second radial scan, and then execute step 5.2.6;

[0037] Step 5.2.6: Let s represent the serial number of the second conical scanning surface, and calculate the one-dimensional azimuth index array B established based on the second conical scanning surface s The subscript u corresponding to the azimuth γ is u=int(γ×10). If B s [u] is -1, which means that the radial scan in the second conical scan surface is no longer searched, V3 and V4 are set to invalid values, and step 5.3 is executed again; if B s If [u] is not -1, the sequence number is B s The radial scan of [u] is defined as the third radial scan, and step 5.2.7 is performed again;

[0038] Step 5.2.7: When the azimuth angle γ is less than or equal to the azimuth angle of the third radial scan, if B s [u]-1 is greater than or equal to 0, then the sequence number is Bs The radial scan of [u]-1 is defined as the fourth radial scan, and then execute step 5.2.8; if B s If [u]-1 is less than 0, the radial scan numbered M-1 is defined as the fourth radial scan, and step 5.2.8 is executed again;

[0039] When the azimuth angle γ is greater than the azimuth angle of the third radial scan, if B s [u]+1 is less than or equal to M-1, then the sequence number is B s The radial scan of [u]+1 is defined as the fourth radial scan, and then step 5.2.8 is executed; if B s If [u]+1 is greater than M-1, the radial scan with sequence number 0 is defined as the fourth radial scan, and step 5.2.8 is executed again;

[0040] Step 5.2.8: Determine whether the absolute value of the difference between the azimuth angle of the third radial scan and the azimuth angle of the fourth radial scan is greater than the missing measurement threshold. If so, set V3 and V4 to invalid values, and then execute step 5.3; otherwise, calculate the index z of the distance library, z is equal to the distance D divided by the length of the distance library on the second conical scanning surface. When z is greater than the number of distance libraries on the second conical scanning surface, set V3 and V4 to invalid values, and then execute step 5.3; when z is less than or equal to the number of distance libraries on the second conical scanning surface, V3 is the data on the distance library with index z on the third radial scan, and V4 is the data on the distance library with index z on the fourth radial scan, and then execute step 5.3;

[0041] Step 5.3: If V1, V2, V3, and V4 are all valid values, use the four valid values ​​to perform bilinear interpolation to obtain the data of the grid point; if V1, V2, V3, and V4 contain both invalid values ​​and valid values, use the valid values ​​to perform inverse distance square interpolation to obtain the data of the grid point; if V1, V2, V3, and V4 are all invalid values, then the grid point is set to an invalid value.

[0042] In step 7, the triangle construction process is as follows: the terrain grid data is an equally spaced two-dimensional grid, and the length and width of the terrain grid data are Width and Height respectively; then let (x, y) represent the two-dimensional coordinates of a grid point in the terrain grid data, the value range of x is [0, Width-2], and the value range of y is [0, Height-2], and for the grid point with two-dimensional coordinates (x, y), three grid points with two-dimensional coordinates (x, y+1), (x+1, y+1), and (x+1, y) are selected; then, according to the diagonal principle, a triangle with three grid points with two-dimensional coordinates (x, y), (x, y+1), and (x+1, y+1) as vertices and a triangle with three grid points with two-dimensional coordinates (x+1, y+1), (x+1, y), and (x, y) as vertices are constructed; and then all the triangles form a terrain grid surface.

[0043] The value range of the missing detection threshold is [3,5] degrees.

[0044] In step 4, the number of horizontal grid points is W, and the W value is dynamically adjusted according to the hardware configuration of the computer, or is a preset fixed value; the longitude and latitude of the starting grid point of the horizontal plane are calculated according to the coordinate position of the upper left corner of the rectangular reconstruction area, and the longitude and latitude of the ending grid point of the horizontal plane are calculated according to the coordinate position of the lower right corner of the rectangular reconstruction area; the total height value is H, and for S-band and C-band radars, the value range of H is 20 to 25 km, and for X-band radars, the value of H is 10 km to 15 km; the number of height layers is L, and the L value is dynamically adjusted according to the hardware configuration of the computer, or is a preset fixed value.

[0045] Compared with the prior art, the advantages of the present invention are:

[0046] 1) Users can select a rectangular reconstruction area to reconstruct the 3D structure, and dynamically adjust the resolution. When analyzing a large area, the resolution is automatically reduced to increase the interpolation speed and reduce the data volume, while when analyzing a small area, the resolution is increased to enhance the detail performance, thereby effectively avoiding the problem of missing details or processing and storage difficulties caused by inappropriate resolution. That is, while ensuring interactive performance, the ability to reconstruct large-scale overall and local high-precision 3D structures is taken into account, allowing users to directly perform interactive 3D structure reconstruction analysis of any area based on weather radar base data.

[0047] 2) A one-dimensional elevation index array is established based on the weather radar base data, and a one-dimensional azimuth index array is established based on each conical scanning surface. They are used for fast interpolation indexing, which reduces the amount of search calculations for specified elevation and azimuth data in interpolation. Compared with the traditional method of aligning data to a fixed angle, it not only retains the integrity of the original data and makes full use of data with higher angular resolution, but also can quickly find data in the original data, thereby improving the efficiency of real-time data processing. Especially in interactive meteorological data analysis, it can quickly respond to user needs, perform real-time calculations and analysis on unspecified areas, and meet the high efficiency requirements of real-time interactive systems.

[0048] 3) A fast terrain mesh construction method is adopted. Through a simplified four-corner mesh splitting strategy, each four-grid mesh surrounded by four grid points is split into two triangles according to a fixed diagonal to form a terrain network surface. Compared with the traditional triangulation method to construct an accurate triangulated mesh surface, while maintaining the basic visual effect, the amount of calculation is reduced and the construction speed of the terrain mesh surface is accelerated, so that DEM data can be efficiently used for terrain analysis in interactive systems, which improves the practicality and flexibility of real-time interactive analysis.

[0049] 4) By combining DEM terrain data and weather radar base data and using a variety of methods to accelerate data processing, users can achieve three-dimensional structure reconstruction and terrain grid surface construction of radar echoes in any area at an acceptable performance, so that the display engine can vividly display the distribution of precipitation echoes in three dimensions. Compared with other software in the industry that mainly uses two-dimensional geographic information, it can better reflect the impact of terrain on precipitation.

[0050] 5) As an independent analysis tool, the method of the present invention does not need to rely on online connection, providing higher flexibility and autonomy. Users can freely conduct interactive analysis of meteorological data in the local environment, and can perform effective data analysis and decision support whether there is a network connection or not. At the same time, data security and privacy protection are enhanced because it allows sensitive data to be processed locally, reducing the risk of data exposure during network transmission. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] Figure 1 It is a general implementation block diagram of the method of the present invention;

[0052] Figure 2 A schematic diagram of the index between two elevation angles;

[0053] Figure 3 Schematic diagram of the index between two azimuths. DETAILED DESCRIPTION

[0054] The present invention is further described in detail below with reference to the accompanying drawings.

[0055] The present invention proposes a method for reconstructing three-dimensional structure with dynamic resolution of weather radar base data combined with terrain, and its overall implementation block diagram is as follows: Figure 1 As shown, it includes the following steps:

[0056] Step 1: Select a DEM (Digital Elevation Model) data source in text form and clip the DEM terrain data. The clipping range is limited to the coverage of the radar of interest. Then, the clipped DEM terrain data is dumped into binary grid data, and the number of grid points, resolution, and latitude and longitude of the starting and ending grid points of the binary grid data are obtained for subsequent processing. Then, the binary grid data is used, and the mapping relationship between elevation and color is combined to generate the corresponding bitmap image, which will be used as the source of the texture image matched with the binary grid data.

[0057] Here, dumping the DEM terrain data into binary network data can be achieved through conventional technical means; the mapping relationship between elevation and color reflects that a section of elevation corresponds to a color.

[0058] Step 2: For the weather radar base data to be processed, it contains N conical scanning surfaces with different elevation angles, and each conical scanning surface contains M radial scans, wherein each conical scanning surface records the elevation angle, the length and number of the distance library on the conical scanning surface, and each radial scan records its azimuth in the conical scanning surface to which it belongs. A serial number is provided for each conical scanning surface in the order of the elevation angles of all conical scanning surfaces from small to large, and the serial numbers are 0 to N-1. A serial number is provided for each radial scan contained in any conical scanning surface in the order of the azimuth angles of these radial scans from small to large, and the serial numbers are 0 to M-1.

[0059] Step 3: Based on the weather radar base data, a one-dimensional elevation index array A of size 1800 is established. Array A quickly maps the serial number of the nearest conical scanning surface to any of the 1800 angle values ​​from -90 degrees to 90 degrees with a resolution increment of 0.1 degrees; based on the conical scanning surface with serial number n, a one-dimensional azimuth index array B of size 3600 is established. n , array B n Rapidly map the radial scan number closest to any of the 3600 angle values ​​from 0 to 360 degrees with a resolution of 0.1 degree, where n = 0, 1, ..., N-1.

[0060] In this specific embodiment, in step 3, the process of establishing the one-dimensional elevation angle index array A is as follows:

[0061] Step 3.1a: The subscripts of all elements in array A are 0 to 1799, which correspond to 1800 angle values ​​from -90 degrees to 90 degrees with a resolution of 0.1 degree. The values ​​of all elements in array A are initialized to -1. When the value of an element in array A is -1, it means that the angle value corresponding to the subscript of the element has no corresponding conical scanning surface.

[0062] Step 3.2a: Traverse all conical scanning surfaces contained in the weather radar base data; set the serial number of the currently traversed conical scanning surface to n, record the elevation angle of the currently traversed conical scanning surface as α, and record the elevation angle of the conical scanning surface with serial number n+1 as β, where n=0,1,…,N-2; then calculate the subscript i corresponding to each angle θ with a resolution of 0.1 degrees in the interval [α,α+(β-α) / 2), i=int(θ×10)+90 0, and then assign the element A[i] with subscript i in the array A to the serial number of the currently traversed conical scanning surface, that is, n, where int() is a rounding function; and calculate the subscript j corresponding to each angle η with a resolution of 0.1 degrees in the interval [α+(β-α) / 2,β], j=int(η×10)+900; and then assign the element A[j] with subscript j in the array A to the serial number of the next conical scanning surface of the currently traversed conical scanning surface, that is, n+1.

[0063] Figure 2 A schematic diagram of the index between two elevation angles is given. Array A can quickly map any angle from -90 degrees to 90 degrees to the serial number of the closest conical scanning surface, so as to quickly locate the nearest radar conical scanning surface in three-dimensional structure reconstruction.

[0064] In this specific embodiment, in step 3, the one-dimensional azimuth index array B n The establishment process is:

[0065] Step 3.1b: Array B n The subscripts of all elements in are 0 to 3599, corresponding to 3600 angle values ​​from 0 to 360 degrees with a resolution of 0.1 degree. n All elements in array B are initialized to -1. n When the value of an element in is -1, it means that the angle value corresponding to the subscript of the element has no corresponding radial scan.

[0066] Step 3.2b: Traverse all radial scans contained in the conical scanning surface numbered n; set the number of the currently traversed radial scan to m, record the azimuth angle of the currently traversed radial scan as a, and record the azimuth angle of the radial scan numbered m+1 as b, where m=0, 1,…, M-2.

[0067] When abs(ba) is greater than the missing threshold, keep array B n The value of the element with subscript int(a×10) to the element with subscript int(b×10) remains unchanged at -1, indicating missing measurement, where abs() is the absolute value function.

[0068] When abs(ba) is less than or equal to the missing threshold, calculate the subscript g corresponding to each angle δ in the interval [a, a+(ba) / 2) with a resolution of 0.1 degrees, g = int(δ×10), and then convert the array B n The element B with subscript g n [g] is assigned the serial number of the currently traversed radial scan, i.e. m; and the subscript h corresponding to each angle σ with a resolution of 0.1 degrees in the interval [a+(ba) / 2,b] is calculated, where h=int(σ×10), and then the array B is n The element B with subscript h n [h] is assigned the serial number of the next radial scan of the currently traversed radial scan, that is, m+1.

[0069] Step 3.3b: Based on the azimuth wraparound characteristic, the azimuth of the radial scan numbered M-1 is recorded as c, and the azimuth of the radial scan numbered 0 is recorded as d; then calculate the subscript e corresponding to each angle ρ with a resolution of 0.1 degrees in the interval [c, c+(d+360-c) / 2). If int(ρ×10) is greater than or equal to 3600, then e=int(ρ×10)-3600; if int(ρ×10) is less than 3600, then e=int(ρ×10), and then replace array B with n The element B with subscript e n [e] is assigned a value of M-1; and calculate the subscript f corresponding to each angle ω in the interval [c+(d+360-c) / 2, d+360] with a resolution of 0.1 degrees. If int(ω×10) is greater than or equal to 3600, then f=int(ω×10)-3600; if int(ω×10) is less than 3600, then f=int(ω×10). Then replace array B n The element B with subscript f n [f] is assigned a value of 0.

[0070] Figure 3 A schematic diagram of the index between two azimuths is given. Array B 0 ,B 1 ,…,B N-1 Used to quickly locate the closest radial scan to a specific azimuth.

[0071] Step 4: Within the coverage of the radar of interest on the electronic map, determine a rectangular reconstruction area by dragging a box, and record the coordinate positions of the upper left corner and the lower right corner of the rectangular reconstruction area; then establish a three-dimensional grid on the rectangular reconstruction area, and determine the parameters of the three-dimensional grid, which include the number of grid points in the horizontal longitude and latitude directions, the longitude and latitude of the starting grid point of the horizontal plane, the longitude and latitude of the ending grid point of the horizontal plane, the total height value, and the number of height layers.

[0072] In this specific embodiment, in step 4, the number of horizontal grid points is W, and the W value is dynamically adjusted according to the hardware configuration of the computer, or is a preset fixed value; the longitude and latitude of the starting grid point of the horizontal plane are calculated according to the coordinate position of the upper left corner of the rectangular reconstruction area, and the longitude and latitude of the ending grid point of the horizontal plane are calculated according to the coordinate position of the lower right corner of the rectangular reconstruction area; the total height value is H, and for S-band and C-band radars, the value range of H is 20 to 25 km, which can be adjusted according to user needs in practice. For X-band radar, the value of H is 10 km to 15 km; the number of height layers is L, and the L value is dynamically adjusted according to the hardware configuration of the computer, or is a preset fixed value.

[0073] Step 5: Traverse all grid points in the three-dimensional grid established in step 4. For the currently traversed grid point, calculate the elevation angle φ, azimuth angle γ and distance D of the grid point relative to the radar position based on the latitude, longitude and altitude values ​​of the grid point and the latitude, longitude and altitude values ​​of the radar of interest, where the accuracy of the elevation angle and azimuth angle is taken to 0.1 degree; then, based on the elevation angle φ, azimuth angle γ and distance D, and combined with A and B n , n=0,1,…,N-1, obtain the data of the distance library closest to the grid point on the four radial scans closest to the grid point, and then perform bilinear interpolation or inverse distance square interpolation to obtain the data of the grid point.

[0074] In this specific embodiment, in step 5, the following process is performed for the currently traversed grid point:

[0075] Step 5.1: Based on the latitude, longitude and altitude values ​​of the grid point, as well as the latitude, longitude and altitude values ​​of the radar of interest, calculate the elevation angle φ, azimuth angle γ and distance D of the grid point relative to the radar position, where the accuracy of the elevation angle and azimuth angle is taken to 0.1 degree.

[0076] Step 5.2: According to the elevation angle φ, the azimuth angle γ and the distance D, obtain the data V1, V2, V3, V4 of the distance library closest to the grid point on the four radial scans closest to the grid point. The specific acquisition process is as follows:

[0077] Step 5.2.1: Calculate the subscript p corresponding to the elevation angle φ in array A, p=int(φ×10)+900, and then define the conical scanning surface numbered A[p] as the first conical scanning surface.

[0078] Step 5.2.2: When the elevation angle φ is lower than or equal to the elevation angle of the first conical scanning surface, if A[p]-1 is greater than or equal to 0, the conical scanning surface numbered A[p]-1 is defined as the second conical scanning surface, and then step 5.2.3 is executed; if A[p]-1 is less than 0, it means that the grid point is in the blind area, and the data of the grid point is set to an invalid value.

[0079] When the elevation angle φ is higher than the elevation angle of the first conical scanning surface, if A[p]+1 is less than or equal to N-1, the conical scanning surface numbered A[p]+1 is defined as the second conical scanning surface, and step 5.2.3 is executed; if A[p]+1 is greater than N-1, it means that the grid point is in the blind area, and the data of the grid point is set to an invalid value.

[0080] Here, the first conical scanning surface is the conical scanning surface that is closest to the grid point and located above, and the second conical scanning surface is the conical scanning surface that is closest to the grid point and located below.

[0081] Step 5.2.3: Calculate the one-dimensional azimuth index array B based on the first conical scanning surface A[p] The subscript q corresponding to the azimuth γ is q = int(γ×10). If B A[p] [q] is -1, which means that the radial scan in the first conical scan surface is no longer searched, V1 and V2 are set to invalid values, and step 5.2.6 is executed again; if B A[p] [q] is not -1, then the sequence number is B A[p] The radial scan of [q] is defined as the first radial scan, and step 5.2.4 is performed again.

[0082] Step 5.2.4: When the azimuth angle γ is less than or equal to the azimuth angle of the first radial scan, if B A[p] [q]-1 is greater than or equal to 0, then the sequence number is B A[p] The radial scan of [q]-1 is defined as the second radial scan, and then step 5.2.5 is executed; if B A[p] If [q]-1 is less than 0, the radial scan numbered M-1 is defined as the second radial scan, and step 5.2.5 is executed again.

[0083] When the azimuth angle γ is greater than the azimuth angle of the first radial scan, if B A[p] [q]+1 is less than or equal to M-1, then the sequence number is B A[p] The radial scan of [q]+1 is defined as the second radial scan, and then step 5.2.5 is executed; if B A[p] If [q]+1 is greater than M-1, the radial scan with sequence number 0 is defined as the second radial scan, and step 5.2.5 is executed again.

[0084] Here, the first radial scan is a radial scan that is closest to the grid point in the counterclockwise direction in the first conical scan surface, and the second radial scan is a radial scan that is closest to the grid point in the clockwise direction in the first conical scan surface.

[0085] Step 5.2.5: Determine whether the absolute value of the difference between the azimuth angle of the first radial scan and the azimuth angle of the second radial scan is greater than the missing measurement threshold. If so, set V1 and V2 to invalid values, and then execute step 5.2.6; otherwise, calculate the index t of the distance library, t is equal to the distance D divided by the length of the distance library on the first conical scanning surface. When t is greater than the number of distance libraries on the first conical scanning surface, set V1 and V2 to invalid values, and then execute step 5.2.6; when t is less than or equal to the number of distance libraries on the first conical scanning surface, V1 is the data on the distance library with index t on the first radial scan, and V2 is the data on the distance library with index t on the second radial scan, and then execute step 5.2.6.

[0086] Step 5.2.6: Let s represent the serial number of the second conical scanning surface, and calculate the one-dimensional azimuth index array B established based on the second conical scanning surface s The subscript u corresponding to the azimuth γ is u=int(γ×10). If B s [u] is -1, which means that the radial scan in the second conical scan surface is no longer searched, V3 and V4 are set to invalid values, and step 5.3 is executed again; if B s If [u] is not -1, the sequence number is B s The radial scan of [u] is defined as the third radial scan, and step 5.2.7 is performed again.

[0087] Step 5.2.7: When the azimuth angle γ is less than or equal to the azimuth angle of the third radial scan, if B s [u]-1 is greater than or equal to 0, then the sequence number is B s The radial scan of [u]-1 is defined as the fourth radial scan, and then execute step 5.2.8; if B s If [u]-1 is less than 0, the radial scan numbered M-1 is defined as the fourth radial scan, and step 5.2.8 is executed again.

[0088] When the azimuth angle γ is greater than the azimuth angle of the third radial scan, if B s [u]+1 is less than or equal to M-1, then the sequence number is B s The radial scan of [u]+1 is defined as the fourth radial scan, and then step 5.2.8 is executed; if B s If [u]+1 is greater than M-1, the radial scan with sequence number 0 is defined as the fourth radial scan, and step 5.2.8 is executed again.

[0089] Here, the third radial scan is a radial scan that is closest to the grid point in the counterclockwise direction in the second conical scan surface, and the fourth radial scan is a radial scan that is closest to the grid point in the clockwise direction in the second conical scan surface.

[0090] Step 5.2.8: Determine whether the absolute value of the difference between the azimuth angle of the third radial scan and the azimuth angle of the fourth radial scan is greater than the missing measurement threshold. If so, set V3 and V4 to invalid values, and then execute step 5.3; otherwise, calculate the index z of the distance library, z is equal to the distance D divided by the length of the distance library on the second conical scanning surface. When z is greater than the number of distance libraries on the second conical scanning surface, set V3 and V4 to invalid values, and then execute step 5.3; when z is less than or equal to the number of distance libraries on the second conical scanning surface, V3 is the data on the distance library with index z on the third radial scan, and V4 is the data on the distance library with index z on the fourth radial scan, and then execute step 5.3.

[0091] Step 5.3: If V1, V2, V3, and V4 are all valid values, use the four valid values ​​to perform bilinear interpolation to obtain the data of the grid point; if V1, V2, V3, and V4 contain both invalid values ​​and valid values, use the valid values ​​to perform inverse distance square interpolation to obtain the data of the grid point; if V1, V2, V3, and V4 are all invalid values, then the grid point is set to an invalid value.

[0092] Step 6: Based on step 5, a smoothing range is preset. For the data of each grid point in the three-dimensional grid, if the data of the grid point is within the smoothing range, the data of the grid point is subjected to median filtering to obtain a three-dimensional grid with smooth edges; if the data of the grid point is not within the smoothing range, the data of the grid point is not processed. Here, the preset smoothing range is different depending on the observed data, such as the reflectivity data. The preset smoothing range can be set according to the target to be achieved, such as edge smoothing.

[0093] Step 7: Set the number of grids of the terrain grid data, and set the area formed by the longitude and latitude of the horizontal plane starting grid points and the horizontal plane ending grid points of the three-dimensional grid as the longitude and latitude area of ​​the terrain grid data, and calculate the longitude and latitude resolution of the terrain grid data; then, according to the longitude and latitude area and longitude and latitude resolution of the terrain grid data, resample the terrain grid data of the longitude and latitude area from the binary grid data; then use the terrain grid data to construct triangles according to the diagonal principle to form a terrain grid surface.

[0094] In this specific embodiment, in step 7, the triangle construction process is as follows: the terrain grid data is an equally spaced two-dimensional grid, and the length and width of the terrain grid data are Width and Height respectively; then let (x, y) represent the two-dimensional coordinates of a grid point in the terrain grid data, the value range of x is [0, Width-2], and the value range of y is [0, Height-2], and for the grid point with two-dimensional coordinates (x, y), three grid points with two-dimensional coordinates (x, y+1), (x+1, y+1), and (x+1, y) are selected; then, according to the diagonal principle, triangles with three grid points with two-dimensional coordinates (x, y), (x, y+1), and (x+1, y+1) as vertices and triangles with three grid points with two-dimensional coordinates (x+1, y+1), (x+1, y), and (x, y) as vertices are constructed; and then all triangles form a terrain grid surface.

[0095] Here, the method of constructing triangles, compared with the traditional triangulation method, only needs to combine vertices, which improves the calculation efficiency and can meet the needs of real-time interactive systems.

[0096] Step 8: According to the longitude and latitude of the starting grid point of the horizontal plane of the three-dimensional grid and the longitude and latitude of the ending grid point of the horizontal plane, the area formed by the longitude and latitude of the starting grid point of the horizontal plane of the three-dimensional grid and the longitude and latitude of the ending grid point of the horizontal plane is cut out from the bitmap image as the mapping image of the terrain grid surface.

[0097] Step 9: Send the three-dimensional grid obtained in step 6, the terrain grid surface obtained in step 7, and the texture image of the terrain grid surface obtained in step 8 to the display engine to realize three-dimensional volume rendering and terrain construction to obtain a three-dimensional image.

[0098] In this embodiment, the value range of the missing detection threshold is [3, 5] degrees. In this embodiment, the missing detection threshold is 5 degrees, that is, the missing detection processing is performed when it is greater than 5 degrees.

[0099] In this embodiment, the user can repeatedly perform steps 4 to 9, select rectangular reconstruction areas of different sizes and positions, and obtain a reconstructed three-dimensional image.

[0100] The contents described in the embodiments of this specification are merely an enumeration of the implementation forms of the inventive concept. The protection scope of the present invention should not be regarded as limited to the specific forms described in the embodiments. The protection scope of the present invention also extends to the equivalent technical means that can be conceived by those skilled in the art based on the inventive concept.

Claims

1. A method for reconstructing three-dimensional structure with dynamic resolution of weather radar base data combined with terrain, characterized in that The following steps are involved: Step 1: Select a data source in the form of DEM text, and clip the DEM terrain data. The clipping range is limited to the coverage of the radar of interest. Then, the clipped DEM terrain data is dumped into binary grid data, and the number of grid points, resolution, and latitude and longitude of the starting and ending grid points of the binary grid data are obtained. Then, the binary grid data is used, and the mapping relationship between elevation and color is combined to generate the corresponding bitmap image. Step 2: For the weather radar base data to be processed, it contains N conical scanning surfaces with different elevation angles, each conical scanning surface contains M radial scans, wherein each conical scanning surface records the elevation angle, the length and number of the range library on the conical scanning surface, and each radial scan records its azimuth in the conical scanning surface to which it belongs. A serial number is provided for each conical scanning surface in the ascending order of the elevation angles of all conical scanning surfaces, and the serial number is 0 to N-1. A serial number is provided for each radial scan contained in any conical scanning surface in the ascending order of the azimuth angles of these radial scans, and the serial number is 0 to M-1; Step 3: Based on the weather radar base data, a one-dimensional elevation index array A of size 1800 is established. Array A quickly maps the serial number of the nearest conical scanning surface to any of the 1800 angle values ​​from -90 degrees to 90 degrees with a resolution increment of 0.1 degrees; based on the conical scanning surface with serial number n, a one-dimensional azimuth index array B of size 3600 is established. n , array B n Rapidly map the number of the radial scan closest to any of the 3600 angle values ​​from 0 to 360 degrees with a resolution of 0.1 degree, where n = 0, 1, ..., N-1; Step 4: Within the coverage of the radar of interest on the electronic map, a rectangular reconstruction area is determined by drawing a frame, and the coordinate positions of the upper left corner and the lower right corner of the rectangular reconstruction area are recorded; then a three-dimensional grid is established on the rectangular reconstruction area, and the parameters of the three-dimensional grid are determined, including the number of grid points in the horizontal longitude and latitude directions, the longitude and latitude of the starting grid point of the horizontal plane, the longitude and latitude of the ending grid point of the horizontal plane, the total height value, and the number of height layers; Step 5: Traverse all grid points in the three-dimensional grid established in step 4. For the currently traversed grid point, calculate the elevation angle φ, azimuth angle γ and distance D of the grid point relative to the radar position based on the latitude, longitude and altitude values ​​of the grid point and the latitude, longitude and altitude values ​​of the radar of interest, where the accuracy of the elevation angle and azimuth angle is taken to 0.1 degree; then, based on the elevation angle φ, azimuth angle γ and distance D, and combined with A and B n , n=0,1,…,N-1, obtain the data of the distance library closest to the grid point on the four radial scans closest to the grid point, and then perform bilinear interpolation or inverse square interpolation to obtain the data of the grid point; Step 6: Based on step 5, a smoothing range is preset. For the data of each grid point in the three-dimensional grid, if the data of the grid point is within the smoothing range, the data of the grid point is subjected to median filtering to obtain a three-dimensional grid with smooth edges; if the data of the grid point is not within the smoothing range, the data of the grid point is not processed; Step 7: Set the number of grids of the terrain grid data, and set the area formed by the latitude and longitude of the horizontal plane starting grid point and the horizontal plane ending grid point of the three-dimensional grid as the latitude and longitude area of ​​the terrain grid data, and calculate the latitude and longitude resolution of the terrain grid data; then, according to the latitude and longitude area and the latitude and longitude resolution of the terrain grid data, resample the terrain grid data of the latitude and longitude area from the binary grid data; then use the terrain grid data to construct triangles according to the diagonal principle, and then form a terrain grid surface; Step 8: According to the longitude and latitude of the starting grid point of the horizontal plane of the three-dimensional grid and the longitude and latitude of the ending grid point of the horizontal plane, a region consisting of the longitude and latitude of the starting grid point of the horizontal plane of the three-dimensional grid and the longitude and latitude of the ending grid point of the horizontal plane is cut out from the bitmap image as a mapping image of the terrain grid surface; Step 9: Send the three-dimensional grid obtained in step 6, the terrain grid surface obtained in step 7, and the texture image of the terrain grid surface obtained in step 8 to the display engine to realize three-dimensional volume rendering and terrain construction to obtain a three-dimensional image.

2. The method for reconstructing a three-dimensional structure with dynamic resolution of weather radar base data combined with terrain according to claim 1 is characterized in that In step 3, the process of establishing the one-dimensional elevation index array A is as follows: Step 3.1a: The subscripts of all elements in array A are 0 to 1799, which correspond to 1800 angle values ​​from -90 degrees to 90 degrees with a resolution of 0.1 degree. The values ​​of all elements in array A are initialized to -1. When the value of an element in array A is -1, it means that the angle value corresponding to the subscript of the element has no corresponding conical scanning surface. Step 3.2a: Traverse all conical scanning surfaces contained in the weather radar base data; set the serial number of the currently traversed conical scanning surface to n, record the elevation angle of the currently traversed conical scanning surface as α, and record the elevation angle of the conical scanning surface with serial number n+1 as β, where n=0,1,…,N-2; then calculate the subscript i corresponding to each angle θ in the interval [α,α+(β-α) / 2] with a resolution of 0.1 degrees, i=int(θ×10)+90 0, and then assign the element A[i] with subscript i in the array A to the serial number of the currently traversed conical scanning surface, that is, n, where int() is a rounding function; and calculate the subscript j corresponding to each angle η with a resolution of 0.1 degrees in the interval [α+(β-α) / 2,β], j=int(η×10)+900; and then assign the element A[j] with subscript j in the array A to the serial number of the next conical scanning surface of the currently traversed conical scanning surface, that is, n+1.

3. The method for reconstructing a three-dimensional structure with dynamic resolution of weather radar base data combined with terrain according to claim 2, characterized in that In step 3, the one-dimensional azimuth index array B n The establishment process is: Step 3.1b: Array B n The subscripts of all elements in are 0 to 3599, corresponding to 3600 angle values ​​from 0 to 360 degrees with a resolution of 0.1 degree. n All elements in array B are initialized to -1. n When the value of an element in is -1, it means that the angle value corresponding to the subscript of the element has no corresponding radial scan; Step 3.2b: traverse all radial scans contained in the conical scanning surface with serial number n; set the serial number of the currently traversed radial scan to m, record the azimuth angle of the currently traversed radial scan as a, and record the azimuth angle of the radial scan with serial number m+1 as b, where m=0, 1, ..., M-2; When abs(ba) is greater than the missing threshold, keep array B n The value of the element with subscript int(a×10) to the element with subscript int(b×10) remains unchanged at -1, indicating missing measurement, where abs() is the absolute value function; When abs(ba) is less than or equal to the missing threshold, calculate the subscript g corresponding to each angle δ in the interval [a, a+(ba) / 2] with a resolution of 0.1 degrees, g = int(δ×10), and then convert the array B n The element B with subscript g n [g] is assigned the serial number of the currently traversed radial scan, i.e. m; and the subscript h corresponding to each angle σ with a resolution of 0.1 degrees in the interval [a+(ba) / 2,b] is calculated, where h=int(σ×10), and then the array B is n The element B with subscript h n [h] is assigned the serial number of the next radial scan of the currently traversed radial scan, i.e. m+1; Step 3.3b: Based on the azimuth wraparound characteristic, the azimuth of the radial scan numbered M-1 is recorded as c, and the azimuth of the radial scan numbered 0 is recorded as d; then calculate the subscript e corresponding to each angle ρ with a resolution of 0.1 degrees in the interval [c, c+(d+360-c) / 2]. If int(ρ×10) is greater than or equal to 3600, then e=int(ρ×10)-3600; if int(ρ×10) is less than 3600, then e=int(ρ×10). Then, replace array B with n The element B with subscript e n [e] is assigned a value of M-1; and calculate the subscript f corresponding to each angle ω in the interval [c+(d+360-c) / 2, d+360] with a resolution of 0.1 degrees. If int(ω×10) is greater than or equal to 3600, then f=int(ω×10)-3600; if int(ω×10) is less than 3600, then f=int(ω×10). Then replace array B n The element B with subscript f n [f] is assigned a value of 0.

4. The method for reconstructing a three-dimensional structure with dynamic resolution of weather radar base data combined with terrain according to claim 3 is characterized in that In step 5, the following process is performed for the currently traversed grid point: Step 5.1: Based on the latitude, longitude and altitude values ​​of the grid point, as well as the latitude, longitude and altitude values ​​of the radar concerned, calculate the elevation angle φ, azimuth angle γ and distance D of the grid point relative to the radar position, where the accuracy of the elevation angle and azimuth angle is taken to 0.1 degree; Step 5.2: According to the elevation angle φ, the azimuth angle γ and the distance D, obtain the data V1, V2, V3, and V4 of the distance library closest to the grid point on the four radial scans closest to the grid point. The specific acquisition process is as follows: Step 5.2.1: Calculate the subscript p corresponding to the elevation angle φ in array A, p = int(φ×10)+900, and then define the conical scanning surface numbered A[p] as the first conical scanning surface; Step 5.2.2: When the elevation angle φ is lower than or equal to the elevation angle of the first conical scanning surface, if A[p]-1 is greater than or equal to 0, the conical scanning surface numbered A[p]-1 is defined as the second conical scanning surface, and then step 5.2.3 is executed; if A[p]-1 is less than 0, it means that the grid point is in the blind area, and the data of the grid point is set to an invalid value; When the elevation angle φ is higher than the elevation angle of the first conical scanning surface, if A[p]+1 is less than or equal to N-1, the conical scanning surface numbered A[p]+1 is defined as the second conical scanning surface, and step 5.2.3 is executed again; if A[p]+1 is greater than N-1, it means that the grid point is in the blind area, and the data of the grid point is set to an invalid value; Step 5.2.3: Calculate the one-dimensional azimuth index array B based on the first conical scanning surface A[p] The subscript q corresponding to the azimuth γ is q=int(γ×10). If B A[p] [q] is -1, which means that the radial scan in the first conical scan surface is no longer searched, V1 and V2 are set to invalid values, and step 5.2.6 is executed again; if B A[p] [q] is not -1, then the sequence number is B A[p] The radial scan of [q] is defined as the first radial scan, and then step 5.2.4 is performed; Step 5.2.4: When the azimuth angle γ is less than or equal to the azimuth angle of the first radial scan, if B A[p] [q]-1 is greater than or equal to 0, then the sequence number is B A[p] The radial scan of [q]-1 is defined as the second radial scan, and then step 5.2.5 is executed; if B A[p] If [q]-1 is less than 0, the radial scan numbered M-1 is defined as the second radial scan, and step 5.2.5 is executed again; When the azimuth angle γ is greater than the azimuth angle of the first radial scan, if B A[p] [q]+1 is less than or equal to M-1, then the sequence number is B A[p] The radial scan of [q]+1 is defined as the second radial scan, and then step 5.2.5 is executed; if B A[p] If [q]+1 is greater than M-1, the radial scan with sequence number 0 is defined as the second radial scan, and step 5.2.5 is executed again; Step 5.2.5: Determine whether the absolute value of the difference between the azimuth angle of the first radial scan and the azimuth angle of the second radial scan is greater than the missing measurement threshold. If so, set V1 and V2 to invalid values, and then execute step 5.2.6; otherwise, calculate the index t of the distance library, t is equal to the distance D divided by the length of the distance library on the first conical scanning surface. When t is greater than the number of distance libraries on the first conical scanning surface, set V1 and V2 to invalid values, and then execute step 5.2.6; when t is less than or equal to the number of distance libraries on the first conical scanning surface, V1 is the data on the distance library with index t on the first radial scan, and V2 is the data on the distance library with index t on the second radial scan, and then execute step 5.2.6; Step 5.2.6: Let s represent the serial number of the second conical scanning surface, and calculate the one-dimensional azimuth index array B established based on the second conical scanning surface s The subscript u corresponding to the azimuth γ is u=int(γ×10). If B s [u] is -1, which means that the radial scan in the second conical scan surface is no longer searched, V3 and V4 are set to invalid values, and step 5.3 is executed again; if B s If [u] is not -1, the sequence number is B s The radial scan of [u] is defined as the third radial scan, and step 5.2.7 is performed again; Step 5.2.7: When the azimuth angle γ is less than or equal to the azimuth angle of the third radial scan, if B s [u]-1 is greater than or equal to 0, then the sequence number is B s The radial scan of [u]-1 is defined as the fourth radial scan, and then execute step 5.2.8; if B s If [u]-1 is less than 0, the radial scan numbered M-1 is defined as the fourth radial scan, and step 5.2.8 is executed again; When the azimuth angle γ is greater than the azimuth angle of the third radial scan, if B s [u]+1 is less than or equal to M-1, then the sequence number is B s The radial scan of [u]+1 is defined as the fourth radial scan, and then step 5.2.8 is executed; if B s If [u]+1 is greater than M-1, the radial scan with sequence number 0 is defined as the fourth radial scan, and step 5.2.8 is executed again; Step 5.2.8: Determine whether the absolute value of the difference between the azimuth angle of the third radial scan and the azimuth angle of the fourth radial scan is greater than the missing measurement threshold. If so, set V3 and V4 to invalid values, and then execute step 5.3; otherwise, calculate the index z of the distance library, z is equal to the distance D divided by the length of the distance library on the second conical scanning surface. When z is greater than the number of distance libraries on the second conical scanning surface, set V3 and V4 to invalid values, and then execute step 5.3; when z is less than or equal to the number of distance libraries on the second conical scanning surface, V3 is the data on the distance library with index z on the third radial scan, and V4 is the data on the distance library with index z on the fourth radial scan, and then execute step 5.3; Step 5.3: If V1, V2, V3, and V4 are all valid values, use the four valid values ​​to perform bilinear interpolation to obtain the data of the grid point; if V1, V2, V3, and V4 contain both invalid values ​​and valid values, use the valid values ​​to perform inverse distance square interpolation to obtain the data of the grid point; if V1, V2, V3, and V4 are all invalid values, then the grid point is set to an invalid value.

5. The method for reconstructing a three-dimensional structure with dynamic resolution of weather radar base data combined with terrain according to claim 1, characterized in that In step 7, the triangle construction process is as follows: the terrain grid data is an equally spaced two-dimensional grid, and the length and width of the terrain grid data are Width and Height respectively; then let (x, y) represent the two-dimensional coordinates of a grid point in the terrain grid data, the value range of x is [0, Width-2], and the value range of y is [0, Height-2], and for the grid point with two-dimensional coordinates (x, y), three grid points with two-dimensional coordinates (x, y+1), (x+1, y+1), and (x+1, y) are selected; then, according to the diagonal principle, a triangle with three grid points with two-dimensional coordinates (x, y), (x, y+1), and (x+1, y+1) as vertices and a triangle with three grid points with two-dimensional coordinates (x+1, y+1), (x+1, y), and (x, y) as vertices are constructed; and then all the triangles form a terrain grid surface.

6. The method for reconstructing a three-dimensional structure with dynamic resolution of weather radar base data combined with terrain according to claim 4, characterized in that The value range of the missing detection threshold is [3,5] degrees.

7. The method for reconstructing a three-dimensional structure with dynamic resolution of weather radar base data combined with terrain according to claim 1, characterized in that In step 4, the number of horizontal grid points is W, and the W value is dynamically adjusted according to the hardware configuration of the computer, or is a preset fixed value; the longitude and latitude of the starting grid point of the horizontal plane are calculated according to the coordinate position of the upper left corner of the rectangular reconstruction area, and the longitude and latitude of the ending grid point of the horizontal plane are calculated according to the coordinate position of the lower right corner of the rectangular reconstruction area; the total height value is H, and for S-band and C-band radars, the value range of H is 20 to 25 km, and for X-band radars, the value of H is 10 km to 15 km; the number of height layers is L, and the L value is dynamically adjusted according to the hardware configuration of the computer, or is a preset fixed value.

Images