A method for electromagnetic spectrum mapping in complex three-dimensional urban environments

Abstract

The present invention discloses a method for mapping electromagnetic spectrum in complex three-dimensional urban environments, comprising: determining the scope, frequency points, and signal recovery threshold of the area to be mapped; representing the spectrum data intensity corresponding to each cube obtained after discretizing the area to be mapped in tensor form; determining the number of samples in the area to be mapped, an optimized set of sampling position indexes, and determining the optimal sampling trajectory of an unmanned aerial vehicle; collecting corresponding spectrum data vectors, recovering the electromagnetic target signals in the area to be mapped, and obtaining a preliminary complete spectrum map; calculating the shadow fading components at all locations based on the sparse electromagnetic target signals and the sampled spectrum data, and correcting the spectrum map; correcting the reconstructed spectrum map based on the shadow fading components to obtain a three-dimensional electromagnetic spectrum map of the area to be mapped. The present invention comprehensively considers the characteristics of the scene and the optimization of sampling positions, and only needs to collect spectrum data from a small amount of key locations to achieve high-precision spectrum mapping.

Description

Technical Field

[0001] The present invention relates to the field of wireless information transmission, and in particular to a method for surveying and mapping electromagnetic spectrum maps for three-dimensional complex urban environments, especially for accurate and efficient spectrum mapping applications under conditions where the time, quantity, and location of spectrum data collection in three-dimensional complex scenes are limited. Background Art

[0002] Electromagnetic spectrum maps quantitatively characterize and visualize specific spectrum information from different dimensions, including received signal strength, channel gain, power spectrum density, and radiation source location. They play a vital role in a variety of fields, including dynamic spectrum intervention, radiation source location, and spectrum sharing. With the rapid development of integrated space-ground information networks and the dramatic increase in various wireless communication devices, the spectrum space has expanded from two-dimensional to three-dimensional. Therefore, the construction of three-dimensional spectrum maps is essential.

[0003] The accuracy of three-dimensional spectrum mapping depends not only on the quantity and location of sampled data but also on a precise description of the electromagnetic wave propagation model. A greater number of sampled data leads to a more accurate reconstruction of the electromagnetic spectrum of the area being measured, but this also increases the acquisition workload. In particular, the location and timing of spectrum data acquisition in complex urban scenarios are subject to numerous constraints. Three-dimensional spectrum mapping faces challenges such as limited sample quantity, limited acquisition time, and a complex and changing acquisition environment.

[0004] Existing electromagnetic spectrum mapping methods are mainly aimed at constructing two-dimensional spectrum maps. Most of them adopt random sampling or fixed sampling methods, which usually require a large amount of spectrum data to achieve sufficient accuracy. Summary of the Invention

[0005] The purpose of this invention is to propose an electromagnetic spectrum mapping method for three-dimensional complex urban environments, which comprehensively considers the scene characteristics and the optimization of sampling positions. It only needs to collect spectrum data of a small amount of key positions to achieve high-precision spectrum mapping.

[0006] In order to achieve the above technical objectives, the technical solution adopted by the present invention is:

[0007] A method for surveying and mapping electromagnetic spectrum maps in a three-dimensional complex urban environment, comprising the following steps:

[0008] S1, determine the scope, frequency and signal recovery threshold of the area to be surveyed;

[0009] S2, discretize the area to be surveyed and mapped, and express the spectrum data intensity corresponding to each cube obtained after discretization of the area to be surveyed and mapped into a tensor form;

[0010] S3, based on the signal recovery threshold set by the user, determines the number of samples in the test area, the optimized sampling position index set, and the optimal sampling trajectory of the drone;

[0011] S4, using the optimal sampling trajectory of the UAV to collect the corresponding spectrum data vector, restore the electromagnetic target signal in the test area, and obtain a preliminary completed spectrum map;

[0012] S5, calculating the shadow fading components at all positions based on the sparse electromagnetic target signal and the sampled spectrum data, correcting the spectrum map and obtaining a reconstructed spectrum map;

[0013] S6, correcting the reconstructed spectrum map according to the shadow fading component to obtain a three-dimensional electromagnetic spectrum map of the area to be mapped.

[0014] As a preferred example, in step S1, the content includes the length L, width W and height H of the region, and the frequency point f to be measured. c , signal recovery threshold λ WCEV .

[0015] As a preferred example, in step S2, the spectrum data intensity corresponding to each cube in the area to be mapped is expressed as a tensor form χ:

