A fast GPU-based method for computing terrain occlusions

Abstract

The application discloses a terrain occlusion fast calculation method based on GPU. The method takes a digital elevation model (DEM) as input, determines an observation point and a region to be calculated, generates a plurality of rays in an azimuth angle range covering the region, and performs ray traversal on a GPU in parallel: threads calculate a direction vector according to a ray number, perform step sampling at a preset step length, and map to a DEM grid, calculate an elevation angle or a slope parameter of a sampling point relative to the observation point, and maintain an extreme value, to complete a visible / occlusion determination, so that an initial occlusion matrix is obtained; further, morphological erosion post-processing is performed on the initial occlusion matrix to suppress holes and boundary discontinuity caused by discrete sampling, and a stable terrain occlusion distribution result is output.

Description

Technical Field

[0001] This invention relates to the field of geographic information processing and terrain analysis technology, and in particular to a method for rapid calculation of terrain occlusion based on a digital elevation model (DEM). Specifically, it utilizes the parallel computing power of a graphics processing unit (GPU) to efficiently determine the visibility / occlusion relationship between observation points and terrain grids and output occlusion distribution results. This method can be used for spatial analysis tasks such as view analysis, visibility assessment, coverage analysis, and simulation environment modeling. Technical Background

[0002] Digital Elevation Models (DEMs) describe the distribution of land surface elevation in a raster format. Terrain occlusion calculation based on DEMs is a typical terrain visibility analysis problem. Its basic objective is to determine, for a given observation point and a specified analysis range, whether the line of sight between the observation point and the surface unit is obstructed by higher terrain, thereby obtaining results such as occlusion distribution, visible range, or horizon outline. These results are often used as basic inputs for tasks such as viewpoint analysis, coverage analysis, spatial accessibility assessment, and simulation environment modeling.

[0003] Among existing methods for calculating terrain occlusion, the most common approach is line-of-sight / ray analysis: This involves constructing terrain profiles or paths in multiple directions and discretely sampling elevations along these paths. Visibility at subsequent locations is determined by comparing relative height relationships along the path or extreme occlusion values ​​(e.g., maximum elevation angle, maximum slope, or horizon angle) along the maintained direction. While this approach is conceptually clear, it has shortcomings in practical applications: Firstly, under conditions of large-scale, high-resolution, or high-angular-resolution analysis, the number of directions and sampling points required increases significantly, leading to a rapid increase in the overall computational scale. Secondly, due to the influence of rasterization and discrete sampling, occlusion determination results are prone to local discontinuities in space, such as jagged edges, small burrs, holes, or isolated misclassified points. This results in unstable representation of the occlusion area, affecting the reliability of subsequent spatial analysis and engineering usability.

[0004] In addition, some methods adopt the idea of ​​constructing the boundary first and then generating the region to obtain the occlusion / visibility mask. However, such methods are more sensitive to the continuity of the boundary and the region generation process. When there are discrete errors or breaks in the boundary, it is easy to cause the region results to be unstable or require additional repair processing, which increases the implementation complexity.

[0005] Therefore, it is still necessary to propose a fast calculation method for terrain shading in DEM, which can better adapt to the needs of large-scale data processing while ensuring the effectiveness of shading determination, and output shading distribution results with better spatial continuity and higher stability. Summary of the Invention

[0006] This invention aims to provide a GPU-based method for fast calculation of terrain occlusion, which can efficiently obtain the occlusion distribution results within the calculation area under the condition of large-scale digital elevation model (DEM) data. Based on this, the continuity and stability of the occlusion region representation are improved through image post-processing to meet the requirements of calculation speed and result stability in engineering applications.

[0007] This invention provides a GPU-based method for fast terrain occlusion calculation, comprising the following steps:

[0008] Step 1: Acquire digital elevation model (DEM) data and organize the DEM data into a two-dimensional raster format; establish a shading result output matrix with the same resolution as the DEM to record the shading status of each raster unit within the area to be calculated;

[0009] Step 2: Determine the observation point And the terrain-shading area Ω to be calculated; wherein, the observation point Located inside or outside the shielded area;

[0010] Step 3: Define the rectangular computational region and construct the ray set;

[0011] Let the four corner points of the rectangular calculation region be... Calculate the azimuth angle of each corner point relative to the observation point. :

