Environment map construction method for arbitrary continuous irregular rough diffuse reflection scene

Abstract

The application belongs to the field of wireless communication and perception, and particularly relates to an environment map construction method for an arbitrary continuous irregular rough diffuse reflection scene. The application comprises: modeling a receiving end receiving signal with a continuous rough diffuse reflection probability model; using a two-dimensional off-network sparse Bayesian algorithm to extract channel parameter information from the receiving signal; perceiving the state of scattering points by using parameter information extracted from different observations; fusing the perceived scattering point information of multiple users and multiple base stations, performing data cleaning and clustering, and realizing perception of an arbitrary shape unknown environment including irregular scattering surfaces; and constructing a high-precision physical map by using a rotating polynomial fitting method. The application can greatly reduce the cost of environment map construction and prevent privacy leakage; can be applied to any unknown shape physical scattering scene, and can effectively extract and fuse a higher order of magnitude of environmental parameter state observations.

Description

Technical Field

[0001] This invention belongs to the field of wireless communication and sensing technology, specifically relating to a method for constructing environmental maps for arbitrary continuous, irregular, rough diffuse reflection scenes. Background Technology

[0002] With the rapid development of the national economy and information technology, emerging applications such as smart homes, autonomous driving, and smart cities are flourishing, and the Internet of Things is gradually being realized. This places higher demands on environmental perception technology, and high-precision positioning and tracking, simultaneous positioning and map building, and other related tasks have become hot research topics in various technical fields.

[0003] Past research has yielded numerous solutions based on single or multi-sensor fusion technologies such as LiDAR, computer vision, and inertial measurement units, achieving promising results in devices like connected cars and robotic vacuum cleaners. However, these devices remain expensive, the potential privacy risks associated with computer vision are still highly controversial, and it is also sensitive to lighting conditions and extreme weather.

[0004] With the continuous upgrading of wireless communication systems and the enhancement of communication capabilities, the application of large-scale user access, multi-antenna arrays, and higher-frequency signals has gradually made it possible to achieve higher-precision environmental perception using wireless communication. Deep integration of large-scale access and sensing-based localization is highly integrated with existing wireless communication systems at the algorithm level. It utilizes the transmission signals from a large number of communication devices over a wide area to provide environmental perception capabilities that are not constrained by time and have a wide coverage, achieving a balance between high-quality communication and high-precision perception. Compared to sensor-based SLAM solutions, this method is unaffected by lighting conditions, weather, etc. It can fully utilize the access signals from low-cost users widely distributed in time and space within ultra-large-scale connections, and achieve communication and cooperation between different devices through base stations without adding expensive equipment. Simultaneously, it can effectively avoid the security risks such as privacy leaks associated with visual solutions.

[0005] Currently, the technology of integrating communication SLAM and communication sensing has deeply empowered cutting-edge industries such as intelligent connected vehicles, smart homes, and smart consumer electronics, becoming an important new driving force for building a world of interconnected things. However, many existing key technologies are still in the early stages of development. A large number of low-cost communication devices with limited transmission power are widely distributed, and their data packet frame lengths are generally not high. Existing communication-driven sensing technologies are mostly based on the specular reflection model that relies on virtual base stations (H. Wymeersch et al., “5G mm wave downlink vehicular positioning,” in Proc. IEEE Global Commun. Conf. (GLOBECOM), Dec. 2018, pp. 206–212.). While assuming an overly ideal environment, the environmental sensing information they can provide is also very limited, making it impossible to achieve continuous environmental sensing. Under the assumption of continuous diffuse reflection, there is also a solution based on tensor modeling (Wen F, Kulmer J, Witrisal K, et al.5GPositioning and Mapping With Diffuse Multipath[J].IEEE Transactions on Wireless Communications,2020,PP(99):1-1.DOI:10.1109 / TWC.2020.3031180.), which provides excellent sensing performance. However, the limited number of factors that can be decomposed by the tensor decomposition step means that its performance needs to be improved when sensing irregular scatterers. Summary of the Invention

[0006] The purpose of this invention is to provide a communication-driven method for constructing environmental maps for arbitrary continuous, irregular, rough diffuse reflection scenes, in order to solve the problem of perceiving arbitrary irregular scattering surfaces in continuous diffuse reflection scenes.

[0007] The communication-driven environmental map construction method for arbitrary continuous, irregular, and rough diffuse reflection scenarios provided by this invention models environmental scattering within wireless communication scenarios as a more realistic continuous and rough diffuse reflection form. Then, an improved off-network sparse Bayesian method is used to extract rich channel and environmental parameters from highly correlated and complex (introduced by the complex environment) channel state information. Finally, massive discrete data from multi-user, multi-base station observations are fused to create a continuous physical environment map. The specific steps are as follows:

[0008] Step (1) Model the received signal at the receiver using a continuous coarse diffuse reflection probability model;

[0009] Step (2) uses a two-dimensional off-network sparse Bayesian algorithm to extract channel parameter information such as angle from the received signal;

[0010] Step (3) Use the parameter information extracted from different observations to perceive the state of the scattering point;

[0011] Step (4) Integrate the sensing scattering point information of multiple users and multiple base stations to achieve data cleaning and clustering;

[0012] Step (5) uses the rotational polynomial fitting method to construct a high-precision physical map.

[0013] Furthermore:

[0014] The specific design for the received signal modeling described in step (1) is as follows:

[0015] Consider an uplink communication scenario with M t M user equipments with M antennas and M r For uplink communication at K base stations with one antenna, consider the transmission of T Orthogonal Frequency Division Multiplexing (OFDM) pilot symbols. Each symbol consists of M subcarriers, orthogonally allocated to M users. Assuming the channel remains unchanged within a frame (within the T pilot symbols), the received signal of the m-th subcarrier at the k-th base station and the channel model form are confirmed as follows;

[0016] y m,k =h m,k x m,k +ω m,k (1)

[0017] in, The channel matrix; The pilot sequence is known at the base station. Each element has a mean of 0 and σ. 2 Additive Gaussian white noise with variance;