[0016]

[0017] Where Δx, Δy and Δz are the set spatial resolutions, and are the position coordinates of each small cube relative to the length of the surveying area, the width of the surveying area, and the height of the surveying area after discretization; the three-dimensional space coordinates of each cube are marked as The surveying area includes N=N x ×N y ×N z =(L / Δx)×(W / Δy)×(H / Δz) discrete cubes, N x 、N y and N z are the number of small cubes in the length, width and height directions of the surveying area respectively.

[0018] As a preferred example, in step S3, according to the signal recovery threshold λ set by the user, WCEV , determine the number of samples M in the area to be tested and the optimized sampling position index set And determine the optimal sampling trajectory, including:

[0019] S31, build channel propagation dictionary The elements is the path loss between the i-th cube and the j-th cube, expressed as:

[0020]

[0021] Where η represents the path loss exponent, d i，j represents the distance between the i-th cube and the j-th cube, f c Represents the frequency point to be measured; for the channel propagation dictionary matrix Perform principal component analysis preprocessing and reduce the dimension to n-dimensional principal component space to obtain

[0022] S32, solve the optimal sampling position set

[0023]

[0024] Where, λ1≥λ2≥...≥λ r is the dual observation matrix H p The eigenvalues ​​of Represents the measurement matrix, where each row has an element of 1, indicating the cube position where the sampled data is located, and the rest of the elements are 0;

[0025] S33, determine the sampling positions one by one. When selecting the t-th sampling position, assume that To determine the dual measurement matrix after the first t-1 sampling positions, according to the different conditions of t, define The minimum feature space and

[0026]

[0027] Where, Represents the normalized eigenvector corresponding to the eigenvalue, and span represents the expansion space

[0028] S34, select row vector In the dual observation matrix The t-th sampling position i is determined by maximizing the projection position on the minimum feature space t :

[0029]

[0030] Where, express The i t row vector, I r represents the identity matrix, is the set of determined sampling locations;

[0031] S35, repeat steps S33 and S34 until the cutoff condition λ is met r ≥λ WCEV , and output the final set of sampling locations

[0032] S36, based on the optimal sampling position, adopt the nearest neighbor strategy, starting from the starting point, each time find the point closest to the current point among the points that have not been traversed, as the point to be traversed next, until all points are traversed and return to the starting point, forming the optimal sampling trajectory.

[0033] As a preferred example, step S4 further includes:

[0034] S41, the unknown electromagnetic target signal ω in the area to be measured is expressed as:

[0035]

[0036] Where, Represents the transmission power of the unknown electromagnetic target; the relative position coordinates of K electromagnetic targets in the test area are

[0037] S42, the spectrum intensity information t collected at the sampling position is expressed as:

[0038]

[0039] Where, Represents the shadow fading value at each sampling position, represents the measurement noise at each sampling location, represents the approximate Gaussian distribution noise, is the vectorized form of the spectrum tensor; represents the measurement matrix, represents the channel propagation dictionary, Φ represents the perception matrix;

[0040] S43, assuming that ω obeys a two-layer sparse prior student-t distribution:

[0041]

[0042] Where α0, b0 represent shape parameters and scale parameters;

[0043] Assume that the sampled data t has a mean of 0 and a variance of Gaussian distribution of:

[0044]

[0045] Where, the hyperparameter α=[α1,α2,...,α N ] T , All obey the Gamma distribution, c0, d0 represent the shape parameter and scale parameter;

[0046] Compute the posterior distribution of the electromagnetic target signal:

[0047]

[0048] in, The mean μ is an estimate of ω;

[0049] S44, calculate the update formula of hyperparameters α and β as follows:

[0050]

[0051] in

[0052] S45, calculate α n The iteration criteria of β and β are:

[0053]

[0054]

[0055] Among them n ≡1-α n ∑ nn , μ n represents the nth element in μ, and M represents the number of sampling points;

[0056] S46, repeating steps S43 to S45, iteratively updating the hyperparameters to obtain an estimated value of the electromagnetic target signal And the preliminary recovered spectrum map

[0057] As a preferred example, step S5 further includes:

[0058] S51, calculate sampling position The shadow fading component of

[0059]

[0060] Where s is the simplified representation of the sampling position, t(s) is the shadow fading component of all sampling positions, Φ m，n is the mth row and nth column element of Φ, for The nth element of ;

