A method for quickly obtaining DEM data based on two-point equidistant projection for calculating radar shadow angle

Abstract

This invention discloses a method for rapidly acquiring DEM data based on two-point equidistant projection for calculating radar shielding angles. First, a unified DEM data coordinate system is established and a VRT virtual dataset is constructed. Then, based on a known point P, the corresponding point P with an azimuth angle of n° and a distance H from the positioning point P is calculated. (n) Calculate P and P (n) The projected coordinates are used to construct the data extraction range, and further raster data within the range is extracted and saved; finally, the maximum radar obstruction angle is calculated. After calculating the maximum radar obstruction angle corresponding to multiple azimuth angles using the above method, the data is drawn according to the radar obstruction angle values. This invention cleverly transforms the traversal of data points into a single row of raster data clipping, greatly accelerating the speed of DEM data acquisition.

Description

Technical Field

[0001] This invention belongs to the field of radar technology, and specifically designs a method for rapidly acquiring DEM data based on two-point equidistant projection to calculate radar shielding angle. Background Technology

[0002] The obscuring angle is the angle between the line of sight and the top of an obstacle and the horizon. When the elevation angle of the radar antenna is less than this angle, the radar will be unable to detect the target due to terrain obstruction. The detection blind zone formed by terrain blocking radar signal propagation can be obtained by calculating the obscuring angle.

[0003] In military simulation systems, to optimize radar deployment and enhance combat capabilities, Digital Elevation Models (DEMs) are typically used to acquire terrain data and calculate radar obstruction angles. This allows simulation engineers to intuitively understand the obstruction conditions of the surrounding environment around the proposed radar deployment area. Currently, obstruction angle calculation faces several challenges: First, the sheer volume of DEM data leads to long loading times and high memory consumption. Second, calculating obstruction angles requires extensive calculations of terrain sampling points, which existing methods are time-consuming. Third, it places high demands on computer performance.

[0004] Solving the problems of long loading time, high memory consumption, and long calculation time for occlusion angles in DEM is a current challenge. Summary of the Invention

[0005] To address the aforementioned issues, this invention proposes a method for rapidly acquiring DEM data based on two-point equidistant projection for calculating radar shielding angles. This method does not rely on GPUs to quickly acquire DEM data and achieve shielding angle calculation.

[0006] This invention relates to a method for rapidly acquiring DEM data based on two-point equidistant projection for calculating radar shielding angles. The specific steps are as follows:

[0007] Step 1: Download the target area DEM data to a directory; then process the target area DEM data into a unified coordinate system.

[0008] Step 2: Construct a VRT virtual dataset based on the data results from Step 1.

[0009] Step 3: Based on the known positioning point P, find the corresponding point P with an azimuth angle of n° and a distance H from the positioning point P. (n) ;

[0010] Step 4: Define the algorithm parameters for the equidistant projection of the two points P and P(n).

[0011] Step 5: Define the scope of the extracted data.

[0012] Step 6: Extract raster data within the range of data extracted in Step 5.

[0013] Step 7: Extract and save the output raster data from Step 6.

[0014] Step 8: Calculate the maximum radar shielding angle in the azimuth direction n.

[0015] Step 9: Iteratively calculate the radar shielding angles corresponding to multiple azimuth angles at intervals of m°.

[0016] Step 10: Draw the radar shielding angle value obtained in Step 9 according to the actual size.

[0017] The advantages of this invention are:

[0018] 1. The method for calculating radar shielding angle in this invention does not rely on GPUs to quickly acquire DEM data and calculate the shielding angle, thereby reducing computer hardware costs. Furthermore, it can be widely deployed on cloud servers without GPUs, broadening its application scope.

[0019] 2. The method for calculating radar shielding angle in this invention supports large-scale data expansion and is simple and easy to use.

[0020] 3. The method for calculating radar shielding angle in this invention uses a grid clipping method to acquire DEMs in batches, reducing unnecessary iterative calculations.

[0021] 4. The method for calculating the radar shielding angle in this invention applies the two-point equidistant projection method to grid data clipping, so that only one row of data for the target area is obtained, making data extraction extremely convenient and fast. Attached Figure Description

[0022] Figure 1 This is a flowchart of the method for rapidly acquiring DEM data based on two-point equidistant projection to calculate the radar shielding angle according to the present invention.

[0023] Figure 2 Point P is a location with an azimuth angle of n° and an interval distance of H. (n) Location diagram;

[0024] Figure 3 This is a schematic diagram of an equidistant projection method;

[0025] Figure 4 This is a schematic diagram of the coordinate range of the rectangular frame;

[0026] Figure 5 The results of the radar shielding angle are plotted. Detailed Implementation

[0027] The present invention will now be described in further detail with reference to the accompanying drawings.

[0028] This invention relates to a method for rapidly acquiring DEM data based on two-point equidistant projection for calculating radar shielding angles, such as... Figure 1 As shown, the specific steps are as follows:

[0029] Step 1: Prepare DEM data and establish coordinate system 1.

