CN113066161A - Modeling method of urban radio wave propagation model - Google Patents

Modeling method of urban radio wave propagation model Download PDF

Info

Publication number
CN113066161A
CN113066161A CN202110269609.9A CN202110269609A CN113066161A CN 113066161 A CN113066161 A CN 113066161A CN 202110269609 A CN202110269609 A CN 202110269609A CN 113066161 A CN113066161 A CN 113066161A
Authority
CN
China
Prior art keywords
point cloud
primitive
point
modeling
model
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
CN202110269609.9A
Other languages
Chinese (zh)
Other versions
CN113066161B (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.)
Wuhan University WHU
Original Assignee
Wuhan University WHU
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 Wuhan University WHU filed Critical Wuhan University WHU
Priority to CN202110269609.9A priority Critical patent/CN113066161B/en
Publication of CN113066161A publication Critical patent/CN113066161A/en
Application granted granted Critical
Publication of CN113066161B publication Critical patent/CN113066161B/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
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/005—Tree description, e.g. octree, quadtree
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/04—Architecture, e.g. interconnection topology
    • G06N3/045—Combinations of networks
    • 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/10—Constructive solid geometry [CSG] using solid primitives, e.g. cylinders, cubes
    • 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/20—Finite element generation, e.g. wire-frame surface description, tesselation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/11—Region-based segmentation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/10—Segmentation; Edge detection
    • G06T7/187—Segmentation; Edge detection involving region growing; involving region merging; involving connected component labelling
    • 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/10028—Range image; Depth image; 3D point clouds

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Software Systems (AREA)
  • Computer Graphics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • General Health & Medical Sciences (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Molecular Biology (AREA)
  • Computing Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Linguistics (AREA)
  • Evolutionary Computation (AREA)
  • Biophysics (AREA)
  • Biomedical Technology (AREA)
  • Artificial Intelligence (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Processing Or Creating Images (AREA)

Abstract

The invention discloses a modeling method of an urban radio wave propagation model, which comprises the following steps: a: acquiring a laser point cloud picture and an optical image of a modeling area, and establishing a dense point cloud picture containing electromagnetic material information according to a laser radar and the optical image; b: dividing the point cloud in the dense point cloud picture according to the corresponding category of the electromagnetic material information of each point in the dense point cloud to obtain an entity formed by the point cloud; c: identifying corresponding object characteristics based on the geometric structure, and fitting a solid model by using the object characteristics; d: a plurality of transmitters and receivers are arranged in a three-dimensional city model formed by a solid model, and the receiving intensity of different positions corresponding to the transmitters is calculated by utilizing a ray tracing algorithm. By applying the method, the very accurate urban radio wave propagation model suitable for electromagnetic communication simulation can be obtained.

Description

Modeling method of urban radio wave propagation model
Technical Field
The invention relates to the technical field of electromagnetic environment simulation, in particular to a modeling method of an urban radio wave propagation model.
Background
The urban electromagnetic radiation environment is an important environmental element, and in recent 20 years, the development and application of electromagnetic technology are increasingly wide, and the influence on the surrounding environment, human health and the like is increasingly serious. In addition, in order to obtain the prediction of the field strength or the signal strength under the condition of some parameters of the known wireless system, such as frequency, an environment model, a transmitting antenna, and the like, and solve the radio wave propagation and power coverage problems, a modeling process needs to be performed on the urban radio wave propagation model, so how to obtain the urban radio wave propagation model is an urgent technical problem to be solved.
In general, it is difficult for an actual propagation environment to have an analytic method for solving the electromagnetic field. In the prior art, the traditional propagation modeling methods are an empirical model and a theoretical model. The empirical model is based on a large number of experimental measurements and mainly comprises a Hata-0kumura model and a COST-Hata model; theoretical models are relatively simple, such as Friis model, Over-rofoop diffraction model. Both of these modeling methods are relatively simple and computationally fast, however the main disadvantages of empirical models are that they are generally not general, regional and that the predicted urban environment is similar to the urban environment from which the model was derived. Moreover, for a MIMO (multiple Input multiple Output) system, the empirical model cannot provide accurate MIMO simulation key parameters, such as space-time and angular delay parameters, and the like.
Disclosure of Invention
The technical problem to be solved by the invention is how to establish an urban electromagnetic environment model suitable for an electromagnetic communication simulation scene.
The invention solves the technical problems through the following technical scheme:
the invention provides a modeling method of an urban radio wave propagation model, which comprises the following steps:
a: acquiring a laser point cloud picture and an optical image of a modeling area, and establishing a dense point cloud picture containing electromagnetic material information according to a laser radar and the optical image;
b: segmenting the dense point cloud graph based on corresponding electromagnetic material information, wherein the electromagnetic material information comprises: plants, glass, concrete, ground;
c: identifying corresponding object features based on the geometric structure, and fitting a solid model by using the object features, wherein the object features comprise: wall, door, roof, protrusion and window, and the geometry includes: size, position, orientation, and topology.
D: a plurality of transmitters and receivers are arranged in a three-dimensional city model formed by a solid model, and the receiving intensity of different positions corresponding to the transmitters is calculated by utilizing a ray tracing algorithm.
Optionally, the specific implementation of step a includes:
acquiring a laser point cloud picture scanned by a laser radar in a modeling area and electromagnetic material information of each point, and adding the electromagnetic material information into a label of the point;
acquiring an optical image of a modeling area based on a full-band multi-view optical camera; identifying object features based on an algorithm of a motion acquisition structure, and constructing a sparse point cloud corresponding to the object features;
acquiring a depth point cloud image corresponding to the optical image;
identifying and calibrating electromagnetic material information of a corresponding region in the depth point cloud picture by using a pre-trained neural network model;
and fusing the calibrated depth point cloud picture, the calibrated sparse point cloud picture and the laser point cloud picture added with the label according to the geometric mapping relation, and adding corresponding object structure characteristics into the label to obtain the dense point cloud picture.
Optionally, the identifying and calibrating the electromagnetic material information of the corresponding region in the depth point cloud chart by using the pre-trained neural network model includes: :
the method comprises the steps of training a pre-built neural network model by using an optical image of an object marked with electromagnetic material information as a sample to obtain a converged neural network model, then segmenting the optical image of a modeling area by using the neural network model, and calibrating a corresponding area in a depth point cloud image by using a segmentation result.
Optionally, the process of acquiring the laser spot cloud image includes:
the modeling area is scanned by using the laser radar, and based on the three-dimensional coordinates of the laser radar, the modeling area is modeled by using a formula,
Figure BDA0002973716130000031
scanning the modeling area by using a laser radar to obtain a laser point cloud picture of the modeling area, wherein,
Xpcoordinates on the x-axis for points in the point cloud; xGCoordinates of the laser radar on an x axis; y ispCoordinates on the y-axis for points in the point cloud; y isGCoordinates of the laser radar on the y axis are obtained; s is a laser radar point-to-point vector; theta is an included angle between the pixel corresponding to the laser ranging point P and the middle pixel in the scanning period; b ═ cos ω · sin α · cos κ + sin κ · sin ω; aerial attitude parameters of laser radar such as roll angle, pitch angle and yaw angle
Figure BDA0002973716130000032
Obtaining photogrammetric attitude angles (alpha, omega, kappa) through coordinate system transformation; zpHeight of a point in the point cloud; zGIs the height of the lidar.
Optionally, the method further includes, before step B, performing outlier removal processing on the dense point cloud graph.
Optionally, the step B of dividing the dense point cloud image includes:
and partitioning the dense point cloud by a point cloud partitioning algorithm based on region growing.
Optionally, the fitting the physical model by using the object features in the step C includes:
determining whether the corresponding entity belongs to undulating terrain, flat ground, simple buildings, complex buildings and trees according to the characteristics of the object;
building modeling is carried out based on a global fitting algorithm under the condition that the entity belongs to a simple building;
under the condition that the entity belongs to a tree, performing data modeling by using a tree skeleton model constructed in advance;
and under the condition that the entity belongs to a complex building or an undulating terrain, carrying out data modeling by utilizing a Poisson reconstruction algorithm.
Under the condition that the entity belongs to a flat ground, modeling by adopting a least square fitting plane method;
and taking the model obtained by modeling as a solid model.
Optionally, the building modeling based on the global fitting algorithm includes:
establishing constraints on the characteristics of the building model object, wherein the constraints comprise: area, position, direction, topology;
the entity is divided by utilizing a random sampling consistency algorithm to obtain a geometric primitive set { chi } composed of local geometric primitivesi}；
According to the geometric primitive set { χiThe parallel and orthogonal relations between every two geometric primitives contained in the method obtain the combined relation between the geometric primitives, and the maximum non-conflict set is selected from the set formed by the combined relation between the geometric primitives
Figure BDA0002973716130000041
Aligning each geometric primitive according to the combination relation in the maximum non-conflict set by utilizing a nonlinear optimization algorithm with constraints;
for aligned maximum non-conflicting set
Figure BDA0002973716130000042
Establishing a plurality of geometric primitive pairs according to the same angle, and establishing a relation by taking the geometric primitive pairs as vertexesDrawing;
calculating the space distance between every two geometric primitive pairs in the relational graph, and deleting the geometric primitive pairs of which the space distances to other geometric primitive pairs are larger than a set distance;
for each geometric primitive pair, connecting the geometric primitive pair with other geometric primitive pairs by using edges, wherein the geometric primitive pair and the other geometric primitive pairs have similar angles; at the same time, by using the formula,
Figure BDA0002973716130000051
a constraint is added to the newly added edge, wherein,
sca confidence score for an edge;
Figure BDA0002973716130000052
the included angle between the direction vector of the geometric primitive i and the direction vector of the geometric primitive j is shown;
Figure BDA0002973716130000053
the included angle between the direction vector of the geometric primitive k and the direction vector of the geometric primitive l is shown; gcA constraint that is an edge; n isiIs a geometric primitive i; n isjIs a geometric primitive j; n iskIs a geometric primitive k; n islIs a geometric primitive l;
taking the newly added edge set as an initial candidate set, and extracting an isogonal relation set consisting of edges with isogonal relation in the initial candidate set;
by means of the formula (I) and (II),
Figure BDA0002973716130000054
the set of angular relationships is optimized, wherein,
min is a minimum evaluation symbol; sigma is a summation symbol; ed(Pi,χi) As a set of points PiAnd primitive xiThe error between; w is apd2(p,χi) Is the accumulated error of the data; w is apIs the weight of the geometric primitive; d2(p,χi) Distance of p points to the geometric primitive;
judging whether each side in the optimized set is the same as the side before optimization;
if yes, deleting the primitive pair corresponding to the edge with the lowest confidence score from the same edges before and after the same-angle relation set optimization;
if not, returning to the step of executing the optimization of the corresponding angle relation set until the optimization times reach the set times;
acquiring a set of primitive pairs corresponding to the edge with the lowest confidence score, and aiming at each primitive pair in the set, utilizing a formula,
Figure BDA0002973716130000061
computing a confidence score for the primitive pair and a constraint, wherein,
spis the confidence score of the primitive pair; p is a radical ofiPoints on each axis of the geometric primitive i; p is a radical ofjPoints on each axis of the geometric primitive j; | | is a modulo symbol; (p)i-pj)TIs the transposition of the vector; gpConstraint of the primitive pair;
carrying out coaxial alignment processing on the primitive pairs in the maximum relation subset according to the direction of a normal vector between each primitive pair in the maximum relation subset, and sorting the primitive pairs in the optimized isometric relation set according to the confidence scores of the primitive pairs to obtain a corresponding maximum relation subset; and obtaining a three-dimensional model of the building according to the relationship between the entities corresponding to the maximum relationship subset.
Optionally, in a case that the entity belongs to a complex building, performing data modeling by using a poisson reconstruction algorithm includes:
taking the point cloud data corresponding to the fitted building entity as input, and determining an indication function of the point cloud data according to the vector direction of each point cloud;
after smooth filtering is carried out on the indicating function, a gradient field of the indicating function is calculated, and the gradient field is discretized by establishing an implicit function of an octree for adjusting the depth of a grid according to the density of point cloud data;
carrying out segmentation sampling on the discretized gradient field, and calculating a vector field of a sampling point by using a cubic linear interpolation method;
the poisson equation is solved and the solution,
Figure BDA0002973716130000062
a scalar indicator function of the point cloud data is obtained, wherein,
Δ is laplace operator; x is an indicator function;
Figure BDA0002973716130000063
is a vector differential operator;
Figure BDA0002973716130000064
a space vector of the point cloud data;
based on a scalar indication function, extracting an isosurface by adopting a mobile cube algorithm, and establishing a three-dimensional model of the building based on the isosurface.
Optionally, in the step D, in the three-dimensional city model formed by the solid model, there are the following three ray propagation modes:
direct ray tracing: judging whether a line segment between the transmitter and the receiver has an intersection point with a building or not, if not, a direct path exists;
tracking by reflected rays: firstly, making a plane-related mirror image point for a transmitter, then connecting the mirror image point with a receiver, connecting a line segment consisting of the two points with a plane equation, solving an intersection point with the plane, if the intersection point exists, then, a reflection point exists, and connecting the transmitter, the reflection point and the receiver to obtain a reflection ray path;
tracing by diffracted rays: firstly, obtaining diffraction points on a straight line where split edges of a transmitter and a receiver are located, and then judging the positions of the diffraction points; if the edge is outside the cleaved edge, no diffraction exists; when the diffraction point is positioned on the splitting edge, the path where the transmitter and the diffraction point are positioned is subjected to shielding judgment; if the shielding does not exist, then shielding judgment is carried out on the path where the diffraction point and the receiver are located; if there is no occlusion, it indicates that a diffracted ray path exists.
Compared with the prior art, the invention has the following advantages:
by applying the embodiment of the invention, when the urban three-dimensional model is established, the urban three-dimensional point cloud data is obtained by selecting the multi-source point cloud data combining optical image matching and laser scanning to carry out matching fusion, and the established urban three-dimensional model also comprises the electromagnetic medium types of all parts of an object, so that the high-precision urban three-dimensional model suitable for electromagnetic calculation is obtained, and the urban radio wave propagation model suitable for electromagnetic communication simulation can be obtained by utilizing a ray tracing method based on the three-dimensional model.
In addition, because the data of the laser radar is used as supplement of the optical image data, the established three-dimensional model of the city has high precision, and therefore, the embodiment of the invention can quickly and accurately establish the three-dimensional model of the city; finally, the optical image acquisition cost is lower, and the acquisition is more convenient, so that the urban three-dimensional model can be established more quickly.
Drawings
Fig. 1 is a schematic flow chart of a modeling method of an urban radio wave propagation model according to an embodiment of the present invention;
fig. 2 is a schematic diagram illustrating a modeling method of a city radio wave propagation model according to an embodiment of the present invention;
FIG. 3 is a schematic diagram of a process for obtaining a dense point cloud according to the present invention;
FIG. 4 is a schematic diagram of a dense point cloud of a modeled region obtained by an embodiment of the present invention;
FIG. 5 is a simplified three-dimensional model of a building reconstructed according to an embodiment of the present invention;
fig. 6 is a three-dimensional model of a complex building reconstructed according to an embodiment of the present invention.
Detailed Description
The following examples are given for the detailed implementation and specific operation of the present invention, but the scope of the present invention is not limited to the following examples.
Fig. 1 is a schematic flow chart of a modeling method of a urban radio wave propagation model according to an embodiment of the present invention, and as shown in fig. 1, the method includes:
s101: and acquiring a laser point cloud picture and an optical image of the modeling area, and establishing a dense point cloud picture containing electromagnetic material information according to the laser radar and the optical image.
Fig. 2 is a schematic diagram illustrating a modeling method of a metro wave propagation model according to an embodiment of the present invention, as shown in fig. 2,
firstly, scanning a modeling area by using an airborne laser radar, and acquiring a high-precision image of the earth surface or a building according to a signal reflected after a laser pulse penetrates by using the characteristic that the laser pulse can partially penetrate through trees to shield; meanwhile, the airborne global positioning unit can obtain the accurate coordinates of the unmanned aerial vehicle, the airborne inertial measurement unit can measure the motion state of the unmanned aerial vehicle, and then laser point cloud pictures of all objects such as trees, the ground, buildings and the like in a modeling area are calculated according to the reflected laser pulse signals. In the embodiment of the invention, the scanning mode of the laser radar is different from the omnidirectional scanning in the prior art for acquiring the point cloud images of all surfaces of the entity, and in the embodiment of the invention, a unidirectional scanning mode or a reciprocating scanning mode is used for acquiring only the point cloud images of one surface or two surfaces and other parts of surfaces of the entity without acquiring all the laser point cloud images of all the surfaces.
When the laser radar is used for measurement, the airborne global positioning system can measure the three-dimensional coordinate (X) of the laser radar in the airG，YG，ZG) The airborne inertia measurement unit can measure the pitch angle, the roll angle and the course angle of the unmanned aerial vehicle
Figure BDA0002973716130000093
The laser pulse emitted by laser radar directly measures the topographic relief condition, and the projection center of laser radar is projected to the unknown point P (X) on the groundG，YG，ZG) The vector S of (a) can be accurately measured, and therefore, the modeling area can be scanned using the lidar, based on the three-dimensional coordinates of the lidar, using a formula,
Figure BDA0002973716130000091
scanning the modeling area by using a laser radar, further measuring the coordinates of each point, and obtaining a laser point cloud picture of the modeling area from each point, wherein,
Xpcoordinates on the x-axis for points in the point cloud; xGCoordinates of the laser radar on an x axis; y ispCoordinates on the y-axis for points in the point cloud; y isGCoordinates of the laser radar on the y axis are obtained; s is a laser radar point-to-point vector; theta is an included angle between the pixel corresponding to the laser ranging point P and the middle pixel in the scanning period;
b ═ cos ω · sin α · cos κ + sin κ · sin ω; aerial attitude parameters of laser radar such as roll angle, pitch angle and yaw angle
Figure BDA0002973716130000094
Obtaining the attitude angle of photogrammetry by coordinate system transformation
Figure BDA0002973716130000092
ZpHeight of a point in the point cloud; zGIs the height of the lidar.
Then, classifying each point in the laser point cloud picture by using a pre-trained neural network model to obtain corresponding electromagnetic material information of each point, and then adding the electromagnetic material information into a label of the point to further obtain the laser point cloud picture containing the electromagnetic material information, wherein the electromagnetic material information refers to the electromagnetic material information of a point corresponding body in the laser point cloud, such as glass, a wall body, the ground, trees and the like. In practical application, the electromagnetic material information of the point cloud in the laser point cloud map of the area with the same type as the modeling area or a part of the modeling area, such as one tenth or one hundredth area, can be marked manually, then the marked laser point cloud map is used for training a neural network model, and then the trained neural network model identifies the electromagnetic material information of each point in the laser point cloud map corresponding to other unmarked modeling areas.
Then, semantic recognition operation of reconstructing sparse point cloud, depth map corresponding to the optical image and optical image is respectively carried out, fig. 3 is a schematic diagram of the process for acquiring dense point cloud provided by the invention, as shown in fig. 3,
firstly, a full-band optical image of a modeling area can be shot by utilizing a five-eye lens carried on an airborne or unmanned aerial vehicle cluster; then, identifying object features based on an algorithm of a motion acquisition structure, and constructing a sparse point cloud corresponding to the object features by utilizing the prior art, wherein the object features comprise: one or a combination of walls, doors, corner points, roofs, protrusions, and windows. The process is the prior art, and comprises the following specific processes: extracting metadata of the image, extracting key features of the image (such as corner points of a building target), matching and tracking the features, and reconstructing the sparse point cloud.
The depth point cloud map corresponding to the whole optical image is obtained, for example, depth information may be calculated for each unmanned aerial vehicle image based on the estimated acquisition point position, and a two-dimensional optical image acquired by an unmanned aerial vehicle may be transformed into a three-dimensional space by using the depth information. Further, can install five-eye optical camera and laser radar on same unmanned aerial vehicle, and then can utilize the unmanned aerial vehicle position that laser radar measured to carry out the generation of degree of depth point cloud picture.
The method comprises the steps of utilizing a pre-trained electromagnetic medium recognition neural network model to segment an optical image of a modeling area, obtaining electromagnetic material information corresponding to objects, such as windows, walls, roofs, the ground and the like, contained in the optical image, and adding the identified electromagnetic material information into labels of corresponding points in a depth point cloud picture.
In practical application, optical images in some urban areas can be shot in advance, electromagnetic materials of different parts of objects in the images are labeled to form training samples, the samples are used for training a preset deep structure neural network model until the neural network model converges to obtain a trained model, and then the neural network model is used for carrying out semantic segmentation on the optical images in the modeling areas.
Then, the sparse point cloud, the depth point cloud image and the calibrated laser point cloud image are fused according to the geometric relation mapping, and the specific process is as follows: firstly, calculating accurate mapping of space geometric relations among different point cloud sets, and solving coordinate conversion parameters; and then carrying out rigid body transformation on the data set to be converted. The embodiment of the invention uses an Iterative Closest Point (ICP) algorithm based on the initial value of local feature matching to perform image fusion. The ICP algorithm is an iterative algorithm for realizing point cloud fusion by global registration, and in order to further improve the precision of the point cloud fusion, the multi-point cloud fusion is realized by using a point cloud rotation matrix constructed by using seven parameters or quaternions obtained by identifying and matching local features such as point-line and the like. And then adding the corresponding labels made of the electromagnetic materials to the corresponding points to obtain the dense point cloud picture added with the labels.
In the process of fusing the laser point cloud picture and the optical image, the depth information of each pixel point is calculated for each unmanned aerial vehicle image based on the estimated acquisition point position by using the multi-view geometry principle, the two-dimensional optical image acquired by the unmanned aerial vehicle can be converted into a three-dimensional space by using the depth information, and then an image loaded with electromagnetic material information, namely an electromagnetic material coefficient, is introduced into the reconstruction. According to the embodiment of the invention, the depth information in the dense point cloud is used as a medium, the idea of depth information sharing is utilized to convert the pixels of the optical picture of the unmanned aerial vehicle loaded with the electromagnetic material information into a three-dimensional space, so that the electromagnetic material information can be loaded into the point cloud to obtain the point cloud containing the electromagnetic material information, and each point in the dense point cloud picture obtained in the way also comprises the electromagnetic material information besides the position and color information.
Fig. 4 is a schematic diagram of the dense point cloud of the modeling area obtained in the embodiment of the present invention, and as shown in fig. 4, the above-described processing is performed on each image, so that the dense 3d point cloud of the entire area based on the optical image reconstruction can be finally obtained, and then the laser point cloud data obtained by the unmanned aerial vehicle cluster is obtained and processed.
Because the dense point cloud includes points on complex ground features such as earth surfaces, trees and the like, the points are distributed in a three-dimensional space in a discrete and irregular manner, and great difficulty is caused to the application of generating a three-dimensional model of an urban area and ground point data. In order to eliminate noise points in the data and further reduce the influence of the noise points on the modeling of the three-dimensional model, before the step S102 is executed, the noise point data needs to be eliminated by performing outlier processing on the acquired dense point cloud map. The processing algorithm needs to select a proper classification threshold according to different terrain conditions, but it is not easy to select a parameter capable of removing the ground features with different sizes for a complex ground surface, so that it is very necessary to individually set the corresponding parameter according to specific ground features and terrain conditions. However, it is not practical if different parameters are artificially set for specific land features and landforms, respectively, and in the embodiment of the present invention, typical characteristics of an urban environment and timeliness of data processing requirements are comprehensively considered, and an outlier removing method based on statistics is selected to process point clouds. The average distance between each point and adjacent points in the point cloud plus the deviation value pre-defined by the user is used as a threshold, namely, if the number of the adjacent points contained in a sphere with the radius of the average distance plus the deviation value of a certain point in the dense point cloud is less than the defined threshold, the point is considered as an outlier and is removed. According to the embodiment of the invention, under the condition of not depending on ground objects and landforms too much, simpler parameters can be set, and ground points in mass point cloud data can be automatically extracted, so that the removal efficiency of noise points is improved.
S102: segmenting the dense point cloud graph based on corresponding electromagnetic material information, wherein the electromagnetic material information comprises: one or a combination of plants, glass, concrete and ground.
In the process of segmenting dense point clouds in urban environments of modeling areas to obtain point clouds of each object in the urban environments, useless information such as vehicles exists near the point clouds of the buildings due to the fact that the environments of the buildings in the urban environments are complex, the point cloud distribution corresponding to the useless information is often located near other electromagnetic media, and therefore the difficulty in removing the useless information is large. The method can be realized by directly segmenting the point cloud by using a point cloud segmentation mode based on the volume element and a point cloud segmentation mode based on a multi-vision picture.
However, the inventors found that the voxel-based method is too resource-consuming, and the way of directly segmenting the point cloud is simpler than the voxel-based segmentation method, but still consumes a large amount of computing resources. Although the reconstruction method based on multi-view reconstruction is simple in calculation, the disadvantages are also obvious: in the process of projecting the 3d structure to 2d, relatively serious structure information loss exists, and in the process of projecting the 3d structure to 2d information, the selection of the projection position, namely the position of the virtual camera, has a great influence on the reconstruction result, so that the consistency of the result based on multi-view reconstruction is poor.
In order to solve the problems, each point in the dense point cloud picture comprises a piece of label information besides position and color information, the label information comprises the characteristics of the electromagnetic material information corresponding to the point, and the electromagnetic material information of different types of objects has larger difference, so point cloud segmentation can be performed by using the label information. And then, according to the distribution rules of buildings, trees, the ground, lakes and the like in the city, processing the point clouds of the same category by using a point cloud segmentation algorithm based on a region growing algorithm to obtain a segmented point cloud picture of the object.
S103: identifying corresponding object features based on the geometric structure, and fitting a solid model by using the object features, wherein the object features comprise: one or a combination of a wall, a window, a roof, a tree, a ground surface, and the geometry comprises: one or a combination of size, position, orientation, and topology.
Firstly, according to the segmented point cloud images of each object, judging which entity belongs to a simple building, a complex building, a tree, or an undulating terrain, a flat ground (in the case that the entity belongs to the flat ground, a least square fitting plane method is adopted for data modeling, the method is the prior art, and the invention is not described in the invention) according to the object characteristics. For example, a pre-trained deep learning model can be used for point cloud semantic segmentation, point cloud is obtained through reconstruction, each point of the point cloud contains semantic category information, and then entities are classified according to the category information.
And under the condition that the entity belongs to a simple building, building modeling is carried out based on a global fitting algorithm.
And under the condition that the entity belongs to the tree, performing data modeling by using a tree skeleton model constructed in advance. The tree skeleton is the basis for building a three-dimensional tree model, and the basic idea is to extract an initial skeleton of a tree from ground laser point cloud, further optimize the initial skeleton and finally realize the extraction of the tree skeleton. Extracting 'similar trunk points' of the trees from the preprocessed original tree point cloud according to the difference of the laser reflection intensity values, organizing the 'similar trunk points' and the rest 'non-trunk points' by adopting a 'minimum spanning tree' algorithm in a graph theory to form an initial skeleton of the trees, and finally further optimizing the initial skeleton of the trees through 'point density adjustment' and 'branch smoothing' to obtain the more accurate and reasonably distributed skeleton of the trees.
And when the number of the planes contained in the fitted building entity is greater than the set number, such as 100, judging that the building belongs to the complex building, and performing building modeling by using a Poisson reconstruction algorithm.
Firstly, the building is modeled in three dimensions, and because the point cloud is obtained in the step S102, the difficulty of directly performing electromagnetic calculation is high. Therefore, surface reconstruction needs to be performed on the monomer point cloud according to the electromagnetic material information, and then electromagnetic calculation is performed on the basis of the reconstructed surface so as to reduce the calculation complexity.
The surface reconstruction methods commonly used in the prior art are the dironi triangulation or von neumoniae method. The reconstruction result can basically reflect the point cloud position and has higher precision, the surface reconstruction mode has wide application and can be suitable for the rapid modeling of almost all point clouds, however, the number of the triangular surfaces contained in the three-dimensional model reconstructed by the method is usually tens of thousands or even tens of thousands, and if the triangular surfaces are used as the direct input of the ray tracing model, the calculation time complexity is very high. At least hundreds of buildings are often distributed in urban areas, which further increases the amount of calculation and even if calculation is possible, the cost is very high.
In order to solve the problems, the inventor finds that the buildings in the modeling area belong to objects with regular structures, and particularly the buildings have larger occupation ratio in the urban environment modeling process; therefore, the building model can be reconstructed by adopting a building facade reconstruction algorithm based on the prior knowledge. The building facade comprises important object characteristics such as electromagnetic reflecting surfaces such as wall surfaces, roofs and protrusions, windows and the like, and each object characteristic has geometric characteristic constraint which is different from other object characteristics, such as one of area, position, direction and topological structure:
size: by area size, the wall surface is easily separated from other features, and some small noise partitions are easily filtered out.
Position: specific features may only be present at specific locations, such as windows and doors on walls and protrusions, while roofs are always above walls.
The direction is as follows: the wall surface is mostly vertical to the ground, the roof is not vertical to the ground, and the wall surface can only be horizontal or inclined: it is also possible to add a floor as an auxiliary feature for better identification of the main wall surface.
Topological structure: the topology between features is often an important cue, such as a wall always intersecting the ceiling, and a roof always intersecting a wall.
Dot density: since the laser beam is reflected to obtain the three-dimensional coordinates of the measured object, the characteristics can be distinguished by the point density.
The method can add the category to each entity in the point cloud through semantic segmentation based on the electromagnetic medium, but in order to more accurately identify the facade elements, each category of semantic features has unique attributes, and the semantic features can be automatically, efficiently and accurately extracted through the attribute regularization into specific feature constraints.
Since artificial engineering elements contain basic geometrical features such as planes, cylinders, spheres and cones. Segmenting the point cloud of the entity by using a random sampling consistency RANSAC algorithm to obtain the basis of the geometrical characteristics contained in the point cloudElements, the elements contained by all entities constituting a set of geometric elements { χ }i}. Then, global constraints such as direction relations, displacement relations, equality relations and the like contained in the geometric primitive set are extracted step by step so as to improve the accuracy of the entity model and the universality of the algorithm. In the process of model reconstruction, the relationship between three types of primitives is mainly concerned, which is also called as alignment relationship: directional relationships, such as parallel and orthogonal, permutation relationships, such as co-planar or coaxial presence between elements; there is an equality relationship between the primitives. Based on the definition of the alignment relationship, the alignment global direction relationship is obtained by optimizing the model parameters from the global consideration, and the alignment global direction relationship can be accurately matched with the original data. Taking into account the sensitivity of different types of relationships to noise, a stepwise extraction strategy is adopted, i.e. the relationship which is easiest to process is prioritized each time. In the embodiment of the invention, the relationship between the element pairs can be well recovered by establishing the relationship graph and the constraint of the data item for the global constraint relationship.
The appearance of the urban artificial building is relatively regular, and objects such as wall surfaces, ground surfaces, roofs and the like are basically in parallel or orthogonal relation. According to the geometric primitive set { χiThe parallel and orthogonal relations between every two geometric primitives contained in the structure are obtained to obtain the combined relation between the geometric primitives, namely a plurality of candidate relation sets C0(ii) a From C0Extracting the largest non-conflict set
Figure BDA0002973716130000161
Associating primitives with relationships using a constrained nonlinear optimization algorithm
Figure BDA0002973716130000162
And data alignment.
Then, alignment of the same direction primitives is performed: in the building of the urban scene, a plurality of regular structures exist, and the regularity often causes equal angles to exist between the elements, and the elements with the same angles form element pairs. While any angular relationship necessarily contains a pair of primitives, i.e., four primitives. Therefore, the relationship graph G is established by taking the geometric primitive pairs as the verticese。
Then, the user can use the device to perform the operation,relationship graph GeEach vertex of (a) represents a primitive pair out of order. In noisy data, spatially distant primitives can produce erroneous relationships. Thus, if the primitive pair { χ }i,χjIf the distance between them exceeds a certain distance, the corresponding node, i.e. vertex, will be deleted. After the above processing, the relationship diagram GeContains only O (m) vertices.
For each geometric primitive pair, connecting the geometric primitive pair with other geometric primitive pairs by using edges, wherein the geometric primitive pair and the other geometric primitive pairs have similar angles; at the same time, by using the formula,
Figure BDA0002973716130000171
a constraint is added to the newly added edge, wherein,
sca confidence score for an edge;
Figure BDA0002973716130000172
the included angle between the direction vector of the geometric primitive i and the direction vector of the geometric primitive j is shown;
Figure BDA0002973716130000173
the included angle between the direction vector of the geometric primitive k and the direction vector of the geometric primitive l is shown; gcA constraint that is an edge; n isiIs a geometric primitive i; n isjIs a geometric primitive j; n iskIs a geometric primitive k; n islIs the geometric primitive/.
According to the above processing rule, the obtained initial candidate set is Ce＝{c1,c2,...}。
Then, alignment of the parallel and orthogonal primitives is performed: similar to the processing method of the direction alignment relationship, the edge C epsilon C is processed step by step according to the order of the scores of the edges from high to lowe. And extracting the primitive pair relation with the equal angle relation to obtain a set of primitive pair relation. Meanwhile, the isometric relationship is transitive, so that the relationship diagram GeIt cannot be a circular graph, and such a process significantly reduces the number of edges in the graph. Further, the extracted equiangular relationshipsThe ensemble may still have conflicting relationships, and the conflicting relationships need to be detected by means of interior point nonlinear programming.
Then the formula is utilized to obtain the final product,
Figure BDA0002973716130000174
the set of angular relationships is optimized, wherein,
min is a minimum evaluation symbol; sigma is a summation symbol; ed(Pi,χi) As a set of points PiAnd primitive xiThe error between; w is apd2(p,χi) Is the accumulated error of the data; w is apIs the weight of the geometric primitive; d2(p,χi) Is the distance of the p point to the geometric primitive. And if the optimization result is not changed, deleting the primitive pair relationship with the lowest score from the optimized pair relationship set. And then optimizing the relation of the rest primitives again, wherein generally, only two times of optimization are needed to obtain a relatively ideal effect.
Then, alignment of the primitives of the permutation relationship is performed: most of urban buildings have coplanar and coaxial parts, and after the directional alignment relationship is processed, the aligned directional relationship is kept unchanged, and then the replacement relationship is processed and aligned. Since the algorithm has already performed a direction-aligned process on the data, if two primitives χiHexix-jIs parallel, then the pair of primitives is likely to be coaxial. By means of the interaxial distance between elements, the fraction s is obtainedpAnd a set C of aligned relationships of the primitives carrying out the permutation relationshippConstraint g of the maximum subset of relationships ofp：
Figure BDA0002973716130000181
Wherein the content of the first and second substances,
spis the confidence score of the primitive pair; p is a radical ofiPoints on each axis of the geometric primitive i; p is a radical ofjPoints on each axis of the geometric primitive j; | | is a modulo symbol; (p)i-pj)TIs the transposition of the vector; gpIs the constraint of the primitive pair. Coplanar surfaceUsually there is a distance in space between pairs of primitives. All cylinders and cones need not be reprocessed due to the extraction via the coaxial relationship. Therefore, the extraction of the coplanar relationship is only to detect the coplanar relationship between the two planar models. For primitive pair χiHexix-jIf their normal vectors are in the same direction, then there is di＝dj(ii) a If the two are reversed then di＝-dj。
And then arranging the extracted relations according to the descending order of confidence scores, and extracting the maximum relation subset of the ordered set
Figure BDA0002973716130000182
And the alignment relationship thereof comprises the three alignment relationships. And finally, carrying out minimum data item error processing to obtain the aligned primitive.
Fig. 5 is a reconstructed simple building three-dimensional model according to an embodiment of the present invention, and as shown in fig. 5, the fitted object features are aligned according to the aligned primitives, so as to obtain a building three-dimensional model.
The building facade (electromagnetic reflecting surface) is fitted using an algorithm of global fitting: given an input point cloud set, the algorithm simultaneously performs local primitive matching and global primitive relationship extraction. The algorithm can finally extract and reconstruct the reflecting surfaces of the buildings in the city, each reflecting surface comprises the electromagnetic material coefficient information of the reflecting surface, and the information is obtained by the electromagnetic medium semantic segmentation based on deep learning. The global fitting algorithm in the embodiment of the invention is the existing GlobFit algorithm.
The three-dimensional reconstruction of the building needs to identify and extract basic geometric structures such as points, lines or surfaces and the like in point cloud data, and the algorithm used by the invention has the following advantages: first, optical image reconstruction and laser point clouds cannot scan a complete building due to the presence of trees or vehicle obscuration. Through the priori semantic knowledge, the defect can be made up, so that the sealing performance of a reconstruction result is ensured; second, depending on the type of semantics, the building model may be given more electromagnetic material information, such as ground, owner, target name or code number, etc., which helps to build a more refined electromagnetic map.
In addition, the embodiment of the invention can directly use the building model and the electromagnetic material information thereof to rapidly acquire the electromagnetic situation of the city, and further can be used for electromagnetic calculation. Furthermore, the influence of absorbers and diffuse scattering objects propagated by electromagnetic waves can be considered in electromagnetic calculation so as to obtain more accurate electromagnetic distribution.
When the fitted building is a complex building or an undulating terrain, the three-dimensional model can be reconstructed by using a poisson reconstruction method, and the reconstruction process of the corresponding three-dimensional model is described below by taking the complex building as an example:
the method of poisson reconstruction is selected for the complex building reconstruction of urban areas. The principle of the method is to convert the surface reconstruction problem of the point cloud into a solving process of a Poisson equation in a space. The poisson equation is used for solving the tone mapping problem of the image in the strong non-static state range at first and also solving the seamless editing problem of the image in the connection area. The typical advantage of poisson's equation is to solve the problem of large scale integrity, which is well-suited for complex architectural point clouds.
The poisson equation under the three-dimensional rectangular coordinate system is as follows:
Figure BDA0002973716130000201
and taking the point cloud data of the urban building as a set S, sampling to obtain a subset S belonging to S, and obtaining an object edge S. The gradient index function is a zero vector field that is nearly ubiquitous, except near the surface, where it is equal to the surface normal. Therefore, the subset of sampling points can be regarded as an index function of a gradient model of the sample, and the problem of calculating the index function is reduced to an inverse gradient operator, i.e. determining the gradient of a scalar function X that most closely approximates the standard vector field. This gradient determines the vector field in dependence on the sample
Figure BDA0002973716130000202
But instead of the other end of the tubeThe gradient and the vector field are not completely fitted, and approximate calculation is needed
Figure BDA0002973716130000203
Adding a non-convergence operator to convert the variation problem into a Poisson problem:
Figure BDA0002973716130000204
obtaining a vector field
Figure BDA0002973716130000205
Then, the Poisson's equation can be solved to obtain a scalar indicator function X, but not an exact solution because of the vector field
Figure BDA0002973716130000206
Usually, the product is not obtained, and a poisson equation needs to be constructed by using a divergence operator to find the optimal approximate least square estimation solution.
Therefore, point cloud data corresponding to the fitted building entity are used as input, and an indicating function of the point cloud data is determined according to the vector direction of each point cloud;
since the indicator function X is a piecewise function, calculating the gradient directly yields infinite values, and, using a formula,
Figure BDA0002973716130000207
after smooth filtering of the indicative function, calculating a gradient field of the indicative function and discretizing the gradient field by establishing an implicit function of an octree adjusting the depth of the mesh according to the density of the point cloud data, wherein,
Figure BDA0002973716130000208
is an object, M is a surface edge, XMIs an object
Figure BDA00029737161300002011
Is used to indicate the function of (a),
Figure BDA0002973716130000209
is a function of a smoothing filter that is,
Figure BDA00029737161300002010
is the inner normal direction of point p. (ii) a
And (3) carrying out segmentation sampling on the discretized gradient field, and calculating a vector field of a sampling point by using a cubic linear interpolation method: dividing discrete points into small curved patches
Figure BDA00029737161300002012
Facet ρSIs approximated by the following equation:
Figure BDA0002973716130000211
thereby obtaining a vector field
Figure BDA0002973716130000212
The poisson equation is solved and the solution,
Figure BDA0002973716130000213
a scalar indicator function of the point cloud data is obtained, wherein,
Δ is laplace operator; x is an indicator function;
Figure BDA0002973716130000214
is a vector differential operator;
Figure BDA0002973716130000215
a space vector of the point cloud data;
based on a scalar indication function, extracting an isosurface by adopting a mobile cube algorithm, and establishing a three-dimensional model of the building based on the isosurface, wherein the established three-dimensional model is shown in FIG. 6.
According to the embodiment of the invention, the dense point cloud result of the whole urban environment is segmented according to the self-adaptive surface reconstruction requirement of the electromagnetic medium, and each ground object is obtained according to the category (buildings, water surface, ground, trees and the like) of the geographic elements; and then, carrying out surface reconstruction on the electromagnetic media of different types to obtain a three-dimensional model for establishing a three-dimensional reconstructed model of a specific region, and further drawing an electromagnetic situation distribution map of the region by using an electromagnetic calculation method.
In addition, according to the characteristics and the timeliness requirements of the urban environment, the multi-source point cloud data fusion mode combining optical image matching and laser scanning is selected to obtain the urban three-dimensional point cloud data.
S104: a plurality of transmitters and receivers are arranged in a three-dimensional city model formed by a solid model, and the receiving intensity of different positions corresponding to the transmitters is calculated by utilizing a ray tracing algorithm.
Illustratively, according to the entity model obtained in step S103, a city three-dimensional model is constructed according to the position of the entity model, and the following three ray propagation modes exist in the city three-dimensional model:
direct ray tracing: it is determined whether there is an intersection between the line segment between the transmitter and the receiver and the building, and if not, a direct path exists.
Tracking by reflected rays: firstly, making a mirror image point related to a plane for a transmitter, then connecting the mirror image point and a receiver, connecting a line segment formed by the two points with a plane equation, solving an intersection point with the plane, if the intersection point exists, then, a reflection point exists, and connecting the transmitter, the reflection point and the receiver to obtain a reflection ray path.
Tracing by diffracted rays: firstly, obtaining diffraction points on a straight line where split edges of a transmitter and a receiver are located, and then judging the positions of the diffraction points; if the edge is outside the cleaved edge, no diffraction exists; when the diffraction point is positioned on the splitting edge, the path where the transmitter and the diffraction point are positioned is subjected to shielding judgment; if the shielding does not exist, then shielding judgment is carried out on the path where the diffraction point and the receiver are located; if there is no occlusion, it indicates that a diffracted ray path exists.
Then, field values are calculated by UTD for the existing direct, reflected and diffracted paths respectively; the modern UTD method mainly comprises the following four modules: geometric modeling, ray tracing, occlusion judgment and field value solving. After the first three steps are completed, a field value solution is performed.
The UTD method calculates the field values, and most importantly, solves for the reflection coefficient and the diffraction coefficient.
For the solution of the reflection coefficient and the diffraction coefficient, the influence of the electrical property of the building surface material is considered, namely the size of the coefficient is related to the surface conductivity of the medium and the relative power saving constant, so that the calculation result is more accurate, and the propagation modeling is closer to the reality.
The relationship between the diffraction coefficient and the reflection coefficients of vertical and parallel polarizations is:
Figure BDA0002973716130000221
Figure BDA0002973716130000222
wherein eta is epsilonr-j×18×109σ/f，εrIs a relative dielectric constant; σ is the conductivity; phi is the diffraction coefficient.
And aiming at the combined distribution mode of each transmitter, utilizing the GPU to calculate the receiving intensity of the receiver at each position in parallel, and further drawing an electromagnetic situation map of the urban environment according to the coordinates of the receiver.
In practical applications, ray tracing can be accelerated by using an optimization technique based on back face rejection and occlusion determination: firstly, the orientation and the back of the building are tested, and the surface in the back ray propagation direction can not participate in the reflection calculation, so that the calculation efficiency can be improved; after the orientation and back test of the building is completed, the shielding test is carried out on the face of the building, and whether the face of the building is visible with the source point or not is judged.
Furthermore, two-dimensional grid division can be performed on the urban environment model, each grid point serves as a receiver, the position of a radiation source is defined, the radiation source transmitting power and the frequency band are specified, and the receiving intensity of each receiver is calculated in a traversing mode. The receivers in the observation area can be averagely distributed to each processor according to the number of CPUs, and then summarized to obtain the electromagnetic situation distribution of the whole urban environment. Because the receivers are independent, a more ideal parallel acceleration ratio can be obtained.
In addition, since ray tracing is the most common deterministic propagation modeling method at present, path loss, arrival angle and time delay can be calculated based on geometric optics, geometric diffraction theory and coherence diffraction theory. However, the ray tracing method is far more complex than an empirical model and a theoretical model, and has no simple calculation loss formula, so that a computer program must be developed to numerically solve the Maxwell equation. Compared with methods such as physical optics based on current magnetic current, a full-wave method moment method, a finite element method and the like, which are also common methods for calculating electromagnetism, ray tracing is based on local characteristics of a high-frequency field, the requirement on computer memory is low, and the calculation speed is high, so that the ray tracing method becomes a widely used radio wave propagation modeling method.
Finally, the research of domestic colleges and universities such as the western-security electronic technology university, the beijing post and telecommunications university and the electronic technology university is earlier, relatively mature ray tracing algorithm exists, but certain limitations still exist, for example, a processed scene model is simple, description of various complex buildings in cities is mostly only a simple cuboid, detail expression is lacked, and the fact that the result of simulation calculation is certainly different from an actual result is also determined. In the embodiment of the invention, the urban three-dimensional model containing the electromagnetic material information of each part of the object is used for ray tracing calculation, so that a higher radio wave propagation model can be obtained.
The three-dimensional city model can restore the three-dimensional structure of the city, so that the three-dimensional city model theoretically has very accurate solving precision, and the calculation is very quick after an accelerated calculation mode is adopted. However, the existing open source map operator or government measured data obtains a three-dimensional city model with a slow updating speed. Modern cities are rapid in development and construction, and the iteration of a three-dimensional structure is frequent, so that the real scene of the current city cannot be reflected in real time by the existing three-dimensional city model. In the embodiment of the invention, the optical image which is more convenient to collect is used, and the laser point cloud picture which is collected by the laser radar is used as an auxiliary, so that the urban three-dimensional model can be established more quickly.
The above description is only for the purpose of illustrating the preferred embodiments of the present invention and is not to be construed as limiting the invention, and any modifications, equivalents and improvements made within the spirit and principle of the present invention are intended to be included within the scope of the present invention.