[0018] The signal propagates to the receiver via a line-of-sight (LoS) path and a non-line-of-sight (NLoS) path through diffuse reflection from N well-separated rough scattering surfaces of unknown location and shape. Estimating the complex rough scattering components makes it possible to create high-precision maps. The LoS path is denoted by n = 0, l = 0, and the l-th NLoS path from the scattering cluster corresponding to the n-th rough surface is denoted by n = 1, 2, ..., N, l = 1, 2, ..., L. m,k,n It means that γ m,k,n,l θ m,k,n,l , Let A, B, C, and D represent the equivalent complex gain, angle of arrival (AOA), and angle of departure (AOD) of the corresponding paths, respectively. The scattering path cluster is characterized using a multipath model.

[0019]

[0020] Where: considering a uniform linear array (ULA) antenna with half-wavelength spacing, the receive array response vector a in the equation is... r (θ m,k,n,l )and for:

[0021]

[0022]

[0023] Among them, the equivalent complex gain γ m,k,n,l :

[0024]

[0025] It is the spacing between adjacent subcarriers, B is the total bandwidth, and τ is the distance between adjacent subcarriers. m,k,n,l To account for the time delay of the corresponding propagation path, It is a random phase with complex gain phase on [0, 2π].

[0026] In step (2), a two-dimensional off-network sparse Bayesian algorithm is used to extract channel parameter information such as angle from the received signal; for the signal y sent by the m-th user received at the k-th base station... m,k The complex angle parameters of diffuse reflection scenes are extracted using an off-grid sparse Bayesian method. The specific design is as follows:

[0027] (2.1) The Least Squares (LS) method is used for channel estimation, where the channel estimation for the uplink communication channel between the m-th user and the k-th base station is:

[0028]

[0029] Where, x m,k H Guide frequency matrix x m,k The conjugate transpose of;

[0030] (2.2) Preliminary determination of the AOD and AOA distributions of several scattering clusters composed of the LoS radius angle and the NLoS radius:

[0031] The angular spectrum of the channel is obtained using a two-dimensional fast Fourier transform (2D-FFT), which will help in understanding the angular spread of the channel. To improve the resolution of the subsequent 2D-FFT, the matrix... Perform zero-padding to expand its dimension to M. f ×M f matrix M f Satisfy M f >>M r M f >>M t The specific operation is as follows:

[0032]

[0033] That is: to use a matrix with smaller dimensions All elements are filled into a larger-dimensional matrix in their original order. In the middle, other elements are set to 0.

[0034] For matrix Perform 2D-FFT to obtain This operation is recorded as:

[0035]

[0036] Where FFT2{·} denotes performing a 2D-FFT transformation on the matrix.

[0037] To extract the region of concentrated large values ​​in the matrix (the clustered distribution region of NLoS paths), it is necessary to process the resulting matrix. Real-valued representation. Two-dimensional matrix. The modulo operation on each element can be represented in the following mathematical form:

[0038]

[0039] in, For the resulting matrix after modulo operation, |·| represents the modulo operation. Representation matrix The (i,j)th element will represent the complex gain of the channel at the i-th arrival angle and the j-th departure angle, and its corresponding discrete angle is:

[0040]

[0041] In the formula, i and j are the indices of the corresponding elements in the two-dimensional FFT result, and they represent the normalized frequencies.

[0042] (2.3) Sliding convolution selects the region of concentrated intensity values ​​in the 2D-FFT result;

[0043] Here, a square matrix of all 1s is chosen as the convolution kernel, denoted as . Where M g <<M f For two-dimensional matrices To perform a convolution operation, the convolution kernel slides across the matrix, multiplying each element value of the matrix with the corresponding value on the convolution kernel. The sum of the multiplied values ​​is then used as a representation of the sum of the intensities of AOD and AOA in that local region.

[0044] According to Formula 7, calculate the angle value corresponding to the boundary value of the high-intensity region, which will be used to establish a dictionary matrix for fine search of the range of the corresponding angle in order to achieve high-precision joint estimation of AOD and AOA with low complexity.

[0045] (2.4) Consider the Loss path component and the Nolos path component in the separated channel matrix. The components are respectively denoted as and The specific implementation is as follows:

[0046] make If some elements are 0, the matrix is ​​denoted as . Since the energy of the coarse scattering NLoS path is significantly lower than that of the LoS path, the angle distributions of DOD and DOA corresponding to the LoS and NLoS paths differ greatly. Furthermore, the padding of zeros during the 2D-FFT operation further increases the distance between paths corresponding to different angles on the matrix. Therefore, this operation has little impact on the subsequent algorithm's search for accurate and reliable NLoS path corresponding angles.

[0047] (2.5) Transform the FFT result inversely into the channel matrix form;

[0048] H m,k NLoS After performing a two-dimensional inverse fast Fourier transform (2D-IFFT) and removing the preceding zero-interpolated components, the matrix is ​​obtained. Defined as an NLoS multi-cluster channel matrix;

[0049] (2.6) The sparse Bayesian search method is used to estimate AOD and AOA;

[0050] In existing joint AOD and AOA estimation methods, traditional subspace methods are widely used, mainly including multi-signal classification (MUSIC) and rotation-invariant (ESPRIT) algorithms for estimating signal parameters. These algorithms rely on a large number of observations to capture the signal or noise subspace, so their performance will be greatly reduced when the number of observations is insufficient or the signal-to-noise ratio is low. Especially in the coarse continuous diffuse reflection scenario of this invention, the NLoS path components in the signal are highly correlated, and the performance of such methods will degrade significantly.

[0051] Sparse Bayesian Learning (SBL) is an emerging parameter estimation method based on the CS (Sparse Bayesian) concept. It exhibits excellent estimation performance under conditions of low signal-to-noise ratio and limited observation count. In cases where NLoS path components are highly correlated, off-grid SBL can overcome the parameter-grid mismatch problem and achieve joint estimation of AOD (Aspect-Oriented Distance) and AOA (Aspect-Oriented Distance). Specifically:

[0052] The joint estimation problem of AOD and AOA is modeled as an SBL problem:

