A Channel Estimation and User Localization Method Based on Spherical Array Smart Metasurface

Abstract

This invention relates to a channel estimation and user localization method based on a spherical array intelligent metasurface. Specifically, it includes: designing a spherical array intelligent metasurface (RIS) by arranging reconfigurable reflective elements in an ordered manner according to a three-dimensional spherical topology; constructing a beam training signal model based on standard multivariate decomposition according to the multipath propagation characteristics of cascaded channels, and recovering the path angle parameters on the base station side using a tensor signal processing algorithm; designing a RIS training reflection coefficient beamforming method based on the spherical phase mode excitation principle, and recovering the equivalent path angle parameters on the RIS side using a spatial spectrum estimation algorithm; and achieving unique decoupling of user path angle and time delay parameters based on the free-space propagation characteristics of broadband multicarrier signals, thus completing user localization. This invention can achieve unambiguous channel estimation, which is unsupported by traditional planar array RIS, with lower hardware costs and training overhead, thereby enabling high-precision user localization applications.

Description

Technical Field

[0001] This invention relates to the field of wireless communication and sensing technology, and in particular to a channel estimation and user localization method based on a spherical array intelligent metasurface. Background Technology

[0002] The technological scope of next-generation mobile communication (6G) will extend vertically from traditional wireless communication to information sensing and big data, encompassing information transmission, information collection, and information computing. Emerging mobile internet and vertical industry businesses centered on information sensing are gradually emerging, and the socio-economic informatization needs are driving a new trend in 6G technology development, focusing on Integrated Sensing and Communication (ISAC).

[0003] The ISAC system for 6G requires enhanced wireless air interface technology, driving the development of millimeter-wave / terahertz communication and multiple-input multiple-output (MIMO) technologies. However, the "scale for gain" technology evolution path faces inherent problems of increased cost, energy consumption, and complexity. To develop innovative, efficient, and green wireless solutions, the industry is seeking disruptive changes in information interaction models. Smart Metasurface (RIS) technology stands out with its unique advantages of low overhead, programmability, and ease of deployment, becoming one of the most promising key technologies for 6G integrated sensing and communication construction.

[0004] The physical essence of RIS (Intelligent Metasurface) is a reconfigurable metasurface array composed of regularly arranged subwavelength electromagnetic units. Through digital coding, it dynamically modulates electromagnetic waves to form an electromagnetic field with controllable amplitude, phase, polarization, and frequency. Unlike traditional wireless systems, the Intelligent Metasurface-based Sensing Integration (RIS-ISAC) system uses RIS configuration and control strategies as its core. By enriching scattering channels, superimposing in-phase signals, and controlling beam direction, it improves multiplexing gain, combats multipath fading, and enhances system coverage. RIS-ISAC technology holds the promise of overcoming the inherent constraint of uncontrollable channels, constructing a programmable intelligent wireless environment, and introducing a new paradigm of sensing integration for future 6G.

[0005] To fully realize the performance gain potential of RIS for ISAC wireless systems, accurate channel state information is essential. Since passive RIS units lack signal reception and processing capabilities, the pilot length of traditional channel estimation methods is proportional to the RIS array size, leading to excessive training overhead for large-scale MIMO systems. Furthermore, conventional planar RIS array structures lose three-dimensional physical space channel information, introducing inherent channel estimation ambiguity and failing to support high-precision user positioning, environmental mapping, and other information-aware services. Summary of the Invention

[0006] Technical Problem: The purpose of this invention is to provide a channel estimation and user localization method based on a spherical array intelligent metasurface, which integrates RIS-assisted communication and sensing into a system, greatly expanding the application scenarios of RIS in wireless systems. It can achieve unambiguous parameter decoupling and channel estimation that existing planar array RIS cannot support with lower hardware costs and training overhead, thereby enabling high-precision user localization applications.

[0007] Technical solution:

[0008] S1, RIS systems and channel modeling

[0009] The spherical array RIS uses N R There are n reflecting elements, regularly distributed around a sphere of radius R, where the spherical coordinate azimuth angle of element n is φ. n The pitch angle is θ n Let the incident azimuth angle of the far-field plane electromagnetic wave be φ, the elevation angle be θ, and the direction vector be... The RIS spherical array manifold is then...

[0010]

[0011] Where, k0=2π / λ c Based on wavelength λ c The electromagnetic wave number. The subscript R refers to the RIS-related parameter, exp(x) = e x Let be an exponential function with base e, (·) T Indicates transpose. It represents the imaginary unit.

[0012] The base station in the RIS-assisted ISAC wireless system is configured with N B A single-antenna array with U users configured with a single antenna, a broadband system with K subcarriers, and a fundamental frequency of f. c The transmission bandwidth is f s The spherical array RIS is deployed near the user side. The user-RIS channel is dominated by line-of-sight paths, while the base station-RIS channel has blocked line-of-sight paths and is composed of non-line-of-sight scattering paths. Let... Let m×n represent the real and complex number fields, then the uplink channel from the base station and user u to RIS of subcarrier k. and Modeled as

[0013]

[0014] In this context, subscripts B, U, D, and A refer to base station BS, user User, departure angle AoD, and arrival angle AoA, respectively. P represents H. BR,k Number of scattering paths, f k =fc +kf s / K represents the subcarrier frequency k. τ BR,p ,α BR,p H respectively BR,k The time delay and fading of the scattering path p, τ RU,u ,α RU,u h respectively RU,u,k Time delay and fading of the line-of-sight path. φ B,p ,φ D,p ,θ D,p H respectively BR,k The azimuth angle of arrival at the base station side, the azimuth angle of departure from the RIS side, and the elevation angle of the scattering path p, φ A,u ,θ A,u h respectively RU,u,k The azimuth and elevation angles of arrival on the RIS side of the line-of-sight path. Based on the free-space loss principle, the fading of the scattering path and the line-of-sight path are modeled as α. BR,p ∝(4πf c τ BR,p ) -1 and α RU,u =(4πf c τ RU,u ) -1 The proportional symbol ∝ indicates that the former includes unknown scattering loss.