[0061] S52, construct shadow fading Gaussian process model represents the covariance function, defined as the Mattern covariance function:

[0062]

[0063] where K g (·) denotes the Bessel function of the second kind, g and ρ denote the order and nonnegative spatial decay parameter, d denotes the distance between sampling locations s, σ denotes the marginal standard deviation, and Γ(g) denotes the gamma function with parameter g;

[0064] S53, calculate the parameters to be determined based on maximum likelihood estimation The optimization goal to be solved is:

[0065]

[0066] in is the model output, is a random variable representing noise, δ m is the mth element in δ, is the noise variance of the Gaussian process model;

[0067] S54, solve and obtain parameters according to the gradient descent method The estimated value of is:

[0068]

[0069] S55, for unsampled position v * The shadow fading component value Make a prediction, which obeys the following distribution:

[0070]

[0071] Where, v is the abbreviated form of all cube coordinates, s is the abbreviated form of the sampling position, i * is the coordinate index, represents the predicted mean, represents the prediction variance;

[0072] S56, calculating the mean value based on the model parameters solved in step S54 Get the shadow fading component of the unsampled position And get the shadow fading components of all positions

[0073] As a preferred example, step S6 further includes:

[0074] According to the shadow fading component ξ v Corrected and reconstructed electromagnetic spectrum map Finally, a three-dimensional electromagnetic spectrum map of the area to be measured is obtained:

[0075]

[0076] in represents the Hadamard product.

[0077] Compared with the prior art, the present invention has the following beneficial effects:

[0078] First, the electromagnetic spectrum mapping method of the present invention for three-dimensional complex urban environments obtains spectrum data by selecting a small number of samples and optimized sampling positions, and then restores the spectrum data of unsampled positions according to the laws of radio wave propagation. It is a compromise solution that takes into account both performance and efficiency.

[0079] Second, the electromagnetic spectrum mapping method of the present invention for three-dimensional complex urban environments takes into account the importance of different spatial positions and the amount of useful information they can provide, optimizes the sampling positions, and designs the optimal sampling trajectory of the drone based on the optimized sampling positions, thereby achieving efficient spectrum mapping.

[0080] Third, the electromagnetic spectrum mapping method for three-dimensional complex urban environments of the present invention takes into account the characteristics of three-dimensional complex scenes and electromagnetic propagation characteristics. Based on the compressed sensing model and Bayesian theory, it first uses the electromagnetic target signals of the area to be measured to preliminarily construct a spectrum map, and then uses the collected spectrum data to estimate the shadow fading in the environment, and uses it to correct the spectrum map to improve the accuracy of three-dimensional spectrum data recovery. BRIEF DESCRIPTION OF THE DRAWINGS

[0081] Figure 1 This is a flow chart of the electromagnetic spectrum mapping method for three-dimensional complex urban environments of the present invention;

[0082] Figure 2 Schematic diagram of the area to be tested for the implementation case.

[0083] Figure 3 Schematic diagram of spectrum data sampled for the implementation case optimization.

[0084] Figure 4 Schematic diagram of the electromagnetic target signal estimation results of the implementation case.

[0085] Figure 5 Schematic diagram of the final output results of the spectrum map implementation case. DETAILED DESCRIPTION

[0086] The embodiments of the present invention are described in further detail below with reference to the accompanying drawings.

[0087] The present invention discloses an electromagnetic spectrum map surveying and mapping method for a three-dimensional complex urban environment, the electromagnetic spectrum map surveying and mapping method comprising the following steps:

[0088] S1, determine the scope, frequency and signal recovery threshold of the area to be surveyed;

[0089] S2, discretize the area to be surveyed and mapped, and express the spectrum data intensity corresponding to each cube obtained after discretization of the area to be surveyed and mapped into a tensor form;

[0090] S3, based on the signal recovery threshold set by the user, determines the number of samples in the test area, the optimized sampling position index set, and the optimal sampling trajectory of the drone;

[0091] S4, using the optimal sampling trajectory of the UAV to collect the corresponding spectrum data vector, restore the electromagnetic target signal in the test area, and obtain a preliminary completed spectrum map;

[0092] S5, calculating the shadow fading components at all positions based on the sparse electromagnetic target signal and the sampled spectrum data, correcting the spectrum map and obtaining a reconstructed spectrum map;

[0093] S6, correcting the reconstructed spectrum map according to the shadow fading component to obtain a three-dimensional electromagnetic spectrum map of the area to be mapped.