[0053] The input is the channel matrix h from which parameters need to be extracted. Let h = h(:), that is, let its columns form a new column vector in order.

[0054] Establish a dictionary matrix: Select M in sequence and uniformly within the AOA and AOD angle range searched in step (3). d A point, generally M d >>M r M t ;

[0055] Establish the following dictionary matrix along the Q diagonals of the rectangular region corresponding to the angular range:

[0056]

[0057]

[0058] in, Refers to the (i, j)th point in the two-dimensional grid;

[0059] If the angles to be estimated are all located on grid points:

[0060]

[0061] The coefficients corresponding to each column of the dictionary matrix, for Then there is (Set to a large number; unreliable estimates can be processed by subsequent noise reduction) All elements are not 0, and E is zero-mean Gaussian noise with unknown noise power.

[0062] Since the angle to be estimated generally does not fall exactly on a grid point, a linear approximation off-grid model is used to handle the mismatch problem:

[0063] If the ξ-th element to be estimated The element that is not in the dictionary matrix and whose nearest element in the dictionary matrix is: Then it can be approximated as:

[0064]

[0065] in:

[0066]

[0067] definition:

[0068]

[0069]

[0070]

[0071]

[0072]

[0073]

[0074] The off-grid model can then be established as follows:

[0075] h=G(ε,ρ)X+E, (14)

[0076]

[0077] diag{·} represents the diagonal matrix operator, and ε and ρ represent the difference between the true angle and the nearest grid cell:

[0078]

[0079] The off-network sparse Bayesian model is established as follows:

[0080]

[0081] Noise accuracy: Prior distribution:

[0082]

[0083] in: Let X denote a complex Gaussian distribution, Γ(·) denote a gamma distribution, and a and b be parameters controlling the distribution. A typical Gaussian prior distribution is assigned to each row of X, with different precision τ. i , I is the identity matrix.

[0084] The optimization estimation problem of the relevant parameters is then determined by the following equation:

[0085]

[0086] in, τ * ,ε * ,ρ * This represents the optimal solution of the parameters under criterion (19).

[0087] The above problem is iteratively solved using the Block Majorization-Minimization algorithm.

[0088] ε is obtained in j+1 iterations. j+1 ,ρ j+1 , And update the grid points:

[0089]

[0090] After each iteration Compare Distance from reality Closer;

[0091] After iteration, based on the constraint that there is a certain distance between the angles, potentially duplicated angle combinations are removed, resulting in the angles that need to be estimated. Group angle parameters.

[0092] The specific design for sensing the scattering point state using parameter information extracted from different observations in step (3) is as follows:

[0093] State perception is performed based on the estimated AOD and AOA parameter information from step 2, and the result is recorded as follows: It is the sum of the observation angles obtained based on observations from all users and base stations. Based on The AOA and AOD of the observed l-th NLoS path are resolved in this step. Estimate the state of the scattering point. This refers to the location of the m-th user device. Let k be the location of the kth base station. This refers to the l-th path (l = 1, 2, ..., L) estimated based on the communication channel between the m-th user and the k-th base station. m,kThe coordinates of the scattering point after scattering by the scatterer are obtained, and the specific operation is as follows:

[0094] The location of the scattering point on the scatterer surface sensed during a single frame of communication is determined based on the estimated base station location and user location. Scattering point The angle between the connection to the base station and the transmission of ULA is The angle between the connection to the user equipment and the receiving ULA is...

[0095] The problem of estimating the location of a single scattering point can be characterized as the intersection of two straight lines. Specifically, the l-th scattering point of the uplink communication between the m-th user and the k-th base station is represented by the straight line. With a straight line The intersection, Let f(·) be a real number variable, which refers to the direction vector function determined by the corresponding AoA and AoD.

[0096] The solution can then be obtained based on the estimated number of angles. Location of each scattering point;

[0097] The least squares solution to this problem is:

[0098]

[0099] The process of fusing sensing scattering point information from multiple users and multiple base stations to achieve data cleaning and clustering in step (4) is specifically designed as follows:

[0100] The estimated scattering point coordinates are clustered using the high-density connected component clustering algorithm DBSCAN, jointly achieving noise reduction and scatterer identification. The estimated scattering point coordinates are then transmitted from the base station to the central processing terminal, and the scattering points are recorded in the order of transmission. use This represents the set of coordinates of all observed scattering points;

[0101] In the DBSCAN algorithm, all point coordinates are either clustered into different clusters (each cluster corresponds to a scatterer) or identified as noise; within a specific cluster C, the two core parameters of the algorithm are the neighborhood. Minimum number of points in the neighborhood

[0102] Define the distance to be measured:

[0103]

[0104] The ε-neighborhood of the p-th scattering point is defined as:

[0105]

[0106] when and At that time, it was believed that from arrive It is density-reachable if an ordered sequence exists. and arrive If the density is achievable, then arrive It has an achievable density. If The number of scattering points contained in the neighborhood of the center is not less than N. min ,Right now Then it can be called This is the core point for clustering.

[0107] Algorithm implementation steps:

[0108] ① Given the total number of estimated scattering points, N min ε;

[0109] ② Select any point from the estimated scattering points that satisfies the core point condition as the initial core point;

[0110] ③ Search for all estimated scattering points that are reachable from the initial core point density. If the intersection of adjacent clusters satisfies the clustering condition, merge them and calculate the total number of clusters Z.

[0111] ④ Output: Set The elements in the data are divided into scattering point clusters S. z (z = 1, 2, Z) and noise set S noise The noise set mainly consists of scattering points;

[0112] The high-precision physical map constructed in step (5) using the rotated polynomial fitting method is specifically designed as follows:

[0113] Cluster z S z The included scattering points can be considered as continuous lines on the z-th scattering surface, and these points will originate from S. z The function expression of the line is obtained by training the system using the coordinates of all scattering points as the training set.

[0114] S z The number of scattering points in the middle is counted as The horizontal and vertical coordinates are denoted as:

[0115]

[0116] The output expression for polynomial fitting is:

[0117]

[0118] Where n is the degree of polynomial fitting, a i (i = 1, 2, ..., n) are the real-valued optimization variables to be estimated.

[0119] However, in actual map construction, since polynomial fitting mainly relies on the loss along the y-axis for optimization, and the position and orientation of objects in real space are random, the way the coordinate axes are set will have a significant impact on the results. Therefore, this method uses a dynamically rotated coordinate system to obtain the optimal polynomial fitting method (hereinafter referred to as the rotated polynomial method), as follows:

[0120] For a fixed coordinate axis, rotating the coordinate axis counterclockwise by ψ results in the scattering point rotating clockwise by ψ (ψ∈

[0121] [0,π]),

[0122] A polynomial fitting optimization equation system is established based on a training dataset with a specific angle ψ.

[0123] The system of equations can then be written as: Ab = f. The problem is then transformed into estimating vector b, specifically the following optimization problem:

[0124]

[0125]

[0126] Because there is The problem has an analytical solution:

[0127]

[0128] The preliminary result of the functional relationship of the z-th scattering surface is obtained from:

[0129] y = g(x) = a0 * +a1 * x+a2 * x 2 +a n * x n (28)

[0130] The loss can be calculated as follows:

[0131]

[0132] Optimization results output: Mesh optimization parameter tuning determines the optimal rotation angle ψ * ;

[0133] Q1 grid points are uniformly selected within the interval [0,π] to form a parameter set. Then we have:

[0134]

[0135] Inverse transform output of the fitting result:

[0136] Rotate ψ clockwise around the coordinate axis * Then rotate the fitted curve counterclockwise by ψ * The output is the scattering surface reconstruction result for all datasets S. z The fitting process is performed sequentially (z = 1, 2, Z) to obtain the final physical map result.

[0137] This invention assumes continuous diffuse reflection of the scatterer and extracts complex and highly correlated channel parameters based on a suitable algorithm. It fuses observation results from multiple users and multiple base stations to achieve perception of unknown environments of arbitrary shapes, including regular smooth scattering surfaces and irregular scattering surfaces, and finally constructs a continuous physical map of the environment.

[0138] Compared to existing sensor solutions, this invention can significantly reduce the cost of environmental map construction and prevent privacy leaks; compared to existing communication-driven environmental map construction solutions, it is applicable to physical scattering scenarios of any unknown shape and can effectively extract and fuse observations of environmental parameter states on a higher order of magnitude. Attached Figure Description

[0139] Figure 1 This is a flowchart of the method of the present invention.

[0140] Figure 2 This is a typical application scenario and modeling diagram of the present invention.

[0141] Figure 3 The figure shows the experimental results of this invention. Detailed Implementation

[0142] A communication-driven map construction method for arbitrary continuous, irregular, coarse diffuse reflection scenes, with the following specific steps:

[0143] Consider an uplink communication scenario with M t M user equipments with M antennas and M r For uplink communication of K base stations with one antenna, consider the transmission of T orthogonal frequency division multiplexing (OFDM) pilot symbols. Each symbol consists of M subcarriers, orthogonally allocated to M users. Assuming the channel remains unchanged within a frame (within the T pilot symbols), the received signal and channel model form are determined as follows:

[0144] The k-th base station receives the signal sent by the m-th user:

[0145] y m,k =h m,k x m,k +ω m,k

[0146] in, The channel matrix; The pilot sequence is known at the base station. Each element has a mean of 0 and σ. 2 Additive Gaussian white noise with variance;

[0147] The signal propagates to the receiver via a line-of-sight path (LoS) and diffuse reflection non-line-of-sight paths (NLoS) from N well-separated rough scattering surfaces of unknown location and shape. Estimating the complex rough scattering components makes it possible to create high-precision maps. The LoS path is denoted by n = 0, l = 0, and the l-th NLoS path from the scattering cluster corresponding to the n-th rough surface is denoted by n = 1, 2, ..., N, l = 1, 2, ..., L. m,k,n It means that γ m,k,n,l θ m,k,n,l , Let A, B, C, and D represent the equivalent complex gain, angle of arrival (AOA), and angle of departure (AOD) of the corresponding paths, respectively. The scattering path cluster is characterized using a multipath model.

[0148]

[0149] Where: considering a uniform linear array (ULA) antenna with half-wavelength spacing, the receive array response vector a in the equation is... r (θ m,k,n,l )and for:

[0150]

[0151]

[0152] Where the equivalent complex gain γ m,k,n,l :

[0153]

[0154] It is the spacing between adjacent subcarriers, B is the total bandwidth, and τ is the distance between adjacent subcarriers. m,k,n,l To account for the time delay of the corresponding propagation path, It is a random phase with complex gain phase on [0, 2π].

[0155] Complex gain The determination was based on existing literature (Kulmer J, Wen F, Garcia N, et al. Impact of Rough Surface Scattering on Stochastic Multipath Component Models[C] / / 2018IEEE 29th Annual International Symposium on Personal, Indoor and Mobile Radio Communications (PIMRC).IEEE,2018.DOI:10.1109 / PIMRC.2018.8580964.).

[0156] Step (2): The sparse Bayesian algorithm extracts channel parameter information such as angle from the received signal. The specific method is as follows:

[0157] First, the least squares (LS) method is used for channel estimation. The channel estimation for the uplink communication channel between the m-th user and the k-th base station is as follows:

[0158]

[0159] Where, x m,k H Guide frequency matrix x m,k The conjugate transpose of .

[0160] (2.1) The distribution of AOD and AOA of several scattering clusters composed of the LoS path angle and NLoS path was initially determined. The angular spectrum of the channel was obtained using a two-dimensional fast Fourier transform (2D-FFT), which will help in understanding the angular spread of the channel. To improve the resolution of the subsequent 2D-FFT, the matrix... Perform zero-padding to expand its dimension to M. f ×M f matrix M f Satisfy M f >>M r M f >>M t The specific operation is as follows:

[0161]

[0162] That is: to use a matrix with smaller dimensions All elements are filled into a larger-dimensional matrix in their original order. In the middle, other elements are set to 0.

[0163] For matrix Perform 2D-FFT to obtain This operation is recorded as:

[0164]

[0165] Where FFT2{·} denotes performing a 2D-FFT transformation on the matrix.

[0166] To extract the region of concentrated large values ​​in the matrix (the clustered distribution region of NLoS paths), it is necessary to process the resulting matrix. Real-valued representation. Two-dimensional matrix. The modulo operation on each element can be represented in the following mathematical form:

[0167]

[0168] in, For the resulting matrix after modulo operation, |·| represents the modulo operation. Representation matrix The (i,j)th element will represent the complex gain of the channel at the i-th arrival angle and the j-th departure angle, and its corresponding discrete angle is:

[0169]

[0170] Where: i and j are the indices of the corresponding elements in the two-dimensional FFT result, and they represent the normalized frequencies.

[0171] (2.2) Sliding convolution selects the region of concentrated intensity values ​​in the 2D-FFT result;

[0172] Here, a square matrix of all 1s is chosen as the convolution kernel, denoted as . Where M g <<M f For two-dimensional matrices The convolution operation involves sliding the convolution kernel across the matrix, multiplying each element value of the matrix by the corresponding value on the convolution kernel, and then summing the multiplied values ​​to represent the sum of the intensity of AOD and AOA in that local region.

[0173] The region with the highest intensity value obtained by solving with a convolution kernel of small matrix dimension is the distribution region of the LoS radial component. Extract the relevant coordinate range information of the region and rewrite the values ​​within the coordinates to 0 to reduce their impact on the estimation of the corresponding parameters of the NLOS path. Then, use the same method to preliminarily search for the coordinate information of the clusters formed by the remaining NLOS scattering paths.

[0174] (2.3) Based on the method in step (2), calculate the angle value corresponding to the boundary value of the above region, and use it as the basis for subsequent establishment of dictionary matrix to finely search the range of the corresponding angle in order to achieve high-precision joint estimation of AOD and AOA with low complexity.

[0175] (2.4) Consider the Loss path component and the Nolos path component in the separated channel matrix. The components are respectively denoted as and The specific implementation is as follows:

[0176] make If some elements are 0, the matrix is ​​denoted as . Because the energy of the coarse scattering NLoS path is significantly lower than that of the LoS path, the angle distributions of the corresponding DOD and DOA of the LoS and NLoS paths differ greatly. Furthermore, the padding with zeros during the 2D-FFT operation further increases the distance between paths corresponding to different angles on the matrix. Therefore, this operation has little impact on the subsequent algorithm's search for accurate and reliable NLoS path corresponding angles.

[0177] (2.5) Transform the FFT result inversely into the channel matrix form;

[0178] H m,k NLoS After performing a two-dimensional inverse fast Fourier transform (2D-IFFT) and removing the preceding zero-interpolated components, the matrix is ​​obtained. Defined as an NLoS multi-cluster channel matrix.

[0179] (2.6) Sparse Bayesian search method for estimating AOD and AOA;

[0180] The input is the channel matrix h from which parameters need to be extracted. Let h = h(:), that is, let its columns form a new column vector in order.

[0181] Establish a dictionary matrix: Select M in sequence and uniformly within the AOA and AOD angle range searched in step (3). d A point, generally M d >>M r M t

[0182] To maintain the same precision as a 1D dictionary matrix in a 2D dictionary matrix, a much larger number of grid points are needed, along with strong correlation between dictionary elements under the condition of the same one-dimensional angle. Therefore, the following dictionary matrix is ​​constructed along the m diagonals of the rectangular region corresponding to the angle range:

[0183]

[0184] in, The i-th (i, j) point in the two-dimensional dictionary matrix grid; if the angles to be estimated are all located on grid points:

[0185] The coefficients corresponding to each column of the dictionary matrix, for h m,k NLoS Then there is (Set to a large number; unreliable estimates can be processed by subsequent noise reduction) All elements are not 0, and E is zero-mean Gaussian noise with unknown noise power.

[0186] Since the angle to be estimated generally does not fall exactly on a grid point, a linear approximation off-grid model is used to handle the mismatch problem:

[0187] If the ξ-th element to be estimated The element that is not in the dictionary matrix and whose nearest element in the dictionary matrix is: Then it can be approximated as:

[0188]

[0189] in:

[0190]

[0191] definition:

[0192]

[0193]

[0194]

[0195]

[0196]

[0197]

[0198] The off-grid model can then be established as follows:

[0199] h=G(ε,ρ)X+E

[0200]

[0201] diag{·} represents the diagonal matrix operator, and ε and ρ represent the difference between the true angle and the nearest grid cell:

[0202]

[0203] The off-network sparse Bayesian model is established as follows:

[0204]

[0205] Noise accuracy: Prior distribution:

[0206]

[0207] in: Let X denote a complex Gaussian distribution, Γ(·) denote a gamma distribution, and a and b be parameters controlling the distribution. A typical Gaussian prior distribution is assigned to each row of X, with different precision τ. i , I is the identity matrix.

[0208] The optimization estimation problem of the relevant parameters is then determined by the following equation:

[0209]

[0210] τ * ,ε * ,ρ * This represents the optimal solution for the parameters under this criterion.

[0211] The above problem is iteratively solved using the Block Majorization-Minimization algorithm.

[0212] The (j+1)th iteration yields ε j+1 ,ρ j+1 , And update the grid points:

[0213]

[0214] After each iteration Compare Distance from reality More recently, after the iteration is complete, based on the constraint that there should be a certain distance between the angles, the angle combinations that may be repeatedly estimated are deleted to obtain the required estimated angles. Group parameters.

[0215] Step (3): Based on the estimated AOD and AOA parameter information in step (2), perform state perception and record... It is the sum of the observation angles obtained based on observations from all users and base stations. Based on The AOA and AOD of the observed l-th NLoS path are resolved in this step. Estimate the state of the scattering point. This refers to the location of the m-th user device. Let k be the location of the kth base station. This refers to the l-th path (l = 1, 2, ..., L) estimated based on the communication channel between the m-th user and the k-th base station. m,k The coordinates of the scattering point after scattering by the scatterer are obtained, and the specific operation is as follows:

[0216] The location of the scattering point on the scatterer surface sensed during a single frame of communication is determined based on the estimated base station location and user location. Scattering point The angle between the connection to the base station and the transmission of ULA is The angle between the connection to the user equipment and the receiving ULA is...

[0217] The problem of estimating the location of a single scattering point can be characterized as the intersection of two straight lines. Specifically, the l-th scattering point of the uplink communication between the m-th user and the k-th base station is represented by the straight line. With a straight line The intersection, Let f(·) be a real number variable, which refers to the direction vector function determined by the corresponding AoA and AoD.

[0218] The solution can then be obtained based on the estimated number of angles. Location of each scattering point;

[0219] The least squares solution to this problem:

[0220]

[0221] Step (4) Construct an environment map based on the results of the state-aware module. →Environment: S s ={S n}

[0222] 4.1 The estimated scattering point coordinates are clustered based on the high-density connected region clustering algorithm DBSCAN, which jointly realizes noise reduction and scattering object identification.

[0223] The estimated scattering point coordinates in 4.2(4) are sent to the central processing terminal via the base station, and the scattering points are recorded in the order of transmission. use This represents the set of coordinates of all observed scattering points;

[0224] 4.3 In the DBSCAN algorithm, all point coordinates are either clustered into different clusters (each cluster corresponds to a scatterer) or identified as noise. Within a specific cluster C, the two core parameters of this algorithm are the neighborhood. Minimum number of points in the neighborhood

[0225] Define the distance to be measured:

[0226]

[0227] The ε-neighborhood of the p-th scattering point is defined as:

[0228]

[0229] Unlike clustering algorithms such as K-means, DBSCAN's clustering feature is density-reachability.

[0230] when and At that time, we believed that from arrive It is density-reachable.

[0231] If an ordered sequence exists and arrive If the density is achievable, then arrive It has a achievable density;

[0232] like The number of scattering points contained in the neighborhood of the center is not less than N. min ,Right now Then it is called This is the core point for clustering.

[0233] 4.4 Algorithm Implementation Steps

[0234] 4.4.1 Given the total number of estimated scattering points N min ε;

[0235] 4.4.2 Select any point from the estimated scattering points that satisfies the core point condition as the initial core point;

[0236] 4.4.3 Search for all estimated scattering points that are reachable from the initial core point density. If the intersection of adjacent clusters satisfies the clustering condition, merge them and calculate the total number of clusters Z.

[0237] 4.4.4 Output: Set The elements in the data are divided into scattering point clusters S. z (z = 1, 2, Z) and noise set S noise The noise set mainly consists of scattering points.

[0238] Step (5) Construct a high-precision physical map based on the rotation polynomial fitting method.

[0239] Cluster z S z The included scattering points can be considered as continuous lines on the z-th scattering surface, and these points will come from S. z The function expression of the line is obtained by training the system using the coordinates of all scattering points as the training set.

[0240] 5.1Sz The number of scattering points in the middle is counted as The horizontal and vertical coordinates are denoted as:

[0241]

[0242] 5.2 The relevant data were processed using a rotated polynomial fitting method;

[0243] The output expression for polynomial fitting is:

[0244]

[0245] Where: n in this expression is the polynomial fitting degree a i (i = 1, 2, ..., n) are the real-valued optimization variables to be estimated.

[0246] However, in actual map construction, since polynomial fitting mainly relies on the loss along the y-axis for optimization, and the position and orientation of objects in actual space are random, the way the coordinate axes are set will have a significant impact on the results. Therefore, this method uses a dynamically rotated coordinate system to obtain the optimal polynomial fitting method (hereinafter referred to as the rotated polynomial method), as follows:

[0247] For a fixed coordinate axis, rotating the coordinate axis counterclockwise by ψ results in the scattering point rotating clockwise by ψ (ψ∈

[0248] [0,π]), specifically represented as:

[0249] x ′ = x·cos(ψ) + f·sin(ψ)

[0250] y ′ = -x·sin(ψ) + f·cos(ψ)

[0251]

[0252] 5.3 Optimization Problem Construction: Based on the training dataset with a specific angle ψ, the following system of equations is established:

[0253]

[0254] make:

[0255]

[0256] The system of equations can then be written as: Ab = f. The problem is then transformed into estimating vector b, specifically the following optimization problem:

[0257]

[0258]

[0259] Because there is The problem has an analytical solution:

[0260]

[0261] The preliminary result of the functional relationship of the z-th scattering surface is obtained from:

[0262] y = g(x) = a0 * +a1 * x+a2 * x 2 +a n * x n

[0263] Calculate the loss:

[0264] 5.4 Optimization Result Output

[0265] 5.4.1 Mesh optimization and parameter tuning to determine the optimal rotation angle ψ * ;

[0266] exist Q grid points are uniformly selected within the interval to form the parameter set S. ψ ={ψ0,ψ1,ψ2…ψ Q-1}, then we have:

[0267]

[0268] 5.4.2 Inverse transformation output of the fitting results;

[0269] Rotate ψ clockwise around the coordinate axis * Then rotate the fitted curve counterclockwise by ψ * The output is the result of the scattering surface reconstruction, specifically:

[0270] Based on ψ * Calculate the corresponding b * (ψ * Preliminary results were obtained:

[0271] y = g(x) = a0 * +a1 * x+a2 * x 2 +a n * x n

[0272] Rotate the function curve counterclockwise by ψ * Specifically, this involves transforming each point (x0, y0) on the curve as follows to map the fitted curve to the original coordinate axes:

[0273] x = x0·cos(ψ) * )-y0·sin(ψ * )

[0274] y = x0·sin(ψ) * )+y0·cos(ψ * )

[0275] For all datasets S z The fitting process is performed sequentially (z = 1, 2, Z) to obtain the final physical map result.

