Large-scale urban electromagnetic environment analysis method based on measurement and calculation fusion

Abstract

The present invention discloses a large-scale urban electromagnetic environment analysis method based on measurement and calculation fusion. It mainly solves the problems of high calculation cost, excessive consumption of computing resources, lack of timeliness and low accuracy in the existing technology. The implementation scheme is: extracting the coordinates of all source points and observation points in the urban environment and radiation source model diagram; using a parallel computing strategy based on a message passing interface to allocate observation points to each CPU and perform ray tracing and occlusion judgment from the observation point to the source point; using a high-order moment method to calculate the electromagnetic current on the surface of the radiation source; calculating the reflection field and diffraction field at the occlusion point based on the electromagnetic current on the surface of the radiation source, and then obtaining the simulated field intensity of the observation point by superposition of field intensity; using a spectrum sensor to measure the selected measured points in a real urban environment and using the precise measurement results to correct the simulated field intensity obtained by fusion simulation to obtain a more accurate fusion field intensity. The present invention greatly shortens the calculation time and improves the calculation accuracy, and can be used for communication spectrum planning and electromagnetic environment assessment.

Description

Technical Field

[0001] The present invention belongs to the technical field of electromagnetic measurement and calculation fusion, and specifically relates to a method for fusing measurement and simulation, which can be used for communication spectrum planning and electromagnetic environment assessment. Background Art

[0002] In recent years, with the advancement and application of electromagnetic technology, an increasing number of radiation sources and wireless communication devices have become part of the daily lives of urban residents. These devices, ranging from traditional radio and television to modern mobile communications devices, as well as the soon-to-be-widely deployed 5G and future 6G networks, have exacerbated the complexity of the urban electromagnetic environment. As the number of wireless devices increases, specific frequency bands become congested, leading to signal interference and reduced communication efficiency. Furthermore, strong electromagnetic radiation poses a potential threat to the normal operation of critical infrastructure such as precision instruments, communication networks, and traffic control systems. Therefore, accurate and effective analysis of the complex urban electromagnetic environment is crucial. This analysis can help reduce electromagnetic pollution, optimize spectrum allocation and utilization, and ensure the sustainable development of the urban electromagnetic environment. Traditional methods for analyzing the electromagnetic environment suffer from numerous challenges, including high computational costs, excessive consumption of computing resources and time, a lack of timeliness, and low accuracy.

[0003] Patent document CN201910012530.0 discloses a method for predicting the situation of an urban outdoor electromagnetic environment based on MoM-UTD. Its implementation steps are as follows: 1) establishing a city simulation model and an antenna simulation model, extracting geometric information and electrical parameter information from the city simulation model, and extracting the operating frequency and source location of the antenna simulation model; 2) setting sampling points of the route to be measured as field points; 3) tracing all rays from the source point to the field point, storing ray information and the point of action; 4) using the moment method MoM to calculate the radiated electric field of the antenna simulation model at the point of action; 5) using the UTD method to calculate the radiated electric field from the point of action to the field point for all ray types; 6) superimposing the radiated electric fields of all ray types at the same field point to obtain the prediction result. Although this method is low-cost, accurate, and efficient, it still suffers from excessive computational complexity, insufficient computational speed, and limited accuracy of the final calculation results. Summary of the Invention

[0004] The purpose of the present invention is to overcome the defects in the above-mentioned prior art and propose a large-scale urban electromagnetic environment analysis method based on measurement and calculation fusion to reduce the amount of calculation, improve the calculation speed and accuracy, and meet the needs of accurate analysis of large-scale urban electromagnetic environment.

[0005] The key technology to achieve the purpose of the present invention is to simulate and calculate the urban electromagnetic environment based on the high-order method of moments-uniform diffraction theory HOMOM-UTD based on the message passing interface MPI, so as to accurately and quickly analyze the electromagnetic environment of large-scale urban areas. The steps include the following:

[0006] (1) Extract the location and height parameters of major buildings and facilities from the city model map, and extract the location coordinates of all radiation sources and observation points;

[0007] (2) Using a parallel computing strategy based on the message passing interface (MPI), all observation points are sorted according to their distance from the radiation source. The observation points are then alternately assigned to each CPU for ray tracing and occlusion judgment based on the sorting. Three paths are obtained: one directly from the radiation source to the observation point, one after one reflection or diffraction, and one after two reflections or diffraction. The coordinates of the reflection point and the diffraction point are then calculated.

[0008] (3) Using the high-order moment method HOMoM to calculate the electromagnetic current J on the surface of the radiation source s (p,s) and M s (p,s);