[0012] ;

[0013] right Angle domain processing is performed to obtain the minimum continuous azimuth angle interval covering the rectangular calculation area. ], and generate rays within that azimuth angle range according to a preset number of rays or angular resolution, where the number of rays is N:

[0014] ;

[0015] For angular resolution: ;

[0016] The preset number of rays is used to ensure that multiple boundary sampling points are uniformly selected on the boundary of the computational region to guarantee the ray coverage area. }, from the observation point A ray set is formed pointing towards each boundary sampling point;

[0017] Step 4: Ray traversal and elevation angle extreme value occlusion determination;

[0018] For each ray direction, traverse the terrain units along that direction in order from nearest to farthest; for any terrain unit on the ray... Define horizontal distance for:

[0019] ;

[0020] Define elevation angle parameters :

[0021] ;

[0022] Maintain the current extreme elevation angle for each ray. During the traversal, occlusion is determined according to the following rules, and the initial occlusion result is output. :

[0023] like The terrain unit is determined to be unobstructed, and updates are performed. ;

[0024] like The terrain unit is determined to be occluded.

[0025] The initial occlusion result For two-dimensional raster results with the same resolution as the DEM data:

[0026] ;

[0027] Since the traversal and extremum maintenance of different rays are independent of each other, the occlusion calculation is performed on multiple rays in parallel by the GPU parallel computing unit.

[0028] Step 5: Apply the initial occlusion results obtained in Step 4 Image post-processing, including at least morphological erosion, is performed to reduce local discontinuities, isolated points, and holes caused by discrete ray and raster sampling, thereby obtaining stable occlusion results. ;

[0029] Step 6: Use the occlusion result M1(x,y) obtained from post-processing as the input constraint for subsequent spatial analysis or simulation modeling: Mask the occluded terrain units in the coverage analysis / visibility analysis or attenuate them according to the preset weight coefficient to obtain the analysis results that take into account the terrain occlusion effect.

[0030] Furthermore, the direction vector of the ray in step 3 is determined by the following formula:

[0031] ;