[0030] Since currently available open-source DEM data is often in partitioned blocks, the target area DEM data is downloaded to a directory without needing to be concatenated. Then, tools such as QGIS or gdal_translate (GDAL (Geospatial Data Abstraction Library)) are used to process the target area DEM data into a unified coordinate system, such as WGS84.

[0031] Step 2: Use the gdalbuildvrt tool to build a VRT (Virtual Raster Tables) dataset based on the data results from Step 1. The VRT data format sequence is a format driver for GDAL, which allows virtual GDAL datasets to be composed of other GDAL datasets and have relocation, possible applied algorithms, and various metadata that have been changed or added. The VRT description of the dataset can be saved in XML format, usually using the vrt extension.

[0032] Therefore, unlike previous methods that involved creating spatial indexes from image data or stitching data together to form a complete image dataset, or expanding data through a spatial database, this invention uses a VRT virtual dataset. This allows for the rapid construction of a complete DEM of the target area and facilitates area expansion.

[0033] Step 3: Based on the known location point P, find the corresponding point P with an azimuth angle of n° and an interval distance of H. (n) .

[0034] like Figure 2 As shown, given the coordinates (x, y) of a point P, we need to find the point P' with an azimuth angle of n and a distance H from P. (n) coordinates (x) n ,y n If the coordinate system of the target area DEM data in step 1 is a plane coordinate system, then P (n) The coordinates can be calculated using the Pythagorean theorem; if the DEM data of the target area in step 1 uses a spherical coordinate system, such as WGS84, then it is necessary to refer to the corresponding spherical coordinate coefficients to calculate P. (n) The corresponding coordinates on the spherical (ellipsoidal) coordinate system.

[0035] Step 4: Using the open-source tool Proj4, define the positioning point P and its corresponding point P with an azimuth angle of n and a distance of H. (n) The algorithm parameters for two-point equidistant projection.

