CN111325791B - OSP space reference line through-view domain analysis method based on regular grid DEM - Google Patents
OSP space reference line through-view domain analysis method based on regular grid DEM Download PDFInfo
- Publication number
- CN111325791B CN111325791B CN202010054396.3A CN202010054396A CN111325791B CN 111325791 B CN111325791 B CN 111325791B CN 202010054396 A CN202010054396 A CN 202010054396A CN 111325791 B CN111325791 B CN 111325791B
- Authority
- CN
- China
- Prior art keywords
- point
- reference line
- osp
- grid
- coordinate system
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/70—Determining position or orientation of objects or cameras
- G06T7/73—Determining position or orientation of objects or cameras using feature-based methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
- G06T17/05—Geographic models
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10016—Video; Image sequence
- G06T2207/10021—Stereoscopic video; Stereoscopic image sequence
Abstract
The invention discloses an OSP space reference line through-view domain analysis method based on a regular grid DEM, which comprises the following steps: partitioning the through-view area analysis area based on the regular grid DEM; constructing a corresponding station center coordinate system in each partition, and constructing a corresponding OSP space coordinate system based on the constructed station center coordinate system; converting grid points corresponding to the terrain lines in each partition into an OSP space coordinate system; and judging the visibility of each partition in an OSP space coordinate system, and combining the visibility judgment results of each partition to obtain the visibility judgment result of the visibility domain analysis area. When the through-view judging area is partitioned according to the grid lines parallel to the horizontal/longitudinal grids by the observation point, the precision of the method is the same as that of an R3 algorithm, and the speed is the same as that of a reference plane algorithm; when the perspective judgment area is partitioned according to diagonal lines by the observation point, the accuracy of the method is far higher than that of an XDraw and reference plane algorithm, and the speed obviously exceeds that of the reference plane algorithm within hundreds of meters of height in a mountain area.
Description
Technical Field
The invention belongs to the technical field of terrain visibility analysis, and particularly relates to an OSP space reference line visibility domain analysis method based on a regular grid DEM.
Background
The main function of the perspective area analysis is to calculate the range of a certain distance from a standing point and the perspective area of an observer based on DEM data. The through-view domain analysis is widely applied to the fields of geographic information service, radar detection, communication infrastructure, military planning and the like, can be used for communication base station site selection, transmitting station site selection of a television station, radar station site selection, road planning and the like in the civil field, and can be used for setting observation stations, laying and erecting communication lines, low-altitude flight path planning, evasion and the like in military matters such as laying battlefields (artillery battle fields and electronic countermeasure battle fields).
For general visual field analysis, a great deal of Research has been done by scholars both at home and abroad, franklin et al propose the famous R2, R3, XDraw algorithm (Franklin R. Geological Algorithms for training of Air defenses media Batteries [ R ]. Columbus Division: A Research Project for Batteries [ c ], 1994.), wu-Brillian (Wu-Brillian. Visual field algorithm [ J ]. Mapping information and engineering based on reference planes, 2001 (1): 19-21, 25.) et al propose the "reference plane algorithm", which has far-reaching impact, on the basis of which many scholars have been researched and improved by the scholars so far. The R3 algorithm considers the influence of each grid line passing by sight on the sight, the calculation result is accurate, but the calculation amount is extremely large, so that the calculation cannot be applied to actual operation; the R2, XDraw and reference plane algorithm is an approximate fast algorithm, the calculation speed can meet the requirement of actual use, but errors are necessarily generated because the influence of each grid line through which the sight line passes on the sight line is not considered. For a long time, most scholars think that the influence of algorithm precision on application is not large due to low data precision of the traditional DEM, however, through research, the problem that the application of a calculation result is influenced due to the fact that the error ratio of approximation algorithms such as XDraw is low but the distribution is uneven, and the regional perspective judgment is possibly wrong is caused, and along with the development of the technology, the data quality of the DEM is greatly improved, the traditional approximation algorithm cannot reflect the data advantage of the modern DEM, and the problem is not substantially improved all the time. Some researchers have tried to improve the accuracy problem, such as Wu (Wu H, pan M, yao L, et al. A Partition-Based sequential Algorithm for Generating Viewswed on Massive Dems [ J ]. International Journal of geographic Information Science,2007,21 (9): 955-964.) to propose the drawback of the Reference plane Algorithm for the computation of topography such as the peak and valley, zhi (Zhi Y, wu L, sui Z, et al. An Improved Algorithm for Computing Viewswed Based Reference places [ C ]// Ieee 11 th) to improve the accuracy of the Algorithm by a factor of 50%, but the conventional Algorithm is not designed to be as fast as possible, but the speed is not as fast as possible, and the error is not as fast as the following is not easy to be reduced by a factor of 50%.
In fact, the algorithm is still to be further improved in terms of the calculation speed, and most of the current basic schemes for improving the algorithm speed of scholars shorten the calculation time through parallelization, but the scheme needs online support of large-scale hardware, and is difficult to obtain a good effect under the condition of portable hardware equipment without a network environment, and particularly difficult to give speed advantages under field conditions such as exploration, archaeology, military application and the like, so that the speed bottleneck of view field analysis is solved from the essence of the calculation amount through researching a new basic algorithm.
Disclosure of Invention
The invention provides an OSP space reference line through-view field analysis method based on a regular grid DEM (digital elevation model) aiming at the problem that the existing through-view field analysis method cannot simultaneously consider speed and precision.
In order to achieve the purpose, the invention adopts the following technical scheme:
an OSP spatial reference line through-view domain analysis method based on a regular grid DEM comprises the following steps:
step 1: partitioning the through-view area analysis area based on the regular grid DEM;
step 2: constructing a corresponding station center coordinate system in each partition, and constructing a corresponding OSP space coordinate system based on the constructed station center coordinate system;
and 3, step 3: converting grid points corresponding to the terrain lines in each partition into an OSP space coordinate system;
and 4, step 4: and judging the visibility of each partition in an OSP space coordinate system, and combining the visibility judgment results of each partition so as to obtain the visibility judgment result of the visibility domain analysis area.
Further, the step 1 comprises:
and dividing the longitudinal topographic line of the through-vision area analysis area into a 1 area and a 2 area and dividing the transverse topographic line of the through-vision area analysis area into a 3 area and a 4 area on the regular grid DEM by taking the observation point as a boundary.
Further, the step 1 comprises:
and dividing the analysis area of the visibility region into a region 1, a region 2, a region 3 and a region 4 on the regular grid DEM by taking a diagonal line of the analysis area of the visibility region passing through an observation point as a boundary, wherein the region 1 and the region 2 are transversely centrosymmetric, and the region 3 and the region 4 are longitudinally centrosymmetric.
Further, the step 2 comprises:
step 2.1: constructing a first station center coordinate system in the area 1: establishing a right-hand space rectangular coordinate system by taking the viewpoint of an observer as an original point, taking the north-offset direction of the longitudinal grid line of the DEM as a longitudinal axis Y, taking the east-offset direction of the transverse grid line of the DEM as a transverse axis X and taking the normal vector direction of the reference ellipsoid where the observer is located as a vertical axis Z;
step 2.2: taking the Z axis of the first station center coordinate system as a rotating axis, respectively rotating by 180 degrees, 270 degrees and 90 degrees clockwise to obtain three coordinate systems which are respectively defined as a second station center coordinate system, a third station center coordinate system and a fourth station center coordinate system, and respectively establishing the second station center coordinate system, the third station center coordinate system and the fourth station center coordinate system in 2, 3 and 4 areas by taking an observation point as a center;
step 2.3: constructing an OSP space coordinate system by taking the direction k as a horizontal axis, the proximity p as a vertical axis and the axis slope g as a vertical axis based on each constructed station center coordinate system; the direction k is a tangent value of a survey point in a station center coordinate system relative to a horizontal axis direction angle alpha; the proximity p is the reciprocal of the value of the horizontal axis x of the investigation point in the station center coordinate system; the axial slope g is the product of the relative elevation z of the point under investigation and the proximity p.
Further, in the step 4, the visibility of each partition is determined as follows:
step 4.1: after grid points corresponding to each topographic line in each partition are respectively converted into respective OSP space coordinate systems from respective station center coordinate systems, initial reference lines are formed by the grid points with the maximum proximity according to the sequence of direction values from small to large, the ending direction k and the parameter a of each section of OSP line segment of the initial reference lines are solved, and the initial point axis slope g of the initial reference lines is solved 0 Recording grid points on the initial reference line as a through view;
and 4.2: moving to a next proximity plane according to the sequence of the proximity of grid points in an OSP space coordinate system from large to small;
step 4.3: projecting the current reference line onto the next proximity plane in step 4.2;
step 4.4: finding out a first grid point on a current proximity plane and an OSP line segment L of a reference line containing the direction of the first grid point in a sequence from small direction to large direction, and solving an axial slope difference delta G between the first point and the L 0 By Δ G 0 Judging the visibility of the first grid point and recording the visibility to a corresponding visibility matrix;
step 4.5: if the first grid point is invisible, executing step 4.7;
step 4.6: if the first grid point is in full view, the first grid point and the grid point which is nearest outside the investigation region and has the same proximity are used for forming an OSP line segment, the ending direction k and the parameter a of the OSP line segment are calculated and used as the first segment of the temporary reference line segment, and g is updated 0 The slope of the starting axis of the first section of the current temporary reference line segment;
step 4.7: recording a next grid point as a current calculation point, checking whether the direction of the current point is greater than the tail direction value of the L, if so, recording the next section connected with the L on the reference line as a current reference line segment L until the direction of the current point is less than the tail direction value of the L;
step 4.8: calculating the axial slope difference delta G between the current grid point and the current reference line segment L by an axial slope difference accumulation formula i By Δ G i Judging the visibility and recording;
step 4.9: if the current grid point and the previous grid point are all in sight, the ending direction k and the parameter a of an OSP line segment formed by the two points are solved, and the OSP line segment is added into a temporary reference line segment;
step 4.10: if the current grid point is not visible and the previous grid point is visible, firstly, the parameter a of the OSP line segment formed by the current point and the previous point is solved, and the intersection point direction k between the OSP line segment and the reference line is solved; then, recording the ending direction of the line segment as k and the parameter as a, and adding a temporary reference line segment; finally, updating the line segment between the original reference line and the direction from the intersection point to the current intersection point by using the temporary reference line segment, if the original reference line does not have the previous intersection point, updating the line segment between the original reference line starting point and the current intersection point by using the temporary reference line segment, and emptying the temporary reference line;
step 4.11: if the current grid point and the previous grid point are invisible, no operation is needed;
step 4.12: if the current grid point is visible and the previous grid point is invisible, the ending direction of an OSP line segment formed by the two points is determined to be k and a parameter a, the intersection point direction k' between the OSP line segment and a reference line is determined, and a line segment from the intersection point to the current grid point forms a first section of a temporary reference line segment;
step 4.13: if the current proximity still has a point of investigation, executing step 4.7;
step 4.14: if the last investigation point of the current proximity is visible and the temporary reference line is not empty, updating the line segment from the first intersection point to the last intersection point on the original reference line by using the temporary reference line segment, if the last intersection point does not exist, updating the line segment from the starting point to the last intersection point of the original reference line, and finally emptying the temporary reference line;
step 4.15: taking the updated reference line as the current reference line;
step 4.16: if there are proximity planes yet to be investigated, step 4.2 is performed until the investigation is completed.
Further, the step 4 further includes:
when the observation points are used as boundaries for partitioning, the visibility matrix of the visibility region analysis area is calculated according to the following mode:
V＝(V 1 +V 2 )·(V 3 +V 4 ) (13)
wherein, V 1 、V 2 、V 3 、V 4 The through-view matrixes are respectively 1,2, 3 and 4 areas;
when the element in the V is 0, grid points with the same row number and the same column number as the element in the visibility analysis area are invisible; and when the element in the V is 1, the grid point with the same row number and column number as the element in the visibility analysis area is subjected to visibility.
Further, the step 4 further includes:
when the through-view area analysis area is partitioned by taking a diagonal line of the through-view area analysis area passing through the observation point as a boundary, calculating a through-view matrix of the through-view area analysis area as follows:
V＝V 1 +V 2 +V 3 +V 4 (14)
wherein, V 1 、V 2 、V 3 、V 4 The through-view matrixes are respectively 1,2, 3 and 4 areas;
when the element in the V is 0, grid points with the same row number and the same column number as the element in the visibility analysis area are invisible; and when the element in the V is 1, the grid point with the same row number and column number as the element in the visibility analysis area is subjected to visibility.
Further, the parameter a is obtained as follows:
a＝u(z B -z A )
wherein u is the reciprocal of the grid width of DEM, and z A 、z B The vertical axis values of the station center coordinates of any two adjacent grid points A and B with the same proximity are respectively, and the direction of B is greater than that of A.
Further, the shaft slope difference accumulation formula is:
ΔG i ＝ΔG i-1 +(a Ci -a Ri )(k i -k i-1 ) (11)
wherein, Δ G i Is the axial slope difference, Δ G, between the current grid point and the current reference line segment L i-1 Is the difference of the axis slope of the last segment end and the last reference segment, a Ci And a Ri The slopes of the line segments on the reference line and the investigation line in the current segmentRate, k i 、k i-1 Respectively in the tail direction of the current section and the last section;
further, in step 4.12, the direction k' of the intersection between the OSP line segment and the reference line is determined by the following formula:
k′＝k i -ΔG i /(a Ci -a Ri ) (12)。
compared with the prior art, the invention has the following beneficial effects:
1. firstly, establishing an OSP space coordinate system, then adopting a partition calculation mode, respectively converting grid points into an OSP space in each area, judging the visibility condition of each grid point in the process of constructing and updating a reference line, wherein the influence of each DEM grid line on visibility is transmitted from the center to the periphery in a function mode in the calculation process, sampling is not needed, interpolation is not needed, redundant calculation is not needed, and the method is suitable for large-area long-distance visibility analysis and also suitable for serving as a basic algorithm for other parallel terrain visibility analysis and calculation;
2. a large number of experiments prove that when the perspective judgment area is partitioned by an observation point in parallel to a transverse/longitudinal grid line, the precision of the method is completely the same as that of an R3 algorithm, but the calculation speed is far higher than that of the R3 algorithm and is in the same order as that of an XDraw and a reference plane algorithm, although the speed is lower than that of the XDraw and the reference plane algorithm, the advantages of a modern DEM product can be embodied, the calculation result is reliable, continuous perspective judgment errors cannot occur, and meanwhile, the requirement of large-area calculation can be met in speed; when the partition is carried out by taking the diagonal line of the through-view area analysis region passing through the observation point as a boundary, the algorithm precision is far higher than that of other quick approximation algorithms, the precision can be improved by more than 5 times compared with the precision of the traditional approximation algorithms (and parallel algorithms based on the traditional approximation algorithms) such as XDraw and reference plane algorithms, the speed is far higher than that of the R3 algorithm, the algorithm is higher than that of the XDraw algorithm in most cases, and the algorithm is obviously higher than that of the reference plane algorithm in the range of hundreds of meters away from the ground in mountainous areas. The invention solves the problem that the landform-based through-vision domain analysis can not realize both speed and precision for a long time.
Drawings
FIG. 1 is a basic flow chart of an OSP spatial reference line through-the-horizon analysis method based on a regular grid DEM according to an embodiment of the present invention;
FIG. 2 is an exemplary diagram of a station center coordinate system of an OSP spatial reference line visibility domain analysis method based on a regular grid DEM according to an embodiment of the present invention;
fig. 3 is a schematic view of a partition of an OSP spatial reference line through-view domain analysis method based on a regular grid DEM according to an embodiment of the present invention;
fig. 4 is a schematic diagram of grid points of a 1-region DEM in an OSP spatial reference line through-the-horizon analysis method based on a regular grid DEM according to an embodiment of the present invention;
fig. 5 is a schematic diagram of a reference line updating process of an OSP spatial reference line through-view domain analysis method based on a regular grid DEM according to an embodiment of the present invention;
FIG. 6 is a flow chart of a sectional visibility decision of an OSP spatial reference line visibility domain analysis method based on a regular grid DEM according to an embodiment of the present invention;
fig. 7 is a schematic data structure diagram of a reference line of an OSP spatial reference line through-view domain analysis method based on a regular grid DEM according to an embodiment of the present invention;
FIG. 8 is a second schematic diagram of a data structure of a reference line of an OSP spatial reference line through-view domain analysis method based on a regular grid DEM according to an embodiment of the present invention;
FIG. 9 is an illustration of the difference of the slope of the axes of an OSP spatial reference line through-view domain analysis method based on a regular grid DEM according to an embodiment of the present invention;
FIG. 10 shows three terrain DEMs for experiments according to an OSP spatial reference line visibility domain analysis method based on a regular grid DEM in an embodiment of the present invention;
fig. 11 is an error distribution situation of each method on different terrains/different heights with reference to an R3 algorithm of an OSP spatial reference line through-view domain analysis method based on a regular grid DEM according to an embodiment of the present invention;
fig. 12 is a diagram illustrating a neighborhood of DEM grid points in an OSP spatial reference line visibility domain analysis method based on a regular grid DEM according to an embodiment of the present invention;
fig. 13 is a comparison graph of calculation results of the visibility areas of the OSP reference line algorithm and the XDraw algorithm of the OSP spatial reference line visibility area analysis method based on the regular grid DEM according to the embodiment of the present invention;
fig. 14 is a statistical chart of the frequency of neighborhood errors occurring in different terrains/different heights by using the R3 algorithm as a reference and the XDraw algorithm in the OSP space reference line through-view domain analysis method based on the regular grid DEM according to the embodiment of the present invention;
FIG. 15 is a statistical graph of the trend of the time consumption along with the calculated radius change of the R3 algorithm and the OSP reference line algorithm of the OSP spatial reference line through-view domain analysis method based on the regular grid DEM in the embodiment of the present invention;
FIG. 16 is a time-consuming trend statistical graph with the number of calculated points for each algorithm of an OSP spatial reference line through-view domain analysis method based on a regular grid DEM according to an embodiment of the present invention;
FIG. 17 is a schematic sectional view of another OSP spatial reference line through-view domain analysis method based on a regular grid DEM according to an embodiment of the present invention;
FIG. 18 is a schematic diagram of grid points of a 1-region DEM according to another OSP spatial reference line through-view domain analysis method based on a regular grid DEM in accordance with an embodiment of the present invention;
FIG. 19 is a diagram showing the error distribution of the OSP spatial reference line through-view domain analysis method based on the regular grid DEM according to the embodiment of the present invention, with reference to the R3 algorithm, on different terrains/different heights;
FIG. 20 is a graph showing a trend of time consumption with calculated radius change for another OSP spatial reference line through-view domain analysis method based on a regular grid DEM according to the present invention and an R3 algorithm;
FIG. 21 is a comparison graph of the time consumption of another OSP spatial reference line through-view domain analysis method based on a regular grid DEM according to the embodiment of the present invention and the reference plane algorithm in different terrains and different heights from the ground;
fig. 22 is a comparison graph of the time consumption of the OSP spatial reference line through-view domain analysis method based on the regular grid DEM according to the embodiment of the present invention and the XDraw algorithm in different terrains and different heights from the ground.
Detailed Description
The invention is further illustrated by the following examples in conjunction with the drawings and the accompanying drawings:
as shown in fig. 1, an OSP spatial reference line through-view domain analysis method based on a regular grid DEM includes:
step S101: partitioning the analysis area of the visibility region based on the regular grid DEM;
step S102: constructing a corresponding station center coordinate system in each partition, and constructing an OSP space coordinate system based on the constructed station center coordinate system;
step S103: converting grid points corresponding to the terrain lines in each partition into an OSP space coordinate system;
step S104: and judging the visibility of each partition in an OSP space coordinate system, and combining the visibility judgment results of each partition to obtain the visibility judgment result of the visibility domain analysis area.
Specifically, the step S101 includes:
for the convenience of calculation, the invention partitions the through-view area analysis area. The continuous broken line formed by connecting certain row/column DEM grid points can express the terrain of the corresponding geographic space of the row/column, and the continuous broken line is called as a transverse/longitudinal terrain line. As shown in fig. 3, the observation points are used as boundaries, the longitudinal topographic lines of the analysis area of the through-vision area in the regular grid DEM are divided into 1 area and 2 areas, and the transverse topographic lines are divided into 3 areas and 4 areas.
Specifically, the step S102 includes:
step S102.1: constructing a first station center Coordinate System (Topocentric Coordinate System) in region 1: establishing a right-hand space rectangular coordinate system by taking the viewpoint of an observer as an origin, taking the north offset direction of the longitudinal grid line of the DEM as a longitudinal axis Y (longitudinal axis), taking the east offset direction of the transverse grid line of the DEM as a transverse axis X (transverse axis), and taking the normal vector direction of a reference ellipsoid where the observer is located as a vertical axis Z (vertical axis); in the center of gravity coordinate system of fig. 2, the center of gravity of point a is (x, y, z), the horizontal distance from the origin is l, the viewing direction angle is α, and the elevation angle is
Step S102.2: taking the Z axis of the first station center coordinate system as a rotating axis, respectively rotating by 180 degrees, 270 degrees and 90 degrees clockwise to obtain three coordinate systems which are respectively defined as a second station center coordinate system, a third station center coordinate system and a fourth station center coordinate system, and respectively establishing the second station center coordinate system, the third station center coordinate system and the fourth station center coordinate system in 2, 3 and 4 areas by taking an observation point as a center;
step S102.3: constructing an OSP space coordinate system by taking the direction k as a horizontal axis, the proximity p as a vertical axis and the axis slope g as a vertical axis on the basis of each constructed station center coordinate system;
whether the target looks through or not is determined by the relative position relationship between the target and the shelter, and the following three shelter conditions are necessarily met simultaneously if the shelter shields the target on the terrain:
the barrier is in front of the target;
(ii) the obstruction and the target are in the same orientation of the observer;
(iii) the obstruction line of sight elevation is greater than the target line of sight elevation.
To quantify the three conditions, the invention defines 3 concepts of proximity, direction and axis slope respectively.
Proximity (denoted p) is defined as: the reciprocal of the value of x on the horizontal axis of the point of interest in the station center coordinate system is defined as formula (1). The meaning of proximity is: when the direction is the same, an object with a large proximity (close to the observer) may occlude an object with a small proximity (far from the observer), whereas if not, the proximity is a determination value of the "occlusion condition (i)".
p＝1/x (1)
The direction (denoted k) is defined as: the tangent value of the point relative to the horizontal axis direction angle alpha in the station center coordinate system is considered, and the formula is defined as formula (2). If the k values of the two points are different, representing that the directions of the observation points are different, and the objects with non-overlapping directions do not form an occlusion relation. The value k can be used as a determination value of the "occlusion condition (ii)".
k＝tan(α)＝yp (2)
The axis slope (labeled g) is defined as: the product of the relative elevation z and the proximity p of the survey point is defined as equation (3). According to this formula, the meaning of g in geographic space is:since the slope of a point (x, Y, z) after translation to (x, 0, z) in the Y-axis direction is considered, g is referred to as "axial slope" in the present invention. Slope (denoted as s) is the elevation of the point of interest in the centroid coordinate system
The larger the inclination, the larger the line of sight elevation angle, and the axis inclination and the inclination satisfy the relation (4), the larger the axis inclination, the larger the line of sight elevation angle of the point of interest, and the axis inclination can be used as the judgment value of the "occlusion condition (iii)", when the direction k is the same.
g＝zp (3)
The OSP spatial coordinate system is defined as: a spatial rectangular coordinate system is formed by taking the direction (k) as a horizontal axis, the proximity (p) as a vertical axis and the axis slope (g) as a vertical axis, and a survey point positioned in a geographic space is converted into the coordinate system to form a spatial relationship of direction-axis slope-proximity.
For any line in any centroid coordinate system that is not in the ZOY plane, there must be appropriate constants c, d, m, n, v, w that enable the line to be expressed as:
dividing the left and right of the equation (5) by x yields equation (6), and equation (6) is the expression after the conversion of the straight line to the OSP space. By combining with relevant mathematical principles, the linear relation among the direction, the axial slope and the proximity of a straight line is known, the property is called as the linear invariance of an OSP space, and the property lays a theoretical foundation for the through-the-horizon algorithm.
Specifically, the step S103 includes:
and converting any horizontal/vertical terrain line in each partition from the corresponding station center coordinate system to the corresponding OSP space coordinate system.
Specifically, the expression of the terrain line in the OSP space includes:
and (3) setting u as the reciprocal of the grid width of the DEM, wherein any small section on the topographic line is a straight line segment constructed by two adjacent grid points, and according to the OSP space property, an appropriate constant a and b is determined to enable the line segment to be consistent with the formula (7) in the OSP space, and the topographic line is positioned on an OSP plane with constant proximity after being converted into the OSP space and is also a continuous broken line. The invention records the line segment converted into OSP space as L (a, b | k) s ,k e ) Wherein k is s ,k e Left and right boundaries of the domain, a, b, k, respectively, for the line segment directions s ,k e A parameter called L. If the station center coordinates of two adjacent grid points are known to be A (x, y) respectively A ,z A )、B(x,y B ,z B ) Where the direction of B is greater than the direction of a, the value of the constants a, B can be expressed as (8):
g＝ak+b (7)
accordingly, it is easy to know that if the axis slope of the starting point of the first segment is known for a certain continuous OSP spatial broken line, only the slope parameter a and the ending direction k of each segment need to be recorded to calculate the arbitrary position configuration of the broken line. The present invention uses the structure shown in fig. 7 to record the polyline after the topographical lines represented by the grid lines are transformed into OSP space. The structure only records the a and ending direction k of each segment in the broken lines and the addresses of the previous segment and the next segment, because the invention needs to sequentially search the intersection point of the two broken lines along the broken lines, the structure can correctly express the structure of the broken lines and is convenient to search, and when the number of the broken lines needs to be modified, only the addresses of the front segment and the back segment of the first segment and the back segment of the last segment of the related segments need to be modified, as shown in figure 8, the sequence is not needed, the space complexity is lower, and the operation speed is higher.
Specifically, taking zone 1 as an example, as shown in FIG. 4, zone C 1 Constructing a terrain line by using the grid points, and converting the terrain line into an OSP space to form an initial reference line L R With C 2 Constructing a terrain line by using the grid points, and converting the terrain line into an OSP space to form a broken line L to be examined 2 From L to L R Translate to L along the P axis 2 In the plane of the plane, with L R And L 2 Dividing two topographic lines into several segments by using the intersection point as boundary, and combining L R And L 2 The section with higher middle shaft slope forms a new reference line L' R And L' R Will be further reacted with C 3 The terrain lines constructed by the column grid points are continuously updated through the steps until the area calculation is finished, and the reference line updating process is shown in fig. 5. Specifically, as shown in fig. 6, the visibility of the 1-region is determined as follows:
step S104.1: after grid points corresponding to all terrain lines in the region 1 are respectively converted into respective OSP space coordinate systems from respective station center coordinate systems, initial reference lines are formed by the grid points with the maximum proximity according to the sequence of direction values from small to large, the ending direction k and the parameter a of each section of OSP line segment of the initial reference lines are solved, and the initial point axis slope g of the initial reference lines is solved 0 Recording grid points on the initial reference line as a through view;
step S104.2: moving to the next proximity plane (the proximity plane is a special OSP plane, namely, the OSP plane with the proximity being a constant, and the OSP plane is a plane in the OSP space coordinate system) according to the sequence of the proximity of grid points in the OSP space coordinate system from large to small;
step S104.3: projecting the current reference line onto the next proximity plane in step S104.2;
step S104.4: finding out a first grid point on a current proximity plane and an OSP line segment L of a reference line containing the direction of the first grid point in a sequence from small direction to large direction, and solving an axial slope difference delta G between the first point and the L 0 By Δ G 0 Judging the visibility of the first grid point and recording the visibility to a corresponding visibility matrix;
specifically, the definition of the axial slope difference and the visibility decision function are:
as shown in FIG. 9, let L (a, b | k) 0 ,k 1 ) With known parameters, adjacentA proximity of p 0 The OSP coordinate of a certain point C in space is (k) C ,g C ,p C ) And has p C <p 0 ，k C ∈[k 0 ,k 1 ]. Note L (a, b | k) 0 ,k 1 ) And OSP plane k = k C Intersects the point D, and the coordinate of the intersection point is set as (k) C ,g D ,p 0 ). Define C Point relative to L (a, b | p) 0 ) The axial gradient difference of (1) is equation (9), corresponding to the black directed line segment in FIG. 9, and is easily known from the definition of the axial gradient, Δ G C >Point C is visible at 0, otherwise not visible. If used, r C If the visibility of the point C is marked, the value 0 represents the non-visibility, and the value 1 represents the visibility, the visibility discrimination function of the point C can be defined as a formula (10).
ΔG C ＝g C -g D (9)
Step S104.5: if the first grid point is invisible, step S104.7 is executed;
step S104.6: if the first grid point is in full view, the first grid point and the grid point which is nearest to the outside of the investigation region and has the same proximity form an OSP line segment, the ending direction k and the parameter a of the OSP line segment are solved and used as the first segment of the temporary reference line segment, and g is updated 0 The initial axis slope of the first section of the current temporary reference line segment;
step S104.7: recording a next grid point as a current calculation point, checking whether the direction of the current point is greater than the tail direction value of the L, if so, recording the next section connected with the L on the reference line as a current reference line segment L until the direction of the current point is less than the tail direction value of the L;
step S104.8: calculating the axial slope difference delta G between the current grid point and the current reference line segment L by an axial slope difference accumulation formula i By Δ G i Judging the visibility and recording;
specifically, the key to the reference line update is to find the direction value of the intersection point between the two terrain lines. Since the terrain lines are broken lines, the intersection points should be segmentedAnd (5) solving. Suppose L R Is a reference line, L C To explore the line (OSP broken line composed of all grid points in the current proximity), take L C And L R The direction values of inflection points of all the broken lines are segmentation points, two topographic lines are divided into a plurality of segments, and the tail direction of the i-th segment is set as k i In this paragraph L C And L R The slope of the line segment on each line segment is a Ci And a Ri The slope difference of the tail axis at the end of the period is recorded as Δ G i Due to the straight line characteristic, the following relation exists:
ΔG i ＝ΔG i-1 +(a Ci -a Ri )(k i -k i-1 ) (11)
equation (11) is called the cumulative equation of axial slope difference, and only L needs to be obtained by the equation C The direction k and the parameter a of each grid point can be cumulatively calculated to obtain L C Each grid point is relative to L R Differential axial slope of (d); judging the visibility of grid points through a formula (10), and recording the visibility into a visibility matrix which is in the same dimension as the whole calculation area; if Δ G i-1 And Δ G i Opposite positive and negative, then L R And L C The i-th segment has an intersection, and the direction of the intersection is k', which can be obtained from equation (12):
k′＝k i -ΔG i /(a Ci -a Ri ) (12)
a grid point must be judged as being in-view in both the 1/2 zone and the 3/4 zone, and the point can be finally determined as being in-view.
It should be noted that, in the process of executing step S104.7 and step S104.8, the formula (12) is used to perform cumulative calculation once each time the current reference line segment needs to be set as a next segment above the reference line or the axial slope difference of the current grid point needs to be obtained.
Step S104.9: if the current grid point and the previous grid point are all in sight, the ending direction k and the parameter a of an OSP line segment formed by the two points are solved, and the OSP line segment is added into a temporary reference line segment;
step S104.10: if the current grid point is invisible and the previous grid point is visible, firstly, the parameter a of an OSP line segment formed by the current point and the previous point is solved, and the intersection point direction k between the OSP line segment and a reference line is solved; then, recording the ending direction of the line segment as k and the parameter as a, and adding a temporary reference line segment; finally, updating the line segment between the original reference line and the direction from the intersection point to the current intersection point by using the temporary reference line segment, if the original reference line does not have the previous intersection point, updating the line segment between the original reference line starting point and the current intersection point by using the temporary reference line segment, and emptying the temporary reference line;
step S104.11: if the current grid point and the previous grid point are not in sight, no operation is needed;
step S104.12: if the current grid point is visible and the previous grid point is invisible, the ending direction k and the parameter a of an OSP line segment formed by the two points are solved, the intersection point direction k' between the OSP line segment and the reference line is solved according to a formula (12), and the line segment from the intersection point to the current grid point forms the first segment of the temporary reference line segment;
step S104.13: if the current proximity still has a point of investigation, go to step S104.7;
step S104.14: if the last investigation point of the current proximity is visible and the temporary reference line is not empty, updating the line segment from the first intersection point to the last intersection point on the original reference line by using the temporary reference line segment, if the last intersection point does not exist, updating the line segment from the starting point to the last intersection point of the original reference line, and finally emptying the temporary reference line;
step S104.15: taking the updated reference line as the current reference line;
step S104.16: if the proximity plane still remains to be investigated, executing step S104.2 until the investigation is completed;
the visibility determination of the other several sections (2, 3, and 4) is performed with reference to the section 1.
Specifically, the step S104 further includes:
according to the method, a 0-1 matrix which is in the same dimension as a DEM (digital elevation model) of a visibility domain is used for recording the analysis result of the visibility domain, the visibility calculation result of each grid point is stored in a corresponding element with the same column number, the value 0 represents the non-visibility, and the value 1 represents the visibility. 1. The calculation results of 2, 3 and 4 areas are respectively recorded in a visibility matrix V with the same dimension as the calculation area 1 、V 2 、V 3 、V 4 In (1), the final see-through matrix should be：
V＝(V 1 +V 2 )·(V 3 +V 4 ) (13)
When the element in the V is 0, grid points with the same row number and the same column number as the element in the visibility analysis area are invisible; and when the element in the V is 1, the grid point with the same row number and column number as the element in the visibility analysis area is subjected to visibility.
To verify the effect of the present invention, the following experiment was performed:
the invention adopts a C # programming mode to realize that the OSP space reference line through-view domain analysis method (for simple description, the OSP reference line algorithm is abbreviated as OSP algorithm in the following) based on DEM, the R3 algorithm, the XDraw algorithm, the reference plane algorithm and other through-view domain analysis algorithms are used for comparative analysis, and simultaneously adopts a front-end interactive platform for development of Cesum to carry out interactive verification on the algorithms.
Experimental data from the http:// www. Gscloud. Cn website, 30 m grid spacing DEM for 3 ASTER GDEM, collected in 2009, in the format TIFF, numbered N34E114, N41E119, N28E097, respectively, as shown in fig. 10, where (a) is N34E114, (b) is N41E119, and (c) is N28E097. The position and terrain characteristics of the three DEMs are as follows:
table 1DEM experimental data overview
The total area of the three DEMs is about 3.6 kilo square kilometers, the three DEMs cover the typical fluctuation forms of various terrains such as plain terrains, hilly terrains, mountain terrains and the like, the three DEMs have wide representativeness, and all statistical calculation and comparison of the experiment are carried out in the three areas so as to ensure that the analysis result has certain universality.
1. Comparison of precision experiment
(1) Algorithm accuracy analysis
The R3 algorithm is used as a traditional classical sight field analysis classical algorithm, the influence of all grid lines on sight direction on sight is considered, the calculation precision is high, and the result is reliable, so that the R3 algorithm result is used as a standard in the experiment to evaluate the precision of other algorithms. In the experiment, the observation positions are randomly selected on 3 DEMs at every 25 m increase from the height of 1 m away from the ground, the analysis radius (the radius range is 270-10 m, 680 m) is randomly set, and the experiment is carried out for 50 times respectively. In each experiment, the same area is calculated once by using an OSP reference line algorithm, an XDraw algorithm, a reference plane algorithm, and an R3 algorithm, and then the algorithm results are compared with the R3 algorithm results to analyze the error of each algorithm, and the statistical result is shown in fig. 11. The curve in fig. 11 is the fitting result of the experimental polynomial, and can be found by comparing the error rate upper limit, the average error rate and the error rate median of each algorithm:
the OSP reference line algorithm has the advantages that the result is completely consistent with the R3 algorithm in all experiments;
(II) the error of the XDraw/reference plane algorithm is maximum near the ground, and gradually decreases along with the increase of the observation height;
thirdly, the error of the XDraw/reference surface algorithm is smaller when the XDraw/reference surface algorithm is close to the ground in the region with relief, and the XDraw/reference surface algorithm is slower to descend along with the observation height;
(IV) the XDraw/reference surface algorithm has larger error when being close to the ground in an area with flatter terrain, and is faster to descend along with the observation height;
in conclusion, compared with the XDraw/reference plane algorithm, the OSP reference line algorithm has greater precision advantage in a terrain relief area and a low-altitude terrain flat area.
(2) Degree of aggregation of error points
In the through-the-horizon analysis, compared with the error rate, the method is more concerned about whether error points are gathered, and even if the error rate of the algorithm is low, if the error points are gathered, the application of the result is greatly influenced. In order to determine whether the situation of error point aggregation exists in the result of the through-view field analysis, the experiment introduces the concept of x-neighborhood errors. An x-neighborhood of a point refers to the area x grid distances around the point, as shown in FIG. 12. Whether x-neighborhood errors occur at a certain grid point means that whether the same type of through-vision judgment errors occur at 90% grid points in the x-neighborhood of the certain grid point on the premise that the through-vision judgment errors occur at the certain grid point. The experiment evaluates the error point aggregation degree of the algorithm by investigating whether neighborhood errors occur in the result of the through-view domain and the size of the neighborhood errors. The DEM resolution used in this experiment was 30 meters and the areas and number of grid points involved in each type of neighborhood error are shown in table 2.
TABLE 2 number of regions and grid points covered by each type of neighborhood
And 3 DEMs are respectively observed for 500 times at intervals of 2 meters and 2 meters between 2 meters and 52 meters from the ground, and the observation area is 120.5604 square kilometers. Experiments show that the OSP reference line algorithm has no errors, correspondingly, the possibility of error point aggregation does not exist, and most experimental results of the XDraw and reference plane algorithms have the condition of error aggregation. Fig. 13 shows the result of the through-view analysis calculated by the OSP reference line algorithm and the XDraw algorithm respectively in one experiment. As shown in fig. 13, (b) a 1-neighborhood error occurred in part of the XDraw calculation results, while (a) part of the OSP reference line algorithm did not have an error.
The error distribution conditions of the neighborhood points of all XDraw algorithms in the experiment are tabulated as figure 14 and table 3 according to height and terrain, and can be seen from the related chart: XDraw has more neighborhood errors on the whole, and 1-neighborhood calculation errors can occur in more than half of calculation results in mountainous (N28E 097) and hilly (N41E 119); the larger the neighborhood range, the lower the occurrence probability, but in the calculation result of the mountain land (N28E 097), the 10-neighborhood errors still have the occurrence frequency of 0.15%, which means that in this experiment, 0.15% of the calculation result of the N28E097 has continuous 396 common view field analysis errors, and the area of the related continuous errors is as high as 39.69 ten thousand square meters.
The experiment shows that the reliability of the calculation result of the XDraw algorithm and the similar LOS approximation algorithm (such as a reference plane algorithm) is low, the XDraw algorithm and the LOS approximation algorithm are difficult to apply to practice, and the improvement of the through-vision field calculation accuracy is very necessary.
Table 3XDraw algorithm neighborhood error rate distribution over three terrains
2. Algorithm speed contrast experiment
The algorithm of the invention has the same precision as the R3 algorithm, but the complexity of the R3 algorithm is O (R) 3 ) The algorithm time consumption increases very fast as the calculation range increases, and the complexity of the algorithm is O (r) 2 ) In the same order as the XDraw and reference plane algorithms. The invention carries out 9000 speed comparison experiments on three DEMs of different terrains at the heights of 1 meter, 100 meters and 1000 meters by randomly selecting the radius, draws the time consumption of an R3 algorithm and an OSP reference line algorithm in a graph 15 according to the radius size sequence, draws the time consumption comparison of the OSP reference line algorithm, an XDraw algorithm and a reference plane algorithm in a graph 16 according to the calculation point quantity sequence (curves in the two graphs are both obtained by polynomial fitting), and lists the calculation speed comparison statistical conditions of the OSP reference line algorithm, the XDraw algorithm and the reference plane algorithm in a table 4 (the head of the table is an algorithm for comparison).
The experiment shows that the method has the speed far higher than that of the R3 algorithm, is close to that of a reference plane algorithm and an XDraw algorithm, and can achieve the calculation speed meeting the actual application on the premise of keeping the same precision as that of the R3 algorithm.
TABLE 4 speed comparison of OSP reference line algorithm to reference plane algorithm, XDraw algorithm
According to the experimental result, the accuracy of the OSP reference line algorithm is completely the same as that of the R3 algorithm, and the OSP reference line algorithm is an accurate algorithm, but the calculation speed is far higher than that of the R3 algorithm and is in the same order as that of the XDraw and reference plane algorithms. Although the speed is not as fast as the XDraw and reference plane algorithm, the advantages of the modern DEM product can be reflected, the calculation result is reliable, the continuous film communication judgment error can not occur, and meanwhile, the speed can also meet the requirement of large-area calculation.
Example 2
Another OSP spatial reference line through-view domain analysis method based on the regular grid DEM of the present invention is substantially the same as the OSP spatial reference line through-view domain analysis method based on the regular grid DEM in embodiment 1, and the difference is that:
when the through-view area analysis area is divided, taking a diagonal line of the through-view area analysis area passing through an observation point as a boundary, dividing the through-view area into an area 1, an area 2, an area 3 and an area 4, wherein the area 1 and the area 2 are transversely and centrally symmetrical, and the area 3 and the area 4 are longitudinally and centrally symmetrical, as shown in fig. 17, the DEM grid points of the area 1 are shown in fig. 18;
calculating a visibility matrix of the visibility region analysis area according to the following mode:
V＝V 1 +V 2 +V 3 +V 4 (14)
wherein, V 1 、V 2 、V 3 、V 4 The through-view matrixes are respectively 1,2, 3 and 4 areas;
when the element in the V is 0, grid points with the same row number and the same column number as the element in the visibility analysis area are invisible; and when the element in the V is 1, the grid point with the same row number and column number as the element in the visibility analysis area is subjected to visibility.
To verify the effect of the present invention, the following experiments were performed:
for convenience of description, another OSP spatial reference line through-the-horizon analysis method based on the regular grid DEM in this embodiment is called a fast OSP reference line algorithm.
1. Precision comparison experiment
The algorithm precision analysis experiment is substantially the same as that of the algorithm in the embodiment 1, except that a fast OSP reference line algorithm is added, and the experimental result is shown in FIG. 19.
2. Speed contrast experiment
Overall, the algorithm speed of the present invention is much higher than that of the R3 algorithm, and is close to that of the XDraw and reference plane algorithms, as shown in fig. 20 to 22. The calculation speed of the invention can be influenced by comprehensive factors such as the breaking degree of the terrain, the visibility and the like. According to 50000 times of experimental results of different heights on three terrains of mountains, hills and plains, the speed of the invention is higher at the position where a survey point is closer to the ground under the same condition, and then the speed is reduced along with the increase of the height of the viewpoint from the ground, and finally the speed is relatively stable. The calculation speed of the invention is compared with the XDraw and reference plane algorithms such as Table 5:
TABLE 5 speed comparison of fast OSP reference line algorithm with reference plane algorithm, XDraw algorithm
According to the experimental result, the fast OSP reference line algorithm can improve the speed by 35% compared with the reference plane algorithm and can improve the speed by 45% compared with the XDraw algorithm in the mountainous region close to the ground. In high altitude, the speed of the fast OSP reference line algorithm is slightly lower than that of the reference plane algorithm and is slightly better than that of the XDraw algorithm.
The above shows only the preferred embodiments of the present invention, and it should be noted that it is obvious to those skilled in the art that various modifications and improvements can be made without departing from the principle of the present invention, and these modifications and improvements should also be considered as the protection scope of the present invention.
Claims (7)
1. An OSP space reference line through-view domain analysis method based on a regular grid DEM is characterized by comprising the following steps:
step 1: partitioning the through-view area analysis area based on the regular grid DEM; the step 1 comprises the following steps:
dividing a longitudinal topographic line of a through-view area analysis area into a 1 area and a 2 area and dividing a transverse topographic line of the through-view area analysis area into a 3 area and a 4 area on a regular grid DEM by taking an observation point as a boundary; the longitudinal topographic line is a continuous broken line formed by connecting column DEM grid points, and the transverse topographic line is a continuous broken line formed by connecting row DEM grid points;
step 2: constructing a corresponding station center coordinate system in each partition, and constructing a corresponding OSP space coordinate system based on the constructed station center coordinate system; the step 2 comprises the following steps:
step 2.1: constructing a first station center coordinate system in the area 1: establishing a right-hand space rectangular coordinate system by taking an observer viewpoint as an original point, taking the north-bias direction of a longitudinal topographic line of the DEM as a longitudinal axis Y, taking the east-bias direction of a transverse topographic line of the DEM as a transverse axis X and taking the normal vector direction of a reference ellipsoid where an observer is located as a vertical axis Z;
step 2.2: taking a Z axis of a first station center coordinate system as a rotating axis, respectively rotating by 180 degrees, 270 degrees and 90 degrees clockwise to obtain three coordinate systems which are respectively defined as a second station center coordinate system, a third station center coordinate system and a fourth station center coordinate system, and respectively establishing the second station center coordinate system, the third station center coordinate system and the fourth station center coordinate system in 2, 3 and 4 areas by taking an observation point as a center;
step 2.3: constructing an OSP space coordinate system by taking the direction k as a horizontal axis, the proximity p as a vertical axis and the axis slope g as a vertical axis on the basis of each constructed station center coordinate system; the direction k is a tangent value of a survey point in a station center coordinate system relative to a horizontal axis direction angle alpha; the proximity p is the reciprocal of the value of the horizontal axis x of the investigation point in the station center coordinate system; the axial slope g is the product of the relative elevation z of the survey point and the proximity p;
and step 3: converting grid points corresponding to the terrain lines in each partition into an OSP space coordinate system;
and 4, step 4: judging the visibility of each partition in an OSP space coordinate system, and combining the visibility judgment results of each partition so as to obtain the visibility judgment result of a visibility domain analysis area;
in the step 4, the visibility of each partition is determined as follows:
step 4.1: after the grid points corresponding to each topographic line in each partition are converted from the respective station center coordinate system to the respective OSP space coordinate system, the maximum proximity is obtainedForming an initial reference line according to the sequence of the direction values from small to large, solving the ending direction k and the parameter a of each OSP line segment of the initial reference line, and solving the initial point axis slope g of the initial reference line 0 Recording grid points on the initial reference line as a through view; the axial slope is the product of the relative elevation and proximity of the survey point;
step 4.2: moving to a next proximity plane according to the sequence of the proximity of grid points in an OSP space coordinate system from large to small;
step 4.3: projecting the current reference line onto the next proximity plane in step 4.2;
step 4.4: finding out a first grid point on a current proximity plane and an OSP line segment L of a reference line containing the direction of the first grid point in a sequence from small direction to large direction, and solving the axial slope difference delta G between the first grid point and the L 0 By Δ G 0 Judging the visibility of the first grid point and recording the visibility to a corresponding visibility matrix;
step 4.5: if the first grid point is invisible, executing step 4.7;
step 4.6: if the first grid point is in full view, the first grid point and the grid point which is nearest outside the investigation region and has the same proximity are used for forming an OSP line segment, the ending direction k and the parameter a of the OSP line segment are calculated and used as the first segment of the temporary reference line segment, and g is updated 0 The initial axis slope of the first section of the current temporary reference line segment;
step 4.7: recording a next grid point as a current calculation point, checking whether the direction of the current point is greater than the tail direction value of the L, if so, recording the next section connected with the L on the reference line as a current reference line segment L until the direction of the current point is less than the tail direction value of the L;
step 4.8: calculating the axial slope difference delta G between the current grid point and the current reference line segment L by an axial slope difference accumulation formula i By Δ G i Judging and recording the visibility;
step 4.9: if the current grid point and the previous grid point are all in sight, the ending direction k and the parameter a of an OSP line segment formed by the two points are solved, and the OSP line segment is added into a temporary reference line segment;
step 4.10: if the current grid point is not visible and the previous grid point is visible, firstly, the parameter a of the OSP line segment formed by the current point and the previous point is solved, and the intersection point direction k between the OSP line segment and the reference line is solved; then, recording the ending direction of the line segment as k and the parameter as a, and adding a temporary reference line segment; finally, updating the line segment between the original reference line and the direction from the intersection point to the current intersection point by using the temporary reference line segment, if the original reference line does not have the previous intersection point, updating the line segment between the original reference line starting point and the current intersection point by using the temporary reference line segment, and emptying the temporary reference line;
step 4.11: if the current grid point and the previous grid point are invisible, no operation is needed;
step 4.12: if the current grid point is visible and the previous grid point is invisible, the ending direction of an OSP line segment formed by the two points is determined to be k and a parameter a, the intersection point direction k' between the OSP line segment and a reference line is determined, and a line segment from the intersection point to the current grid point forms a first section of a temporary reference line segment;
step 4.13: if the current proximity still has a point of investigation, executing step 4.7;
step 4.14: if the last investigation point of the current proximity is visible and the temporary reference line is not empty, updating the line segment from the first intersection point to the last intersection point on the original reference line by using the temporary reference line segment, if the last intersection point does not exist, updating the line segment from the starting point to the last intersection point of the original reference line, and finally emptying the temporary reference line;
step 4.15: taking the updated reference line as the current reference line;
step 4.16: if there are proximity planes yet to be investigated, step 4.2 is performed until the investigation is completed.
2. The OSP spatial reference line through-the-horizon analysis method based on a regular grid DEM of claim 1, wherein step 1 comprises:
and dividing the analysis area of the visibility region into a region 1, a region 2, a region 3 and a region 4 on the regular grid DEM by taking a diagonal line of the analysis area of the visibility region passing through an observation point as a boundary, wherein the region 1 and the region 2 are transversely centrosymmetric, and the region 3 and the region 4 are longitudinally centrosymmetric.
3. The OSP spatial reference line through-the-horizon analysis method based on regular grid DEM as claimed in claim 1, wherein said step 4 further comprises:
when the observation points are used as boundaries for partitioning, the visibility matrix of the visibility region analysis area is calculated according to the following mode:
V＝(V 1 +V 2 )·(V 3 +V 4 ) (13)
wherein, V 1 、V 2 、V 3 、V 4 The perspective matrixes are respectively 1,2, 3 and 4 areas;
when the element in the V is 0, grid points with the same row number and the same column number as the element in the visibility analysis area are invisible; and when the element in the V is 1, the grid point with the same row number and column number as the element in the visibility analysis area is subjected to visibility.
4. The OSP spatial reference line through-the-horizon analysis method based on a regular grid DEM of claim 1, wherein step 4 further comprises:
when the partition is carried out by taking the diagonal line of the through-view area analysis area passing through the observation point as a boundary, the through-view matrix of the through-view area analysis area is calculated as follows:
V＝V 1 +V 2 +V 3 +V 4 (14)
wherein, V 1 、V 2 、V 3 、V 4 The perspective matrixes are respectively 1,2, 3 and 4 areas;
when the element in the V is 0, the grid points with the same row number and the same column number as the element in the visibility analysis area are invisible; and when the element in the V is 1, the grid point with the same row number and column number as the element in the visibility analysis area is subjected to visibility.
5. The OSP spatial reference line through-the-horizon analyzing method based on regular grid DEM of claim 1, wherein the parameter a is derived as follows:
a＝u(z B -z A )
wherein u is the reciprocal of the grid width of DEM, and z A 、z B The vertical axis values of the station center coordinates of any two adjacent grid points A and B with the same proximity are respectively, and the direction of B is greater than that of A.
6. The OSP spatial reference line through-the-horizon analysis method based on regular grid DEM of claim 1 wherein the axis slope difference cumulative formula is:
ΔG i ＝ΔG i-1 +(a Ci -a Ri )(k i -k i-1 ) (11)
wherein, Δ G i Is the axial slope difference, Δ G, between the current grid point and the current reference line segment L i-1 Is the difference in the axial slope of the last segment end and the last reference segment, a Ci And a Ri The slope, k, of the line segment on the reference line and the investigation line in the current segment i 、k i-1 Respectively the current segment and the last segment.
7. The OSP spatial reference line through-the-horizon analysis method based on a regular grid DEM as claimed in claim 6, wherein the direction k' of intersection between the OSP line segments and the reference line is found in step 4.12 by:
k′＝k i -ΔG i /(a Ci -a Ri ) (12)。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010054396.3A CN111325791B (en) | 2020-01-17 | 2020-01-17 | OSP space reference line through-view domain analysis method based on regular grid DEM |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010054396.3A CN111325791B (en) | 2020-01-17 | 2020-01-17 | OSP space reference line through-view domain analysis method based on regular grid DEM |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111325791A CN111325791A (en) | 2020-06-23 |
| CN111325791B true CN111325791B (en) | 2023-01-06 |
Family
ID=71168645
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010054396.3A Active CN111325791B (en) | 2020-01-17 | 2020-01-17 | OSP space reference line through-view domain analysis method based on regular grid DEM |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111325791B (en) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112162248A (en) * | 2020-08-21 | 2021-01-01 | 中国人民解放军93114部队 | Rapid calculation method and device for radar terrain shielding blind area |
| CN113313826B (en) * | 2021-05-27 | 2022-08-30 | 南京邮电大学 | Scale-adaptive visual domain analysis method |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR20110085107A (en) * | 2010-01-19 | 2011-07-27 | 경북대학교 산학협력단 | Geographic information analysis method using of viewshed frequency analysis |
| CN106530398A (en) * | 2016-12-01 | 2017-03-22 | 南京师范大学 | Terrain visibility analysis-oriented visibility graph network construction method |
| CN107886570A (en) * | 2017-09-22 | 2018-04-06 | 中国矿业大学 | A kind of visible range computational methods for taking into account precision and efficiency |
| CN108919319A (en) * | 2018-05-15 | 2018-11-30 | 中国人民解放军战略支援部队信息工程大学 | Sea island reef satellite image Pillarless caving localization method and system |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104717462A (en) * | 2014-01-03 | 2015-06-17 | 杭州海康威视系统技术有限公司 | Supervision video extraction method and device |
-
2020
- 2020-01-17 CN CN202010054396.3A patent/CN111325791B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR20110085107A (en) * | 2010-01-19 | 2011-07-27 | 경북대학교 산학협력단 | Geographic information analysis method using of viewshed frequency analysis |
| CN106530398A (en) * | 2016-12-01 | 2017-03-22 | 南京师范大学 | Terrain visibility analysis-oriented visibility graph network construction method |
| CN107886570A (en) * | 2017-09-22 | 2018-04-06 | 中国矿业大学 | A kind of visible range computational methods for taking into account precision and efficiency |
| CN108919319A (en) * | 2018-05-15 | 2018-11-30 | 中国人民解放军战略支援部队信息工程大学 | Sea island reef satellite image Pillarless caving localization method and system |
Non-Patent Citations (6)
| Title |
|---|
| An Improved Algorithm for Computing Viewshed Based on Reference Planes;Ye Zhi etal.;《An Improved Algorithm for Computing Viewshed Based on Reference Planes》;20110811;全文 * |
| Geometric Algorithms for Siting of Air Defense Missile Batteries;Prof Wm Randolph Franklin etal.;《https://www.researchgate.net/publication/2808061》;19980730;全文 * |
| 基于DEM实现城市高斯投影变形可视化分析;尚金光等;《城市勘测》;20181031(第05期);全文 * |
| 基于DEM的地形可视域分析关键技术;张斌等;《计算机科学》;20130915(第09期);全文 * |
| 基于大地坐标的通视性检查算法;余文广等;《计算机仿真》;20081115(第11期);全文 * |
| 雷达探测范围地形可视域分析方法;赵永刚等;《雷达探测范围地形可视域分析方法》;20190630;第44卷(第3期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111325791A (en) | 2020-06-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US7646916B2 (en) | Linear analyst | |
| CN106649821B (en) | Spatial target index construction, collision early warning, region and nearest neighbor query method | |
| CN111325791B (en) | OSP space reference line through-view domain analysis method based on regular grid DEM | |
| CN103324916B (en) | Vehicle-mounted and aviation LiDAR data method for registering based on building profile | |
| CN111540005B (en) | Loop detection method based on two-dimensional grid map | |
| CN109345617B (en) | Chain type high-precision splicing and adjustment method based on long-strip multi-station point cloud | |
| CN111398918B (en) | Radar detection capability analysis method under complex mountain environment | |
| CN105868326A (en) | Pipeline data storage method | |
| CN111899332A (en) | Overhead transmission line three-dimensional design method based on oblique photogrammetry technology | |
| CN111553292A (en) | Rock mass structural plane identification and occurrence classification method based on point cloud data | |
| CN107119657A (en) | A kind of view-based access control model measures foundation ditch monitoring method | |
| CN110827405A (en) | Digital remote sensing geological mapping method and system | |
| CN108153979A (en) | Deformation information extraction method based on InSAR, terminal and storage medium | |
| CN111765869A (en) | Different-gradient road earthwork measurement method based on oblique photography technology | |
| CN113238228B (en) | Three-dimensional earth surface deformation obtaining method, system and device based on level constraint | |
| CN103116183A (en) | Method of oil earthquake collection surface element covering degree property body slicing mapping | |
| CN111832582A (en) | Method for classifying and segmenting sparse point cloud by using point cloud density and rotation information | |
| CN111681313B (en) | Space vision analysis method based on digital topography and electronic equipment | |
| CN103777229A (en) | VSP observation system design method facing objective layer | |
| CN114529688B (en) | Method for calculating global topography correction based on multi-resolution polyhedron under hexagonal grid index | |
| CN115930800A (en) | Tunnel face displacement field monitoring method based on three-dimensional laser point cloud | |
| CN110909448B (en) | High-frequency sky wave return scattering ionization diagram inversion method | |
| CN114877860A (en) | Long tunnel multi-station combined measurement combination resolving method and device and storage medium | |
| Guan et al. | Fast approximate viewshed analysis based on the regular-grid digital elevation model: X-type partition proximity-direction-elevation spatial reference line algorithm | |
| Zhou et al. | An optimized fuzzy K-means clustering method for automated rock discontinuities extraction from point clouds |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |