HDS architecture-based DOA estimation method in THz super-large scale MIMO system

Abstract

The invention relates to a DOA estimation method based on an HDS architecture in a THz super-large-scale MIMO system, and provides a novel hybrid dynamic sub-array (HDS) architecture, an antenna unit is strategically divided into sub-arrays, and a unique one-dimensional DOA related to phase difference between the sub-arrays is ensured; a dimensionality reduction MUSIC (IMRD-MUSIC) algorithm is provided by utilizing the framework innovation, and on the basis of the HDS framework design, by utilizing the characteristics of a signal subspace matrix, a characteristic value phase is calculated through pseudo-inverse operation, and initial elevation angle estimation is obtained; small-range one-dimensional search is carried out near the initial elevation angle, so that the estimation precision is further improved; obtaining azimuth angle estimation through least square optimization by using the extracted direction vector phase information; and automatic pairing of the elevation angle and the azimuth angle is realized. A large number of simulation results show that the solution realizes superior estimation precision and calculation efficiency, so that the solution is particularly suitable for an actual terahertz UMIMO system.

Description

Technical Field

[0001] The present invention relates to the field of wireless communication technologies, and in particular, to a method for estimating the two-dimensional Direction-of-Arrival (DOA) based on a Hybrid Dynamic Subarray (HDS) architecture in terahertz (THz) communication and ultra-large-scale MIMO systems, and proposes a DOA estimation method. Background Art

[0002] With the development of 6G wireless communication systems, terahertz communication has become an important technology to support ultra-high-speed data transmission due to its extremely wide bandwidth characteristics. However, due to the relatively high frequency band of terahertz communication, the signal propagation loss is significant, which greatly limits its effective communication distance. To overcome this problem, the ultra-large-scale Multiple-Input Multiple-Output (UM-MIMO) system provides high gain through beamforming technology, thus effectively compensating for the signal propagation loss.

[0003] The application of UM-MIMO systems in terahertz communication not only increases the transmission distance but also enables high-precision direction positioning and beam tracking of signals. However, the direct application of traditional DOA estimation methods in UM-MIMO systems faces the following problems:

[0004] (1) High computational complexity: Traditional two-dimensional spectral search algorithms, such as the MUSIC algorithm, although effective in a fully digital (FD) architecture, have an exponentially increasing computational complexity and are not suitable for real-time processing required by terahertz communication. In THz UM-MIMO, the scale of the antenna array is usually extremely large, and the computational overhead of traditional methods is too high to meet the actual needs of the system;

[0005] (2) Difficulty in adapting to the Hybrid Analog-Digital (HAD) architecture: Due to hardware complexity and cost limitations, UM-MIMO systems mostly adopt the HAD architecture to reduce the number of Radio Frequency (RF) chains. The analog / digital signal processing methods in the HAD architecture bring new problems in signal synthesis and decomposition, especially facing challenges in obtaining accurate spatial information;

[0006] (3) Estimation error of the noise subspace: Under a large-scale antenna array, the estimation of the eigenvalues and eigenvectors of the noise subspace matrix is easily affected by the array scale and the number of observation samples, resulting in significant deviations, thus affecting the DOA estimation accuracy. Especially in high-dimensional scenarios, the estimation deviation of traditional methods cannot be ignored. Summary of the Invention

[0007] The present invention overcomes the deficiencies in the existing methods. To solve the above problems, the present invention combines a closed-form solution and an improved MUSIC algorithm (IMRD-MUSIC) and proposes a DOA estimation method based on the HDS architecture in a THz ultra-large-scale MIMO system.

[0008] The technical solution proposed by the embodiment of the present invention is a DOA estimation method based on the HDS architecture in a THz ultra-large-scale MIMO system. In the HDS architecture, it is assumed that L far-field narrowband signals are incident on the antenna array from different angles, and the entire antenna array is divided into N RF sub-arrays, and each sub-array corresponds to a radio frequency chain, so the number of radio frequency chains is N RF , each sub-array includes N s = N r / N RF antennas, where N r represents the number of receiving antennas in the antenna array, and N RF << N r . The switching network in the HDS architecture allows each radio frequency chain to be dynamically connected to any one of the sub-arrays; the method includes the following steps:

[0009] Step 1: Based on T observation matrices, aggregate the observation data of each radio frequency chain. The observation data of N RF radio frequency chains form multi-dimensional observation data. According to the multi-dimensional observation data, construct an extended array output matrix and calculate the corresponding covariance matrix R y . Perform eigenvalue decomposition on the sample covariance matrix R y to obtain the noise subspace matrix and the signal subspace matrix;

[0010] Step 2: Extract two signal subspaces U s and U s1 from all the signal subspaces U s2 . The signal subspace U s1 corresponds to the signal subspace of the observation data of the first radio frequency chain, and the signal subspace U s2 corresponds to the signal subspace of the observation data of the second radio frequency chain. The first radio frequency chain is any one of the N RF radio frequency chains, and the second radio frequency chain is the radio frequency chain adjacent to any one of the N RF radio frequency chains; calculate the pseudo-inverse product of the signal subspaces U s1 and U s2 . Perform eigenvalue decomposition on the pseudo-inverse product to obtain a diagonal matrix Ψ. The diagonal matrix Ψ represents the phase difference between the observation data of the first radio frequency chain and the observation data of the second radio frequency chain, which contains the elevation angle information of each incident signal . Wherein, ​Let \(L\) denote the number of signal paths of the incident signal, and the initial elevation angle direction values of each signal path are estimated using a closed-form solution.

[0011] Step 3: Based on the initial elevation angle direction values described above Perform a one-dimensional search within a small range on the initial elevation angle direction values to optimize them. For each of the signal paths, find the corresponding optimized elevation angle estimation values to obtain the optimized elevation angle estimation values for all paths

[0012] Step 4: According to the optimized elevation angle estimation values for all paths obtained in Step 3, extract their phases and combine them with the linear least squares method to calculate the azimuth angle of each of the signal paths Finally, obtain the two-dimensional DOA estimation values

[0013] The beneficial effects of the present invention are as follows: Using a closed-form solution to quickly obtain the initial estimation values, avoiding complex full-range searches. The elevation angle and azimuth angle further reduce the computational complexity through independent optimization processes. The HDS architecture supports different-scale array configurations and can adjust the hardware design according to actual needs. The HDS architecture provides significant advantages for effective DOA estimation. By simplifying the 2-D problem and utilizing the phase differences between sub-arrays, the accuracy of DOA estimation is improved. This structure achieves the best balance between hardware complexity and DOA estimation efficiency, which is beneficial for accurate beamforming and beam tracking in terahertz UM-MIMO systems.

[0014] Preferably, in Step 1, the specific process of constructing an extended array output matrix is as follows: Based on \(T\) observation matrices, randomly generate multiple simulated combination matrices \(W\), aggregate the observation data of each RF chain, and then collect the observation data of \(N\) RF RF chains as multi-dimensional observation data to construct the extended array output matrix

[0015] wherein,

[0016]

[0017] represents a sub-observation matrix, which contains \(T\) pilot observation data from the \(q\)-th RF chain; \(A\) r represents the array manifold matrix of the entire receiving array, \(S\) represents the signal information of \(L\) signal paths; \(N\) a represents the length of the observation data; is the combination of the simulated combination matrices of all RF chains; Denotes the extended noise matrix, which is the combination of all noise matrices .

[0018] Preferably, in step 1, the sample covariance matrix R y has the following expression: where N a denotes the length of the observed data, denotes the extended observation matrix; the covariance matrix R y is subjected to eigenvalue decomposition to obtain: where Σ n denotes the diagonal eigenvalue matrix of noise obtained after eigenvalue decomposition, thus obtaining the noise subspace matrix U n and the signal subspace matrix U s , where

[0019]

[0020] Preferably, in step 2, the diagonal matrix obtained by eigenvalue decomposition of the pseudo-inverse product is specifically represented as: where T is a non-singular matrix, which represents the relationship between the signal subspace and the array manifold matrix; denotes the augmented array matrix relative to the first RF chain, the initial elevation direction value of each signal path is specifically represented as: where denotes the phase of the eigenvalues {β1, β2,..., β } obtained by eigenvalue decomposition of the pseudo-inverse product L}, which contains the elevation information of each signal path; λ denotes the signal wavelength, d denotes the element spacing, and N z denotes the number of array antennas in the vertical direction.

[0021] Preferably, the specific process of step 3 includes the following steps:

[0022] S3.1. For each initial elevation direction value calculate the direction vector and the objective function, and the direction vector is represented as: The objective function is represented as: where d1 denotes the set constraint vector, and an auxiliary matrix Q(u) is generated according to the direction vector. The auxiliary matrix Q(u) is specifically represented as: where En represents the noise subspace matrix, a z (u) represents the array steering vector in the vertical direction, d represents the element spacing, λ represents the signal wavelength, N x represents the number of array antennas in the vertical direction, u = sinφ, u represents the sine value corresponding to the current elevation angle, I Nx represents the identity matrix of dimension N x ;

[0023] S3.2. Perform a one-dimensional search within a small range according to the objective function to find the minimum of the objective function and simultaneously find the u value that minimizes the objective function: According to the described For each signal path, find the corresponding optimized elevation angle estimate where obtain all

[0024] optimized elevation angle estimates for the paths

[0025] Preferably, the specific process of step 4 includes the following steps:

[0026] S4.1. For each optimized elevation angle estimate Utilize the steering vector obtained in step S3.1 to extract the phase information of its steering vector

[0027] S4.2. Construct a least squares optimization problem according to the phase information where P represents an auxiliary matrix, represents the estimated parameter error, represents the key parameter to be estimated; ;

[0028] S4.3. Solve the least squares optimization problem to obtain a closed-form solution: The closed-form solution contains two components. The first component is the offset The second component is the offset corresponding azimuth angle Use linear least squares method to estimate the second component and estimate the azimuth angle

[0029] Retain the estimated value of each azimuth angle and the azimuth angle and elevation angle estimate Through indexing Automatically pair to obtain the output two-dimensional DOA estimation value Description of the drawings

[0030] Figure 1 It is a schematic structural diagram of the UM-MIMO system model under the HDS architecture in the present invention;

[0031] Figure 2 It is a flowchart of a DOA estimation method based on the HDS architecture in a THz ultra-large-scale MIMO system of the present invention;

[0032] Figure 3 It is a comparison chart of performance curves when the method of the present invention is compared with other algorithms;

[0033] Figure 4 It is a comparison chart of performance curves when the method of the present invention is used under conditions of different numbers of pilots and RF chains. Detailed implementation manners

[0034] The following further describes the invention with reference to the accompanying drawings and in combination with specific implementation manners, so that those skilled in the art can implement it according to the text of the specification. The protection scope of the present invention is not limited to this specific implementation manner.

[0035] The present invention relates to a DOA estimation method based on the HDS architecture in a THz ultra-large-scale MIMO system. In a THz ultra-large-scale MIMO (UM-MIMO) system, due to the extremely short wavelength and high attenuation characteristics of THz signals, in order to achieve high gain and precise beamforming, the system needs to deploy a large number of antennas; as Figure 1 shown in the HDS architecture, in the described HDS architecture, the basic unit of the antenna array is extended from a single antenna to a modular antenna group, that is, the entire antenna array is divided into N RF sub-arrays, and each sub-array corresponds to an RF chain, so the number of RF chains is N RF , each sub-array includes N s = N r / N RF antennas, where N r and N RF respectively represent the number of receiving antennas and RF chains, and N RF << N r, Different partitioning methods can be dynamically adjusted according to actual needs. The switching network in the HDS architecture allows each RF chain to be dynamically connected to any one of the sub-arrays. Under the HDS architecture, by dividing the antenna array into multiple sub-arrays and using a flexible switching network and phase modulation, signal acquisition can be achieved with fewer RF chains. The entire receiving antenna array is divided into multiple sub-arrays, each sub-array consists of a part of the antennas, and each sub-array is connected to an RF chain through the switching network; during signal transmission, the transmitting end transmits signals through pre-designed pilot signals; when each sub-array at the receiving end receives signals, due to their different physical positions, different phase delays will be generated. In order to capture these phase differences caused by the direction of arrival (DOA) of the signals, we introduce a specific phase compensation β in each sub-array l,p =(p - 1)(N z / N RF )sinφ l,r , which enables the phase information of the p-th sub-array on the l-th signal to reflect the elevation angle information of the signal.

[0036] The method of the present invention includes the following steps:

[0037] Step 1, based on T observation matrices, aggregate the observation data of each RF chain. The observation data of N RF RF chains form multi-dimensional observation data, and construct an extended array output matrix according to the multi-dimensional observation data and calculate the corresponding covariance matrix R y , perform eigenvalue decomposition on the sample covariance matrix R y to obtain the noise subspace matrix and the signal subspace matrix;

[0038] In a specific embodiment, the specific process of constructing an extended array output matrix is as follows: perform doa estimation in a classical time slot, and a total of T pilot symbols are used. At each pilot duration, the receiving end adjusts the phase of the phase shifter to capture the omnidirectional space information to obtain the observation data. For the convenience of signal processing, we reconstruct the observation data. First, aggregate the observation data under each RF chain, and then aggregate the observation data of N RF chains. For the q-th RF chain, the set of observation data under T pilots is the observation matrix of the sub-array: where, N a represents the data length, is the combined vector matrix formed by stacking the analog combined vectors in each pilot. S contains the signal information of L signal paths, and A r =A r (θ r ,φ r) represents the array manifold matrix of the entire receiving array, represents the q-th RF chain and the corresponding noise matrix. Then, by stacking the observation matrices of all RF chains (q = 1, 2, 3..., N RF ), we obtain the final extended observation matrix which represents the observation matrix obtained by reconstructing using all the observation data. Among them, is the analog combination matrix of all RF chains merged, represents the merger of all noise matrices result.

[0039] In this specific embodiment, the sample covariance matrix R y has the following expression: where N a represents the length of the observation data, represents the extended observation matrix; performing eigenvalue decomposition on the covariance matrix R y gives: where Σ n represents the diagonal eigenvalue matrix of noise obtained after eigenvalue decomposition, the diagonal eigenvalue matrix of the signal obtained after eigenvalue decomposition thus obtaining the noise subspace matrix and the signal subspace matrix

[0040] Step 2: Extract U s and U s1 and U s2 from the signal subspace matrix U two parts, and calculate the pseudo-inverse product of the two Performing eigenvalue decomposition on the pseudo-inverse product gives the diagonal matrix Ψ, and the diagonal matrix Ψ contains the elevation phase of the signal path of each signal incident on the antenna array where

[0041] In this specific embodiment, the diagonal matrix obtained by performing eigenvalue decomposition on the pseudo-inverse product is specifically expressed as: where T represents a non-singular matrix, represents the augmented array covariance matrix relative to the first RF chain, The initial elevation direction value of each signal path Specifically expressed as: Wherein, Denotes the phase of the eigenvalues {β1, β2, …, β } obtained by performing eigenvalue decomposition on the pseudo-inverse product L}, which contains the elevation information of each signal path; λ represents the signal wavelength, d represents the element spacing, and N z Represents the number of array antennas in the vertical direction; assuming the eigenvalues are {β1, β2,..., β L}, the phase of the eigenvalues Contains the elevation information of the signal path. The initial elevation direction value of the elevation is obtained using the following closed-form solution. For each signal path k, the initial elevation direction value of the elevation is given by the following formula:

[0042] Step 3. Based on the initial elevation direction value Optimize the initial elevation direction value through a small-range one-dimensional search Generate an auxiliary matrix Q(u), and for each of the signal paths, find the corresponding optimized elevation estimate value Obtain the optimized elevation estimate values for all paths

[0043] In this specific embodiment, the specific process of step 3 includes the following steps:

[0044] S3.1. For each initial elevation direction value Calculate the direction vector and the objective function. The direction vector is expressed as: The objective function is expressed as: Wherein, d1 represents a set constraint vector. Generate an auxiliary matrix Q(u) according to the direction vector. The auxiliary matrix Q(u) is specifically expressed as: Wherein, E n Represents the noise subspace matrix, and a z (u) represents the array direction vector in the vertical direction, d represents the element spacing, λ represents the signal wavelength, and N x Represents the number of array antennas in the vertical direction, u = sinφ, and u represents the sine value corresponding to the current elevation, and I Nx Represents the identity matrix of dimension N x ;

[0045] S3.2. According to the objective function, perform a one-dimensional search within a small range to find the minimum of the objective function and simultaneously find the u value that minimizes the objective function: According to the u value obtained, for each signal path, find the corresponding optimized elevation angle estimation value Wherein, Obtain the optimized elevation angle estimation values of all paths

[0046] Step 4: According to the optimized elevation angle estimation values of all paths obtained in Step 3, extract their phases and calculate the azimuth direction by combining the linear least squares method Finally, obtain the two-dimensional DOA estimation value

[0047] In this specific embodiment, the specific process of Step 4 includes the following steps:

[0048] S4.1: For each optimized elevation angle estimation value Extract the phase information of its direction vector Since the actually obtained phase information is not linear like the theoretical phase information g = [0, 2πdv / λ,..., 2π(N x -1)dv / λ] T Therefore, the least squares method needs to be used for fitting to reduce the error;

[0049] S4.2: According to the phase information described above Construct a least squares optimization problem Where P represents an auxiliary matrix, Represents the estimated parameter error, Represents the key parameter to be estimated;

[0050] S4.3: Solve the least squares optimization problem to obtain a closed-form solution: The closed-form solution contains two components. The first component is the offset The second component is the corresponding azimuth angle Use the linear least squares method for the second component Estimate and estimate the azimuth angle

[0051]

[0052] Retain the estimated value of each azimuth angle Of the azimuth angle And the elevation angle estimation value Automatically pair through the index To obtain the output two-dimensional DOA estimation value

[0053] The performance and computational efficiency of a DOA estimation method based on the HDS architecture in a THz ultra-large-scale MIMO system proposed by the present invention are analyzed through simulation experiments below.

[0054] The simulation process is carried out through MATLAB software. In the following two simulation experiments, the simulation configurations are set as follows: the center frequency is f = 1.0 THz, the total number of array antennas is N r = 1024, where the number of antennas on the x and z axes is N x = N z = 32, and the number of RF chains N RF varies between 4, 8, 16, and 32, while the number of pilot signals T varies between 8, 12, 16, and 20. In addition, three signal propagation paths are considered, consisting of one line-of-sight (LOS) path and two non-line-of-sight (NLOS) paths. The azimuth and elevation angles of the NLOS paths, where the channel attenuation is randomly generated between 5 dB and 10 dB, are randomly generated, respectively, within [-180°, 180°] and [-90°, 90°]. The root mean square error (RMSE) of the DOA estimation obtained through 1000 Monte-Carlo trials is used to evaluate the estimation performance.

[0055] Simulation Experiment 1: Compared with the existing technology using the method proposed by the present invention, the DOA estimation performance of different algorithms at different signal-to-noise ratios when N RF = 8 and T = 12 is evaluated. The simulation results are as Figure 3 shown. It can be observed that the RMSE of all algorithms decreases as the SNR increases. In addition, the IMRD-MUSIC algorithm proposed by the present invention (i.e., the method of the present invention) has an estimation performance very close to that of the HDS-based 2D-MUSIC and RD-MUSIC algorithms when SNR ≥ 0 dB, effectively demonstrating the effectiveness of the algorithm proposed by the present invention. On the other hand, it can be further observed that IMRD-MUSIC is superior to RD-MUSIC at low SNR, which can be attributed to the special design of IMRD-MUSIC. The application of the adopted HDS architecture and the LS solution in it. Combining the computational efficiency of IMRD-MUSIC, we can consider IMRD-MUSIC to be a practical solution very suitable for practical applications.

[0056] Simulation Experiment 2: Figure 2 shows the DOA estimation performance of the proposed IMRD-MUSIC at different RF chains N RF and the number of pilot signals T, with SNR and N RF fixed at 10 dB and 16, respectively. From Figure 4As can be seen, with the increase in the number of pilots, the RMSE of DOA estimation decreases significantly, highlighting the significant impact of array aperture on the estimation accuracy. However, it is worth noting that when T = 8, the RMSE of 1280 antennas is actually higher than that of 1024 antennas. This is because the number of RF chains remains unchanged, meaning that the dimension of the observed signal remains unchanged, but the increase in the number of antennas will introduce phase ambiguity, thus affecting the final estimation accuracy. In this case, increasing the number of pilots used in the estimation process can help mitigate this impact to a certain extent.

Claims

1. A DOA estimation method based on the HDS architecture in a THz ultra-large-scale MIMO system, characterized in that: In the HDS architecture, L far-field narrowband signals are assumed to be incident on the antenna array from different angles, and the entire antenna array is divided into N RF subarrays, each subarray corresponds to a RF chain, then the number of RF chains is N RF , each subarray consists of N s =N r / N RF antennas, of which N r represents the number of receiving antennas in the antenna array, and N RF <<N r , the switch network in the HDS architecture allows each RF chain to be dynamically connected to any sub-array; the method comprises the following steps: Step 1: Aggregate the observation data of each RF chain based on T observation matrices. The observation data of N RF RF chains form multi-dimensional observation data. Construct an extended array output matrix according to the multi-dimensional observation data and calculate the corresponding covariance matrix R y . Perform eigenvalue decomposition on the sample covariance matrix R y to obtain the noise subspace matrix and the signal subspace matrix; Step 2: From all the signal subspaces U s Extract two signal subspaces U s1 and U s2 , the signal subspace U s1 The signal subspace corresponding to the observation data of the first radio frequency chain, the signal subspace U s2 The signal subspace of the observed data corresponding to the second radio frequency chain, the first radio frequency chain is the N RF Any one of the N radio frequency chains, the second radio frequency chain is the N RF The radio frequency chain adjacent to any one of the radio frequency chains in the radio frequency chains; calculate the signal subspace U s1 and U s2 The pseudo-inverse product of For the pseudo-inverse product The eigenvalue decomposition is performed to obtain a diagonal matrix Ψ, which represents the phase difference between the observation data of the first radio frequency chain and the observation data of the second radio frequency chain, and contains the elevation angle information β of each incident signal l , where l = 1, ..., L, L represents the number of signal paths of the incident signal, and the closed-form solution is used to estimate the initial elevation direction value of each signal path Step 3: Based on the initial elevation angle direction value Optimize the initial elevation angle direction value through one-dimensional search within a small range For each of the signal paths, find the corresponding optimized elevation angle estimate respectively Obtain the optimized elevation angle estimates for all paths Step 4: Extract the phase from the optimized elevation angle estimates of all the paths obtained in Step 3 and calculate the azimuth angle of each of the signal paths in combination with the linear least squares method Finally, obtain the two-dimensional DOA estimate 2. The DOA estimation method based on the HDS architecture in a THz ultra-large-scale MIMO system according to claim 1, characterized in that: In step 1, the specific process of constructing an extended array output matrix is as follows: Based on T observation matrices, multiple simulated combination matrices W are randomly generated, the observation data of each RF chain are aggregated, and then the observation data of N RF RF chains are collected as multi-dimensional observation data to construct an extended array output matrix Among them, represents a sub-observation matrix, which contains T pilot observation data from the q-th RF chain; A r represents the array manifold matrix of the entire receiving array, S represents the signal information of L signal paths; N a represents the length of the observation data; is the combination of the analog combination matrices of all RF chains; represents the extended noise matrix, which is the combination of all noise matrices combined.

3. The DOA estimation method based on the HDS architecture in a THz ultra-large-scale MIMO system according to claim 2, wherein: In step 1, the sample covariance matrix R y has the following expression: where N a represents the length of the observed data, represents the extended observation matrix; for the covariance matrix R y perform eigenvalue decomposition to obtain: where Σ n represents the diagonal eigenvalue matrix of noise obtained after eigenvalue decomposition, thus obtaining the noise subspace matrix U n and the signal subspace matrix U s , where 4. The DOA estimation method based on the HDS architecture in a THz ultra-large-scale MIMO system according to claim 3, characterized in that: In step 2, for the pseudo-inverse product The diagonal matrix obtained by performing eigenvalue decomposition is specifically expressed as: where T is a non-singular matrix, which represents the relationship between the signal subspace and the array manifold matrix; represents the augmented array matrix relative to the first RF chain, The initial elevation direction value of each signal path is specifically expressed as: where ∠[β l represents the phase of the eigenvalues {β1, β2, …, β } obtained by performing eigenvalue decomposition on the pseudo-inverse product L}, which contains the elevation information of each signal path; λ represents the signal wavelength, d represents the element spacing, and N z represents the number of array antennas in the vertical direction.

5. The DOA estimation method based on the HDS architecture in a THz ultra-large-scale MIMO system according to claim 4, characterized in that: The specific process of Step 3 includes the following steps: S3.

1. For each initial elevation angle value Calculate the direction vector and the objective function. The direction vector is expressed as: The objective function is expressed as: Where d1 represents the set constraint vector. Generate the auxiliary matrix Q(u) according to the direction vector. The auxiliary matrix Q(u) is specifically expressed as: Where E n represents the noise subspace matrix, a z (u) represents the array direction vector in the vertical direction, d represents the element spacing, λ represents the signal wavelength, N x represents the number of array antennas in the vertical direction, u = sinφ, u represents the sine value corresponding to the current elevation angle, I Nx represents the identity matrix of dimension N x ; S3.

2. Perform a one-dimensional search within a small range according to the objective function to find the minimum of the objective function and simultaneously find the value of u that minimizes the objective function: According to the For each signal path, find the corresponding optimized elevation angle estimate where l = 1, …, L, to obtain the optimized elevation angle estimates for all paths 6. The DOA estimation method based on the HDS architecture in a THz ultra-large-scale MIMO system according to claim 5, characterized in that: The specific process of Step 4 includes the following steps: S4.

1. For each optimized elevation angle estimate value Using the direction vector obtained in step S3.1, extract the phase information of the direction vector S4.

2. Construct a least-squares optimization problem according to the phase information where P represents an auxiliary matrix l denotes the estimated parameter error, and v l denotes the key parameter to be estimated; S4.

3. Solve the least - squares optimization problem to obtain a closed - form solution: The closed - form solution contains two components. The first component is the offset The second component is the offset The corresponding azimuth angle Use linear least - squares method to estimate the second component And estimate the azimuth angle Retain the estimated value of each azimuth angle , and pair the said azimuth angle and elevation angle estimate automatically in pairs through index l to obtain the output two-dimensional DOA estimate

Images