[0094] See also Figure 1 , the electromagnetic spectrum mapping method comprises the following steps:

[0095] Step 1: Determine the range of the area to be tested, the frequency point and the signal recovery threshold, including the area length L, width W and height H, the frequency point to be tested f c , signal recovery threshold λ WCEV .

[0096] Step 2: Discretize the area to be measured and set the spatial resolution Δx, Δy and Δz. The position coordinates of each small cube relative to the length of the survey area after discretization Position coordinates relative to the width of the survey area and the position coordinates relative to the height of the survey area The spectrum data intensity corresponding to each cube in the area to be mapped is expressed as a tensor:

[0097]

[0098] At this time, the surveying area includes N=N x ×N y ×N z =(L / Δx)×(W / Δy)×(H / Δz) discrete cubes.

[0099] Step 3: Set the parameters according to the user, that is, the signal recovery threshold λ WCEV , determine the number of samples M in the area to be tested and the optimized sampling position index set And determine the optimal sampling trajectory. The specific implementation steps are as follows:

[0100] 3.1) Constructing channel propagation dictionary The elements is the path loss between the i-th cube and the j-th cube, which can be expressed as:

[0101]

[0102] Where η represents the path loss exponent, d i，j Represents the distance between the i-th cube and the j-th cube. At the same time, the channel propagation dictionary matrix Perform principal component analysis preprocessing and reduce the dimension to r-dimensional principal component space to obtain

[0103] 3.2) Solve the optimal sampling location set The specific optimization methods are as follows:

[0104]

[0105] where λ1≥λ2≥…≥λ r H p The eigenvalues ​​of represents the dual observation matrix, Represents a measurement matrix, where each row has an element of 1, indicating the cube position where the sampled data is located, and the rest of the elements are 0.

[0106] 3.3) Determine the sampling positions one by one. When selecting the t-th sampling position, assume that To determine the dual measurement matrix after the first t-1 sampling positions, according to the different conditions of t, define The minimum feature space of :

[0107]

[0108] in, Represents the normalized eigenvector corresponding to the eigenvalue.

[0109] 3.4) Select row vector In the dual observation matrix The position where the projection on the minimum feature space is maximized is determined by t , the specific optimization methods are as follows:

[0110]

[0111] in, express The i t row vector,

[0112] 3.5) Repeat steps 3.3) and 3.4) until the cutoff condition λ is met r ≥λ WCEV , and output the final set of sampling locations

[0113] 3.6) Based on the optimal sampling location, the optimal flight trajectory of the drone is planned. Using the nearest neighbor strategy, starting from the starting point, the nearest point to the current point is found each time, which is the next point to be traversed. This continues until all points are traversed and the drone returns to the starting point, ultimately forming the optimal sampling trajectory.

[0114] Step 4: Use the flight trajectory to collect the corresponding spectrum data vector Restore the electromagnetic target signal in the area to be measured and obtain a preliminary complete spectrum map. The specific implementation steps are as follows:

[0115] 4.1) The unknown electromagnetic target signal in the area to be measured is expressed as:

[0116]

[0117] in Represents the transmission power of the unknown electromagnetic target. The relative position coordinates of K electromagnetic targets in the test area are

[0118] 4.2) The spectrum intensity information t collected at the corresponding position is expressed as:

[0119]

[0120] in, Represents the shadow fading value at each sampling position, represents the measurement noise at each sampling location, is the vectorized form of the spectrum tensor.

[0121] 4.3) Since the electromagnetic target signal ω is a sparse signal compared to the entire spectrum space, it is assumed that ω obeys a two-layer sparse prior student-t distribution, that is:

[0122]

[0123] Assume that the sampled data t has a mean of 0 and a variance of Gaussian distribution, that is:

[0124]

[0125] where the hyperparameter α=[α1,α2,…,α N ] T , All obey the Gamma distribution.

[0126] 4.4) Calculate the posterior distribution of the electromagnetic target signal as follows:

[0127]

[0128] in, The mean μ is an estimate of ω.

[0129] 4.5) Calculate the update formula for the hyperparameters α and β as follows:

[0130]

[0131] in

[0132] 4.6) Calculate α n The iteration criteria of β and β are:

[0133]

[0134]

[0135] Among them n ≡1-α n ∑ nn .

[0136] 4.7) Repeat steps 4.4) to 4.6) to iteratively update the hyperparameters to obtain the electromagnetic target signal estimate And the preliminary recovered spectrum map

[0137] Step 5: Based on sparse electromagnetic target signals And the sampled spectrum data, calculate the shadow fading component ξ at all locations v , correct the spectrum map and obtain the final reconstructed spectrum map x. The specific implementation steps are as follows:

[0138] 5.1) Calculate sampling position The shadow fading component of the s is:

[0139]

[0140] 5.2) Constructing a shadow fading Gaussian process model represents the covariance function, defined as the Mattern covariance function:

[0141]

[0142] where K g (·) denotes the Bessel function of the second kind, g and ρ denote the order and nonnegative spatial decay parameter, and d denotes the distance between sampling positions s.

[0143] 5.3) Calculate the parameters to be determined based on maximum likelihood estimation The optimization goal to be solved is:

[0144]

[0145] in is the model output,

[0146] 5.4) Solve according to the gradient descent method and obtain the parameters The estimated value of is:

[0147]

[0148] 5.5) For unsampled positions The shadow fading component value Make a prediction, which obeys the following distribution:

[0149]

[0150] 5.6) Calculate the mean value based on the model parameters solved in step 5.4) Get the shadow fading component of the unsampled position And get the shadow fading components of all positions

[0151] Step 6: According to the shadow fading component ξ v Correction of the electromagnetic spectrum map initially restored in step 5 Finally, a three-dimensional electromagnetic spectrum map of the area to be measured is obtained as follows:

[0152]

[0153] in represents the Hadamard product.

[0154] Specific implementation cases

[0155] There are 8 transmitters in the test area of ​​this embodiment, and the location parameters, transmission frequency, and transmission power of the transmitters are shown in Table 1.

[0156] Table 1 Transmitter parameters

[0157]

[0158] Step 1: The user sets the length of the test area L = 1250m, the width W = 1250m and the height H = 60m, and the center frequency point f c =2450MHz, signal recovery threshold λ WCEV =60, the scene of the area to be tested is as follows Figure 2 shown.

[0159] Step 2: The user sets the spatial resolution to 5m×5m×10m and determines the total number of cubes after discretization N=N x ×N y ×N z =250×250×6.

[0160] Step 3: Input the position set of all cubes as the candidate sampling position set and determine the optimized sampling position and number. The specific implementation steps are as follows:

[0161] 3.1) Use formula (2) to construct the channel dictionary The path loss exponent is set to η=2.

[0162] 3.2) Perform principal component analysis on the channel dictionary to obtain Dual observation matrix H p =(Φ p ) T Φ p .

[0163] 3.3) Use formula (4) to calculate the minimum eigenspace of the dual observation matrix under different conditions.

[0164] 3.4) Use formula (5) to select the sampling position index i with the largest projection on the minimum feature space t .

[0165] 3.5) Repeat steps 3.3) to 3.4) and iteratively select sampling points until the minimum eigenvalue λ of the current dual observation matrix is r Satisfy λ r <λ WCEV , end the iteration and output the final set of sampling positions like Figure 3 shown.

[0166] 3.6) According to the optimal sampling position Enter the coordinates of the location to be sampled. Use the nearest neighbor strategy to obtain the optimal sampling trajectory.

[0167] Step 4: According to the flight trajectory designed in step 3, obtain the corresponding spectrum data vector Restore the position of the electromagnetic target signal in the area to be measured. The specific implementation steps are as follows:

[0168] 4.1) Model the sampled spectrum data according to formula (7), and initialize the parameter a0=10 -6 , b0=10 -6 , c0=10 -6 , d0=10 -6 .

[0169] 4.2) Calculate the estimated mean μ and variance ∑ of the electromagnetic target signal according to formula (10): ω .

[0170] 4.3) Update parameters α and β according to formulas (12) and (13).

[0171] 4.4) Repeat steps 4.2) to 4.3) until the cutoff condition is met and the electromagnetic target signal estimation value is obtained like Figure 4 As shown, the preliminary reconstructed electromagnetic spectrum map is calculated

[0172] Step 5: Based on electromagnetic target signals And the sampled spectrum data t, estimate the shadow fading component ξ at all locations v , correct The final reconstructed spectrum map x is obtained. The specific implementation steps are as follows:

[0173] 5.1) Calculate the shadow fading component at the sampling position according to formula (14).

[0174] 5.2) Construct the covariance matrix of the Gaussian process according to formula (15), where g = 3 / 2.