[0009] (4) Calculate the simulated electromagnetic field E at the observation point based on the results of steps (2) and (3) s (r′):

[0010] (4a) For the path from the radiation source directly to the observation point, the direct field E0(r′) is calculated based on the electromagnetic current on the surface of the radiation source;

[0011] (4b) For the path that reaches the observation point after one reflection or diffraction, the reflection field of the first-order ray is calculated by the uniform diffraction theory UTD algorithm based on the electromagnetic current on the surface of the radiation source. or diffraction field and as the electromagnetic field E1(r′) generated by the first-order ray;

[0012] (4c) For the path that reaches the observation point after two reflections or diffractions, the final reflection field of the second-order ray at the observation point is obtained by the UTD algorithm based on the electromagnetic current on the surface of the radiation source using the double uniform diffraction theory or diffraction field and as the electromagnetic field E2(r′) generated by the second-order rays;

[0013] (4d) The vector sum of the direct field E0(r), the electromagnetic field E1(r′) generated by the first-order ray, and the electromagnetic field E2(r′) generated by the second-order ray is obtained to obtain the simulated electromagnetic field E at the observation point. s (r′);

[0014] (5) Deploy radiation sources and a limited number of spectrum sensors throughout the city, use these devices to measure the electromagnetic field at selected measurement points in the real environment to obtain accurate electromagnetic field data at the measurement points, and use the accurate data to calibrate the simulated electromagnetic field E s (r′), and the corrected electromagnetic field is obtained

[0015]

[0016] in, Represents the observation point x p The corrected electromagnetic field, E s (x p ) and E s (x j ) represent the observation points x p and x j The simulated electromagnetic field value at m (x j ) represents the observation point x j The actual measured value at ω(x j )=ω p (x j )·ω d (x j ) is the final weighting factor, ω P (x j ) represents the distance from the radiation source to the measured point x j The power weighting factor, ω d (x j ) represents the measured point x j and the predicted point x p The distance weighting factor between them.

[0017] Compared with the prior art, the present invention has the following advantages:

[0018] First, the present invention uses the high-order moment method HOMoM to calculate the electromagnetic current J on the surface of the radiation source s (p,s) and M s (p,s), which can handle complex models more effectively, with higher calculation accuracy and faster calculation speed while occupying less memory space.

[0019] Secondly, the present invention adopts a parallel computing strategy based on the message passing interface MPI to assign observation points to each central processing unit CPU for ray tracing and occlusion judgment, which not only greatly accelerates the calculation speed and shortens the processing time, but also expands the calculation scale and can quickly analyze the light path between the source point and the observation point.

[0020] Third, the present invention adopts a solution that integrates measurement and simulation to correct the static electromagnetic field simulated by the MPI-based HOMoM-UTD hybrid algorithm, that is, the precise data obtained by measurement with professional equipment is used to correct the simulated static electromagnetic field, so that the corrected data has higher accuracy and real-time performance. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] Figure 1 It is an implementation flow chart of the present invention;

[0022] Figure 2 It is a simulation model diagram and a partial enlarged diagram used in the embodiment of the present invention;

[0023] Figure 3 This is the location diagram of the two radiation sources of 409.5MHz and 2501MHz in the simulation model;

[0024] Figure 4 is a simulated electromagnetic field distribution diagram calculated for all observation points in an embodiment of the present invention;

[0025] Figure 5 This is a diagram of a radiation source placed for measuring the true value of the electromagnetic field in an embodiment of the present invention;

[0026] Figure 6 is a distribution diagram of spectrum sensors placed to measure the true value of the electromagnetic field in an embodiment of the present invention;

[0027] Figure 7 This is an electromagnetic field distribution diagram after the actual values ​​measured in two frequency bands and the simulation values ​​are fused in an embodiment of the present invention;

[0028] Figure 8 This is a local magnified view of the electromagnetic field distribution after the two frequency bands are fused;

[0029] Figure 9 This is a comparison chart of the sampling point distribution and fusion value, simulation value and measurement value in 409.5MHz scenario A;

[0030] Figure 10 The sampling point distribution in 409.5MHz scenario B and the comparison of fusion value, simulation value and measurement value;

[0031] Figure 11 This is a comparison chart of the sampling point distribution and fusion value, simulation value and measurement value in 409.5MHz scenario C;

[0032] Figure 12 This is a comparison chart of the sampling point distribution and fusion value, simulation value and measurement value in 2501MHz scenario D;

[0033] Figure 13This is a comparison chart of the sampling point distribution and fusion value, simulation value and measurement value in 2501MHz scenario E. DETAILED DESCRIPTION