[0032] Indicates the first Ray in a plane coordinate system Components in direction, Indicates the first A ray in a plane coordinate system The component in the direction, and ( ( ) represents the corresponding azimuth angle The unit direction vector in the plane.

[0033] Furthermore, in step 3, the ray azimuth angle covers the entire angular domain [0, 2π), and satisfies:

[0034] .

[0035] Furthermore, the method for setting the number of rays N in step 3 includes:

[0036] Method 1: The number of rays N is set proportionally to the discrete perimeter P of the region to be calculated, satisfying N=kP, where k>0 is the scaling factor.

[0037] Method 2: Determined by angular resolution. If the angular resolution is determined... So there are a total of / A ray.

[0038] Furthermore, in step 4, the mapping from the sampling points to the DEM data raster index uses a rounding method:

[0039] ;

[0040] in, This represents the DEM raster column index after the sampling points are mapped. This represents the DEM raster row index after the sampling points are mapped. This indicates the floor function.

[0041] Furthermore, the morphological etching process in step 5 satisfies:

[0042] ;

[0043] Where S represents the structural element used in the morphological processing, which defines the shape and scale of the local neighborhood of the erosion operation, preferably... Square structural element; This represents a morphological erosion operation used to reduce isolated points, local holes, and boundary discontinuities in the initial occlusion result, thereby obtaining an occlusion result with better spatial continuity. .

[0044] Furthermore, steps 1 to 6 can be executed in series according to the following process:

[0045] First, the host computer reads the digital elevation model (DEM) raster data and linearizes it using either row-major or column-major order to obtain a one-dimensional DEM array suitable for device access. Simultaneously, observation points are established on the host computer. Region to be calculated DEM dimensional parameters and Raster resolution Total number of rays Step length and safety redundancy coefficient The setup, where the process corresponds to the input preparation and parameter initialization in steps 1, 2 and 3;

[0046] Subsequently, a DEM array and a masking output matrix are requested on the device side. The storage space copies the linearized DEM data from the host to the device; the masking output matrix Initialize the values ​​on the device side to preset values, preferably to all 1s, to implement a conservative decision-making strategy for default occlusion. This process corresponds to the establishment and initialization of the occlusion result matrix in step 1.

[0047] After data transfer and video memory allocation are completed, the host side operates according to the corner range determined in step 3. and total number of rays Configure and start the parallel computing task on the device side to establish a correspondence between each parallel execution unit and each traversable ray; wherein, each parallel execution unit is responsible for sampling traversal, maintenance of elevation angle extreme values ​​and writing of occlusion results in the direction of a ray, thereby completing the ray set construction in step 3 and the parallel task allocation in step 4.

[0048] During parallel computing on the device side, each parallel execution unit executes its corresponding ray number. Calculate the azimuth angle of the ray. And further obtain the corresponding planar unit direction vector. ; by step length By advancing step by step along the ray direction, the coordinates of each sampling point are obtained in sequence. And map the sampling points to the DEM raster index. When the sampling points exceed the valid boundaries of the DEM or the current traversal step count reaches the maximum step count limit; When the sampling traversal corresponding to the ray is terminated, the ray step sampling and boundary termination control in step 4 are completed.

[0049] For each valid sampling point, the device reads the terrain elevation of the corresponding raster from the DEM data. And combined with the height of the observation point and the horizontal distance between sampling points and observation points Calculate the elevation angle or slope parameter corresponding to this sampling point. Simultaneously, the maximum elevation angle extreme value parameter in that direction is continuously maintained during the traversal of the corresponding ray. When the current sampling point satisfies When the sampling point is in a visible state, the terrain unit corresponding to the sampling point is determined to be visible, and the maximum elevation angle extreme value parameter is updated; when the current sampling point satisfies When the sampling point is located, it is determined that the terrain unit corresponding to the sampling point is occluded by the terrain in front; based on the determination result, the occlusion status is written into the initial occlusion result matrix. This completes the visibility determination and initial occlusion result generation in step 4;

[0050] After completing the parallel traversal of all rays, the initial occlusion result matrix is ​​further refined. Post-processing is performed to reduce the impact of isolated misjudgment points and local boundary burrs on the continuity of the occlusion result; morphological erosion is performed on the initial occlusion result matrix to obtain the final occlusion result matrix. This completes the occlusion result optimization process in step 5.

[0051] Finally, the final occlusion result matrix obtained from the device side is... Or the initial occlusion result matrix Returning to the host, it serves as the input for subsequent terrain visibility analysis, radar coverage analysis, terrain constraint modeling, or clutter simulation processing. Areas identified as occluded are shielded or attenuated to obtain subsequent analysis results that take into account terrain occlusion effects, thus completing the result output and application in step 6.

[0052] Furthermore, in step 4, when traversing the terrain units in this direction in order from near to far, the maximum number of steps is set to:

[0053] ;

[0054] in Where C is the diagonal length of the DEM region, and C is the safety redundancy factor. Indicates DEM raster in Width in the direction, Indicates DEM raster in Height in direction Indicates the DEM raster resolution. Indicates the ray step length. This represents the safety redundancy coefficient, used to compensate for boundary rounding errors and step margins caused by the initial offset of observation points.

[0055] Furthermore, a terrain occlusion fast calculation system includes a processor, a memory, and a parallel computing device; the memory stores program instructions that execute the above-described method on the processor or the parallel computing device.

[0056] Compared with the prior art, the present invention has at least the following beneficial effects:

[0057] (1) By “computational region corner domain clipping + ray set construction”, the occlusion calculation is limited to the orientation range related to the simulation area, reducing the amount of calculation in invalid directions;

[0058] (2) By using the ray traversal criterion of "elevation angle extreme value maintenance + threshold comparison", the occlusion determination can be updated quickly, reducing the computational complexity and making it suitable for parallel execution to improve the computational efficiency under large-scale DEM conditions.

[0059] (3) By performing erosion calculations on the initial shading results, the local discontinuities, isolated points and holes caused by discrete rays and grid sampling are effectively eliminated, making the shading area expression more continuous and stable, which is more in line with the engineering application requirements for result stability.

[0060] (4) The output occlusion distribution results can be directly integrated into subsequent spatial analysis or simulation modeling, improving the overall efficiency of analysis / simulation and engineering applicability. Attached Figure Description

[0061] Figure 1 This is a flowchart illustrating the overall process of the method of the present invention.

[0062] Figure 2 A schematic diagram showing the determination of observation points, computational regions, and angular domains.

[0063] Figure 3 This is a schematic diagram of ray angle discretization and step sampling.

[0064] Figure 4 This is a schematic diagram of the shading determination principle based on the maintenance of elevation angle extreme values. Detailed Implementation

[0065] The present invention will be further described below with reference to specific embodiments. Those skilled in the art should understand that equivalent substitutions can be made to the parameters, data structures, and implementation platforms in the embodiments without departing from the concept of the present invention, and all such substitutions fall within the scope of the present invention.

[0066] Step 1: DEM data preparation and output initialization;

[0067] Acquire digital elevation model (DEM) data and organize it into a two-dimensional raster array. ,in , This represents the corresponding raster terrain elevation. In this embodiment, the raster resolution of the DEM is r=1 m, meaning that the distance between the centers of adjacent raster cells in the east-west and north-south directions is 1 m.

[0068] Establish a binary occlusion result matrix with the same resolution as the DEM. This is used to store the initial occlusion determination results, and the following encoding conventions apply:

[0069] ;

[0070] In this embodiment, initially Initialize to all 1s (default occlusion), that is:

[0071] ;

[0072] This initialization strategy corresponds to the "conservative occlusion" principle: a grid point is marked as unoccluded only when it is confirmed that the grid point is visible in a certain direction.

[0073] Step 2: Determining the observation point and calculation area

[0074] Determine the observation point (starting point of the line of sight): P0 = (x0, y0, h0), where (x0, y0) are the projected coordinates of the observation point in the DEM plane coordinate system, and h0 is the height of the observation point. The observation point can be located inside or outside the area to be calculated.

[0075] Define the terrain shading calculation area Ω. Ω is a rectangular area, and its four corner points in the DEM coordinate system are:

[0076] ;

[0077] in ;

[0078] Step 3: Calculate and determine the angular domain of the calculation area.

[0079] To reduce computational overhead in irrelevant directions, the azimuth angles relative to the observation point are calculated based on the four corner points of the rectangular region:

[0080] ;

[0081] For the set of angles { } Perform angle domain processing to obtain the minimum continuous azimuth angle interval covering the calculation area. ].like Figure 2 As shown, observation point It can be located inside or outside the calculation region Ω, through the four corner points of the calculation region. Compared to The azimuth angle determines the angular domain range. During subsequent ray construction, rays are only deployed within this angular domain to reduce the scale of parallel computing tasks and improve efficiency. Note: When the angular domain crosses the −π and π boundaries, the “minimum continuous interval” is determined using an angular expansion method to ensure the continuity of the angular domain.

[0082] Step 4: Boundary parameterization sampling and ray endpoint determination

[0083] Within the azimuth range of the calculation region Ω [ To ensure sufficient ray coverage while maintaining computational efficiency, this implementation establishes a correspondence between the number of rays N and the perimeter of the computational region boundary.

[0084] Given a DEM raster resolution of r=1m, the calculation region Ω is rectangular, with a width in pixels of:

[0085] ;

[0086] The height in pixels is:

[0087] ;

[0088] The discrete perimeter pixel length of this region is:

[0089] ;

[0090] In this embodiment, the number of rays is set to twice the circumference:

[0091] ;

[0092] Therefore, the number of rays increases linearly with the scale of the computational region, which can maintain a relatively stable angular sampling density when the region expands, thereby improving the continuity of occlusion determination.

[0093] Furthermore, the azimuth interval [ The angular resolution is defined as follows: The rays are uniformly discretized according to their number. ;

[0094] The azimuth angle of the i-th ray is:

[0095] ;

[0096] In other embodiments, the number of rays N is set to K*P, where K is a proportionality coefficient greater than zero.

[0097] Step 5: Calculate the ray direction vector and perform step traversal;

[0098] For the i-th ray obtained in step 4, its azimuth angle Calculate its unit direction vector in the plane:

[0099] ;

[0100] Set the step length to In this example, we take: Considering the DEM resolution r=1 m, this setting is equivalent to "sampling once every 0.5 m", which enhances the stability of the determination at the occlusion boundary.

[0101] The step increment for the i-th ray is:

[0102] ;

[0103] Starting from the observation point, the sampling point is progressively advanced along the ray. The sampling point after the m-th step is represented as:

[0104] ;

[0105] in This represents the horizontal distance from the sampling point to the observation point. The sampling point is mapped to the DEM raster index (rounding is used in this embodiment): The ray traversal terminates when the index exceeds the DEM range. The traversal terminates upon reaching the observation point. To avoid numerical instability near the observation point, points that are extremely close can be skipped. When, skip this traversal, where 0.1m can be taken. For example... Figure 3 As shown, the i-th ray originates from the azimuth angle. It is determined that its direction vector is ( ), and adopt step size Sampling is performed stepwise along the ray direction and mapped to the DEM raster index (ix, iy).

[0106] To prevent a thread from entering an infinite loop, one implementation sets a maximum step count limit as follows:

[0107] ;

[0108] in Where C is the diagonal length (meters) of the DEM region, and C is the safety redundancy factor. In one implementation, C is 100, used for coverage boundary rounding and starting point offset.

[0109] Step 6: Shading Determination for Extreme Value Maintenance of Slope (Elevation Angle Parameter)

[0110] For the terrain elevation H(ix,iy) corresponding to the ray sampling point, the slope / elevation coefficient is defined as:

[0111] ;

[0112] And maintain the extreme value for each ray. Initialize it to a sufficiently small value.

[0113] When traversing along the ray from near to far, the following rules apply to determination and updating:

[0114] like This indicates the emergence of a new "maximum elevation angle" point, which is unobstructed and visible. Update: and will Set it to zero; otherwise, continue iterating. Using the above method of "extreme value maintenance + threshold comparison," each ray only requires one forward traversal to complete the occlusion determination.

[0115] Step 7: Parallel Implementation

[0116] Since the rays are independent of each other, this embodiment uses a GPU parallel approach: each ray is assigned to a parallel execution unit (thread / computation unit) for independent execution; DEM data is preloaded to device-side storage, and the occlusion matrix... It is also generated on the device side, and after calculation, it can be sent back to the host for subsequent simulation. Note: The above parallel implementation can be implemented through CUDA or other parallel computing frameworks; this invention is not limited to a specific framework. Figure 4 As shown, the elevation angle parameters for terrain points are calculated along the same ray direction. and maintain extreme values When a new maximum elevation angle point appears, it is determined to be unobstructed and the status is updated. After that, it satisfies The point is judged to be occluded.

[0117] In one specific implementation, GPU-side data organization includes: reading the DEM elevation raster on the host side and storing it linearly; calling cudaMalloc to allocate storage space for the DEM array and occlusion output matrix on the device side; and copying the DEM data to the device side using cudaMemcpy. The occlusion output matrix is ​​initialized to a preset value on the device side to reduce host-side involvement.

[0118] In one kernel function mapping method, the thread index `ray_id` corresponds one-to-one with the ray number. The thread calculates the azimuth angle θ and direction vector (cosθ, sinθ) based on `ray_id`, performs step sampling with `step_len`, and maintains the extreme parameters `max_slope` (or `max_angle`) of the elevation / slope direction for that ray, thereby completing the visibility / occlusion determination and writing the occlusion result to the device. The thread block size can be 128 or 256, and the mesh size is configured by rounding up `num_rays`. Each thread terminates traversal when it traverses outside the DEM area or reaches the maximum traversal step count.

[0119] To further improve throughput, the occlusion matrix preferably uses an 8-bit or 16-bit compact data type to store the occlusion markers; the post-erosion processing can be executed in parallel on the GPU as an independent kernel function, and after processing, the final occlusion matrix is ​​sent back to the host for display, archiving, or as input for subsequent spatial analysis processes.

[0120] Step 8: Image Post-processing - Erosion to Expand Unmasked Areas

[0121] To correct the gaps and discontinuities in the occlusion boundary caused by ray discretization and grid sampling, this embodiment modifies the initial occlusion matrix. Perform morphological erosion to expand the unmasked area, resulting in the final occlusion matrix: ,in For structural elements, considering a DEM resolution of 1m, a 3x3 square structural element is preferred;

[0122] Step 9: Output and Application of Results

[0123] The final occlusion matrix M1(x,y) is output as an occlusion distribution result, which can be used as an input constraint in processes such as view analysis, coverage analysis, or simulation environment modeling. For example, it can mask occlusion grids that satisfy M1(x,y)=1, or reduce their contribution according to a preset weight coefficient, thereby reflecting the terrain occlusion effect and improving the rationality of the analysis.

Claims

1. A GPU-based method for fast terrain occlusion calculation, comprising the following steps: Step 1: Acquire digital elevation model (DEM) data and organize the DEM data into a two-dimensional raster format; establish a shading result output matrix with the same resolution as the DEM to record the shading status of each raster unit within the area to be calculated; Step 2: Determine the observation point And the terrain-shading area Ω to be calculated; where, The observation point Located inside or outside the shielded area; Step 3: Define the rectangular computational region and construct the ray set; Let the four corner points of the rectangular calculation region be... Calculate the azimuth angle of each corner point relative to the observation point. : ； right Angle domain processing is performed to obtain the minimum continuous azimuth angle interval covering the rectangular calculation area. ], and generate rays within that azimuth angle range according to a preset number of rays or angular resolution, where the number of rays is N: ； For angular resolution: ; The preset number of rays is used to ensure that multiple boundary sampling points are uniformly selected on the boundary of the computational region to guarantee the ray coverage area. }, from the observation point A ray set is formed pointing towards each boundary sampling point; Step 4: Ray traversal and elevation angle extreme value occlusion determination; For each ray direction, traverse the terrain units along that direction in order from nearest to farthest; for any terrain unit on the ray... Define horizontal distance for: ； Define elevation angle parameters : ； Maintain the current extreme elevation angle for each ray. During the traversal, occlusion is determined according to the following rules, and the initial occlusion result is output. : like The terrain unit is determined to be unobstructed, and updates are performed. ; like The terrain unit is determined to be occluded. The initial occlusion result For two-dimensional raster results with the same resolution as the DEM data: ； Since the traversal and extremum maintenance of different rays are independent of each other, the occlusion calculation is performed on multiple rays in parallel by the GPU parallel computing unit. Step 5: Apply the initial occlusion results obtained in Step 4 Image post-processing, including at least morphological erosion, is performed to reduce local discontinuities, isolated points, and holes caused by discrete ray and raster sampling, thereby obtaining stable occlusion results. ; Step 6: Use the occlusion result M1(x,y) obtained from post-processing as the input constraint for subsequent spatial analysis or simulation modeling: Mask the occluded terrain units in the coverage analysis / visibility analysis or attenuate them according to the preset weight coefficient to obtain the analysis results that take into account the terrain occlusion effect.

2. The GPU-based fast terrain occlusion calculation method as described in claim 1, characterized in that, In step 3, the direction vector of the ray is determined by the following formula: ； Indicates the first Ray in a plane coordinate system Components in direction, Indicates the first A ray in a plane coordinate system The component in the direction, and ( ( ) represents the corresponding azimuth angle The unit direction vector in the plane.

3. The GPU-based fast terrain occlusion calculation method as described in claim 1, characterized in that, In step 3, the ray azimuth angle covers the entire angular domain [0, 2π) and satisfies: 。 4. The GPU-based fast terrain occlusion calculation method as described in claim 1, characterized in that, The method for setting the number of rays N in step 3 includes: Method 1: The number of rays N is set proportionally to the discrete perimeter P of the region to be calculated, satisfying N=kP, where k>0 is the scaling factor. Method 2: Determined by angular resolution. If the angular resolution is determined... So there are a total of / A ray.

5. The GPU-based fast terrain occlusion calculation method as described in claim 1, characterized in that, In step 4, the mapping from the sampling points to the DEM data raster index uses a rounding method: ； in, This represents the DEM raster column index after the sampling points are mapped. This represents the DEM raster row index after the sampling points are mapped; This indicates the floor function.

6. The GPU-based fast terrain occlusion calculation method as described in claim 1, characterized in that, The morphological etching process in step 5 satisfies the following: ； Where S represents the structural element used in the morphological processing, which defines the shape and scale of the local neighborhood of the erosion operation, preferably... Square structural element; This represents a morphological erosion operation used to reduce isolated points, local holes, and boundary discontinuities in the initial occlusion result, thereby obtaining an occlusion result with better spatial continuity. .

7. The GPU-based fast terrain occlusion calculation method as described in claim 1, characterized in that, In step 4, when traversing the terrain units in this direction in order from near to far, the maximum number of steps is set to: ； in Where C is the diagonal length of the DEM region, and C is the safety redundancy factor. Indicates DEM raster in Width in the direction, Indicates DEM raster in Height in direction Indicates the DEM raster resolution. Indicates the ray step length. This represents the safety redundancy coefficient, used to compensate for boundary rounding errors and step margins caused by the initial offset of observation points.

8. The GPU-based fast terrain occlusion calculation method as described in claim 1, characterized in that, Steps 1 to 6 can be executed in sequence according to the following process: First, the host computer reads the digital elevation model (DEM) raster data and linearizes it using either row-major or column-major order to obtain a one-dimensional DEM array suitable for device access. Simultaneously, observation points are established on the host computer. Region to be calculated DEM dimensional parameters and Raster resolution Total number of rays Step length and safety redundancy coefficient The setup, where the process corresponds to the input preparation and parameter initialization in steps 1, 2 and 3; Subsequently, a DEM array and a masking output matrix are requested on the device side. The storage space copies the linearized DEM data from the host to the device; the masking output matrix Initialize the values ​​on the device side to preset values, preferably to all 1s, to implement a conservative decision-making strategy for default occlusion. This process corresponds to the establishment and initialization of the occlusion result matrix in step 1. After data transfer and video memory allocation are completed, the host side operates according to the corner range determined in step 3. and total number of rays Configure and start the parallel computing task on the device side to establish a correspondence between each parallel execution unit and each traversable ray; wherein, each parallel execution unit is responsible for sampling traversal, maintenance of elevation angle extreme values ​​and writing of occlusion results in the direction of a ray, thereby completing the ray set construction in step 3 and the parallel task allocation in step 4. During parallel computing on the device side, each parallel execution unit executes its corresponding ray number. Calculate the azimuth angle of the ray. And further obtain the corresponding planar unit direction vector. ; by step length By advancing step by step along the ray direction, the coordinates of each sampling point are obtained in sequence. And map the sampling points to the DEM raster index. When the sampling points exceed the valid boundaries of the DEM or the current traversal step count reaches the maximum step count limit; When the sampling traversal corresponding to the ray is terminated, the ray step sampling and boundary termination control in step 4 are completed. For each valid sampling point, the device reads the terrain elevation of the corresponding raster from the DEM data. And combined with the height of the observation point and the horizontal distance between sampling points and observation points Calculate the elevation angle or slope parameter corresponding to this sampling point. Simultaneously, the maximum elevation angle extreme value parameter in that direction is continuously maintained during the traversal of the corresponding ray. When the current sampling point satisfies When the sampling point is in a visible state, the terrain unit corresponding to the sampling point is determined to be visible, and the maximum elevation angle extreme value parameter is updated; when the current sampling point satisfies When the sampling point is located, it is determined that the terrain unit corresponding to the sampling point is occluded by the terrain in front; based on the determination result, the occlusion status is written into the initial occlusion result matrix. This completes the visibility determination and initial occlusion result generation in step 4; After completing the parallel traversal of all rays, the initial occlusion result matrix is ​​further refined. Post-processing is performed to reduce the impact of isolated misjudgment points and local boundary burrs on the continuity of the occlusion result; morphological erosion is performed on the initial occlusion result matrix to obtain the final occlusion result matrix. This completes the occlusion result optimization process in step 5. Finally, the final occlusion result matrix obtained from the device side is... Or the initial occlusion result matrix Returning to the host, it serves as the input for subsequent terrain visibility analysis, radar coverage analysis, terrain constraint modeling, or clutter simulation processing. Areas identified as occluded are shielded or attenuated to obtain subsequent analysis results that take into account terrain occlusion effects, thus completing the result output and application in step 6.

9. A computing system for performing the GPU-based fast terrain occlusion calculation method as described in claim 1, comprising a processor, a memory, and a parallel computing device; wherein the memory stores program instructions that execute the method described above on the processor or the parallel computing device.

Images