CN111325791A - 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 PDF

Info

Publication number
CN111325791A
CN111325791A CN202010054396.3A CN202010054396A CN111325791A CN 111325791 A CN111325791 A CN 111325791A CN 202010054396 A CN202010054396 A CN 202010054396A CN 111325791 A CN111325791 A CN 111325791A
Authority
CN
China
Prior art keywords
osp
point
reference line
grid
visibility
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.)
Granted
Application number
CN202010054396.3A
Other languages
Chinese (zh)
Other versions
CN111325791B (en
Inventor
吴传均
管凌霄
夏青
陈刚
甄学忠
沈宝鸿
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Information Engineering University of PLA Strategic Support Force
Original Assignee
Information Engineering University of PLA Strategic Support Force
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Information Engineering University of PLA Strategic Support Force filed Critical Information Engineering University of PLA Strategic Support Force
Priority to CN202010054396.3A priority Critical patent/CN111325791B/en
Publication of CN111325791A publication Critical patent/CN111325791A/en
Application granted granted Critical
Publication of CN111325791B publication Critical patent/CN111325791B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • Software Systems (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Remote Sensing (AREA)
  • Computer Graphics (AREA)
  • Processing Or Creating Images (AREA)

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 a through-view area analysis area based on a 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 so as to obtain the visibility judgment result of the visibility domain analysis area. When the through-vision 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 in the same order as that of a reference plane algorithm; when the perspective judgment area is partitioned according to diagonal lines by an observation point, the precision 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

OSP space reference line through-view domain analysis method based on regular grid DEM
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 a certain distance range of a standing point and the area through which an observer looks. 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 well-known R2, R3, XDraw algorithm (Franklin R. geological Algorithms for training of Air Defense Missile batteries [ R ]. Columbus Division: A Research Project for batteries [ c ],1994.), Wu Brillian blue (Wu Brillian. reference plane-based visual field algorithm [ J ] mapping information and engineering, 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 studied and improved. 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 method cannot be applied to actual operation; the R2, XDraw, and reference plane algorithm is an approximate fast algorithm, the calculation speed of which can meet the requirement of practical use, but errors are inevitably 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 the algorithm precision on the application is not large due to low data precision of the traditional DEM, however, the study finds that although the error ratio of the approximation algorithms such as XDraw is low, the distribution is not uniform, the regional perspective judgment is possibly wrong, so that the application of the calculation result is influenced, along with the development of the technology, the DEM data quality 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 precision problem, such as Wu (Wu H, Pan M, Yao L, et al. A Partition-Based Serial Algorithm for generating Viewswed on Massive Dems [ J ]. International Journal of Geographanimation Science,2007,21(9): 955. 964.) to improve the 2011 problem by proposing the calculated defects of the Reference plane Algorithm for the topography such as the peaks and valleys, Zhi (Zhi Y, Wu L, Sui Z, et al. Improved Algorithm for computing Viewswed Based Reference Planes [ C ]// Ieee 19 International Conference Reference gels, Shanghai, Shannao, 1-5, etc.), but the result error is only reduced by less than 50%, and the time consumption of the algorithm is increased by 4 to 10 times, which shows that the traditional approximate fast algorithm is difficult to combine the advantages of speed and precision from the basis design.
In fact, in terms of calculation speed, the algorithm is still to be further improved at present, 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 good effects under the condition of portable hardware equipment without network environment, and particularly difficult to exert speed advantages under field conditions such as exploration, archaeology, military application and the like, so that the speed bottleneck of analysis through a visual field is solved from the essence of the calculated 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 a through-view area analysis area based on a 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 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 taking the observation points as boundaries, dividing the longitudinal topographic line of the analysis area of the visibility region into a region 1 and a region 2, and dividing the transverse topographic line of the analysis area of the visibility region into a region 3 and a region 4 on the regular grid DEM.
Further, the step 1 comprises:
and dividing the through-view area analysis area into a 1 area, a 2 area, a 3 area and a 4 area on the regular grid DEM by taking a diagonal line of the through-view area analysis area passing through an observation point as a boundary, wherein the 1 area and the 2 area are transversely centrosymmetric, and the 3 area and the 4 area 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;
and 2.3, constructing an OSP space coordinate system by taking a direction k as a horizontal axis, a proximity p as a vertical axis and an axis slope g as a vertical axis on the basis of each constructed station center coordinate system, wherein the direction k is a tangent value of an investigation point in the station center coordinate system relative to a direction α of the horizontal axis, the proximity p is an inverse number of a horizontal axis x value of the investigation point in the station center coordinate system, and the axis slope g is a product of a relative elevation z of the investigation point 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 terrain line in each partition are respectively converted into each OSP space coordinate system from each station center coordinate system, the grid points with the maximum proximity are formed into an initial reference line 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 line are solved, and the initial point axial slope g of the initial reference line is solved0Recording grid points on the initial reference line as a through view;
step 4.2: moving to a next proximity plane according to the sequence of the proximity of the grid points in the 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 L0By Δ G0Judging 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 sight, the first grid point and the grid point which is nearest outside 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 calculated to be used as the first segment of the temporary reference line segment, and g is updated0The 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: through the formula of shaft slope difference accumulationCalculating the axial slope difference delta G between the current grid point and the current reference line segment LiBy Δ GiJudging 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 intersection point by using the temporary reference line segment, if the original reference line does not have the intersection point, updating the line segment between the original reference line starting point and the 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 not in sight, 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 a line segment from a first intersection point to the tail of the original reference line by using the temporary reference line segment, if the last intersection point does not exist, updating a line segment from the starting point to the tail 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＝(V1+V2)·(V3+V4) (13)
wherein, V1、V2、V3、V4The 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 perspective analysis area are not in perspective; 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 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＝V1+V2+V3+V4(14)
wherein, V1、V2、V3、V4The 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 perspective analysis area are not in perspective; 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(zB-zA)
wherein u is the reciprocal of the grid width of DEM, and zA、zBThe vertical axis values of the station coordinates of any two adjacent grid points A, B of the same proximity, respectively, and the direction of B is greater than the direction of a.
Further, the shaft slope difference accumulation formula is:
ΔGi＝ΔGi-1+(aCi-aRi)(ki-ki-1) (11)
wherein, Δ GiIs the axial slope difference, Δ G, between the current grid point and the current reference line segment Li-1Is the difference in the axial slope of the last segment end and the last reference segment, aCiAnd aRiThe slope, k, of the line segments on the reference line and the investigation line in the current segmenti、ki-1Respectively 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′＝ki-ΔGi/(aCi-aRi) (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, and the method is not required to sample, interpolate and calculate redundantly, is suitable for large-area long-distance visibility analysis and is also suitable for serving as a basic algorithm for other parallel topographic 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 an 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 modern DEM products 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 diagonals of the observation points of the through-view area analysis area are used as boundaries for partitioning, the algorithm precision is far higher than that of other rapid approximation algorithms, the precision can be improved by more than 5 times compared with that of traditional approximation algorithms (such as XDraw and reference plane algorithms) and the like (and parallel algorithms based on the traditional approximation algorithms), the speed is far higher than that of an 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 within the height 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 sectional view 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 embodiment of the present invention based on a regular grid DEM OSP spatial reference line through vision field analysis method;
FIG. 11 is an error distribution of the OSP spatial reference line through vision field 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 altitudes;
FIG. 12 is a diagram illustrating an exemplary neighborhood of DEM grid points for an OSP spatial reference line through view analysis method based on a regular grid DEM according to an embodiment of the present invention;
fig. 13 is a comparison graph of the 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 graph of the frequency of occurrence of neighborhood errors in different terrains/different heights by using the R3 algorithm as a reference and the XDraw algorithm in the OSP spatial reference line through vision field 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 with the calculated radius of the R3 algorithm and the OSP reference line algorithm of the OSP spatial reference line through vision field analysis method based on the regular grid DEM according to the embodiment of the present invention;
FIG. 16 is a time consumption and calculated point amount variation trend statistical chart of each algorithm of an OSP space reference line through vision field analysis method based on a regular grid DEM in the embodiment of the present invention;
FIG. 17 is a schematic sectional view of another OSP spatial reference line through vision field 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 for another OSP spatial reference line-through-view field analysis method based on a regular grid DEM in accordance with an embodiment of the present invention;
FIG. 19 is a graph of the error distribution of another method for OSP spatial reference line through vision field analysis based on regular grid DEM according to the present invention, with reference to the R3 algorithm, over different terrains/different altitudes;
FIG. 20 is a time consuming trend statistical plot of calculated radius variation for another method for OSP spatial reference line through view analysis based on a regular grid DEM in accordance with the present invention and the R3 algorithm;
FIG. 21 is a comparison graph of the time consumption of another OSP spatial reference line through vision field analysis method based on a regular grid DEM according to the present invention and a reference plane algorithm at different terrains and different terrain heights;
fig. 22 is a comparison graph of the time consumption of the method for analyzing the OSP spatial reference line-through-view area based on the regular grid DEM according to the embodiment of the present invention and the XDraw algorithm in different terrains and different ground heights.
Detailed Description
The invention is further illustrated by the following examples in conjunction with 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 so as 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. If the connecting line of two adjacent grid points on the DEM can express the terrain between the two points, 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, with the observation points as boundaries, the vertical topographic lines of the through-view area analysis area in the regular grid DEM are divided into 1 area and 2 area, and the horizontal topographic lines are divided into 3 area and 4 area.
Specifically, the step S102 includes:
step S102.1: constructing a first station center Coordinate System (Topocentric Coordinate System) in region 1: the observer viewpoint is taken as the origin, the north direction of the longitudinal grid line of the DEM is taken as the longitudinal axis Y (longitudinal axis), the east direction of the transverse grid line of the DEM is taken as the transverse axis X (transverse axis), andthe normal vector direction of the reference ellipsoid where the observer is located is a vertical axis Z (vertical axis) to establish a right-hand space rectangular coordinate system, in the station center coordinate system of FIG. 2, the station center coordinate of point A is (x, y, Z), the horizontal distance from the origin is l, the sight line direction angle is α, and the elevation angle is α
Figure BDA0002372309010000081
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 sheltering object, and the following three sheltering conditions are necessarily met simultaneously if the sheltering object shelters 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 block an object with a small proximity (far from the observer), whereas otherwise, the proximity is a determination value of the "blocking condition (i)".
p＝1/x (1)
The direction (denoted as k) is defined as the tangent value of a point under consideration to the direction angle α of the horizontal axis in the station center coordinate system, and is defined as formula (2). if the values of k at the two points are different, the object representing that the directions at the observation point are different and the directions are not overlapped does not form an occlusion relation.
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
Figure BDA0002372309010000091
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)
Figure BDA0002372309010000092
The OSP spatial coordinate system is defined as: a space 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 space relation 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 the following equation:
Figure BDA0002372309010000101
dividing the left and right of equation (5) by x yields equation (6), and equation (6) is the expression after the straight line conversion to the OSP space. By combining with the relevant mathematical principle, the direction, the axis slope and the proximity of the straight line are known to be in linear relation, the property is called as the linear invariance of the OSP space, and the property lays a theoretical foundation for the view-through algorithm of the invention.
Figure BDA0002372309010000102
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 terrain line is a straight line segment constructed by two adjacent grid points, and according to the OSP space property, appropriate constants a and b are determined to enable the line segment to be consistent with the formula (7) in the OSP space, and the terrain 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,ke) Wherein k iss,keLeft and right boundaries of the domain, a, b, k, respectively, for the line segment directionss,keA parameter called L. If the station center coordinates of two adjacent grid points are known to be A (x, y) respectivelyA,zA)、B(x,yB,zB) 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)
Figure BDA0002372309010000103
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 section in the broken lines and the addresses of the sections before and after the section, and the invention needs to sequentially search the intersection point of the two broken lines along the broken lines, so the structure has the advantages of correctly expressing the structure of the broken lines and being convenient for searching, and when determining that a plurality of sections of broken lines need to be modified, only the addresses of the sections before and after the first and last sections of the relevant sections 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 C1Constructing a terrain line by using the grid points, and converting the terrain line into an OSP space to form an initial reference line LRWith C2Constructing a terrain line by the grid points, and converting the terrain line into an OSP space to form a broken line L to be inspected2Is prepared by mixing LRTranslate to L along the P axis2In the plane of the plane, with LRAnd L2Dividing two topographic lines into several segments by using their cross points as boundary, combining LRAnd L2The section with higher middle axle slope forms a new reference line L'RAnd L'RWill be further reacted with C3The terrain lines constructed by the column grid points are continuously updated through the steps until the region 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 all OSP space coordinate systems from respective station center coordinate systems, the grid points with the maximum proximity are formed into initial reference lines 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 solved0Recording 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 the 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 the first grid point on the current proximity plane and the reference line containing the direction of the point in the order of small direction to large directionOSP line segment L, and calculating the axial slope difference delta G between the first point and L0By Δ G0Judging 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,k1) With a known parameter, the proximity is p0The OSP coordinate of a certain point C in space is (k)C,gC,pC) And has pC＜p0，kC∈[k0,k1]. Note L (a, b | k)0,k1) And OSP plane k ═ kCIntersects the point D, and the coordinate of the intersection point is set as (k)C,gD,p0). Define C Point relative to L (a, b | p)0) Has an axial gradient difference of formula (9), corresponding to the black directed line segment in FIG. 9, and is easily defined from the definition of the axial gradient, Δ GCIf the vision is larger than 0, the vision is obstructed by C points, otherwise, the vision is obstructed. If used, rCIf 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).
ΔGC＝gC-gD(9)
Figure BDA0002372309010000121
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 outside 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 calculated to be used as the first segment of the temporary reference line segment, and g is updated0The 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 axis slope difference delta G between the current grid point and the current reference line segment L by an axis slope difference accumulation formulaiBy Δ GiJudging 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 geodesic is a polyline, the intersection points should be solved piecewise. Suppose LRIs a reference line, LCTo explore the line (OSP broken line composed of all grid points in the current proximity), take LCAnd LRThe 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 kiIn this paragraph LCAnd LRThe slope of the line segment on the upper line is aCiAnd aRiThe slope difference of the tail axis at the end of the period is recorded as Δ GiDue to the straight line characteristic, the following relation exists:
ΔGi＝ΔGi-1+(aCi-aRi)(ki-ki-1) (11)
equation (11) is called the cumulative equation of the axial slope difference, and only L needs to be obtained by the equationCThe direction k and the parameter a of each grid point can be accumulated to calculate LCEach grid point is relative to LRDifferential 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 Δ Gi-1And Δ GiOpposite positive and negative, then LRAnd LCThe i-th segment has an intersection point, the direction of which is k', and the direction can be obtained by the formula (12):
k′＝ki-ΔGi/(aCi-aRi) (12)
a grid point must be determined to be in-view in both regions 1/2 and 3/4 before the point can be finally determined to be 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 the cumulative calculation once each time the current reference line segment needs to be set as a next segment of 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 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 intersection point by using the temporary reference line segment, if the original reference line does not have the intersection point, updating the line segment between the original reference line starting point and the 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 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 according to a formula (12), and a line segment from the intersection point to the current grid point forms a first segment of a temporary reference line segment;
step S104.13: if the current proximity still has a point of investigation, executing 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 invention, a 0-1 matrix which is in the same dimension with the DEM of the visibility domain is used for recording the analysis result of the visibility domain, the visibility calculation result of each grid point is stored in the 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 area1、V2、V3、V4In (3), the final see-through matrix should be:
V＝(V1+V2)·(V3+V4) (13)
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 perspective analysis area are not in perspective; 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 in the following) based on DEM, R3 algorithm, XDraw algorithm, reference plane algorithm and other through-view domain analysis algorithms are used for comparative analysis, and simultaneously adopts a Cesium development front-end interaction platform to carry out interactive verification on the algorithms.
Experimental data from the http:// www.gscloud.cn website was acquired in 2009 at a 30-meter grid spacing DEM of 3 frames ASTER GDEM, in TIFF format, numbered N34E114, N41E119, and N28E097, respectively, as shown in FIG. 10, where (a) is N34E114, (b) is N41E119, and (c) is N28E 097. The position and terrain characteristics of the three DEMs are as follows:
table 1 DEM experimental data overview
Figure BDA0002372309010000141
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. Precision comparison 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 from the height of 1 m from the ground, and the analysis radius (the radius range is 270-10 m and 680 m) is randomly set every 25 m, and the experiment is carried out for 50 times. 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 results of the algorithms are compared with the results of the R3 algorithm to analyze the errors of the algorithms, wherein 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 result of the OSP reference line algorithm in all the experiments is completely consistent with that of the R3 algorithm;
(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;
(III) the error of the XDraw/reference surface algorithm is smaller near the ground in the region with relief topography, and the XDraw/reference surface algorithm is slower to decline 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 decline 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 analysis of the general view field, we are concerned about whether error points are gathered or not compared with the error rate, 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 judge whether the situation of error point aggregation exists in the result of the through-view field analysis, the concept of x-neighborhood error is introduced in the experiment. The x-neighborhood of a point refers to the area that is x grid distances around the point, as shown in FIG. 12. Whether x-neighborhood errors occur to a certain grid point means whether the same type of through-view judgment errors occur to 90% grid points in the x-neighborhood under the premise that the through-view judgment errors occur to the certain grid point. The experiment evaluates the error point aggregation degree of the algorithm by investigating whether neighborhood errors occur in the through-view domain result and the neighborhood size of the neighborhood errors. The DEM resolution used in this experiment was 30 meters, and the areas and grid point numbers involved in various types of neighborhood errors are shown in Table 2.
TABLE 2 number of regions and grid points covered by each type of neighborhood
Figure BDA0002372309010000151
Figure BDA0002372309010000161
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, and 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. Shown in fig. 13 are the results 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, while (a) part of the OSP reference line algorithm did not have an error.
The error distribution conditions of the neighborhood points of all the XDraw algorithms in the experiment are shown in FIG. 14 and Table 3 according to the height and terrain statistics, and can be seen from the related chart: as the XDraw has more neighborhood errors on the whole, more than half of the calculation results in mountainous (N28E097) and hilly (N41E119) regions have 1-neighborhood calculation errors; the larger the neighborhood range, the lower the occurrence probability, but in the calculation result of the mountain land (N28E097), 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 396 continuous point panoramic 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 the reference plane algorithm) is low, the XDraw algorithm and the LOS approximation algorithm are difficult to be applied to practice, and the improvement of the through-view calculation precision is very necessary.
Table 3 XDraw algorithm neighborhood error rate distribution over three terrains
Figure BDA0002372309010000162
Figure BDA0002372309010000171
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 along with the enlargement of the calculation range, and the complexity of the algorithm is O (r)2) In the same order as XDraw and reference plane algorithms. The invention carries out 9000 speed comparison experiments on three DEMs of different terrains at heights of 1 meter, 100 meters and 1000 meters at randomly selected radiuses, draws time consumption of an R3 algorithm and an OSP reference line algorithm in a radius size sequence in a graph 15, draws time consumption comparison of the OSP reference line algorithm and an XDraw and reference plane algorithm in a graph 16 (curves in the two graphs are both obtained by polynomial fitting) according to a calculation point quantity sequence, and lists calculation speed comparison statistics of the OSP reference line algorithm, the XDraw algorithm and the reference plane algorithm in a table 4 (a table head 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
Figure BDA0002372309010000172
According to experimental results, the accuracy of the OSP reference line algorithm is completely the same as that of the R3 algorithm, 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 visibility judgment error cannot occur, and meanwhile, the requirement of large-area calculation can be met in speed.
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 analysis area of the through-view area is partitioned, the through-view area is divided into an area 1, an area 2, an area 3 and an area 4 by taking a diagonal line of the analysis area of the through-view area passing an observation point as a boundary, the area 1 and the area 2 are transversely and centrally symmetrical, the area 3 and the area 4 are longitudinally and centrally symmetrical, as shown in fig. 17, a DEM grid point of the area 1 is shown in fig. 18;
calculating a visibility matrix of the visibility region analysis area according to the following mode:
V＝V1+V2+V3+V4(14)
wherein, V1、V2、V3、V4The 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 perspective analysis area are not in perspective; 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:
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 embodiment 1, except that a fast OSP reference line algorithm is added, and the experiment result is shown in FIG. 19.
2. Speed contrast experiment
On the whole, the algorithm speed of the 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 FIGS. 20-22. The calculation speed of the invention can be influenced by comprehensive factors such as the degree of terrain breakage, the visibility and the like. According to the results of 50000 times of experiments on three terrains of mountains, hills and plains with different heights, the speed of the invention is higher at the position where the observation 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
Figure BDA0002372309010000181
Figure BDA0002372309010000191
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 slightly superior to 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 several 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 (10)

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 a through-view area analysis area based on a 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 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.
2. The OSP spatial reference line through-the-horizon analysis method based on regular grid DEM as claimed in claim 1, wherein said step 1 comprises:
and dividing the longitudinal topographic line of the through-view area analysis area into a 1 area and a 2 area and dividing the transverse topographic line of the through-view area analysis area into a 3 area and a 4 area on the regular grid DEM by taking the observation point as a boundary.
3. The OSP spatial reference line through-the-horizon analysis method based on regular grid DEM as claimed in claim 1, wherein said 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.
4. The OSP spatial reference line through-the-horizon analysis method based on a regular grid DEM according to claim 2 or 3, wherein 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 a 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;
and 2.3, constructing an OSP space coordinate system by taking a direction k as a horizontal axis, a proximity p as a vertical axis and an axis slope g as a vertical axis on the basis of each constructed station center coordinate system, wherein the direction k is a tangent value of an investigation point in the station center coordinate system relative to a direction α of the horizontal axis, the proximity p is an inverse number of a horizontal axis x value of the investigation point in the station center coordinate system, and the axis slope g is a product of a relative elevation z of the investigation point and the proximity p.
5. The OSP spatial reference line visibility domain analysis method based on regular grid DEM as claimed in claim 4, wherein in 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 solved0Recording grid points on the initial reference line as a through view;
step 4.2: moving to a next proximity plane according to the sequence of the proximity of the grid points in the 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 L0By Δ G0Judging 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 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 updated0The 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 formulaiBy Δ GiJudging 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 not in sight, 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.
6. The OSP spatial reference line through-the-horizon analysis method based on regular grid DEM of claim 5, wherein the 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＝(V1+V2)·(V3+V4) (13)
wherein, V1、V2、V3、V4The 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.
7. The OSP spatial reference line through-the-horizon analysis method based on regular grid DEM of claim 5, wherein the 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＝V1+V2+V3+V4(14)
wherein, V1、V2、V3、V4The 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.
8. The OSP spatial reference line through-the-horizon analyzing method based on regular grid DEM of claim 5, wherein the parameter a is derived as follows:
a＝u(zB-zA)
wherein u is the reciprocal of the grid width of DEM, and zA、zBThe vertical axis values of the station coordinates of any two adjacent grid points A, B of the same proximity, respectively, and the direction of B is greater than the direction of a.
9. The OSP spatial reference line through field analysis method based on a regular grid DEM of claim 5, wherein the axis slope difference accumulation formula is:
ΔGi＝ΔGi-1+(aCi-aRi)(ki-ki-1) (11)
wherein, Δ GiIs the axial slope difference, Δ G, between the current grid point and the current reference line segment Li-1Is the difference in the axial slope of the last segment end and the last reference segment, aCiAnd aRiThe slope, k, of the line segment on the reference line and the investigation line in the current segmenti、ki-1Respectively the current segment and the last segment.
10. The OSP spatial reference line through-the-horizon analysis method based on a regular grid DEM as claimed in claim 9, wherein the direction k' of intersection between the OSP line segments and the reference line is found in step 4.12 by:
k′＝ki-ΔGi/(aCi-aRi) (12)。
CN202010054396.3A 2020-01-17 2020-01-17 OSP space reference line through-view domain analysis method based on regular grid DEM Active CN111325791B (en)

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 true CN111325791A (en) 2020-06-23
CN111325791B 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)

Cited By (2)

* Cited by examiner, † Cited by third party
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
CN113313826A (en) * 2021-05-27 2021-08-27 南京邮电大学 Scale-adaptive visual domain analysis method

Citations (5)

* Cited by examiner, † Cited by third party
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
US20160323535A1 (en) * 2014-01-03 2016-11-03 Hangzhou Hikvision Digital Technology Co., Ltd. Method and apparatus for extracting surveillance recording videos
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

Patent Citations (5)

* Cited by examiner, † Cited by third party
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
US20160323535A1 (en) * 2014-01-03 2016-11-03 Hangzhou Hikvision Digital Technology Co., Ltd. Method and apparatus for extracting surveillance recording videos
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)

* Cited by examiner, † Cited by third party
Title
PROF WM RANDOLPH FRANKLIN ETAL.: "Geometric Algorithms for Siting of Air Defense Missile Batteries", 《HTTPS://WWW.RESEARCHGATE.NET/PUBLICATION/2808061》 *
YE ZHI ETAL.: "An Improved Algorithm for Computing Viewshed Based on Reference Planes", 《AN IMPROVED ALGORITHM FOR COMPUTING VIEWSHED BASED ON REFERENCE PLANES》 *
余文广等: "基于大地坐标的通视性检查算法", 《计算机仿真》 *
尚金光等: "基于DEM实现城市高斯投影变形可视化分析", 《城市勘测》 *
张斌等: "基于DEM的地形可视域分析关键技术", 《计算机科学》 *
赵永刚等: "雷达探测范围地形可视域分析方法", 《雷达探测范围地形可视域分析方法》 *

Cited By (3)

* Cited by examiner, † Cited by third party
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
CN113313826A (en) * 2021-05-27 2021-08-27 南京邮电大学 Scale-adaptive visual domain analysis method
CN113313826B (en) * 2021-05-27 2022-08-30 南京邮电大学 Scale-adaptive visual domain analysis method

Also Published As

Publication number Publication date
CN111325791B (en) 2023-01-06

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
CN103324916B (en) Vehicle-mounted and aviation LiDAR data method for registering based on building profile
CN111398918B (en) Radar detection capability analysis method under complex mountain environment
CN110111274B (en) Method for calibrating exterior orientation elements of satellite-borne push-broom optical sensor
CN111325791B (en) OSP space reference line through-view domain analysis method based on regular grid DEM
AU2013237642A1 (en) Multiple Viewshed Analysis
CN109345617B (en) Chain type high-precision splicing and adjustment method based on long-strip multi-station point cloud
CN113514041B (en) Engineering construction project multi-measurement-in-one data acquisition and library building method
CN111899332A (en) Overhead transmission line three-dimensional design method based on oblique photogrammetry technology
CN110827405A (en) Digital remote sensing geological mapping method and system
Liu et al. An improved line-of-sight method for visibility analysis in 3D complex landscapes
CN111832582A (en) Method for classifying and segmenting sparse point cloud by using point cloud density and rotation information
CN116011291A (en) Slope stability digital evaluation platform for tunnel portal multisource information fusion
Sun et al. Building displacement measurement and analysis based on UAV images
CN106291671A (en) A kind of automatic troubleshooting method of stereo observing system based on satellite image data
Schunert et al. Grouping of persistent scatterers in high-resolution SAR data of urban scenes
Sun et al. An adaptive cross-section extraction algorithm for deformation analysis
CN111681313B (en) Space vision analysis method based on digital topography and electronic equipment
CN109166084A (en) A kind of SAR geometric distortion quantitative simulation method based on neighbor point gradient relation
US11361502B2 (en) Methods and systems for obtaining aerial imagery for use in geospatial surveying
Harshit et al. Analysis of survey approach using UAV images and lidar for a Chimney study
CN115018973A (en) Low-altitude unmanned-machine point cloud modeling precision target-free evaluation method
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
CN110909448B (en) High-frequency sky wave return scattering ionization diagram inversion method

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