[0175] 5.3) According to formula (16), a Gaussian process model is established for the shadow component at the sampling position to obtain the likelihood function.

[0176] 5.4) Solve the Gaussian process regression parameters according to formula (17)

[0177] 5.5) Estimate the shadow fading component value at the unsampled position based on formula (18) and the model parameters solved above Get the shadow fading components at all locations

[0178] Step 6: Exploiting Electromagnetic Target Signals and the shadow fading component ξ v , calculate the three-dimensional spectrum data of the entire test area according to formula (19), and reconstruct the spectrum map as follows Figure 5 shown.

[0179] Those skilled in the art will appreciate that the embodiments of the present application can be provided as methods, systems, or computer program products. Therefore, the application can adopt the form of a complete hardware embodiment, a complete software embodiment, or an embodiment in combination with software and hardware. Moreover, the application can adopt the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) that contain computer-usable program code. The scheme in the embodiment of the present application can be implemented in various computer languages, for example, object-oriented programming language Java and literal translation scripting language JavaScript, etc.

[0180] The present application is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or box in the flowchart and / or block diagram, as well as the combination of the processes and / or boxes in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the steps in the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0181] These computer program instructions may also be stored in a computer readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 a process or multiple processes and / or boxes Figure 1 The function specified in one or more boxes.

[0182] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing the instructions for executing the steps in the process. Figure 1 a process or multiple processes and / or boxes Figure 1 A step that specifies a function in one or more boxes.

[0183] Although the preferred embodiments of the present application have been described, those skilled in the art may make additional changes and modifications to these embodiments once they have learned the basic creative concept. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments and all changes and modifications that fall within the scope of the present application.

[0184] Obviously, those skilled in the art may make various changes and modifications to this application without departing from the spirit and scope of this application. Thus, if these modifications and variations of this application fall within the scope of the claims of this application and their equivalents, this application is intended to include these modifications and variations.

Claims

1. A method for electromagnetic spectrum mapping in a three-dimensional complex urban environment, characterized in that: The electromagnetic spectrum mapping method comprises the following steps: S1, determine the scope, frequency and signal recovery threshold of the area to be surveyed; S2, discretize the area to be surveyed and mapped, and express the spectrum data intensity corresponding to each cube obtained after discretization of the area to be surveyed and mapped into a tensor form; S3, based on the signal recovery threshold set by the user, determines the number of samples in the test area, the optimized sampling position index set, and the optimal sampling trajectory of the drone; S4, using the optimal sampling trajectory of the UAV to collect the corresponding spectrum data vector, restore the electromagnetic target signal in the test area, and obtain a preliminary completed spectrum map; S5, calculating the shadow fading components at all positions based on the sparse electromagnetic target signal and the sampled spectrum data, correcting the spectrum map and obtaining a reconstructed spectrum map; S6, correcting the reconstructed spectrum map according to the shadow fading component to obtain a three-dimensional electromagnetic spectrum map of the area to be mapped.

2. The electromagnetic spectrum mapping method for a three-dimensional complex urban environment according to claim 1, characterized in that: In step S1, the content includes the length L, width W and height H of the area, the frequency point f to be measured c , signal recovery threshold λ WCEV .

3. The electromagnetic spectrum mapping method for a three-dimensional complex urban environment according to claim 1, characterized in that: In step S2, the spectrum data intensity corresponding to each cube in the area to be mapped is expressed as a tensor form χ: Where Δx, Δy and Δz are the set spatial resolutions, and are the position coordinates of each small cube relative to the length of the surveying area, the width of the surveying area, and the height of the surveying area after discretization; the three-dimensional space coordinates of each cube are marked as The surveying area includes N=N x ×N y ×N z =(L / Δx)×(W / Δy)×(H / Δz) discrete cubes, N x 、N y and N z are the number of small cubes in the length, width and height directions of the surveying area respectively.