[0276] Experiment: Consider multi-user access scenarios (e.g.) Figure 3 As shown in (a), the carrier frequency is considered to be 28 GHz, the pilot bandwidth is 25 MHz, the number of subcarriers M = 50, and the number of OFDM pilot symbols T = 400. There are two scatterers in the scene. One scatterer is represented by a straight line, and the surfaces of both scatterers are described as rough surfaces without specular components. Each surface is approximated by 60 scattering points. The number of users K = 50, the number of base stations M = 4, and both users and base stations are set as uniform linear arrays (ULAs) with 16 antennas. The antenna element spacing is set to half a wavelength, and the transmission signal-to-noise ratio is set to 10 dB. The environmental reconstruction performance of the method of this invention and the existing method (CP decomposition-ESPRIT parameter estimation) in this scene is compared in the experiment. The results are as follows: Figure 3 As shown in (b)(c)(d).

Claims

1. A method for constructing an environmental map for an arbitrary continuous irregular rough diffuse reflection scene, characterized by, The environment scattering in a wireless communication scenario is modeled as a more realistic continuous rough diffuse reflection form, then rich channel and environment parameters are extracted from the strong correlation and complex component channel state information using an improved off-grid sparse Bayesian method, and finally a continuous physical environment map is drawn by fusing massive discrete data observed by multiple users and multiple base stations; The specific steps are as follows: Step (1) modeling the received signal at the receiving end with a continuous rough diffuse reflection probability model; Step (2) using a two-dimensional off-grid sparse Bayesian algorithm to extract channel parameter information from the received signal; Step (3) using parameter information extracted by different observations to perceive the state of scattering points; Step (4) fusing the perceived scattering point information of multiple users and multiple base stations to realize data cleaning and clustering; Step (5) using a rotating polynomial fitting method to construct a high-precision physical map; Wherein: In step (4), the perceived scattering point information of multiple users and multiple base stations is fused to realize data cleaning and clustering, as follows: High-density connected region clustering algorithm The estimated scattering point coordinates are clustered to jointly achieve noise reduction and scatterer identification. The estimated scattering point coordinates are sent to a central processing terminal by a base station. The scattering points are denoted as The set of all observed scattering point coordinates is denoted as ​ In high-density connected region clustering algorithms, all point coordinates are clustered into different clusters or identified as noise; in a particular cluster , two core parameters of the algorithm are neighborhood and minimum number of points within the neighborhood ; Defining the measured distance: , (22) Then the neighborhood of the scatter point is defined as: ， （23） when and At that time, it was believed that from arrive It is density-reachable if an ordered sequence exists. ,and arrive If the density is achievable, then arrive It is achievable in terms of density; if The number of scattering points contained in the neighborhood of the center is not less than ,Right now Then it is called As the core point of clustering; The specific steps of the algorithm are as follows:

1. The total number of estimated scattering points is given by reference , ; ② Select any point that meets the core point condition from the estimated scattering points as the initial core point; ③Search all estimated scattering points which satisfy the initial core point density reachable, if the intersection of adjacent clusters satisfies the clustering condition, merge, calculate the total number of clusters ; ④ Output: Set The elements in the data are divided into scattering point clusters S. z and noise set S z In the range z=1, 2, Z, the noise set mainly consists of scattering points; In step (5), the high-precision physical map is constructed by using a rotating polynomial fitting method, as follows: Cluster z S z The included scattering points can be considered as the first On a scattering surface, these points are considered as continuous lines in two-dimensional space, and these points will come from S. z The coordinates of all scattering points are used as the training set to train and obtain the functional relationship expression of the line; S z The number of scattering points in the image is counted as The horizontal and vertical coordinates are denoted as: ， （24） The output expression of polynomial fitting is: ， （25） wherein is the polynomial fit degree, is a real optimization variable to be estimated, ; A dynamic rotating coordinate system is used to obtain the optimal polynomial fitting method, as follows: For fixed axes, rotate the rotating axes counterclockwise then the scatter point rotates clockwise counterclockwise ; Based on a specific angle of the training dataset to establish a polynomial fitting optimization equation set; The system of equations is written as: The problem then reduces to estimating the vector which is specifically the following optimization problem: ， （26） ， Because of , the problem has an analytical solution: = ， （27） Then the function relationship of the first scatter surface is the preliminary result: ， （28） The loss is calculated as: = ； （29） Optimization result output: grid optimization parameter tuning determines the best rotation angle ; In Uniformly select grid points in the interval, to form a parameter set , then: ， （30） The fitting result is output by inverse transformation: Rotate the coordinate axes clockwise Rotate the fitted curve counterclockwise Output is the scatter plane recovery result, for all data sets S z Sequentially perform fitting to obtain the final physical map result.