Claims (10)

1. A modeling method of a city radio wave propagation model is characterized by comprising the following steps:
a: acquiring a laser point cloud picture and an optical image of a modeling area, and establishing a dense point cloud picture containing electromagnetic material information according to a laser radar and the optical image;
b: segmenting the dense point cloud graph based on corresponding electromagnetic material information, wherein the electromagnetic material information comprises: plants, glass, concrete, ground;
c: identifying corresponding object features based on the geometric structure, and fitting a solid model by using the object features, wherein the object features comprise: wall, door, roof, protrusion and window, and the geometry includes: size, position, orientation, and topology;
d: a plurality of transmitters and receivers are arranged in a three-dimensional city model formed by a solid model, and the receiving intensity of different positions corresponding to the transmitters is calculated by utilizing a ray tracing algorithm.
2. The modeling method of a urban electric wave propagation model according to claim 1, characterized in that: the specific implementation of the step A comprises the following steps:
acquiring a laser point cloud picture scanned by a laser radar in a modeling area and electromagnetic material information of each point, and adding the electromagnetic material information into a label of the point;
acquiring an optical image of a modeling area based on a full-band multi-view optical camera; identifying object features based on an algorithm of a motion acquisition structure, and constructing a sparse point cloud corresponding to the object features;
acquiring a depth point cloud image corresponding to the optical image;
identifying and calibrating electromagnetic material information of a corresponding region in the depth point cloud picture by using a pre-trained neural network model;
and fusing the calibrated depth point cloud picture, the calibrated sparse point cloud picture and the laser point cloud picture added with the label according to the geometric mapping relation, and adding corresponding object structure characteristics into the label to obtain the dense point cloud picture.
3. The modeling method of a urban electric wave propagation model according to claim 2, characterized in that: the method for identifying and calibrating the electromagnetic material information of the corresponding region in the depth point cloud picture by utilizing the pre-trained neural network model comprises the following steps:
the method comprises the steps of training a pre-built neural network model by using an optical image of an object marked with electromagnetic material information as a sample to obtain a converged neural network model, then segmenting the optical image of a modeling area by using the neural network model, and calibrating a corresponding area in a depth point cloud image by using a segmentation result.
4. The modeling method of a urban electric wave propagation model according to claim 2, characterized in that: the laser point cloud picture acquiring process comprises the following steps:
the modeling area is scanned by using the laser radar, and based on the three-dimensional coordinates of the laser radar, the modeling area is modeled by using a formula,
Figure FDA0002973716120000021
a laser spot cloud map of the modeled region is obtained, wherein,
Xpcoordinates on the x-axis for points in the point cloud; xGFor laser minesTo coordinates on the x-axis; y ispCoordinates on the y-axis for points in the point cloud; y isGCoordinates of the laser radar on the y axis are obtained; s is a laser radar point-to-point vector; theta is an included angle between the pixel corresponding to the laser ranging point P and the middle pixel in the scanning period; b ═ cos ω · sin α · cos κ + sin κ · sin ω; aerial attitude parameters of laser radar such as roll angle, pitch angle and yaw angle
Figure FDA0002973716120000022
Obtaining photogrammetric attitude angles (alpha, omega, kappa) through coordinate system transformation; zpHeight of a point in the point cloud; zGIs the height of the lidar.
5. The modeling method of a urban electric wave propagation model according to claim 1, characterized in that: the method further comprises the following steps: and before the step B, performing outlier removal processing on the dense point cloud picture.
6. The modeling method of a urban electric wave propagation model according to claim 1, characterized in that: in the step B, the step of dividing the dense point cloud picture includes:
and partitioning the dense point cloud by a point cloud partitioning algorithm based on region growing.
7. The modeling method of a urban electric wave propagation model according to claim 1, characterized in that: and C, fitting a solid model by using the object characteristics, wherein the step C comprises the following steps:
determining whether the corresponding entity belongs to undulating terrain, flat ground, simple buildings, complex buildings and trees according to the characteristics of the object;
building modeling is carried out based on a global fitting algorithm under the condition that the entity belongs to a simple building;
under the condition that the entity belongs to a tree, performing data modeling by using a tree skeleton model constructed in advance;
and under the condition that the entity belongs to a complex building or an undulating terrain, carrying out data modeling by utilizing a Poisson reconstruction algorithm.
Under the condition that the entity belongs to a flat ground, performing data modeling by adopting a least square fitting plane method;
and taking the model obtained by modeling as a solid model.
8. The modeling method of a urban electric wave propagation model according to claim 7, characterized in that: the building modeling based on the global fitting algorithm comprises the following steps:
establishing constraints on the characteristics of the building model object, wherein the constraints comprise: area, position, orientation, topology;
the entity is divided by utilizing a random sampling consistency algorithm to obtain a geometric primitive set { chi } composed of local geometric primitivesi}；
According to the geometric primitive set { χiThe parallel and orthogonal relations between every two geometric primitives contained in the method obtain the combined relation between the geometric primitives, and the maximum non-conflict set is selected from the set formed by the combined relation between the geometric primitives
Figure FDA0002973716120000041
Aligning each geometric primitive according to the combination relation in the maximum non-conflict set by utilizing a nonlinear optimization algorithm with constraints;
for aligned maximum non-conflicting set
Figure FDA0002973716120000042
Establishing a plurality of geometric primitive pairs according to the same angle, and establishing a relational graph by taking the geometric primitive pairs as vertexes;
calculating the space distance between every two geometric primitive pairs in the relational graph, and deleting the geometric primitive pairs of which the space distances to other geometric primitive pairs are larger than a set distance;
for each geometric primitive pair, connecting the geometric primitive pair with other geometric primitive pairs by using edges, wherein the geometric primitive pair and the other geometric primitive pairs have similar angles; at the same time, by using the formula,
Figure FDA0002973716120000043
a constraint is added to the newly added edge, wherein,
sca confidence score for an edge;
Figure FDA0002973716120000044
the included angle between the direction vector of the geometric primitive i and the direction vector of the geometric primitive j is shown;
Figure FDA0002973716120000045
the included angle between the direction vector of the geometric primitive k and the direction vector of the geometric primitive l is shown; gcA constraint that is an edge; n isiIs a geometric primitive i; n isjIs a geometric primitive j; n iskIs a geometric primitive k; n islIs a geometric primitive l;
taking the newly added edge set as an initial candidate set, and extracting an isogonal relation set consisting of edges with isogonal relation in the initial candidate set;
by means of the formula (I) and (II),
Figure FDA0002973716120000046
the set of angular relationships is optimized, wherein,
min is a minimum evaluation symbol; sigma is a summation symbol; ed(Pi,χi) As a set of points PiAnd primitive xiThe error between; w is apd2(p,χi) Is the accumulated error of the data; w is apIs the weight of the geometric primitive; d2(p,χi) Distance of p points to the geometric primitive;
judging whether each side in the optimized set is the same as the side before optimization;
if yes, deleting the primitive pair corresponding to the edge with the lowest confidence score from the same edges before and after the same-angle relation set optimization;
if not, returning to the step of executing the optimization of the corresponding angle relation set until the optimization times reach the set times;
acquiring a set of primitive pairs corresponding to the edge with the lowest confidence score, and aiming at each primitive pair in the set, utilizing a formula,
Figure FDA0002973716120000051
computing a confidence score for the primitive pair and a constraint, wherein,
spis the confidence score of the primitive pair; p is a radical ofiPoints on each axis of the geometric primitive i; p is a radical ofjPoints on each axis of the geometric primitive j; | | is a modulo symbol; (p)i-pj)TIs the transposition of the vector; gpConstraint of the primitive pair;
carrying out coaxial alignment processing on the primitive pairs in the maximum relation subset according to the direction of a normal vector between each primitive pair in the maximum relation subset, and sorting the primitive pairs in the optimized isometric relation set according to the confidence scores of the primitive pairs to obtain a corresponding maximum relation subset; and obtaining a three-dimensional model of the building according to the relationship between the entities corresponding to the maximum relationship subset.
9. The modeling method of a urban electric wave propagation model according to claim 7, characterized in that: under the condition that the entity belongs to a complex building, data modeling is carried out by utilizing a Poisson reconstruction algorithm, and the method comprises the following steps:
taking the point cloud data corresponding to the fitted building entity as input, and determining an indication function of the point cloud data according to the vector direction of each point cloud;
after smooth filtering is carried out on the indicating function, a gradient field of the indicating function is calculated, and the gradient field is discretized by establishing an implicit function of an octree for adjusting the depth of a grid according to the density of point cloud data;
carrying out segmentation sampling on the discretized gradient field, and calculating a vector field of a sampling point by using a cubic linear interpolation method;
the poisson equation is solved and the solution,
Figure FDA0002973716120000061
obtaining scalar quantities of point cloud dataA function is indicated in which, among other things,
Δ is laplace operator; x is an indicator function;
Figure FDA0002973716120000062
is a vector differential operator;
Figure FDA0002973716120000063
a space vector of the point cloud data;
based on a scalar indication function, extracting an isosurface by adopting a mobile cube algorithm, and establishing a three-dimensional model of the building based on the isosurface.
10. The modeling method of a urban electric wave propagation model according to claim 1, characterized in that: in the step D, the following three ray propagation modes exist in the three-dimensional city model formed by the solid model:
direct ray tracing: judging whether a line segment between the transmitter and the receiver has an intersection point with a building or not, if not, a direct path exists;
tracking by reflected rays: firstly, making a plane-related mirror image point for a transmitter, then connecting the mirror image point with a receiver, connecting a line segment consisting of the two points with a plane equation, solving an intersection point with the plane, if the intersection point exists, then, a reflection point exists, and connecting the transmitter, the reflection point and the receiver to obtain a reflection ray path;
tracing by diffracted rays: firstly, obtaining diffraction points on a straight line where split edges of a transmitter and a receiver are located, and then judging the positions of the diffraction points; if the edge is outside the cleaved edge, no diffraction exists; when the diffraction point is positioned on the splitting edge, the path where the transmitter and the diffraction point are positioned is subjected to shielding judgment; if the shielding does not exist, then shielding judgment is carried out on the path where the diffraction point and the receiver are located; if there is no occlusion, it indicates that a diffracted ray path exists.
CN202110269609.9A 2021-03-12 2021-03-12 Modeling method of urban radio wave propagation model Active CN113066161B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110269609.9A CN113066161B (en) 2021-03-12 2021-03-12 Modeling method of urban radio wave propagation model

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110269609.9A CN113066161B (en) 2021-03-12 2021-03-12 Modeling method of urban radio wave propagation model