[0015] S2, Beam Training Method Design

[0016] Channel estimation beam training uses S tr One data stream, T tr Each time frame and K tr Each frequency domain subcarrier. Base station beamforming vectors are configured on the data stream s. Configure the RIS reflection coefficient vector within training frame t User u sends pilot symbol x u,k Then the uplink received training signal of subcarrier k is modeled as follows:

[0017]

[0018] in, For base station beamformers, For equivalent concatenated channels, It is Gaussian white noise. ⊙ represents the Khatri-Rao product, diag(x) represents a matrix with x as its diagonal elements, and eq refers to the equivalent channel correlation parameters. H eq,u,k According to equation (2), it can be rewritten as follows:

[0019]

[0020] in, This represents the Hadamard product. Substituting equation (4) into equation (3) yields...

[0021]

[0022] Where I = PU represents the number of cascaded paths, and the two-dimensional index is mapped to a one-dimensional index. and Let be the equivalent fading and equivalent delay of path i, respectively. and These are the sets of azimuth and elevation angle parameters for path i, respectively. and χ i,k These are the equivalent angle of arrival and equivalent pilot signal at the base station side for path i, respectively. Define the mapping relationship. Then there is and

[0023] Merge K tr T on each subcarrier tr Received training signals within each training frame to construct a third-order data tensor The pattern-1 fiber vector with index (t,k) Store signal vector The effective signal part follows the standard multivariate (CP) tensor model.

[0024]

[0025] in, This is a RIS reflection coefficient sequence. This is the base station equivalent array response vector. This is the RIS array response vector. For path equivalent gain, where [g eq,i ] k =χ i,k exp(-j2πf k τ eq,i ). For the kernel tensor, where All other elements are 0. It is a Gaussian noise tensor. × represents the cross product of vectors. n Tensor-matrix product representing pattern-n.

[0026] Assume the base station antenna array manifold is a B (μ)=[exp(-jμ(N B -1) / 2),...,1,...,exp(jμ(N B -1) / 2)] T ,in Antenna spacing d B The spatial frequency corresponding to the lower azimuth angle φ. The base station beamformer is designed using Discrete Fourier Transform (DFT).

[0027]

[0028] in, For a search accuracy of 2π / N B Integer search codewords.

[0029] To utilize the symmetry of the spherical array, the RIS reflection coefficient is designed based on the phase mode excitation principle of the array element domain manifold. Specifically, the weighted beamforming vector of the (l,m)-order spherical harmonics is expressed as...

[0030]

[0031] Where l∈{0,1,...},m∈{-l,...,l}, the mapping relationship between the order (l,m) and the training frame sequence t is as follows: Y represents the weighting of the signal response of array element n. l,m (φ,θ) is the (l,m)-order spherical harmonic function.

[0032]

[0033] Among them, P l m (cosθ) represents the associated Legendre function, defined as follows:

[0034]

[0035] Among them, P l (cosθ) represents the l-th order Legendre function.

[0036] S3, Channel Angle Parameter Estimation

[0037] Based on the properties of the CP tensor, the training signal Mode-1 expansion Represented as

[0038]

[0039] in, This is the actual array response vector of the base station. For noise tensor The mode-1 expansion. Equation (7) shows that the beamforming vector can be used to derive the real-valued beam domain manifold. Design search code {ω s} constitutes D BFor each search domain, the beamformer W B The phase shift rotation invariance property of the base station antenna array is preserved as follows:

[0040]

[0041] in, To select the matrix, the unitary rotation-invariant signal parameter estimation algorithm (U-ESPRIT) is used to obtain the estimated value of the actual angle of arrival of the base station. Reconstructing the base station array response vector calculate in Represents the pseudo-inverse of a matrix. New tensor Pattern-(2,3) slice matrix It can be represented as

[0042]

[0043] Among them, i p (u) = (p-1)U+u is the path index based on the mapping relationship. The equivalent noise. Solving the channel parameters using equation (13) is equivalent to a compressed sensing subproblem with P observation matrices of type Ψ and a component number of U.

[0044] Based on the RIS topology of the spherical array, the path The electromagnetic response at RIS element n is:

[0045]

[0046] Where, r eq,i , These represent the magnitude, azimuth, and elevation angles of the equivalent direction vector of the cascaded path, respectively, and the angular parameters φ that make up the path. eq,i ,θ eq,i It has the following equation relationship

[0047]

[0048] Let d(φ) D,p ,θ D,p )+d(φ A,u ,θ A,u The element value of ) is x eq,i ,y eq,i ,z eq,i Then there is mod N (·) represents the modulo-N operation, and arctan2(·,·) and arccos(·) represent the arctangent and arccosine operations in the four quadrants, respectively.

[0049] For simplicity, the subscript eq,i of the relevant equivalent angle parameters is temporarily ignored. The spherical function of equation (8) is decomposed into the sum of the spherical Fourier modes, which can be rewritten as:

[0050]

[0051] Where, j l (·) denotes the first kind of l-order spherical Bessel function, (·) * This represents the conjugate operation. Combining equations (8) and (16), the excitation of the spherical phase mode on the pairwise manifold can be obtained as follows:

[0052]

[0053] The discrete excitation objective of spherical phase modes is to select the corresponding spherical Fourier modes for the array element domain manifold. Due to the exponential decay of the spherical Bessel function, a finite-sized RIS spherical array aperture can only excite a finite number of modes. Let L be the highest order of the exciteable modes, and the value of L is determined empirically.

[0054]

[0055] in, This indicates rounding up. If the modulus r of the equivalent direction vector of the cascaded path is ≤ 2, then the upper bound of the highest order L is... The total number of effectively excited spherical phase modes is T tr =(L+1) 2 .

[0056] Design RIS reflector unit parameters {φ n ,θ n} follows a t-design distribution, for The following spherical harmonic orthogonality holds true

[0057]

[0058] Where δ(x) represents the discrete delta function. The beamforming weights are designed as follows: Equation (17) can be simplified to

[0059]

[0060] design and The beam domain array response vector of RIS beamforming is

[0061]

[0062] The (l,m)th order spherical harmonics have the following recursive properties.

[0063]

[0064] in, To utilize this property, define Let be a subvector of the beam-domain manifold, where To select the matrix, extract sub-vectors of each order. The first, middle, and last (2l+1) elements. The recursive relation of equation (22) is expressed by the following relation: The elements are connected

[0065]

[0066] in, diagonal array Defined respectively

[0067]

[0068] Calculate the autocorrelation matrix of formula (13) in(·) H This represents the conjugate transpose. Eigenvalue decomposition. We can obtain the eigenvector set corresponding to U largest eigenvalues. To characterize the signal subspace. There exists a non-singular transformation matrix. satisfy in for The corresponding beam domain manifold vector is equation (21). Equation (23) can then be rewritten in subspace form.

[0069]

[0070] in, c = -1, 0, 1. in When L 2 When ≥2U, equation (25) is an overdetermined system of equations, and the least squares solution is: Eigenvalue decomposition yields Repeating the ESPRIT subspace algorithm described above can solve the P-group compressed sensing subproblems.

[0071] S4, Parameter Decoupling and User Location

[0072] get Subsequently, the two sets of angle estimates for the equivalent direction vector of the cascaded path are...

[0073]

[0074] Where |·| and arg(·) represent the modulus and phase operations respectively, and arctan(·) represents the arctangent operation. This is used to filter out the correct angle estimates and estimate the modulus {r} of the equivalent direction vector. eq,i},right Eigenvalue decomposition is performed to obtain (T) tr -U) eigenvector groups corresponding to the smallest eigenvalues To characterize the noise subspace. Based on the orthogonality between the signal and noise subspaces, a beam domain spectrum for the Multi-Signal Classification (MUSIC) algorithm is constructed.

[0075]

[0076] Substitute into each group And perform a one-dimensional search; only the correct estimate will... A distinct spectral peak appears, and its location is the estimated value of the equivalent direction vector magnitude.

[0077] according to recover for It is possible to construct about {φ D,p ,θ D,p ,φ A,u ,θ A,u There are a total of 3U effective constraint equations with 2(U+1) unknown variables. When U≥2, it is an overdetermined system of equations. Solving this nonlinear system of equations can achieve unique decoupling of the channel angle parameters. By statistically analyzing the decoupling results of the solutions to the P groups of compressed sensing subproblems, more accurate estimates can be obtained.

[0078] Will Substitute back into the RIS array manifold of equation (6) and calculate in Includes the estimated parameters corresponding to the RIS and the base station array response vector. for The mode-3 expansion. Channel estimation beam training uses continuous K... tr If ≤K adjacent subcarriers transmit pilot symbols, then the equivalent time delay has a least-squares solution.

[0079]

[0080] Equivalent path fading can be recovered using the following formula:

[0081]

[0082] Select the equivalent path corresponding to the base station to RIS channel path p. According to the free space loss model, we can obtain... The estimated propagation delay from the user to the RIS channel can be obtained by the following formula.

[0083]

[0084] In addition, there are Based on the delay parameter estimation results, the estimated value of the corresponding path fading can be easily derived. Thus, the spherical array RIS-assisted channel estimation method achieves unique decoupling of the multipath angle, delay, and fading parameters of the base station and user-to-RIS channel.

[0085] Based on the angle and time delay parameters of the line-of-sight path from the user to the RIS, user positioning applications can be directly implemented using the orientation-distance positioning method.

[0086]

[0087] in, For the prior physical location information of RIS, v c This represents the speed of electromagnetic wave propagation.

[0088] The channel estimation and user localization method based on the spherical array RIS of this invention is summarized as follows: This algorithm mainly uses algebraic operations and does not involve random initialization, loop iteration, etc. Therefore, it can obtain high-precision channel estimation results with low computational complexity and robust operation.

[0089]

[0090] Beneficial effects: Compared with the prior art, the advantages of the present invention are:

[0091] 1) This invention designs a new topological paradigm for spherical array RIS, which supports electromagnetic wave signal reflection and propagation in the three-dimensional full spatial domain, and greatly expands the half-space domain supported by planar array RIS.

[0092] 2) This invention designs a channel estimation method based on spherical array RIS, which supports the accurate recovery and uniqueness decoupling of three-dimensional spatial path parameters. Compared with planar array RIS, which has fuzzy parameter estimation, it achieves a qualitative breakthrough.

[0093] 3) This invention designs a RIS reflection coefficient training pattern based on the phase mode excitation principle, which can significantly reduce the need for large-scale training.

[0094] Pilot overhead for channel estimation in RIS-assisted wireless systems.

[0095] 4) This invention designs a user localization method based on RIS deep channel estimation, realizing RIS-based communication and sensing integration, and effectively expanding the application scenarios of RIS in wireless systems. Attached Figure Description

[0096] Figure 1 This is a flowchart illustrating the implementation of the channel estimation and user localization method based on the spherical array RIS in this invention.

[0097] Figure 2 This is a schematic diagram of a wireless communication and sensing system based on a spherical array RIS assisted according to an embodiment of the present invention.

[0098] Figure 3 This is a diagram showing the relationship between the channel angle parameter decoupling accuracy and signal-to-noise ratio based on the spherical array RIS in an embodiment of the present invention.

[0099] Figure 4 This is a graph showing the relationship between the channel delay parameter decoupling accuracy and signal-to-noise ratio based on a spherical array RIS in an embodiment of the present invention.

[0100] Figure 5 This is a graph showing the relationship between the channel fading parameter decoupling accuracy and signal-to-noise ratio based on a spherical array RIS in an embodiment of the present invention.

[0101] Figure 6 This is a graph showing the relationship between the accuracy and signal-to-noise ratio of user positioning applications based on a spherical array RIS in an embodiment of the present invention. Detailed Implementation

[0102] The present invention will now be described in detail with reference to the accompanying drawings and embodiments, so as to more clearly and thoroughly illustrate the purpose, technical solution and outstanding advantages of the present invention.

[0103] like Figure 1 As shown, the channel estimation and user localization method based on a spherical array RIS provided in this embodiment of the invention includes the following design steps. It should be emphasized that, provided the operational requirements of the method of this invention are met, the specific numerical values ​​in the embodiments do not limit the scope of application of this invention.

[0104] S1, RIS systems and channel modeling

[0105] This invention designs a system consisting of N R =1024 passive reflective elements forming a spherical array RIS with a radius of R = 2λ c / 3, where λ c =0.01m is the system reference wavelength. The spherical coordinate azimuth angle of array element n is φ. n The pitch angle is θ n To avoid spatial aliasing, this embodiment of the invention is designed to satisfy a t-design distribution with parameter t = 31. The direction vector of the far-field plane wave with azimuth and elevation angles (φ, θ) is... Excited RIS spherical array manifold

[0106]

[0107] Where, k0=2π / λ c =200πrad / m is the electromagnetic wave number. The subscript R refers to the RIS-related parameters, exp(x) = e x Let be an exponential function with base e, (·) T Indicates transpose. It represents the imaginary unit.

[0108] like Figure 2 As shown, this embodiment of the invention designs a RIS-assisted ISAC wireless system, wherein the base station is configured with N B =16 antenna elements, forming a uniform linear array with half-wavelength spacing, U=2 user-configured single antennas, the broadband system has K=128 subcarriers, and the fundamental frequency is f. c =30GHz, transmission bandwidth is f s =0.32GHz. The spherical array RIS is deployed near the user side. The base station and users are randomly distributed within a spherical domain with a radius of 15m centered on the RIS. Each user's channel to the RIS is dominated by one line-of-sight path, while the line-of-sight path to the RIS from the base station is blocked, containing P = 2 non-line-of-sight scattering paths. Let Let m×n represent the real and complex number fields, then the uplink channel from the base station and user u to RIS of subcarrier k. and Modeled as

[0109]

[0110] In this context, subscripts B, U, D, and A refer to base station BS, user User, departure angle AoD, and arrival angle AoA, respectively. P represents H. BR,k Number of scattering paths, f k =f c +kf s / K represents the subcarrier frequency k. τ BR,p ,α BR,p H respectively BR,k The time delay and fading of the scattering path p, τ RU,u ,α RU,u h respectively RU,u,k Time delay and fading of the line-of-sight path. φ B,p ,φ D,p ,θ D,p H respectively BR,k The azimuth angle of arrival at the base station side, the azimuth angle of departure from the RIS side, and the elevation angle of the scattering path p, φ A,u ,θ A,u h respectively RU,u,k The azimuth and elevation angles of arrival on the RIS side of the line-of-sight path. Based on the free-space loss principle, the fading of the scattering path and the line-of-sight path are modeled as α.BR,p ∝(4πf c τ BR,p ) -1 and α RU,u =(4πf c τ RU,u ) -1 The proportional symbol ∝ indicates that the former includes unknown scattering loss. According to an embodiment of the present invention, the wireless device distribution is such that the path space angle {φ}... B,p}、{φ D,p ,φ A,u} and {θ D,p ,θ A,u The values ​​are randomly generated within the ranges [60°, 120°], [-180°, 180°], and [0°, 180°], respectively, with a propagation delay of [missing information]. Randomly generated within (0ns, 50ns) (ns: nanosecond), the unknown scattering loss is randomly generated within (0.5, 1.0).

[0111] S2, Beam Training Method Design

[0112] Channel estimation beam training uses S tr = 9 data streams, T tr = 36 time frames and K tr ={64,128} subcarriers. Configure base station beamforming vectors on data stream s. Configure RIS reflection coefficient training vectors within training frame t User u sends pilot symbol x u,k Then the uplink received training signal of subcarrier k is modeled as follows:

[0113]

[0114] in, For base station beamformers, For equivalent concatenated channels, It is Gaussian white noise. ⊙ represents the Khatri-Rao product, diag(x) represents a matrix with x as its diagonal elements, and eq refers to the equivalent channel correlation parameters. H eq,u,k According to equation (2), it can be rewritten as follows:

[0115]

[0116] in, This represents the Hadamard product. Substituting equation (4) into equation (3) yields...

[0117]

[0118] Where I = PU represents the number of cascaded paths, and the two-dimensional index is mapped to a one-dimensional index. and Let be the equivalent fading and equivalent delay of path i, respectively. and These are the sets of azimuth and elevation angle parameters for path i, respectively. and χ i,k These are the equivalent angle of arrival and equivalent pilot signal at the base station side for path i, respectively. Define the mapping relationship. Then there is and

[0119] Merge K tr T on each subcarrier tr Received training signals within each training frame to construct a third-order data tensor The pattern-1 fiber vector with index (t,k) Store signal vector The effective signal part follows the standard multivariate (CP) tensor model.

[0120]

[0121] in, This is a RIS reflection coefficient sequence. This is the base station equivalent array response vector. This is the RIS array response vector. For path equivalent gain, where [g eq,i ] k =χ i,k exp(-j2πf k τ eq,i ). For the kernel tensor, where All other elements are 0. It is a Gaussian noise tensor. × represents the cross product of vectors. n Tensor-matrix product representing pattern-n.

[0122] Assume the base station antenna array manifold is a B (μ)=[exp(-jμ(N B -1) / 2),...,1,...,exp(jμ(N B -1) / 2)] T ,in Antenna spacing d B =λ c / 2=0.5cm corresponds to the spatial frequency of the azimuth angle φ. The base station beamformer is designed using Discrete Fourier Transform (DFT).

[0123]

[0124] in, For a search accuracy of 2π / N B The integer search codeword is given below. In this embodiment of the invention, the search codebook is designed as ω. s ∈{0,...,(S tr -1) / 2,N B -(S tr -1) / 2,...,N B -1}.

[0125] To utilize the symmetry of the spherical array, the RIS reflection coefficient is designed based on the phase mode excitation principle of the array element domain manifold. Specifically, the weighted beamforming vector of the (l,m)-order spherical harmonics is expressed as...

[0126]

[0127] Where l∈{0,1,...},m∈{-l,...,l}, the mapping relationship between the order (l,m) and the training frame sequence t is as follows: Y is the weighting weight of the signal response of array element n. l,m (φ,θ) is the (l,m)-order spherical harmonic function.

[0128]

[0129] Among them, P l m (cosθ) represents the associated Legendre function, defined as follows:

[0130]

[0131] Among them, P l (cosθ) represents the l-th order Legendre function.

[0132] S3, Channel Angle Parameter Estimation

[0133] Based on the properties of the CP tensor, the training signal Mode-1 expansion Represented as

[0134]

[0135] in, This is the actual array response vector of the base station. For noise tensor The mode-1 expansion. Equation (7) shows that the beamforming vector can be used to derive the real-valued beam domain manifold. This invention provides an embodiment for searching the codeword {ω} s} constitutes D B =1 search domain, then beamformer W B The phase shift rotation invariance property of the base station antenna array is preserved as follows:

[0136]

[0137] in, To select the matrix, the unitary rotation-invariant signal parameter estimation algorithm (U-ESPRIT) is used to obtain the estimated value of the actual angle of arrival of the base station. Reconstructing the base station array response vector calculate in Represents the pseudo-inverse of a matrix. New tensor Pattern-(2,3) slice matrix It can be represented as

[0138]

[0139] Among them, i p (u) = (p-1)U+u is the path index based on the mapping relationship. The equivalent noise. Solving the channel parameters using equation (13) is equivalent to a compressed sensing subproblem with P observation matrices of type Ψ and a component number of U.

[0140] Based on the RIS topology of the spherical array, the path The electromagnetic response at RIS element n is:

[0141]

[0142] Where, r eq,i , These represent the magnitude, azimuth, and elevation angles of the equivalent direction vector of the cascaded path, respectively, and the angular parameters φ that make up the path. eq,i ,θ eq,i It has the following equation relationship

[0143]

[0144] Let d(φ) D,p ,θ D,p )+d(φ A,u ,θ A,u The element value of ) is x eq,i ,y eq,i ,z eq,i Then there is mod N(·) represents the modulo-N operation, and arctan2(·,·) and arccos(·) represent the arctangent and arccosine operations in the four quadrants, respectively.

[0145] For simplicity, the subscript eq,i of the relevant equivalent angle parameters is temporarily ignored. The spherical function of equation (8) is decomposed into the sum of the spherical Fourier modes, which can be rewritten as:

[0146]

[0147] Where, j l (·) denotes the first kind of l-order spherical Bessel function, (·) * This represents the conjugate operation. Combining equations (8) and (16), the excitation of the spherical phase mode on the pairwise manifold can be obtained as follows:

[0148]

[0149] The discrete excitation objective of spherical phase modes is to select the corresponding spherical Fourier modes for the array element domain manifold. Due to the exponential decay of the spherical Bessel function, a finite-sized RIS spherical array aperture can only excite a finite number of modes. Let L be the highest order of the exciteable modes, and the value of L is determined empirically.

[0150]

[0151] in, This indicates rounding up. If the modulus r of the equivalent direction vector of the cascaded path is ≤ 2, then the upper bound of the highest order L is... In this embodiment of the invention, the highest order L = 5, corresponding to a total number of effective excitation spherical phase modes of T. tr =(L+1) 2 =36.

[0152] The azimuth and elevation angle parameters of the RIS array elements in this embodiment of the invention are {φ} n ,θ n} follows a t-design distribution with parameter t=31, for The following spherical harmonic orthogonality holds true

[0153]

[0154] Where δ(x) represents the discrete delta function. In this embodiment of the invention, the beamforming weights are set to... Equation (17) can be simplified to

[0155]

[0156] design and The beam domain array response vector of RIS beamforming is

[0157]

[0158] The (l,m)th order spherical harmonics have the following recursive properties.

[0159]

[0160] in, To utilize this property, define Let be a subvector of the beam-domain manifold, where To select the matrix, extract sub-vectors of each order. The first, middle, and last (2l+1) elements. The recursive relation of equation (22) is expressed by the following relation: The elements are connected

[0161]

[0162] in, diagonal array Defined respectively

[0163]

[0164] Calculate the autocorrelation matrix of formula (13) in(·) H This represents the conjugate transpose. Eigenvalue decomposition. We can obtain the eigenvector set corresponding to U largest eigenvalues. To characterize the signal subspace. There exists a non-singular transformation matrix. satisfy in for The corresponding beam domain manifold vector is equation (21). Equation (23) can then be rewritten in subspace form.

[0165]

[0166] in, c = -1, 0, 1. in When L 2 When ≥2U, equation (25) is an overdetermined system of equations, and the least squares solution is: Eigenvalue decomposition yields Repeating the ESPRIT subspace algorithm described above can solve the P-group compressed sensing subproblems.

