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

Abstract

The invention discloses a large-scale urban electromagnetic environment analysis method based on measurement and calculation fusion. The problems that in the prior art, calculation cost is too high, calculation resource consumption is too large, timeliness does not exist, and precision is low are mainly solved. The implementation scheme is as follows: extracting coordinates of all source points and observation points in an urban environment and radiation source model graph; an observation point is distributed for each CPU by adopting a parallel computing strategy based on a message passing interface, and ray tracing and shielding judgment from the observation point to a source point are carried out; calculating a radiation source surface electromagnetic current by using a high-order moment method; calculating a reflection field and a diffraction field at the shielded part according to the electromagnetic current on the surface of the radiation source, and then obtaining the simulation field intensity of the observation point through field intensity superposition; in a real city environment, a frequency spectrum sensor is used to measure a selected actual measurement point, and an accurate measurement result is used to correct a simulation field intensity obtained by fusion simulation so as to obtain a more accurate fusion field intensity. The method greatly shortens the calculation time, improves the calculation precision, and can be used for communication spectrum planning and electromagnetic environment evaluation.

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, more and more radiation sources and wireless communication devices have been integrated into the daily lives of urban residents, including traditional radio and television to various modern mobile communication devices, as well as the 5G and future 6G networks that will be widely deployed. These devices have increased the complexity of the urban electromagnetic environment. With the increase in the number of wireless devices, specific frequency bands will become crowded, resulting in signal interference and reduced communication efficiency. In addition, the interference of strong electromagnetic radiation also poses a potential threat to the normal operation of key infrastructure such as precision instruments, communication networks and traffic control systems. Therefore, it is crucial to accurately and effectively analyze the complex electromagnetic environment in the city. This analysis helps to 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 have many problems such as high computational cost, excessive consumption of computing resources and time, lack of timeliness and low accuracy.

[0003] The patent document with application number: CN201910012530.0 discloses a method for predicting the situation of an outdoor electromagnetic environment in an urban area based on MoM-UTD, and its implementation steps are as follows: 1) Establishing an urban simulation model and an antenna simulation model, extracting geometric information and electrical parameter information in the urban simulation model, and extracting the operating frequency and source point position of the antenna simulation model; 2) Setting the sampling points of the route to be measured as field points; 3) Tracking all rays from the source point to the field point, storing ray information and action points; 4) Using the moment method MoM to calculate the radiated electric field of the antenna simulation model at the action point; 5) Using the UTD method to calculate the radiated electric field from the action point to the field point of all ray types; 6) Superimposing the radiated electric fields of all ray types at the same field point to obtain the prediction results. Although this method has low cost, accurate prediction and high efficiency, it still has the problems of excessive calculation amount, insufficient calculation speed and limited accuracy of the final calculation result. 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 requirements 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 moment method-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, and 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, and the observation points are alternately assigned to each CPU for ray tracing and occlusion judgment according to the sorting. Three paths are obtained: reaching the observation point directly from the radiation source, reaching the observation point after one reflection or diffraction, and reaching the observation point after two reflections or diffractions. The coordinates of the reflection point and the diffraction point are 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) 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′);

[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 at 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 allocate observation points to each central processing unit CPU for ray tracing and occlusion judgment, which not only greatly speeds up 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 a flow chart of the implementation 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 It 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 of all observation points calculated in an embodiment of the present invention;

[0025] Figure 5 It is a physical picture of the radiation source arranged for measuring the real value of the electromagnetic field in the embodiment of the present invention;

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

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

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

[0029] Fig. 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] Fig.10 The sampling point distribution in 409.5MHz scenario B and the comparison of fusion value, simulation value and measurement value;

[0031] Fig.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] Fig.12 This is a comparison chart of the sampling point distribution and fusion value, simulation value and measurement value in 2501MHz scenario D;

[0033] Fig.13It 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 scheme of the present invention, the technical scheme in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments 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 creative work should all 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 in an accurate model of the urban environment and radiation sources.

[0038] 1.1) Extract the location and height parameters of main buildings and facilities in the city model map:

[0039] The city model map uses a CAD model provided by a surveying company, and extracts information about the location, floor level and height of major buildings and facilities, such as Figure 2 As shown, among which Figure 2 (a) is the 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] In order 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, and uses an omnidirectional antenna with a gain of -2 dBi and a power of 45 dBm to establish radiation sources of these frequencies. The operating frequency range of the antenna is extended from 20 MHz to 6 GHz.

[0042] Step 2, calculate the simulated electromagnetic field based on the high-order moment method-uniform geometric diffraction theory hybrid algorithm of the message passing interface.

[0043] 2.1) Perform ray tracing and occlusion judgment according to 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 the tracing needs to be re-performed for different observation points, and the number of observation points in this example is huge, if the serial algorithm is used for tracing, there will be problems such as large amount of calculation and long calculation time. 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, so as to identify each valid ray path that can reach the source point from the observation point and calculate the coordinates of all reflection points and diffraction points.

[0045] This example uses 4096 cores with frequencies of 409.5 MHz and 2501 MHz. The relative dielectric constant 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: reaching the observation point directly from the radiation source, reaching the observation point after one reflection or diffraction, and reaching the observation point after two reflections or diffractions.

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

[0047] 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. High-order geometric patches can flexibly approximate any surface, and 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 current M s (p,s), which is calculated as follows:

[0048]

[0049] in, are polynomial basis functions, is a scaling factor, p i s j is the specific term of 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 surface current distribution Js(p,s) of the radiation source:

[0051]

[0052] Where k is the beam, η is the free space wave impedance, 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 by 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 point 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 propagating 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 incident angle and the transmission angle 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 relevant 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 associated with the incident angle α0,

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

[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) According to the two cases of reflection or diffraction at the obstruction, the electromagnetic field E1(r′) generated by the first-order ray at the observation point is obtained:

[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 UTD algorithm of the double uniform diffraction theory:

[0080] 2.5.1) 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

[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 propagating 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 factor of the reflection field or diffraction field at the first reflection point or diffraction point, respectively.

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

[0087] 2.5.3) According to whether reflection or diffraction occurs at the first obstruction, the incident field at the final reflection point or diffraction point is obtained.

[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) According to the two situations of reflection or diffraction at the second obstruction, the electromagnetic field E1(r′) generated by the second-order ray at the observation point is obtained:

[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] In this example, there are 173,591 observation points in total. The simulation electromagnetic field solution calculation is performed on all observation points to obtain the static electromagnetic field distribution of a 60 square kilometer 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 values ​​to correct and fuse the simulated values.

[0101] 3.1) Use spectrum sensors to measure accurate electromagnetic field data at actual measurement points in the urban electromagnetic 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 is inevitably different from the true value. To this end, these differences are corrected by the fusion of measurement and simulation.

[0103] In order to realize the rapid collection, transmission and analysis of real electromagnetic data in urban blocks, this example uses a spectrum sensor with ultra-wideband electromagnetic signal detection capability to measure the real urban electromagnetic environment. At the same time, to ensure the consistency between the measurement and simulation scenes, the radiation sources with center frequencies of 409.5MHz and 2501MHz are installed at the positions corresponding to the radiation sources in the simulation scenes, 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, Figure 2 The simulation model of the urban area is divided into 60 sub-areas, each of which is about 1 square kilometer in area. A total of 885 measurement points are established in all sub-areas. Spectrum sensors are distributed in each sub-area and continuously measured at fixed measurement points, collecting data once a second. Each point is measured for 3 hours to obtain the measured real electromagnetic field data. Figure 6 The distribution of measurement points and spectrum sensors in a sub-area in the 409.5MHz measurement scenario is shown.

[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′) is fused and corrected:

[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 betweenp,j , calculate the measured point x j With the predicted point x p The distance weighting factor ω 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 ) to calculate 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 simulation 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), respectively. p and x j The simulated electromagnetic field value at the observation point x is represented by Em(xj). j The actual measured value at

[0118] 3.2.6) Change point x p Corrected fusion value at is considered as the measured 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 corrected and fused observation points, we can get Figure 2 The fused 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 in each area. The sampling points of scene A are as follows: Fig. 9 As shown in (a), the sampling points of scene B are as follows Fig.10 As shown in (a), the sampling points of scene C are as follows Fig.11 (a) shown.

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

[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 according to the difference between the fusion value, simulation value and the actual measurement value:

[0125] Comparison 1: Yes Fig. 9 (a) Verify that the fusion value, simulation value, and measurement value of the sampling point in scenario A are plotted as line graphs of the sampling point and electric field value. The results are shown in the figure. Fig. 9 (b) as shown.

[0126] from Fig. 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 Fig.10 (a) Verify that the fusion value, simulation value, and measurement value of the sampling point in scene B are plotted as line graphs of the sampling point and the electric field value. The results are shown in the figure. Fig.10 (b) as shown.

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

[0129] Comparison 3: Yes Fig.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 shown in the figure. Fig.11 (b) as shown.

[0130] from Fig.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: Yes Fig.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 shown in the figure. Fig.12 (b) as shown.

[0132] from Fig.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 Fig.13 (a) The fusion value, simulation value and measurement value of the sampling point in the verification scenario E are plotted as line graphs of the sampling point and the electric field value. The results are shown in the figure. Fig.13 (b) as shown.

[0134] from Fig.13 (b) It can be seen that the simulation value line graph is obviously very different from the measured value line graph, but the fusion value at each sampling point can accurately simulate the measured value in the real environment.

[0135] The root mean square error (RMSE) between the fusion value and simulation value of the sampling points in the above five scenarios and the actual measured value 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] It can be seen from Table 1 that in different scenarios, the root mean square error RMSE between the fusion value and the actual measurement value is lower than the root mean square error RMSE between the simulation value and the actual measurement value, and is maintained 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 new measurement and simulation fusion technology proposed in the present invention can accurately simulate the electromagnetic environment of a large-scale city 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 implementation scheme of the present invention to facilitate understanding, and the order of the numbers is not limited.

Claims

1. A large-scale urban electromagnetic environment analysis technology based on measurement and calculation fusion, characterized in that: 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, and the observation points are alternately assigned to each CPU for ray tracing and occlusion judgment according to the sorting. Three paths are obtained: reaching the observation point directly from the radiation source, reaching the observation point after one reflection or diffraction, and reaching the observation point after two reflections or diffractions. The coordinates of the reflection point and the diffraction point are calculated. (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); (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), respectively. p and x j The simulated electromagnetic field value at the observation point x is represented by Em(xj). j The actual measured 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 the measured point x is ωd(xj). 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. A high-precision reverse ray tracing method is used. The image method and the principle of geometric optics are used to identify each ray from the observation point to the source point, 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 current M s (p,s), the formulas are as follows: in, are polynomial basis functions, is a scaling factor, p i s j is the specific term of 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 directly 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 wave impedance, 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 positions of the source point and the observation point respectively, is the unit vector from the source point to the observation point.

5. The method according to claim 1, characterized in that In step (4b), the reflection field of the first-order ray is calculated by the uniform diffraction theory UTD algorithm according to the electromagnetic current on the radiation source surface. 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 incident angle and the transmission angle 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 relevant 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 associated with the incident angle α0, is the incident angle α n The relevant hard diffraction coefficient is 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 radiation source surface s (p,s) Calculate the incident field at the reflection point or diffraction point Where k is the beam, η is the free space wave impedance, r and r1′ represent the positions of the radiation source point and the reflection point 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 propagating 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], the diffraction coefficient matrix [D] and the 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 UTD algorithm based on the electromagnetic current on the radiation source surface. or diffraction field The implementation steps include the following: 4c1) 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 Where k is the beam, η is the free space wave impedance, r and r2′1 represent the positions of the radiation source point and the first reflection point 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 propagating 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], the diffraction coefficient matrix [D] and the 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 factor of the reflection field or diffraction field at the first reflection point or diffraction point, respectively. and Respectively represent the phase factor of the reflected field or diffracted field at the first reflection point or diffraction point; 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) Based on 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′), 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 rays and the electromagnetic field generated by the second-order rays, 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 between p,j , calculate the measured point x j and the predicted point x p The distance weighting factor ω 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 ) to calculate 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), respectively. p and x j The simulated electromagnetic field value at the observation point x is represented by Em(xj). j The actual measured value at (5f) point x p Corrected fusion value at is considered as the measured 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