[0034] In order to enable those skilled in the art to better understand the solutions of the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings of the embodiments of the present invention. Obviously, the embodiments described are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, other embodiments obtained by ordinary technicians in this field without making creative efforts should fall within the scope of protection of the present invention.

[0035] The implementation scenario of the fusion of actual measurement and simulation in this example is based on an area of ​​60 square kilometers in Xi'an, Shaanxi Province.

[0036] Reference Figure 1 The implementation steps of this example include the following:

[0037] Step 1: Extract key parameters from an accurate model of the urban environment and radiation sources.

[0038] 1.1) Extract the location and height parameters of major buildings and facilities from the city model:

[0039] The city model map uses a CAD model provided by a surveying and mapping company, and extracts information about the location, floor and height of major buildings and facilities, such as Figure 2 As shown, among which Figure 2 (a) is a simulation model diagram, which contains 29,170 buildings and 173,591 observation points. Figure 2 (b) Yes Figure 2 (a) A partial enlarged view;

[0040] 1.2) Extract the position coordinates of the observation point and the radiation source point.

[0041] To avoid interference from other radiation sources, this example selects, but is not limited to, two idle frequency bands centered at 409.5 MHz and 2501 MHz. An omnidirectional antenna with a gain of -2 dBi and a power of 45 dBm is used to establish radiation sources at these frequencies. The antenna's operating frequency range is extended from 20 MHz to 6 GHz.

[0042] Step 2: Calculate the simulated electromagnetic field using a hybrid algorithm of the high-order method of moments and uniform geometric diffraction theory based on the message passing interface.

[0043] 2.1) Ray tracing and occlusion judgment are performed based on the parallel computing strategy based on the message passing interface MPI:

[0044] The reverse ray tracing is based on the principle of geometric optics and searches for all valid ray paths that can reach the source point from the observation point. Compared with the existing forward ray tracing, it has higher calculation accuracy. The reverse ray tracing may include serial algorithm tracing and parallel algorithm tracing. Since tracing needs to be performed again for different observation points, and the number of observation points in this example is huge, if the serial algorithm is used for tracing, problems such as large amount of calculation and long calculation time will arise. Therefore, this example adopts a parallel tracing calculation strategy based on the message passing interface MPI, and allocates observation points to each central processing unit CPU on the supercomputer cluster for parallel ray tracing and occlusion judgment, thereby identifying each valid ray path that can reach the source point from the observation point and calculating the coordinates of all reflection points and diffraction points.

[0045] This example uses 4096 cores at 409.5 MHz and 2501 MHz. The relative permittivity of the building in the model is ε. r =15, conductivity is σ = 0.5S / m. The position of the radiation source in the simulation model is as follows Figure 3 Because the influence of high-order rays on the observation point is very small, only three paths are retained: directly from the radiation source to the observation point, after one reflection or diffraction to the observation point, and after two reflections or diffraction to the observation point.

[0046] 2.2) Calculate the electromagnetic current on the surface of the radiation source according to the high-order moment method HOMoM:

[0047] The Moment Method (MoM) is a numerical method that expands functional equations through basis functions, and converts functional equations into linear equations through a verification process while ensuring that the error is very small. The High-Order Moment Method (HOMoM) consists of two parts: bilinear patches and high-order basis functions. The high-order geometric patches can flexibly approximate any surface, while the high-order basis functions can flexibly approximate any current distribution. Compared with the Moment Method, the High-Order Moment Method has the advantages of stronger ability to handle complex models and higher calculation accuracy. Therefore, this example uses the High-Order Moment Method (HOMoM) to calculate the surface current J of the radiation source. s (p,s) and surface magnetic flux M s (p,s), which is calculated as follows:

[0048]

[0049] in, are polynomial basis functions, is a scale factor, p i s j is a specific term in the polynomial expansion; a ij and b ij is the coefficient to be determined, N p and Ns are the expansion orders in the p-direction and the s-direction, respectively.

[0050] 2.3) Calculate the direct field E0(r′) from the radiation source to the observation point based on the current distribution Js(p,s) on the radiation source surface:

[0051]

[0052] Where k is the beam, η is the free space impedance, and exp(-jk|r′-r) represents the phase change of the electromagnetic wave propagating from the source point to the observation point. Indicates that the source current element is The phase change in the direction, r and r′ are the positions of the source and observation points respectively, is the unit vector from the source point to the observation point;

[0053] 2.4) According to the surface current distribution J of the radiation source s (p,s) The electromagnetic field E1(r′) generated by the first-order ray at the observation point is calculated using the uniform diffraction theory UTD algorithm:

[0054] 2.4.1) According to the electromagnetic current J on the radiation source surface s (p,s) Calculate the incident field at the reflection or diffraction point

[0055]

[0056] Where r1′ represents the position of the reflection point or diffraction point, represents the unit vector from the source point to the reflection point or diffraction point, exp(-jk|r1′-r) represents the phase change of the electromagnetic wave from the source point to the reflection point or diffraction point, Indicates that the source current element is Phase changes in direction;

[0057] 2.4.2) According to the vertical polarization reflection coefficient R (⊥) and the horizontal polarization reflection coefficient R (||) Obtain the reflection coefficient matrix [R]:

[0058]

[0059] in, Represents the wave impedance of medium m, m=1,2, μ m , σ m and ε m are the magnetic permeability, electrical conductivity and dielectric constant of medium m, θ i and θ t denote the angle of incidence and the angle of transmission respectively;

[0060] 2.4.3) According to the soft diffraction coefficient D s and hard diffraction coefficient D h Obtain the diffraction coefficient matrix [D]:

[0061]

[0062] in

[0063]

[0064] is the soft diffraction coefficient related to the incident angle α0,

[0065] is the incident angle α n The related soft diffraction coefficient,

[0066] is the propagation correlation factor of the 0th soft diffraction wave

[0067] is the propagation correlation factor of the nth soft diffraction wave

[0068] is a constant term,

[0069] is the hard diffraction coefficient related to the incident angle α0,

[0070] is the incident angle α n The relevant hard diffraction coefficient,

[0071] is the propagation correlation factor of the 0th hard diffraction wave,

[0072] is the propagation correlation factor of the nth hard diffraction wave;

[0073] 2.4.4) According to the reflection coefficient matrix [R], diffraction coefficient matrix [D] and incident field Obtain the reflection field of the first-order ray or diffraction field

[0074]

[0075] Among them, A(s r ) and A(s d ) represent the amplitude factors of the reflected field and the diffracted field, and They represent the phase factors of the reflected field and the diffracted field, respectively. [R] and [D] represent the reflection coefficient matrix and the diffraction coefficient matrix, respectively.

[0076] 2.4.5) Determine the electromagnetic field E1(r′) generated by the first-order ray at the observation point based on whether reflection or diffraction occurs at the obstruction:

[0077] When reflection occurs at the occlusion,

[0078] When diffraction occurs at the obstruction,

[0079] 2.5) According to the surface current distribution J of the radiation source s (p,s) The electromagnetic field E2(r′) generated by the second-order ray at the observation point is obtained by the double uniform diffraction theory UTD algorithm:

[0080] 2.5.1) According to the electromagnetic current J on the radiation source surface s (p,s) Calculate the incident field at the first reflection or diffraction point

[0081]

[0082] Where r2′1 represents the position of the first reflection point or diffraction point, represents the unit vector from the source point to the first reflection point or diffraction point, exp(-jk|r2′1-r) represents the phase change of the electromagnetic wave from the source point to the first reflection point or diffraction point, Indicates that the source current element is Phase changes in direction;

[0083] 2.5.2) According to the incident field Obtain the reflection field at the first reflection point or diffraction point or diffraction field

[0084]

[0085] Among them, A(s r ) and A(s d ) represent the amplitude factors of the reflection field or diffraction field at the first reflection point or diffraction point, respectively.

[0086] and Represents the phase factor of the reflected field or diffracted field at the first reflection point or diffraction point respectively;

[0087] 2.5.3) Determine the incident field at the final reflection point or diffraction point based on whether reflection or diffraction occurs at the first obstruction.

[0088] When the ray is reflected at the first occlusion,

[0089] When the ray is diffracted at the first obstruction,

[0090] 2.5.4) According to the incident field at the final reflection point or diffraction point Find the final reflection field or diffraction field

[0091]

[0092] Among them, A(s r ) and A(s d ) represent the amplitude factors of the final reflected field and diffracted field, and Represent the phase factors of the final reflected field and diffracted field respectively;

[0093] 2.5.5) Determine the electromagnetic field E1(r′) generated by the second-order ray at the observation point based on whether reflection or diffraction occurs at the second obstruction:

[0094] When reflection occurs at the second occlusion,

[0095] When diffraction occurs at the second obstruction,

[0096] 2.6) Add the vectors of the direct field E0(r), the electromagnetic field E1(r′) generated by the first-order ray, and the electromagnetic field E2(r′) generated by the second-order ray to obtain the simulated electromagnetic field E at the observation point. s (r′):

[0097] E s (r′)=E0(r′)+E1(r′)+E2(r′);

[0098] 2.7) The static electromagnetic field distribution of the overall model is obtained based on the simulated electromagnetic fields of all observation points:

[0099] This example has a total of 173,591 observation points. The simulation electromagnetic field calculation is performed on all observation points to obtain the static electromagnetic field distribution of the 60 square kilometers urban area at 409.5MHz and 2501MHz in the simulation model diagram, as shown in the figure. Figure 4 As shown, Figure 4 (a) is the static electromagnetic field distribution diagram at 409.5MHz, Figure 4 (b) is the static electromagnetic field distribution diagram at 2501MHz.

[0100] Step 3: Use the actual measured value to correct and fuse the simulated value.

[0101] 3.1) Use spectrum sensors to measure accurate electromagnetic field data at actual measurement points in a real urban environment:

[0102] Taking into account the modeling errors introduced by factors such as vegetation, moving vehicles, and the shape and material of buildings, the simulated static electromagnetic field distribution inevitably differs from the real value. To this end, these differences are corrected by fusing measurement and simulation.

[0103] In order to achieve rapid collection, transmission, and analysis of real electromagnetic data in urban blocks, this example uses a spectrum sensor with ultra-wideband electromagnetic signal detection capabilities to measure the real urban electromagnetic environment. At the same time, to ensure consistency between the measurement and simulation scenarios, radiation sources with center frequencies of 409.5 MHz and 2501 MHz are installed in positions corresponding to the radiation sources in the simulation scenario, such as Figure 5 As shown, Figure 5 (a) is a physical picture of the 409.5MHz radiation source. Figure 5 (b) is a physical picture of the 2501MHz radiation source.

[0104] During the test, the Figure 2 The simulation model shows that the urban area is divided into 60 sub-regions, each approximately 1 square kilometer in size. A total of 885 measurement points were established across all sub-regions. Spectrum sensors were distributed throughout each sub-region and continuously measured at fixed points, collecting data once per second for three hours at each point to obtain real electromagnetic field data. Figure 6 Shows the distribution of measurement points and spectrum sensors in a sub-area in a 409.5 MHz measurement scenario.

[0105] 3.2) After the measurement, the data of 885 measured points were collected and sorted, and compared with the electromagnetic field E of each observation point obtained by simulation. s (r′) performs fusion correction:

[0106] 3.2.1) The position of the known simulation observation point and the corresponding simulation value are expressed as (x i ,E s (x i )), the position of the measured point and the corresponding true measurement value are expressed as (x j ,E m (x j )), where i = 1, 2, ..., n, j = 1, 2, ..., m, n is the number of observation points, and m is the number of measured points;

[0107] 3.2.2) According to the measured point x j and observation point x p The distance dp,j , calculate the measured point x j and the predicted point x p The distance weighting factor ω between d (x j ):

[0108]

[0109] Where β is the path loss exponent of the distance weighting factor;

[0110] 3.2.3) According to the radiation source to the measured point x j The distance D j , calculate the radiation source to the measured point x j The power weighting factor ω P (x j ):

[0111]

[0112] Where α is the path loss exponent of the power weighting factor;

[0113] 3.2.4) According to the power weighting factor ω P (x j ) and the distance weighting factor ω d (x j ) calculates the correction weighting factor ω(x j ):

[0114] ω(x j )=ω p (x j )·ω d (x j )

[0115] 3.2.5) For Figure 2 Observation point x in the simulation model diagram p , using the correction weighting factor ω(x j ) to correct the analog value, the formula is as follows:

[0116]

[0117] in, Represents the observation point x p The corrected electromagnetic field at the observation point x is represented by Es(xp) and Es(xj). p and x j The simulated electromagnetic field value at the observation point x is represented by Em(xj). j The actual measurement value at

[0118] 3.2.6) Change point x p Corrected fusion value at is considered as a measurement value and the next uncorrected observation point x p+1 Repeat the same process to calculate the correction values ​​for all observation points;

[0119] 3.2.7) According to the electromagnetic field values ​​of all the observation points after correction and fusion, we can get Figure 2 The fusion electromagnetic field distribution diagram of the simulation model diagram, such as Figure 7 As shown, Figure 7 (a) is the fusion electromagnetic field distribution diagram at 409.5MHz, Figure 7 (b) is the fusion electromagnetic field distribution diagram at 2501MHz, and its local details are as follows Figure 8 As shown, Figure 8 (a) Yes Figure 7 (a) is a partial enlarged view. Figure 8 (b) Yes Figure 7 (b) A partial enlarged view.

[0120] The effect of the present invention can be further illustrated by comparing parameters in different scenarios:

[0121] 1. Scene Setting

[0122] Under the premise that the radiation source frequency is 409.5MHz, Figure 2 In the simulation model diagram, three different verification areas are selected for observation point sampling. These three verification areas are called scene A, scene B and scene C respectively. Different sampling points are selected for each area. The sampling points of scene A are as follows: Figure 9 As shown in (a), the sampling points of scene B are as follows Figure 10 As shown in (a), the sampling points of scene C are as follows Figure 11 As shown in (a).

[0123] Under the premise that the radiation source frequency is 2501MHz, Figure 2 In the simulation model, two different verification areas are selected for observation point sampling. These two verification areas are called scene D and scene E. Different sampling points are selected for each area. The sampling points of scene D are as follows: Figure 12 As shown in (a), the sampling points of scene E are as follows Figure 13 As shown in (a).

[0124] 2. Compare the fusion value, simulation value and measurement value of different scene sampling points, and judge whether the fusion value is more accurate based on the difference between the fusion value, simulation value and actual measurement value:

[0125] Comparison 1: Yes Figure 9 (a) Verify the fusion value, simulation value and measurement value of the sampling point in scene A and draw the sampling point-electric field value line graph respectively. The results are as follows: Figure 9 (b) shown.

[0126] from Figure 9 (b) It can be seen that the simulation value line graph fluctuates greatly, but the fusion value line graph after fusion is smoother than the simulation value line graph, and at some observation points, the fusion true value is closer to the measured value than the simulation value.

[0127] Comparison 2: Yes Figure 10 (a) Verify the fusion value, simulation value and measurement value of the sampling point in scene B and draw the line graph of the sampling point-electric field value. The results are as follows: Figure 10 (b) shown.

[0128] from Figure 10 (b) It can be seen that the error between the simulation value line graph and the measurement value line graph is large, but the difference between the fusion value line graph and the measurement value line graph after fusion is obviously much smaller than the error between the simulation value line graph and the measurement value line graph.

[0129] Comparison 3: Yes Figure 11 (a) Verify that the fusion value, simulation value, and measurement value of the sampling point in scenario C are plotted as line graphs of the sampling point and the electric field value. The results are as follows: Figure 11 (b) shown.

[0130] from Figure 11 (b) It can be seen that although the simulation value line graph and the fusion value line graph are consistent with the overall fluctuation trend of the measured value, it is obvious that the fusion value line graph is more stable and closer to the measured value line graph.

[0131] Comparison 4: Pair Figure 12 (a) Verify that the fusion value, simulation value, and measurement value of the sampling point in scene D are plotted as line graphs of the sampling point and the electric field value. The results are as follows: Figure 12 (b) shown.

[0132] from Figure 12 (b) It can be seen that the amplitude of the simulation value line graph fluctuates greatly. In addition, there is a certain gap between the simulation value and the measured value at many sampling points. However, the fluctuation amplitude of the fusion value line graph is much smaller and the gap with the measured value is also reduced.

[0133] Comparison 5: Yes Figure 13 (a) Verify that the fusion value, simulation value, and measurement value of the sampling point in scenario E are plotted as line graphs of the sampling point and the electric field value. The results are as follows: Figure 13 (b) shown.

[0134] from Figure 13 (b) It can be seen that the line graph of the simulation value has a large error compared with the line graph of the measurement value, but the fusion value at each sampling point can accurately simulate the measurement value in the real environment.

[0135] The root mean square error (RMSE) between the fusion value, simulation value and actual measurement value of the sampling points in the above five scenarios is calculated, and the results are shown in Table 1.

[0136] Table 1 Comparison of fusion values, simulation values ​​and measured values ​​in different scenarios

[0137]

[0138]

[0139] As can be seen from Table 1, the root mean square error (RMSE) between the fusion value and the actual measurement value in different scenarios is lower than the root mean square error (RMSE) between the simulation value and the actual measurement value, and remains below 8 dB, which proves the accuracy of the measurement and simulation fusion technology proposed in the present invention.

[0140] The comparison of the above parameters shows that the novel measurement and simulation fusion technology proposed in the present invention can accurately simulate large-scale urban electromagnetic environments with high efficiency.

[0141] It should be noted that the step numbers in the specification and claims of the present invention are only for a clear description of the embodiments of the present invention and for ease of understanding, and the order of the step numbers is not limited.

Claims

1. A large-scale urban electromagnetic environment analysis technology based on measurement and calculation fusion, characterized by: The steps include: (1) Extract the location and height parameters of major buildings and facilities from the city model map, and extract the location coordinates of all radiation sources and observation points; (2) Using a parallel computing strategy based on the message passing interface (MPI), all observation points are sorted according to their distance from the radiation source. The observation points are then alternately assigned to each CPU for ray tracing and occlusion judgment based on the sorting. Three paths are obtained: one directly from the radiation source to the observation point, one after one reflection or diffraction, and one after two reflections or diffraction. The coordinates of the reflection point and the diffraction point are then calculated. (3) Using the high-order moment method HOMoM to calculate the electromagnetic current J on the radiation source surface s (p,s) and M s (p,s); (4) Calculate the simulated electromagnetic field E at the observation point based on the results of steps (2) and (3) s (r′): (4a) For the path from the radiation source directly to the observation point, the direct field E0(r′) is calculated based on the electromagnetic current on the surface of the radiation source; (4b) For the path that reaches the observation point after one reflection or diffraction, the reflection field of the first-order ray is calculated by the uniform diffraction theory UTD algorithm based on the electromagnetic current on the surface of the radiation source. or diffraction field and as the electromagnetic field E1(r′) generated by the first-order ray; (4c) For the path that reaches the observation point after two reflections or diffractions, the final reflection field of the second-order ray at the observation point is obtained by the UTD algorithm based on the electromagnetic current on the surface of the radiation source using the double uniform diffraction theory or diffraction field and as the electromagnetic field E2(r′) generated by the second-order rays; (4d) The vector sum of the direct field E0(r), the electromagnetic field E1(r′) generated by the first-order ray, and the electromagnetic field E2(r′) generated by the second-order ray is obtained to obtain the simulated electromagnetic field E at the observation point. s (r′); (5) Deploy radiation sources and a limited number of spectrum sensors throughout the city, use these devices to measure the electromagnetic field at selected measurement points in the real environment to obtain accurate electromagnetic field data at the measurement points, and use the accurate data to calibrate the simulated electromagnetic field E s (r′), and the corrected electromagnetic field is obtained in, Represents the observation point x p The corrected electromagnetic field at the observation point x is represented by Es(xp) and Es(xj). p and x j The simulated electromagnetic field value at the observation point x is represented by Em(xj). j The actual measurement value at ω(xj)=ωp(xj)·ωd(xj) is the final weighting factor, ω P (xj) represents the distance from the radiation source to the measured point x j The power weighting factor of ωd(xj) represents the measured point x j and observation point x p The distance weighting factor between them.

2. The method according to claim 1, characterized in that In step (2), the observation points are alternately assigned to each CPU for ray tracing and occlusion judgment according to the sorting. This is done by using a high-precision reverse ray tracing method. By using the image method and the principle of geometric optics, each ray from the observation point to the source point is identified, and the path that cannot directly reach the source point from the observation point is judged as an occlusion path.

3. The method according to claim 1, characterized in that In step (3), the high-order moment method HOMoM is used to calculate the surface current J of the radiation source s (p,s) and surface magnetic flux M s (p,s), the formulas are as follows: in, are polynomial basis functions, is a scale factor, p i s j is a specific term in the polynomial expansion; a ij and b ij is the coefficient to be determined, N p and N s are the expansion orders in the p-direction and the s-direction, respectively.

4. The method according to claim 1, characterized in that In step (4a), the direct field E0(r′) of the path from the radiation source to the observation point is calculated based on the electromagnetic current on the radiation source surface. The formula is as follows: Where k is the beam, η is the free space impedance, and exp(-jk|r′-r) represents the phase change of the electromagnetic wave propagating from the source point to the observation point. Indicates that the source current element is Phase change in direction, J s (p, s) represents the current distribution at the radiation source, r and r′ are the source and observation points respectively, is the unit vector from the source point to the observation point.

5. The method according to claim 1, wherein In step (4b), the reflection field of the first-order ray is calculated by the uniform diffraction theory UTD algorithm based on the electromagnetic current on the surface of the radiation source. or diffraction field The implementation steps include the following: 4b1) According to the vertical polarization reflection coefficient R (⊥) and the horizontal polarization reflection coefficient R (||) Obtain the reflection coefficient matrix [R]: in, Represents the wave impedance of medium m, m=1,2, μ m , σ m and ε m are the magnetic permeability, electrical conductivity and dielectric constant of medium m, θ i and θ t denote the angle of incidence and the angle of transmission respectively; 4b2) According to the soft diffraction coefficient D s and hard diffraction coefficient D h Obtain the diffraction coefficient matrix [D]: in is the soft diffraction coefficient related to the incident angle α0, is the incident angle α n The related soft diffraction coefficient, is the propagation correlation factor of the 0th soft diffraction wave is the propagation correlation factor of the nth soft diffraction wave is a constant term, is the hard diffraction coefficient related to the incident angle α0, is the incident angle α n The relevant hard diffraction coefficient, is the propagation correlation factor of the 0th hard diffraction wave, is the propagation correlation factor of the nth hard diffraction wave; 4b3) According to the electromagnetic current J on the surface of the radiation source s (p,s) Calculate the incident field at the reflection or diffraction point Where k is the beam, η is the free space impedance, r and r1′ represent the positions of the radiation source and reflection or diffraction point, respectively. represents the unit vector from the source point to the reflection point or diffraction point, exp(-jk|r1′-r) represents the phase change of the electromagnetic wave from the source point to the reflection point or diffraction point, Indicates that the source current element is Phase changes in direction; 4b4) According to the reflection coefficient matrix [R], diffraction coefficient matrix [D] and incident field Obtain the reflection field of the first-order ray or diffraction field Among them, A(s r ) and A(s d ) represent the amplitude factors of the reflected field and the diffracted field, and represent the phase factors of the reflected field and the diffracted field respectively.

6. The method according to claim 1, characterized in that In step (4c), the final reflection field of the second-order ray at the observation point is obtained by the double uniform diffraction theory UTD algorithm based on the electromagnetic flow on the radiation source surface. or diffraction field The implementation steps include the following: 4c1) According to the electromagnetic current J on the surface of the radiation source s (p,s) Calculate the incident field at the first reflection or diffraction point Where k is the beam, η is the free space impedance, r and r2′1 represent the positions of the radiation source and the first reflection or diffraction point, respectively. represents the unit vector from the source point to the first reflection point or diffraction point, exp(-jk|r2′1-r) represents the phase change of the electromagnetic wave from the source point to the first reflection point or diffraction point, Indicates that the source current element is Phase changes in direction; 4c2) According to the reflection coefficient matrix [R], diffraction coefficient matrix [D] and incident field Obtain the reflection field at the first reflection point or diffraction point or diffraction field Among them, A(s r ) and A(s d ) represent the amplitude factors of the reflection field or diffraction field at the first reflection point or diffraction point, respectively. and Represents the phase factor of the reflected field or diffracted field at the first reflection point or diffraction point respectively; 4c3) According to the reflection field at the first reflection point or diffraction point or diffraction field Obtain the incident field at the final reflection point or diffraction point When the ray is reflected at the first occlusion, When the ray is diffracted at the first obstruction, 4c4) According to the reflection coefficient matrix [R], the diffraction coefficient matrix [D] and the incident field at the final reflection point or diffraction point The final reflection field or diffraction field Among them, A(s r ) and A(s d ) represent the amplitude factors of the final reflected field and diffracted field, and represent the phase factors of the final reflected field and diffracted field respectively.

7. The method according to claim 1, characterized in that The simulated electromagnetic field E at the observation point obtained in step (4d) s (r′), which is expressed as follows: E s (r′)=E0(r′)+E1(r′)+E2(r′) E0(r), E1(r′) and E2(r′) represent the direct field, the electromagnetic field generated by the first-order ray and the electromagnetic field generated by the second-order ray, respectively.

8. The method according to claim 1, characterized in that In step (5), accurate data is used to calibrate the simulated electromagnetic field E s (r′), and the corrected electromagnetic field is obtained The implementation steps include the following: (5a) The position of the known simulation observation point and the corresponding simulation value are expressed as (x i ,E s (x i )), the position of the measured point and the corresponding true measurement value are expressed as (x j ,E m (x j )), where i = 1, 2, ..., n, j = 1, 2, ..., m, n is the number of observation points, and m is the number of measured points; (5b) According to the measured point x j and observation point x p The distance d p,j , calculate the measured point x j and the predicted point x p The distance weighting factor ω between d (x j ): Where β is the path loss exponent of the distance weighting factor; (5c) According to the radiation source to the measured point x j The distance D j , radiation source to the measured point x j The power weighting factor ω P (x j ) is calculated as follows: Where α is the path loss exponent of the power weighting factor; (5d) According to the power weighting factor ω P (x j ) and the distance weighting factor ω d (x j ) calculates the correction weighting factor ω(x j ): ω(x j )=ω p (x j )·ω d (x j ) (5e) For the observation point x p , using the correction weighting factor ω(x j ) to correct the simulation value, the formula is: in, Represents the observation point x p The corrected electromagnetic field at the observation point x is represented by Es(xp) and Es(xj). p and x j The simulated electromagnetic field value at the observation point x is represented by Em(xj). j The actual measurement value at (5f) point x p Corrected fusion value at is considered as a measurement value and the next uncorrected observation point x p+1 Repeat the same process to calculate the corrected fusion values ​​for all observation points.

Images