Publications (2)

Publication Number Publication Date
CN113066161A true CN113066161A (en) 2021-07-02
CN113066161B CN113066161B (en) 2022-04-29

Family

ID=76560182

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110269609.9A Active CN113066161B (en) 2021-03-12 2021-03-12 Modeling method of urban radio wave propagation model

Country Status (1)

Country Link
CN (1) CN113066161B (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114612536A (en) * 2022-03-22 2022-06-10 北京诺亦腾科技有限公司 Method, device and equipment for identifying three-dimensional model of object and readable storage medium

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR19990080905A (en) * 1998-04-23 1999-11-15 윤종용 Prediction Method of Propagation Characteristics of Radio Wave Considering Polarization Effect in Urban Canyon Model
CN103942420A (en) * 2014-04-08 2014-07-23 北京大学 Rapid estimation method for solar energy in construction size
CN104180793A (en) * 2014-08-27 2014-12-03 北京建筑大学 Device and method for obtaining mobile spatial information for digital city construction
CN104680581A (en) * 2015-03-05 2015-06-03 西安电子科技大学 Architecture object selection method in three-dimensional environment
CN111815776A (en) * 2020-02-04 2020-10-23 山东水利技师学院 Three-dimensional building fine geometric reconstruction method integrating airborne and vehicle-mounted three-dimensional laser point clouds and streetscape images
CN112348867A (en) * 2020-11-18 2021-02-09 南通市测绘院有限公司 Method and system for constructing city high-precision three-dimensional terrain based on LiDAR point cloud data
CN112418245A (en) * 2020-11-04 2021-02-26 武汉大学 Electromagnetic emission point positioning method based on urban environment physical model

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR19990080905A (en) * 1998-04-23 1999-11-15 윤종용 Prediction Method of Propagation Characteristics of Radio Wave Considering Polarization Effect in Urban Canyon Model
CN103942420A (en) * 2014-04-08 2014-07-23 北京大学 Rapid estimation method for solar energy in construction size
CN104180793A (en) * 2014-08-27 2014-12-03 北京建筑大学 Device and method for obtaining mobile spatial information for digital city construction
CN104680581A (en) * 2015-03-05 2015-06-03 西安电子科技大学 Architecture object selection method in three-dimensional environment
CN111815776A (en) * 2020-02-04 2020-10-23 山东水利技师学院 Three-dimensional building fine geometric reconstruction method integrating airborne and vehicle-mounted three-dimensional laser point clouds and streetscape images
CN112418245A (en) * 2020-11-04 2021-02-26 武汉大学 Electromagnetic emission point positioning method based on urban environment physical model
CN112348867A (en) * 2020-11-18 2021-02-09 南通市测绘院有限公司 Method and system for constructing city high-precision three-dimensional terrain based on LiDAR point cloud data

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
张富彬等: "基于深度学习的全球电离层TEC预测", 《电波科学学报》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114612536A (en) * 2022-03-22 2022-06-10 北京诺亦腾科技有限公司 Method, device and equipment for identifying three-dimensional model of object and readable storage medium

Also Published As

Publication number Publication date
CN113066161B (en) 2022-04-29

Similar Documents

Publication Publication Date Title
CN113066162B (en) Urban environment rapid modeling method for electromagnetic calculation
CN106127857B (en) The on-board LiDAR data modeling method of integrated data driving and model-driven
Overby et al. Automatic 3D building reconstruction from airborne laser scanning and cadastral data using Hough transform
CN112418245B (en) Electromagnetic emission point positioning method based on urban environment physical model
EP2849117B1 (en) Methods, apparatuses and computer program products for automatic, non-parametric, non-iterative three dimensional geographic modeling
WO2006027339A2 (en) Method and system for 3d scene change detection
Budroni et al. Automatic 3D modelling of indoor manhattan-world scenes from laser data
US20080123961A1 (en) Scalable method for rapidly detecting potential ground vehicle under cover using visualization of total occlusion footprint in point cloud population
CN114066960B (en) Three-dimensional reconstruction method, point cloud fusion method, device, equipment and storage medium
Jebur et al. Assessing the performance of commercial Agisoft PhotoScan software to deliver reliable data for accurate3D modelling
Vasquez-Gomez et al. Hierarchical ray tracing for fast volumetric next-best-view planning
CN110838115A (en) Ancient cultural relic three-dimensional model change detection method by contour line extraction and four-dimensional surface fitting
CN112305559A (en) Power transmission line distance measuring method, device and system based on ground fixed-point laser radar scanning and electronic equipment
Jin et al. An indoor location-based positioning system using stereo vision with the drone camera
Stal et al. Test case on the quality analysis of structure from motion in airborne applications
CN116543117B (en) High-precision large-scene three-dimensional modeling method for unmanned aerial vehicle images
WO2023124676A1 (en) 3d model construction method, apparatus, and electronic device
Yang et al. Linking persistent scatterers to the built environment using ray tracing on urban models
Özdemir et al. A multi-purpose benchmark for photogrammetric urban 3D reconstruction in a controlled environment
CN116563466A (en) Deep learning-based three-dimensional Shan Mudian cloud completion method
CN113066161B (en) Modeling method of urban radio wave propagation model
Ma et al. Rapid reconstruction of a three-dimensional mesh model based on oblique images in the Internet of Things
CN113345072A (en) Multi-view remote sensing topographic image point cloud reconstruction method and system
Wan et al. Enhance accuracy: Sensitivity and uncertainty theory in LiDAR odometry and mapping
Fiocco et al. Multisensor fusion for volumetric reconstruction of large outdoor areas

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