[0036] Two-point equidistant projection is a modified azimuth projection that preserves the distances from two selected points on a map. If two points (P and P') are equidistant from each other... (n) If the coordinates (x, y) are the same, then the generated projection is an equidistant azimuth projection. For example: the location point P(x, y) and the azimuth point P... (n) (x (n) ,y (n) The algorithm parameters for two-point equidistant projection are:

[0037] +proj=tpeqd+lon_1=x+lat_1=y+lon_2=x (n) +lat_2=y (n) +datum=WGS84+units=m+no_defs+type=crs

[0038] Step 5: Define the scope of the extracted data.

[0039] Using the projection parameters defined in step 4, P and P can be calculated using the porj4 tool. (n) The projected coordinates. For example... Figure 3 As shown, according to the principle of equidistant projection of two points, P and P (n) The midpoint is the origin after projection.

[0040] like Figure 4 As shown, calculate the coordinate range of the rectangle for use in step 6.

[0041] Xmin = x' – PixelSize / 2

[0042] Ymin = y' – PixelSize / 2

[0043] Xmax = x (n) +PixelSize / 2

[0044] Ymax = y (n) +PixelSize / 2

[0045] Where x', y' are the x and y coordinates of point P after projection, x (n) ',y (n) 'For P (n) After projection of point P (n)The x and y coordinates of the point. PixelSize is the resolution of the current DEM raster data, which is a known parameter; Xmin and Xmax are the minimum and maximum x coordinates of the rectangle, respectively; Ymin and Ymax are the minimum and maximum y coordinates of the rectangle, respectively.

[0046] Step 6: Use the gdal_warp method to extract the raster data within the rectangular frame calculated in Step 5.

[0047] like:

[0048] gdalwarp dem_data.vrt output.tif-overwrite -teXmin, Ymin, XMax, YMax -t_srs″+proj=tpeqd+lon_1=x+lat_1=y+lon_2=x(n)+lat_2=y(n)+datum=WGS84+ units=m+no_defs+type=crs″-ts(Xmax-Xmin) / PixelSize1 -r near

[0049] The parameters indicated by underscores have the following meanings:

[0050] dem_data.vrt is the virtual dataset constructed in step 2;

[0051] output.tif is the address of the output raster file, which will be used in subsequent step 7;

[0052] -te Xmin,Ymin,XMax,YMax are the rectangular ranges calculated in step 5;

[0053] -t_srs″+proj=tpeqd+lon_1=x+lat_1=y+lon_2=x(n)+lat_2=y(n)+datum=WGS84+units=m+no_defs+type=crs″ are the projection method parameters defined in step 4;

[0054] -ts(Xmax-Xmin) / PixelSize 1 is the metadata parameter for the output raster file;

[0055] The parameter (Xmax-Xmin) / PixelSize 1 indicates the output raster data. (Xmax-Xmin) / PixelSize List However, there is only one line.

[0056] Step 7: Extract and save the output raster data from Step 6.

[0057] Using GDAL or other raster data reading tools, read the raster data output.tif from step 6 and obtain the result:

[0058] num index =[num1,num2,num3...num k ]

[0059] The result is a single row with multiple columns, and this data is the extracted DEM data required for radar shielding angle calculation.

[0060] Where k = (Xmax - Xmin) / PixelSize.

[0061] For example: [1397.099 1360.8303 1339.6602...76.00187 96.683014 131.35663]

[0063] Step 8: Calculate the maximum radar shielding angle in the azimuth direction n.

[0064] There are various methods for calculating radar shielding angles, which are not the focus of this study. This method selects one of these formulas, which takes into account standard atmospheric refraction and the equivalent Earth radius (8500 km), for calculation. The specific method is as follows:

[0065]

[0066] Among them, h a is the vertical distance (m) from the radar antenna height (above sea level), h is the height of the ground obstacle (m), and d is the horizontal distance (m) from the ground obstacle to the radar.

[0067] h a =num1 + radar antenna height (known)

[0068] d = PixelSize * index(0 <index<k+1)

[0069] h = num index

[0070] When index = 1, 2, 3...k, the maximum value of α can be obtained by iterative calculation, which is the maximum radar shielding angle of that azimuth angle n.

[0071] Step 9: Steps 3-8 provide detailed steps for calculating the radar shielding angle in the direction of azimuth n. Since the radar shielding angle often requires calculation across the entire 360° direction, it is necessary to iteratively calculate the corresponding radar shielding angle when the azimuth n is 1°, 2°, 3°, ..., 360°.

[0072] If higher accuracy is required for radar shielding angle calculation, such as a 0.5° sampling interval, then the azimuth angle n needs to be calculated for 720 cycles of 0.5°, 1°, 1.5°, ..., 360°.

[0073] Step 10: Step 9 will return a calculated radar shielding angle value corresponding to 360 (or 720) azimuth angles. These values ​​will be plotted according to their actual sizes. The final result is as follows: Figure 5 As shown.

[0074] This invention presents a method for rapidly acquiring DEM data for calculating radar shielding angles based on two-point equidistant projection. By cleverly transforming the traversal of data points into a single row of raster data, the speed of DEM data acquisition is significantly accelerated. Without two-point equidistant projection, the acquired DEM data would be multi-row and multi-column, making it impossible to quickly load target data into memory. The only solution is to use the currently known invention patent (patent number 202010156881.1), "A Radar Shielding Angle Calculation Method and System Based on DEM," which upgrades hardware computing power to increase the speed of DEM data acquisition. Extensive experiments have shown that acquiring the DEM data required for calculating the shielding angle is the most time-consuming module. This patent does not provide a clear implementation method for addressing this issue; furthermore, this patent relies on a GPU to achieve high-performance computing results, making its applicability and cost significantly higher than this solution. Additionally, the process of constructing the DEM index file in this patent is overly cumbersome, requiring detailed information such as the data range of the DEM. In summary, the method of this invention cleverly utilizes the VRT dataset of GDAL, which greatly simplifies the data indexing process.

Claims

1. A method for rapidly acquiring DEM data based on two-point equidistant projection for calculating radar shielding angle, characterized in that: The specific steps are as follows: Step 1: Download the target area DEM data to a directory; further process the target area DEM data into a unified coordinate system; Step 2: Construct a VRT virtual dataset based on the data results from Step 1; Step 3: Based on the known positioning point P, find the corresponding point P with an azimuth angle of n° and a distance H from the positioning point P. (n) ; Step 4: Define the algorithm parameters for the equidistant projection of the two points P and P(n); Step 5: Define the scope of the data to be extracted; Step 6: Extract raster data within the range of data extracted in Step 5; Step 7: Extract and save the output raster data from Step 6; Extract the raster data output in step 6 and obtain the following result: index =[ 1, 2, 3... k ] The result is a single row with multiple columns, and this data is the extracted DEM data required for radar shielding angle calculation; Where k = (Xmax - Xmin) / PixelSize; PixelSize is the resolution of the current DEM raster data; Xmin and Xmax are the minimum and maximum x-coordinates of the rectangle, respectively; Step 8: Calculate the maximum radar shielding angle in the azimuth n direction; Step 9: Iteratively calculate the radar shielding angles corresponding to multiple azimuth angles at intervals of m°; Step 10: Draw the radar shielding angle value obtained in Step 9 according to the actual size.

2. The method for rapidly acquiring DEM data based on two-point equidistant projection for calculating radar shielding angle as described in claim 1, characterized in that: The method for constructing the range of data extracted in step 5 is as follows: First, calculate P and P (n) The projected coordinates are used to extract the data range of the rectangle's coordinates. Xmin=x'–PixelSize / 2 Ymin=y'–PixelSize / 2 Xmax=x (n) '+PixelSize / 2 Ymax=y (n) '+PixelSize / 2 Where x', y' are the x and y coordinates of point P after projection, x (n) ',y (n) 'For P (n) After projection of point P (n) The x and y coordinates of the point; PixelSize is the resolution of the current DEM raster data; Xmin and Xmax are the minimum and maximum x coordinates of the rectangle, respectively; Ymin and Ymax are the minimum and maximum y coordinates of the rectangle, respectively.

3. The method for rapidly acquiring DEM data based on two-point equidistant projection for calculating radar shielding angle as described in claim 1, characterized in that: In step 9, the interval m° is 1° or 0.5°.

Images