[0167] S4, Parameter Decoupling and User Location

[0168] get Subsequently, the two sets of angle estimates for the equivalent direction vector of the cascaded path are...

[0169]

[0170] Where |·| and arg(·) represent the modulus and phase operations respectively, and arctan(·) represents the arctangent operation. This is used to filter out the correct angle estimates and estimate the modulus {r} of the equivalent direction vector. eq,i},right Eigenvalue decomposition is performed to obtain (T) tr -U) eigenvector groups corresponding to the smallest eigenvalues To characterize the noise subspace. Based on the orthogonality between the signal and noise subspaces, a beam domain spectrum for the Multi-Signal Classification (MUSIC) algorithm is constructed.

[0171]

[0172] Substitute into each group And perform a one-dimensional search; only the correct estimate will... A distinct spectral peak appears, and its location is the estimated value of the equivalent direction vector magnitude.

[0173] according to recover for It is possible to construct about {φ D,p ,θ D,p ,φ A,u ,θ A,u The system comprises 3U effective constraint equations with a total of 2(U+1) unknown variables. When U≥2, it is an overdetermined system of equations. This invention employs the interior-point method or the Levenberg-Marquardt algorithm to solve this nonlinear system of equations, achieving unique decoupling of the channel angle parameters. By statistically analyzing the decoupling results of P groups of compressed sensing subproblems, more accurate estimates can be obtained.

[0174] Will Substitute back into the RIS array manifold of equation (6) and calculate in Includes the estimated parameters corresponding to the RIS and the base station array response vector. for The mode-3 expansion. This embodiment of the invention employs continuous K... tr If ≤K adjacent subcarriers transmit training pilot symbols, then the equivalent time delay has a least-squares solution.

[0175]

[0176] Equivalent path fading can be recovered using the following formula:

[0177]

[0178] Select the equivalent path corresponding to the base station to RIS channel path p. According to the free space loss model, we can obtain... The estimated propagation delay from the user to the RIS channel can be obtained by the following formula.

[0179]

[0180] In addition, there are Based on the delay parameter estimation results, the estimated value of the corresponding path fading can be easily derived. Thus, the spherical array RIS-assisted channel estimation method achieves unique decoupling of the multipath angle, delay, and fading parameters of the base station and user-to-RIS channel.

[0181] Based on the angle and time delay parameters of the line-of-sight path from the user to the RIS, user positioning applications can be directly implemented using the orientation-distance positioning method.

[0182]

[0183] in, For the prior physical location information of RIS, v c =3×10 8 m / s is the propagation speed of electromagnetic waves. In this embodiment of the invention, RIS is set at the origin p of the Cartesian coordinate system. R =[0,0,0] T .

[0184] Based on the above system parameters and scheme design, this embodiment of the invention is repeatedly run for 10 cycles. 4 The Monte Carlo simulation was used to test the performance of the channel estimation and user localization method based on the spherical array RIS. The relationships between the recovery and decoupling accuracy of the path angle, delay, and fading parameters obtained by the channel estimation method and the system signal-to-noise ratio (SNR) are as follows: Figure 3 , Figure 4 and Figure 5 As shown in the figure. The relationship between user positioning accuracy and received SNR is as follows. Figure 6 As shown, the distance between the actual location and the estimated location is the absolute error, and the ratio of the absolute error to the distance from the user to the RIS is the relative error. Channel estimation and user positioning performance are measured by the root mean square error (RMSE). Simulation results show that the channel estimation method provided in this embodiment can achieve accurate channel parameter recovery and unique decoupling performance, in T... tr <<N R Under these conditions, it is particularly capable of achieving angle estimation accuracy on the order of 0.1° and 10... -2 The delay estimation accuracy is on the order of nanoseconds, thus supporting user positioning accuracy on the order of 1.0 cm.

[0185] In summary, the channel estimation and user localization method based on spherical array RIS provided in this invention fully explores the potential of the RIS topology to "sense" the spatial propagation environment by modeling and analyzing the channel parameter information contained in the RIS array manifold, thereby overcoming the major challenge of unambiguous decoupling of passive RIS cascade channel path parameters. Furthermore, this invention utilizes the electromagnetic reconfigurability of the RIS reflector units and designs a special RIS beamforming pattern based on the phase mode excitation principle, significantly reducing training pilot overhead.

[0186] In contrast, generally (N) R,h ×N R,v The RIS uniform planar array manifold is in These are uniform linear array manifolds in the horizontal and vertical directions, respectively. Signal processing can recover at most the spatial angular parameters of the apertures in two orthogonal directions, i.e., the equivalent direction vector d(φ). D,p ,θ D,p )+d(φ A,u ,θ A,u Two sets of nonlinear constraints with three element values ​​are used to construct at most 2(P+U-1) sets of constraints for 2(P+U) unknown variables. The solution to this underdetermined system of equations is not unique, therefore it cannot support the unique decoupling and recovery of angle, delay, and fading parameters. Although wireless communication applications such as beamforming and capacity optimization only require cascaded channels H... BR,k diag(h RU,u,k While wireless sensing applications such as user positioning and environmental mapping can obtain information, they require real and unique channel propagation information. The channel estimation and user positioning method based on a spherical array RIS provided in this invention overcomes the bottlenecks of existing technologies through low-cost, low-power, programmable, and easily deployable RIS technology, achieving the organic integration of wireless communication and sensing systems.

Claims

1. A channel estimation and user localization method based on a spherical array RIS, characterized in that, The method includes the following steps: Step S1, RIS System and Channel Modeling: Design the topology and cell layout of the spherical array RIS, and establish the RIS-assisted communication and sensing system and multipath channel model; Step S2, beam training method design: Establish a multi-carrier beam training signal model, and design a base station beamforming and RIS reflection coefficient training pattern assignment method based on beam domain manifold transformation and spherical phase mode excitation; Step S3, Channel equivalent parameter estimation: Establish the compressed sensing subproblem and design a path equivalent angle parameter estimation method based on the characteristics of spherical array manifold and spatial spectrum estimation algorithm; Step S4, Parameter Decoupling and User Localization: Design an angle, delay, and fading decoupling method based on RIS array topology and free space propagation characteristics to realize user localization application; Step S1, which involves designing the topology and cell layout of the spherical array RIS, and establishing a RIS-assisted communication and sensing system and a multipath channel model, specifically includes: The spherical array RIS uses N R There are n reflecting elements, regularly distributed around a sphere of radius R, where the spherical coordinate azimuth angle of element n is φ. n The pitch angle is θ n Let the incident azimuth angle of the far-field plane electromagnetic wave be φ, the elevation angle be θ, and the direction vector be... The RIS spherical array manifold is then... Where, k0=2π / λ c Based on wavelength λ c The electromagnetic wave number, where the subscript R refers to the RIS-related parameter, exp(x) = e x Let be an exponential function with base e, (·) T Indicates transpose. Represents the imaginary unit; The base station in the RIS-assisted communication and sensing system is configured with N B A single-antenna array with U users configured with a single antenna, a broadband system with K subcarriers, and a fundamental frequency of f. c The transmission bandwidth is f s The spherical array RIS is deployed near the user side. The user-to-RIS channel is dominated by the line-of-sight path, while the base station-to-RIS channel has an obstructed line-of-sight path and is composed of non-line-of-sight scattering paths. Let m×n represent the real and complex number fields, then the uplink channel from the base station and user u to RIS of subcarrier k. and Modeled as Wherein, subscripts B, U, D, and A refer to the base station BS, user User, departure angle AoD, and arrival angle AoA, respectively; P represents H BR,k Number of scattering paths, f k =f c +kf s / K is the subcarrier frequency k, τ BR,p ,α BR,p H respectively BR,k The time delay and fading of the scattering path p, τ RU,u ,α RU,u h respectively RU,u,k The time delay and fading of the line-of-sight path, φ B,p ,φ D,p ,θ D,p H respectively BR,k The azimuth angle of arrival at the base station side, the azimuth angle of departure from the RIS side, and the elevation angle of the scattering path p, φ A,u ,θ A,u h respectively RU,u,k The azimuth and elevation angles of arrival on the RIS side of the line-of-sight path; based on the free-space loss principle, the fading of the scattering path and the line-of-sight path are modeled as α respectively. BR,p ∝(4πf c τ BR,p ) -1 and α RU,u =(4πf c τ RU,u ) -1 The proportional symbol ∝ indicates that the former includes unknown scattering loss; Step S2 involves establishing a multi-carrier beam training signal model and designing a base station beamforming and RIS reflection coefficient training pattern assignment method based on beam domain manifold transformation and spherical phase mode excitation. Specifically, this includes: Channel estimation beam training uses S tr One data stream, T tr Each time frame and K tr Each frequency domain subcarrier; base station beamforming vector configured on data stream s. Configure the RIS reflection coefficient vector within training frame t User u sends pilot symbol x u,k Then the uplink received training signal of subcarrier k is modeled as follows: in, For base station beamformers, For equivalent concatenated channels, For Gaussian white noise, ⊙ represents the Khatri-Rao product, diag(x) represents a matrix with x as its diagonal element, and eq refers to the equivalent channel correlation parameter; H eq,u,k Rewritten according to equation (2) as follows: in, Representing the Hadamard product; substituting equation (4) into equation (3) yields Where I = PU represents the number of cascaded paths, and the two-dimensional index is mapped to a one-dimensional index. and Let be the equivalent fading and equivalent delay of path i, respectively. and These are the sets of azimuth and elevation angle parameters for path i, respectively. and χ i,k These represent the base station-side equivalent angle of arrival and equivalent pilot signal for path i, respectively; define the mapping relationship. Then there is and Merge K tr T on each subcarrier tr Received training signals within each training frame to construct a third-order data tensor The pattern-1 fiber vector with index (t,k) Store signal vector The effective signal part obeys the standard multivariate tensor model. in, This is a RIS reflection coefficient sequence. This is the base station equivalent array response vector. This is the RIS array response vector. For path equivalent gain, where [g eq,i ] k =χ i,k exp(-j2πf k τ eq,i ); For the kernel tensor, where The remaining elements are 0. Let be the Gaussian noise tensor; o denotes the vector outer product, × n Tensor-matrix product representing pattern-n; The base station antenna array manifold is in Antenna spacing d B The spatial frequency corresponding to the lower azimuth angle φ; the base station beamformer is designed using discrete Fourier transform. in, For a search accuracy of 2π / N B Integer search codewords below; To utilize the symmetry of the spherical array, based on the phase mode excitation principle of the array element domain manifold, the spherical Fourier transform is used to design the RIS reflection coefficient; specifically, the weighted beamforming vector of the (l,m)th order spherical harmonics is expressed as: Where l∈{0,1,...},m∈{-l,...,l}, the mapping relationship between the order (l,m) and the training frame sequence t is as follows: Y is the weighting of the signal response of array element n; l,m (φ,θ) is the (l,m)-order spherical harmonic function. Among them, P l m (cosθ) represents the associated Legendre function, defined as follows: Among them, P l (cosθ) represents the l-th order Legendre function. Step S3 involves establishing a compressed sensing subproblem and designing a path equivalent angle parameter estimation method based on the characteristics of a spherical array manifold and a spatial spectrum estimation algorithm. Specifically, this includes: Based on the properties of the CP tensor, the training signal Mode-1 expansion Represented as in, This is the actual array response vector of the base station. For noise tensor The pattern-1 expansion; Equation (7) derives the real-valued beam domain manifold from the beamforming vector. Design search code {ω s } constitutes D B For each search domain, the beamformer W B The phase shift rotation invariance property of the base station antenna array is preserved as follows: in, To select the matrix, a unitary rotation-invariant signal parameter estimation algorithm is used to obtain an estimate of the actual angle of arrival at the base station. Reconstructing the base station array response vector calculate in Represents the pseudo-inverse of a matrix; new tensor Pattern-(2,3) slice matrix Represented as Among them, i p (u) = (p-1)U+u is the path index based on the mapping relationship. The equivalent noise; solving the channel parameters by equation (13) is equivalent to a compressed sensing subproblem with P observation matrices of type Ψ and number of components U; According to the RIS topology of the spherical array, path i∈L p ={i p (1),...,i p The electromagnetic response of (u)} at the RIS array element n is: Where, r eq,i , Represent the magnitude, azimuth, and elevation angles of the cascaded path equivalent direction vector, respectively, and the channel path angle parameter φ. eq,i ,θ eq,i It has the following equation relationship Let d(φ) D,p ,θ D,p )+d(φ A,u ,θ A,u The element value of ) is x eq,i ,y eq,i ,z eq,i Then there is mod N (·) represents the modulo-N operation, and arctan2(·,·) and arccos(·) represent the arctangent and arccosine operations in the four quadrants, respectively. Ignoring the subscripts eq,i of the relevant equivalent angle parameters for the time being, the spherical function of equation (8) is decomposed into a sum of spherical Fourier modes, and rewritten as Where, j l (·) denotes the first kind of l-order spherical Bessel function, (·) * Represents the conjugate operation; combining equations (8) and (16), the excitation of the spherical phase mode on the pairwise manifold is obtained as follows: The discrete excitation objective of spherical phase modes is to select the corresponding spherical Fourier modes for the array element domain manifold. Due to the exponential decay of the spherical Bessel function, a finite-sized RIS spherical array aperture can only excite a finite number of modes. Let L be the highest order of the exciteable modes, and determine the value of L according to empirical rules. in, This indicates rounding up; if the magnitude of the equivalent direction vector of the cascaded path is r≤2, then the upper bound of the highest order L is... The total number of effectively excited spherical phase modes is T tr =(L+1) 2 ; Design RIS reflector unit parameters {φ n ,θ n } follows a t-design distribution, for For |m|≤l,|m′|≤l′, the following spherical harmonic orthogonality holds. Where δ(x) represents the discrete delta function; the designed beamforming weights are... Equation (17) simplifies to design and The beam domain array response vector of RIS beamforming is The (l,m)th order spherical harmonics have the following recursive properties. in, To utilize this property, define c = -1, 0, 1 are subvectors of the beam domain manifold, where To select the matrix, extract sub-vectors of each order. The first, middle, and last (2l+1) elements; the recursive relation of equation (22) is expressed by the following relation: The elements are connected in, diagonal array Defined respectively Calculate the autocorrelation matrix of formula (13) in(·) H Represents conjugate transpose; eigenvalue decomposition Obtain the eigenvector group corresponding to the U largest eigenvalues To characterize the signal subspace; there exists a non-singular transformation matrix. satisfy in for The corresponding beam domain manifold vector is equation (21); at this time, equation (23) is rewritten in subspace form. in, in When L 2 When ≥2U, equation (25) is an overdetermined system of equations, and the least squares solution is: Eigenvalue decomposition yields Repeat the ESPRIT subspace algorithm described above to solve the P groups of compressed sensing subproblems; Step S4 involves designing a decoupling method for angle, delay, and fading based on the RIS array topology and free-space propagation characteristics to implement user positioning applications. Specifically, this includes: get Subsequently, the two sets of angle estimates for the equivalent direction vector of the cascaded path are... Where |·| and arg(·) represent the modulus and phase operations respectively, and arctan(·) represents the arctangent operation; to filter out the correct angle estimates and estimate the modulus {r} of the equivalent direction vector. eq,i },right Eigenvalue decomposition is performed to obtain (T) tr -U) eigenvector groups corresponding to the smallest eigenvalues To characterize the noise subspace; based on the orthogonality between the signal and noise subspaces, a beam domain spectrum for a multi-signal classification algorithm is constructed. Substitute into each group And perform a one-dimensional search; only the correct estimate will... A distinct spectral peak appears, and its location is the estimated value of the equivalent direction vector magnitude. according to recover for Construct about {φ D,p ,θ D,p ,φ A,u ,θ A,u A total of 3U effective constraint equations with 2(U+1) unknown variables are used. When U≥2, it is an overdetermined system of equations. Solving this nonlinear system of equations achieves unique decoupling of the channel angle parameters. The decoupling results of the solutions to the P groups of compressed sensing subproblems are statistically analyzed to obtain more accurate estimates. Will Substitute back into the RIS array manifold of equation (6) and calculate in Includes the estimated parameters corresponding to the RIS and the base station array response vector. for Mode-3 expansion; channel estimation beam training uses continuous K tr If ≤K adjacent subcarriers transmit pilot symbols, then the equivalent time delay has a least-squares solution. Equivalent path fading is recovered by the following formula Select the equivalent path i corresponding to the base station to RIS channel path p. p (u)∈L p According to the free space loss model, The estimated propagation delay from the user to the RIS channel is obtained by the following formula. In addition, there are Based on the time delay parameter estimation results, it is easy to deduce the estimated value of the corresponding path fading. Thus, the spherical array RIS-assisted channel estimation method has achieved unique decoupling of multipath angle, delay, and fading parameters of the base station and user-to-RIS channel; Based on the angle and time delay parameters of the line-of-sight path from the user to the RIS, the user positioning application directly achieves positioning through the orientation-distance method. in, For the prior physical location information of RIS, v c This represents the speed of electromagnetic wave propagation.

Images