2. The environmental map construction method according to claim 1, characterized by, In step (1), the received signal is modeled as follows: Consider an uplink communication scenario with M user equipments with M root antennas communicating with M base stations with M orthogonal frequency division multiplexing pilot symbols, each symbol consisting of M subcarriers, orthogonally assigned to M users; assume that the channel remains constant over a frame, the Mth subcarrier is received at the Mth base station, the received signal and channel model form ， （1） wherein, is a channel matrix; is a pilot sequence known by the base station end; is an additive white Gaussian noise with mean and variance and variance signals propagate to the receiver through line-of-sight (LoS) paths and non-line-of-sight (NLoS) paths with rough scattering surfaces, which are well separated from each other and whose shape is unknown. The estimation of the complex scattering components also brings the possibility to depict high-precision maps. The LoS path is denoted as , the first NLoS path from the scattering cluster corresponding to the first rough surface is denoted as , and the th NLoS path from the scattering cluster is denoted as , , where denotes the number of diffuse multipaths scattered by the nth rough surface, and denotes the equivalent complex gain, the angle of arrival, and the angle of departure of the corresponding path, respectively. The scattering path cluster is characterized by a multipath model. ； （2） where: considering the uniform linear array antenna with half-wavelength interval, the received array response vector in the formula is and is: ， ， ； （3） ， ， ； （4） wherein the equivalent complex gain : is the spacing of adjacent subcarriers, B is the total bandwidth, is the delay of the corresponding propagation path, is the real complex gain, the complex gain phase is a random phase on 3. The environmental map construction method according to claim 2, characterized by, In step (2), a two-dimensional off-network sparse Bayesian algorithm is used to extract channel parameter information from the received signal; for the first... The first base station received the first Signal sent by a user The complex angle parameters of the diffuse scene are extracted using the off-network sparse Bayesian method, as follows: (2.1) Channel estimation is performed using the least squares method. The user and the first Channel estimation of uplink communication channels for each base station for: ， （5） wherein guide frequency matrix the conjugate transpose of (2.2) preliminary determination angle and distribution of the angles of departure, angles of arrival of the several scattering clusters of the composition: The two-dimensional fast Fourier transform is used to obtain an angle spectrum of the channel, which will help to understand the angle dispersion of the channel; in order to improve the resolution of the subsequent two-dimensional fast Fourier transform, zero padding is performed on the matrix to expand its dimension to , , , , , and the specific operation is represented as: ， （6） That is: fill all elements of a matrix of smaller dimension into a matrix of larger dimension in the original order with other elements set to 0. on the matrix a two-dimensional fast Fourier transform is performed, resulting in This operation is denoted by ， wherein denotes a two-dimensional fast Fourier transform on a matrix; In order to extract the region of the larger value in the matrix, i.e. the region of the cluster of the non-line-of-sight path, the result matrix Real number representation; two-dimensional matrix The modulo operation on each element of the real number representation; two-dimensional matrix is expressed in the following mathematical form: wherein, is the modulo result matrix, represents the modulo operation, denotes the matrix , the element of the matrix , the complex gain of the representative channel at the th angle of arrival and the th angle of departure position, the corresponding discrete angle is: ， （7） Where i and j are the indices of the corresponding elements in the two-dimensional fast Fourier transform result, which represent the normalized frequency; (2.3) Select the high-intensity value area of the two-dimensional fast Fourier transform result by sliding convolution; Here, we select the full. A square matrix is ​​used as the convolution kernel, denoted as . ,in For two-dimensional matrices To perform a convolution operation, the convolution kernel slides across the matrix, multiplying each element value of the matrix with the corresponding value on the convolution kernel. The sum of the multiplied values ​​is then used as a representation of the sum of the intensity of the departure angle and arrival angle of the local region. According to formula 7, the angle value corresponding to the high-intensity area boundary value is calculated as the range of the corresponding angle for subsequent dictionary matrix fine search to realize high-precision departure angle and arrival angle joint estimation with low complexity; (2.4) the separation channel matrix is part of the composition of part in the components, respectively ; (2.5) Inverse transform the fast Fourier transform result into a channel matrix form; Will After performing a two-dimensional inverse fast Fourier transform and removing the preceding zero-interpolated components, the matrix is ​​obtained. Defined as Multi-cluster channel matrix; (2.6) Estimate the departure angle and arrival angle by using the sparse Bayesian search method.

4. The environmental map construction method according to claim 3, characterized by, In step (2.6), the departure angle and arrival angle are estimated by using the sparse Bayesian search method, as follows: The joint estimation problem of departure angle and arrival angle is modeled as a sparse Bayesian search problem: Input channel matrix from which parameters are to be extracted , let , i.e. let its columns form in order new column vectors; Establishing dictionary matrix: in the angle range of the arrival angle and the departure angle searched in step (3), select points in order and uniformly , ; The following dictionary matrix is established on the Q diagonal lines of the rectangular area corresponding to the angle range: （8） wherein , denotes the entry in the two-dimensional dictionary matrix grid at point; If the estimated angles are all on the grid points: ； （9） corresponding to each column of the dictionary matrix, for then elements are not zero, is a 0-mean Gaussian noise with unknown noise power. Since the estimated angles are not exactly on the grid points, a linear approximation off-grid model is used to handle the mismatch problem: If the first One element to be estimated The element that is not in the dictionary matrix and is the closest element in the dictionary matrix to that point is denoted as . Then it can be approximated as: ， （10） Where: ， （11） Definition: （12） ， ， ， （13） Then the off-grid model is established as: ， （14） ， （15） , With represents the difference between the real angle and the nearest grid: ， （16） The off-grid sparse Bayesian model is established as: ， （17） Noise precision: Prior distribution: ; ； （18） where: represents a complex Gaussian distribution, represents a Gamma distribution, is a parameter of the control distribution, assigning a typical Gaussian prior distribution to each row of X with different precision , , is the identity matrix; Then the optimization estimation problem of related parameters is determined by the following formula: ； （19） represents the optimal solution of the parameter under the (19) criterion; The iteration of the above problem is performed using the block max-min algorithm, the first step iteration gives and the grid points are updated: ， ， （20） after each iteration than closer to the real closer; After the iteration is completed, according to the constraint that there is a certain distance between the angles, the angle combinations that are possibly repeatedly estimated are deleted to obtain the angle parameters that need to be estimated group angles.

5. The environmental map construction method according to claim 4, characterized by, In step (3), the state of the scattering points is perceived by using the parameter information extracted by different observations, as follows: Based on the estimated departure angle and arrival angle parameters from step 2, state perception is performed, and the following is recorded: It is the sum of the observation angles obtained based on observations from all users and base stations. based on The first observation The arrival and departure angles of a non-line-of-sight path are determined in this step. Estimate the state of the scattering point. Refers to the first The location of each user device. For the first The location of each base station Refers to the first The user and the first The estimated number of base station communication channels in the nth base station communication channel The coordinates of the scattering point of the stripe diameter after being scattered by the scatterer. The specific steps are as follows: Solving perceived scatterer surface scatterer point locations from base station location and user location estimates for single frame communication scatterer point Angle of base station connection and transmit uniform linear array is Angle of user equipment connection and receive uniform linear array is ; The problem of estimating the location of a single scattering point can be characterized as the intersection of two straight lines, specifically the first line. The user and the first The first base station uplink communication The scattering points are a straight line. , , For real number variables, Refers to the direction vector function determined by the corresponding arrival and departure angles; Then the estimated number of angles can be solved based on scatter point positions; The least squares solution of this problem is: ，（21）。

Images