4. The electromagnetic spectrum mapping method for a three-dimensional complex urban environment according to claim 1, characterized in that: In step S3, the signal recovery threshold λ set by the user is WCEV , determine the number of samples M in the area to be tested and the optimized sampling position index set And determine the optimal sampling trajectory, including: S31, build channel propagation dictionary The elements is the path loss between the i-th cube and the j-th cube, expressed as: Where η represents the path loss exponent, d i，j represents the distance between the i-th cube and the j-th cube, f c Represents the frequency point to be measured; for the channel propagation dictionary matrix Perform principal component analysis preprocessing and reduce the dimension to n-dimensional principal component space to obtain S32, solve the optimal sampling position set Where, λ1≥λ2≥...≥λ r is the dual observation matrix H p The eigenvalues ​​of Represents the measurement matrix, where each row has an element of 1, indicating the cube position where the sampled data is located, and the rest of the elements are 0; S33, determine the sampling positions one by one. When selecting the t-th sampling position, assume that To determine the dual measurement matrix after the first t-1 sampling positions, according to the different conditions of t, define The minimum feature space and Where, Represents the normalized eigenvector corresponding to the eigenvalue, and span represents the expansion space; S34, select row vector In the dual observation matrix The t-th sampling position i is determined by maximizing the projection position on the minimum feature space t : Where, express The i t row vector, I r represents the identity matrix, is the set of determined sampling locations; S35, repeat steps S33 and S34 until the cutoff condition λ is met r ≥λ WCEV , and output the final set of sampling locations S36, based on the optimal sampling position, adopt the nearest neighbor strategy, starting from the starting point, each time find the point closest to the current point among the points that have not been traversed, as the point to be traversed next, until all points are traversed and return to the starting point, forming the optimal sampling trajectory.

5. The electromagnetic spectrum mapping method for a three-dimensional complex urban environment according to claim 1, characterized in that: Step S4 further includes: S41, the unknown electromagnetic target signal ω in the area to be measured is expressed as: Where, Represents the transmission power of the unknown electromagnetic target; the relative position coordinates of K electromagnetic targets in the test area are S42, the spectrum intensity information t collected at the sampling position is expressed as: Where, Represents the shadow fading value at each sampling position, represents the measurement noise at each sampling location, represents the approximate Gaussian distribution noise, is the vectorized form of the spectrum tensor; represents the measurement matrix, represents the channel propagation dictionary, Φ represents the perception matrix; S43, assuming that ω obeys a two-layer sparse prior student-t distribution: Where a0, b0 represent shape parameters and scale parameters; Assume that the sampled data t has a mean of 0 and a variance of Gaussian distribution of: Where, the hyperparameter α=[α1,α2,...,α N ] T , All obey the Gamma distribution, c0, d0 represent the shape parameter and scale parameter; Compute the posterior distribution of the electromagnetic target signal: in, The mean μ is an estimate of ω; S44, calculate the update formula of hyperparameters α and β as follows: in S45, calculate α n The iteration criteria of β and β are: where γ n ≡1-α n ∑ nn , μ n represents the nth element in μ, and M represents the number of sampling points; S46, repeating steps S43 to S45, iteratively updating the hyperparameters to obtain an estimated value of the electromagnetic target signal And the preliminary recovered spectrum map 6. The electromagnetic spectrum mapping method for a three-dimensional complex urban environment according to claim 1, characterized in that: Step S5 further includes: S51, calculate sampling position The shadow fading component of Where s is the simplified representation of the sampling position, t(s) is the shadow fading component of all sampling positions, Φ m，n is the mth row and nth column element of Φ, for The nth element of ; S52, construct shadow fading Gaussian process model represents the covariance function, defined as the Mattern covariance function: where K g (·) denotes the Bessel function of the second kind, g and ρ denote the order and nonnegative spatial decay parameter, d denotes the distance between sampling locations s, σ denotes the marginal standard deviation, and Γ(g) denotes the gamma function with parameter g; S53, calculate the parameters to be determined based on maximum likelihood estimation The optimization goal to be solved is: in is the model output, m=1,2,...,M is a random variable representing noise, δ m is the mth element in δ, is the noise variance of the Gaussian process model; S54, solve and obtain parameters according to the gradient descent method The estimated value of is: S55, for unsampled position v * The shadow fading component value Make a prediction, which obeys the following distribution: Where, v is the abbreviated form of all cube coordinates, s is the abbreviated form of the sampling position, i * is the coordinate index, represents the predicted mean, represents the prediction variance; S56, calculating the mean value based on the model parameters solved in step S54 Get the shadow fading component of the unsampled position And get the shadow fading components of all positions 7. The electromagnetic spectrum mapping method for a three-dimensional complex urban environment according to claim 1, characterized in that: Step S6 further includes: Correction and reconstruction of electromagnetic spectrum map based on shadow fading component ξv Finally, a three-dimensional electromagnetic spectrum map of the area to be measured is obtained: in represents the Hadamard product.

Images