A Joint Estimation Method for the Number of Signal Sources and Directions of Arrival in a Large-Scale MIMO Square Array

Abstract

The present invention discloses a method for jointly estimating the number of signal sources and direction of arrival in a massive MIMO array, which includes: establishing a massive MIMO array model in an impulsive noise environment; constructing a weighted signal subspace fitting equation based on Laplacian kernel correntropy; using the piecewise idea to simplify and obtain the objective function; initializing the individual quantum positions to obtain the global optimal quantum position; initializing the quantum ice crystal energy value to determine the position of the temporary lake center; updating the energy value and the historical quantum position space; updating the quantum positions; generating a new generation of quantum ice crystals according to roulette wheel selection and updating the global optimal quantum position; determining whether the maximum number of iterations is reached. If not, return to step five; determining whether the signal source exists. If it exists, return to step three, otherwise output the number of signal sources and the corresponding direction of arrival. The present invention has the characteristics of robustness, high precision and a wider application range.

Description

Technical Field

[0001] The present invention belongs to the field of array signal processing, and relates to a method for jointly estimating the number of signal sources and direction of arrival (DOA) of a large-scale MIMO square array, in particular to a method for jointly estimating the number of signal sources and DOA of a Massive MIMO square array based on a quantum ice crystal optimization mechanism in an impulsive noise environment. Background Art

[0002] DOA estimation technology, as one of the key research directions in array signal processing, has very wide applications in multiple fields such as wireless communication, radar, and positioning. Especially in the fifth-generation mobile communication system, large-scale MIMO technology is considered one of the key technologies of 5G and has received extensive attention from industry insiders. In a large-scale MIMO system, two-dimensional DOA estimation can simultaneously estimate the elevation angle and azimuth angle, can provide a more accurate position estimate for the mobile station, and can make the directivity of beamforming stronger. It is the key to the implementation of 3D-MIMO technology. Therefore, it is of great significance to conduct in-depth research on DOA estimation algorithms. However, in current two-dimensional DOA estimation algorithms, most are extended from the MUSIC algorithm and ESPRIT algorithm based on linear arrays, and mainly obtain the DOA estimation value through joint search in the two-dimensional angle domain of the azimuth angle and elevation angle. These methods not only require the prior knowledge of the number of signal sources, but also often require a large number of snapshots to ensure the estimation accuracy and accuracy. In the face of the huge antenna array of a large-scale MIMO system, some high-precision spectral peak search methods have very high computational complexity. And for coherent signal sources, it is necessary to sacrifice the array aperture to perform the decoherence operation. Therefore, for the planar array system of large-scale MIMO, it is very necessary to design a two-dimensional DOA estimation method that is resistant to impulsive noise, fast, efficient, and has low computational complexity in the case of unknown number of signal sources.

[0003] After a search of the existing literature, it was found that Wu W et al. proposed a traditional Capon DOA estimation algorithm in "A low copmlexity 2D DOAestimation algorithm for massive MIMO systems" published in the 《IEEE Intemational Conference onSignal,Information and DataProcessing》. First, the DOA of the signal is initially estimated by discrete Fourier transform, and then the reduced-dimension Capon algorithm is used to further estimate it within a small range. This algorithm can appropriately reduce the computational complexity, but the resolution of this algorithm is limited and the accuracy still needs to be further improved. Zheng Z et al. proposed using a generalized array manifold matrix for beamspace transformation of distributed sources and constructing a shift-invariant structure in the beamspace to avoid the operation of spectral peak search in "Efficient beam space-based algorithm for two-dimensional DOA estimation ofincoherently distributed sources in massive MIMO systems" published in the 《IEEE Transactions on Vehicular Technology》, so as to reduce the dimension of matrix operations and have great advantages in terms of computational complexity. However, its accuracy performance is average, and the algorithm performance deteriorates seriously in the case of impulsive noise. Zhu Y et al. proposed a two-dimensional DOA angle estimation method based on the MUSIC algorithm in "DoA estimation and capacity analysis for 2D activemassive MIMO systems" published in the 《IEEE International Conference onCommunications》. Since the MUSIC algorithm needs to perform spectral peak search in the entire domain, although the literature uses a segmented method for angle estimation, in the context of a large-scale antenna array system, the proposed algorithm still has a high computational complexity and inevitably has quantization errors.Wang A et al. published "Low complexity direction of arrival(DoA)estimation for 2D massive MIMO systems" in 《IEEE GLOBECOM Workshops》, proposing to use the Unitary ESPRIT algorithm for DoA angle estimation of large-scale antenna arrays. The Unitary ESPRIT algorithm converts the direction matrix and covariance matrix into real values, cleverly simplifying the imaginary part, thus reducing the computational complexity. However, the article does not sufficiently analyze the joint analysis of the elevation angle and azimuth angle, the estimation accuracy is limited, and a large number of snapshots are required. It can be seen that most of the existing research focuses on reducing the computational complexity of the method at the cost of sacrificing the estimation accuracy. However, improving the accuracy and success probability of the direction finding method under low signal-to-noise ratio is also an important research direction. Moreover, the above methods all assume that the number of signal sources is known, but in actual applications, the number of signal sources is not known. Therefore, additional signal source number estimation is required first, and the signal source number estimation lags behind the direction of arrival estimation, affecting the practical process of high-precision direction of arrival estimation, and none of them consider the impact of impulsive noise on the direction finding results.

[0004] In addition, in the existing literature, most of these direction finding methods assume that the background noise is Gaussian noise, and ideal results can be obtained by using second-order or higher-order cumulants for analysis. However, in actual signal transmission, the noise more often does not follow a Gaussian distribution, such as sea clutter noise, atmospheric noise, and wireless channel noise, etc. These noises can be modeled by the Alpha stable distribution, which is mismatched with the Gaussian noise model, making the traditional algorithms based on second-order or higher-order cumulants ineffective. Therefore, it is of great significance and value to study the joint estimation of the number of signal sources and two-dimensional direction of arrival with decoherence, small snapshots, high precision, and high robustness in a large-scale MIMO system under impulsive noise. Summary of the Invention

[0005] Aiming at the above-mentioned existing technologies, the technical problem to be solved by the present invention is to provide a joint estimation method for the number of signal sources and direction of arrival of a large-scale MIMO square matrix based on a quantum ice crystal optimization mechanism with higher effectiveness and robustness, realizing the joint estimation of the number of signal sources and two-dimensional direction of arrival in a Massive MIMO system in the presence of impulsive noise, small snapshots, and coherent sources.

[0006] To solve the above technical problems, a joint estimation method for the number of signal sources and direction of arrival of a large-scale MIMO square matrix of the present invention includes:

[0007] Step 1: Establish a Massive MIMO square matrix model in an impulsive noise environment and obtain the snapshot data received by the array;

[0008] Step 2: Use the received data to construct a low-order covariance matrix based on the Laplacian kernel correlation entropy, and obtain a general weighted signal subspace fitting equation based on the Laplacian kernel correlation entropy;

[0009] Step 3: Use the piecewise-by-piece solution method to obtain the final objective function of the incoming wave to be searched;

[0010] Step 4: Initialize and estimate the number of ice crystal individuals of the quantum ice crystal mechanism for the l-th incoming wave, the quantum positions of each quantum ice crystal and the corresponding mapped state positions, and obtain the global optimal quantum position;

[0011] Step 5: Initialize the quantum ice crystal energy value, determine a temporary lake center position, and exchange energy with the outside world;

[0012] Step 6: Update the quantum ice crystal energy. The N q quantum ice crystals with the lowest energy values start to precipitate and freeze, are added to the existing ice crystal shell, and the historical quantum position space is updated

[0013] Step 7: According to the quantum evolution rule, use the simulated quantum rotation gate to update the quantum positions of the unprecipitated quantum ice crystals and the corresponding mapped state positions;

[0014] Step 8: Calculate the fitness value of each quantum ice crystal, select quantum ice crystals according to the roulette wheel selection strategy to generate the final positions of the new generation of quantum ice crystals, and then update the global optimal quantum position;

[0015] Step 9: Determine whether the maximum iteration number T max is reached. If not, let t = t + 1, and return to Step 5 to continue the iteration; if the maximum iteration number is reached, select the mapped state position of the current global optimal quantum position as the final result, and output the mapped state position at this time

[0016] Step 10: Determine whether the signal source l exists: Set a threshold If then it is judged that there are no unknown signal sources in the space, and then output the number of the measured unknown signal sources in the space and the estimated values of the corresponding two-dimensional directions of arrival. The number of signal sources is l - 1; if then select the current mapped state position as the estimated value of the two-dimensional direction of arrival of this signal source, then let l = l + 1, return to Step 3, and substitute the obtained direction of arrival estimated value as the prior information obtained into the simplification of the objective function for the next solution, and then search again.

[0017] Furthermore, the snapshot data received by the array in Step 1 satisfies:

[0018] Assume there are B narrowband far-field signal sources in space, with azimuth angles θ b and elevation angles incident on a Massive MIMO horizontal uniform square array composed of M×N array elements. For b = 1, 2,..., B, the element spacing is d, and the incident wavelength is λ. Then the mathematical model of the k-th snapshot sampling data received by the array is:

[0019]

[0020] where k = 1, 2,..., K, K is the maximum number of snapshots, z(k)=[z1(k),z2(k),...,z MN (k)] T is the k-th snapshot data vector received by the MN×1-dimensional array. The superscript T represents transpose, is the MN×B-dimensional manifold matrix, where θ = [θ1,θ2,...,θ B and are the azimuth angle vector and elevation angle vector of the signal sources respectively, represents the steering vector of the b-th signal source, where is the steering vector of the b-th signal source of the manifold matrix in the x-axis direction, is the steering vector of the b-th signal source of the manifold matrix in the y-axis direction, is the Kronecker product, s(k) is the B×1-dimensional signal vector, and n(k) is the MN×1-dimensional impulse noise vector subject to the SαS stable distribution.

[0021] Furthermore, the low-order covariance matrix based on the Laplacian kernel correlation entropy and the general weighted signal subspace fitting equation described in Step 2 are specifically:

[0022] The element in the row and column of the low-order covariance matrix R based on the Laplacian kernel correlation entropy is expressed as:

[0023]

[0024] In the formula, and respectively represent the [[ID=[]52]]-th dimension and the -th dimension of the k-th snapshot signal data vector received, (·) * represents the conjugate operation, and η is the kernel length of the kernel function;

[0025] Given the number B of unknown signal sources in space, then perform eigenvalue decomposition on the low-order covariance matrix R:

[0026]

[0027] Among them, U S is the signal subspace spanned by the eigenvectors corresponding to the B larger eigenvalues, and Σ S is the diagonal matrix composed of the B larger eigenvalues, and U N is the noise subspace spanned by the eigenvectors corresponding to the remaining smaller eigenvalues, and Σ N is the diagonal matrix composed of the remaining smaller eigenvalues, and (·) H represents the conjugate transpose operation, and the orthogonal projection matrix of the uniform square matrix is (·) -1 represents the inverse operation;

[0028] Then the general formula of the weighted signal subspace fitting equation based on the Laplacian kernel correlation entropy is:

[0029]

[0030] Among them, is the optimal weight matrix of the signal subspace, where σ 2 represents the noise power, and tr(·) is the matrix trace function.

[0031] Furthermore, the final objective function of the incoming wave to be searched in step three includes:

[0032]

[0033] Among them, the label of the l-th incoming wave unknown in space is l, and initially the unknown incoming wave label l = 1, is the defined normalized vector, satisfying:

[0034]

[0035] Among them, and respectively represent the azimuth angle parameter vector and the elevation angle parameter vector of the (l - 1)-dimensional that have been searched, and the low-order covariance matrix R is decomposed into: Among them is the signal subspace spanned by the eigenvectors corresponding to the first l larger eigenvalues, is the diagonal matrix composed of the first l larger eigenvalues, and σ 2 represents the noise power.

[0036] Furthermore, initializing and estimating the number of ice crystal individuals of the quantum ice crystal mechanism of the l-th incoming wave, the quantum position of each quantum ice crystal, and the corresponding mapped state position and obtaining the global optimal quantum position in step four include:

[0037] Set the number of individual quantum ice crystals as N p , and the spatial dimension of each quantum ice crystal is The maximum number of iterations is T max , where t represents the number of iterations. For N p quantum ice crystals, their quantum positions are randomly initialized within the quantum domain [0, 1]. Then, when estimating the l-th source, the quantum position of the i-th quantum ice crystal in the t-th generation is defined as where represents the -th dimension of the quantum position of the i-th quantum ice crystal in the t-th generation, and the corresponding mapped state position is The mapping relationship is where and represent the lower and upper boundaries of the -th dimensional variable respectively, and i = 1, 2,..., N p , Calculate the fitness value of the mapped state position of each quantum ice crystal. When estimating the l-th source, substitute the mapped state position of the i-th quantum ice crystal in the t-th generation into the fitness function expression to obtain Record the quantum position with the maximum fitness value among the quantum ice crystals guiding the t-th generation as the global optimal quantum ice crystal position The corresponding mapped state position is denoted as

[0038] Furthermore, initialize the energy value of the quantum ice crystal as described in step five, determine a temporary lake center position, and the energy exchange with the outside world includes:

[0039] For the i-th quantum ice crystal in the t-th generation, its initial energy is Then determine a temporary lake center, and use the method of combining angle and distance to dynamically adjust the lake center position. Define the diversity metric equation of the i-th quantum ice crystal in the t-th generation as represents the i-th quantum ice crystal in the t-th generation The minimum angle between the i-th quantum ice crystal and the j-th crystal, where i = 1, 2,..., N p , j = 1, 2,..., N p and i ≠ j, The larger the value, the more obvious the separation of the crystal from other crystals, and the better its diversity. Among them ||·||2 represents the second-order norm of the vector, and (·) T represents the transpose; then define the dynamic adjustment equation of the angle combined distance of the i-th quantum ice crystal in the t-th generation as where Z(ρ) is defined as the angle influence factor, τ is the control coefficient, is the energy value of the \(i\)-th quantum ice crystal in the \(t\)-th generation, and finally the temporary center position of the lake water is obtained as

[0040] The freezing process of the lake water is summarized into two stages: precipitation and freezing. First, it is judged whether the iteration number \(t\) is less than or equal to \(T\) max / 3. If so, the quantum ice crystal enters the precipitation stage and updates the energy change value of each crystal; otherwise, the quantum ice crystal enters the freezing stage and updates the energy change value of each crystal;

[0041] For the \(i\)-th quantum ice crystal in the \(t\)-th generation, its energy change equation is defined as which is mainly divided into three parts. The first part is the energy obtained by the \(i\)-th quantum ice crystal in the \(t\)-th generation from the center of the lake water, is the distance from the \(i\)-th quantum ice crystal in the \(t\)-th generation to the center of the lake water, and \(q\) is an energy constant; the second part represents the energy exchange between the \(i\)-th quantum ice crystal in the \(t\)-th generation and the ice crystals that have completed freezing, and \(O\) is an energy constant; the third part \(F\) i t =C represents the energy carried away by the wind of the \(i\)-th quantum ice crystal in the \(t\)-th generation, \(C\) is an energy constant, and \(\lambda_1\), \(\lambda_2\), and \(\lambda_3\) are all proportionality coefficients. By adjusting the three proportionality coefficients \(\lambda_1\), \(\lambda_2\), and \(\lambda_3\) to distinguish between the precipitation stage or the freezing stage.

[0042] Furthermore, as described in step six, when updating the energy of the quantum ice crystal, the \(N\) q quantum ice crystals with the lowest energy values start to precipitate and freeze, are added to the existing ice crystal shell, and the historical quantum position space is updated including:

[0043] Update the energy of the quantum ice crystal. The energy of the \(i\)-th quantum ice crystal in the \((t + 1)\)-th generation is After the energy is updated, the \(N\) q quantum ice crystals with the lowest energy values precipitate and freeze. At the same time, record the quantum positions corresponding to the first q of the largest fitness values among these \(N\) mapping positions, which form the historical quantum position space where the space represented by its corresponding mapped state position is where As the iteration progresses, the quantum positions in the historical quantum position space are continuously updated and replaced in it.

[0044] Furthermore, as described in step seven, according to the quantum evolution rule, using the simulated quantum rotation gate to update the quantum positions and the corresponding mapped state positions of the non-precipitated quantum ice crystals includes:

[0045] Define the -dimensional quantum rotation angle of the th quantum ice crystal in the th generation as where is the transfer factor, ν t is the inertia weight coefficient in the ν max is the maximum value of t ν min v t is the minimum value of is used to find the -dimensional quantum position in the th generation that is closest to Finally, update the -dimensional quantum position of each quantum ice crystal through a simplified simulated quantum rotation gate. Then, the -dimensional quantum position of the th quantum ice crystal in the

[0046] th generation is updated as

[0047] First, calculate the probability that the th quantum ice crystal individual is selected. For the th quantum ice crystal in the th generation, the probability that it is selected is If is greater than or equal to the cumulative probability and less than then the th quantum ice crystal is selected to guide the generation of a new quantum ice crystal. The formula for generating its -dimensional quantum position needs to be constrained between [0, 1]. If it exceeds the boundary value, take the boundary value. is the weight factor, is a random number uniformly distributed between [0, 1], and then calculate its corresponding mapped state position Repeat this process to generate N q quantum ice crystals, supplement the quantum ice crystals that precipitate and freeze during the update process to ensure that the number of quantum ice crystal individuals remains unchanged; finally, substitute the mapped state positions of the generated ice crystals into the fitness value formula, calculate the corresponding fitness values, and put them together with the fitness values of the original updated quantum ice crystals. Select the maximum value among them and compare it with the fitness value of the global optimal quantum position of the previous generation. Finally, select the quantum position corresponding to the larger fitness value as the global optimal quantum position guiding the (t + 1)-th generation of quantum ice crystals, denoted as Its corresponding mapped state position is denoted as

[0048] Furthermore, the setting of the threshold includes:

[0049] Set where δ is the average fitness of all mapped positions during the calculation of the iteration process, is a set multiple, and at the same time determine a lower threshold limit Take and The larger value of is used as the threshold.

[0050] Advantages of the present invention: The present invention combines the quantum optimization theory and the ice crystal optimization mechanism to obtain the quantum ice crystal optimization mechanism, further improving the global convergence performance of the ice crystal algorithm. Among them, the information data received by the Massive MIMO antenna array is used to construct a low-order covariance matrix based on the Laplacian kernel correlation entropy, eliminating the deficiency of the weak anti-impulse noise ability of the second and higher moments, and designing a weighted signal subspace fitting equation based on the Laplacian kernel correlation entropy. Moreover, the weighted signal subspace fitting method has excellent direction finding performance in the case of small snapshots and low signal-to-noise ratio, and can also process coherent sources without additional de-coherence processing. The only drawback is that as a non-linear multi-dimensional search problem, it has extremely high complexity, especially when applied to the Massive MIMO system. Therefore, by using the idea of solving one by one in a segmented manner, only one unknown source is solved in each search, while the parameters of other already obtained sources remain unchanged. This not only converts the search area from high-dimensional to low-dimensional, overcomes the problem of large computational amount of large-scale antenna arrays, greatly reduces the computational complexity, but also can jointly estimate the number of unknown sources and the direction of arrival in the case where the number of sources is unknown. And based on quantum coding and the simulated quantum rotation gate, the quantum rotation angle and the quantum evolution equation are designed, and then the quantum ice crystal optimization mechanism with a faster convergence speed is designed, further improving the search efficiency, search accuracy and anti-interference ability of the existing DOA estimation method in the search interval, and can obtain the global optimal solution of the number of sources and the two-dimensional direction of arrival more quickly. Compared with the prior art, the present invention has the following characteristics:

[0051] (1) When using traditional classical methods, such as the MUSIC algorithm, the ESPRIT algorithm or other extended algorithms therefrom, when facing the complex electromagnetic environment in practical applications, especially when performing two-dimensional direction of arrival estimation of the Massive MIMO square matrix under low signal-to-noise ratio, small snapshots and impulse noise background, its robustness and accuracy are often difficult to meet the requirements of practical applications. By adopting the weighted signal subspace fitting algorithm, the present invention designs a weighted signal subspace fitting equation based on the Laplacian kernel correlation entropy, which can effectively eliminate the influence of impulse noise and still perform effective direction of arrival estimation in the case of low signal-to-noise ratio, small snapshots and coherent sources.

[0052] (2) When using traditional methods to estimate the direction of arrival (DOA) in a Massive MIMO system, it is inevitable to face a problem of huge computational amount and very high computational complexity brought about by a large-scale antenna array. Therefore, it often sacrifices a part of the estimation accuracy of the algorithm to reduce the computational complexity of the algorithm. However, improving the estimation accuracy of the algorithm is also very important. Therefore, in the present invention, the idea of solving one by one in segments is used to further transform the weighted signal subspace fitting equation, and the high-dimensional search problem is cleverly solved through a low-dimensional search method, greatly reducing the computational complexity. At the same time, it also maintains the high-precision characteristics of the weighted signal subspace fitting algorithm. And in the present invention, the quantum optimization theory is combined with the ice crystal optimization mechanism. The quantum evolution equation and the quantum ice crystal optimization mechanism based on quantum coding and simulated quantum rotation gates have a faster convergence speed, can quickly and accurately solve the designed objective function, can further reduce the amount of operation, improve the accuracy and robustness of the solution, and make the DOA estimation result have better robustness and convergence.

[0053] (3) Traditional methods often require the number of signal sources as prior knowledge when estimating the DOA, which is a prerequisite for the implementation of the algorithm. Therefore, it is often necessary to first estimate the number of unknown signal sources additionally. However, by using the idea of solving one by one in segments, the present invention not only does not require the prior knowledge of the number of unknown signal sources, but also can complete the joint estimation of the number of unknown signal sources and the two-dimensional DOA. Therefore, the present invention has less restrictions and a wider applicable range. Description of the Drawings

[0054] Figure 1 is a flow chart of a joint estimation method for the number of signal sources and the direction of arrival in a large-scale MIMO square matrix

[0055] Figure 2 is a schematic diagram of the joint estimation of the number of signal sources and the two-dimensional direction of arrival

[0056] Figure 3 is a schematic diagram of the DOA estimation of two signal sources 30 times under impulse noise

[0057] Figure 4 is a curve of the root mean square error of two DOA estimation methods versus the number of snapshots

[0058] Figure 5 is a curve of the root mean square error of two DOA estimation methods versus the generalized signal-to-noise ratio Detailed Embodiment

[0059] The present invention will be further described below with reference to the accompanying drawings of the specification and embodiments.

[0060] The overall process of the method for jointly estimating the number of signal sources and direction of arrival in a Massive MIMO square array based on the quantum ice crystal optimization mechanism of the present invention is as follows Figure 1 as shown, the present invention includes the following steps:

[0061] Step 1: Establish a Massive MIMO square array model in an impulsive noise environment and obtain the snapshot data received by the array.

[0062] Assume that there are B narrowband far-field signal sources in space, with azimuth angles θ b and elevation angles incident on a Massive MIMO horizontal uniform square array composed of M×N array elements. b = 1, 2,..., B, the element spacing is d, and the incident wavelength is λ. Then the mathematical model of the k-th snapshot sampling data received by the array is where k = 1, 2,..., K, K is the maximum number of snapshots, z(k)=[z1(k),z2(k),…,z MN (k)] T is the k-th snapshot data vector received by the MN×1-dimensional array. The superscript T represents the transpose. is the MN×B-dimensional manifold matrix, where θ = [θ1,θ2,...,θ B and are the azimuth angle vector and elevation angle vector of the signal source respectively. represents the steering vector of the b-th signal source, where is the steering vector of the b-th signal source of the manifold matrix in the x-axis direction. is the steering vector of the b-th signal source of the manifold matrix in the y-axis direction. is the Kronecker product, s(k) is the B×1-dimensional signal vector, and n(k) is the MN×1-dimensional impulsive noise vector subject to the SαS stable distribution.

[0063] Step 2: Use the received data to construct a low-order covariance matrix based on the Laplacian kernel correlation entropy, and obtain a general weighted signal subspace fitting equation based on the Laplacian kernel correlation entropy.

[0064] Use the received data to construct a low-order covariance matrix R based on the Laplacian kernel correlation entropy. Its row column element is expressed as In the formula and respectively represent the dimension and the dimension of the k-th snapshot signal data vector received. (·) *Denotes the conjugate operation, and η is the kernel length of the kernel function.

[0065] For the direction-of-arrival (DOA) estimation using the traditional weighted signal subspace fitting method, first, the number of unknown signal sources B in the space needs to be known. Then, the low-order covariance matrix R constructed from the received data is eigen-decomposed. where U S is the signal subspace spanned by the eigenvectors corresponding to B larger eigenvalues, and Σ S is the diagonal matrix composed of B larger eigenvalues. U N is the noise subspace spanned by the eigenvectors corresponding to the remaining smaller eigenvalues, and Σ N is the diagonal matrix composed of the remaining smaller eigenvalues. (·) H denotes the conjugate transpose operation, and the orthogonal projection matrix of the uniform square matrix is (·) -1 denotes the inverse operation. Then, based on the two-dimensional DOA estimation problem of the MIMO square matrix, the general weighted signal subspace fitting equation based on the Laplacian kernel correntropy is designed as is the optimal weight matrix of the signal subspace, where σ 2 denotes the noise power, and tr(·) is the matrix trace function.

[0066] Step 3: Using the idea of piecewise and sequential solution, obtain the final objective function of the incoming wave to be searched.

[0067] Based on the weighted signal subspace fitting method, the present invention uses the idea of piecewise and sequential solution. Each search only solves for one unknown signal source, while the parameters of other already obtained signal sources remain unchanged, and no prior information about the number of signal sources is required.

[0068] First, let the label of the l-th incoming wave in the space be l. Initially, the label of the unknown incoming wave l = 1. Then, for the general weighted signal subspace fitting equation perform decomposition and transformation. Then, the weighted signal subspace fitting equation for the estimated value of the l-th unknown signal source after transformation can be expressed as where respectively represent the azimuth angle parameter vector and elevation angle parameter vector of the already searched (l - 1) dimensions. where is the signal subspace spanned by the eigenvectors corresponding to the first l larger eigenvalues, is the diagonal matrix composed of the first l larger eigenvalues, is the noise subspace spanned by the eigenvectors corresponding to the remaining smaller eigenvalues, is the diagonal matrix composed of the remaining smaller eigenvalues, Then, decompose and simplify the projection matrix in the fitting equation to obtain The first term is independent of the elevation angle and azimuth angle estimation of the l-th signal source, so it is ignored. Then, the weighted signal subspace fitting equation of the l-th signal source after further simplification is obtained. Where Finally, a normalized vector is defined. Obtain Therefore, the weighted signal subspace fitting objective equation for searching the l-th unknown incoming wave assumed to exist in the space using the quantum ice crystal optimization mechanism is expressed as

[0069] Step 4: Initialize and estimate the number of ice crystal individuals of the quantum ice crystal mechanism for the l-th incoming wave, the quantum positions and the corresponding mapped state positions of each quantum ice crystal, and obtain the global optimal quantum position.

[0070] Set the number of individuals of the quantum ice crystal to N p , and the spatial dimension of each quantum ice crystal is The maximum number of iterations is T max , where t represents the number of iterations. Randomly initialize the quantum positions of N p quantum ice crystals within the quantum domain [0, 1]. Then, the quantum position of the i-th quantum ice crystal in the t-th generation when estimating the l-th signal source is defined as Where represents the -th dimension of the quantum position of the i-th quantum ice crystal in the t-th generation, and the corresponding mapped state position is The mapping relationship is Where and represent the lower and upper boundaries of the -th dimensional variable respectively, i = 1, 2,..., N p , Calculate the fitness value of the mapped state position of each quantum ice crystal. When estimating the l-th signal source, substitute the mapped state position of the i-th quantum ice crystal in the t-th generation into the fitness function expression to obtain Record the quantum position with the maximum fitness value among the quantum ice crystals guiding the t-th generation as the global optimal quantum ice crystal position The corresponding mapped state position is denoted as

[0071] Step 5: Initialize the quantum ice crystal energy value, determine a temporary lake center position, and exchange energy with the outside world.

[0072] For the i-th quantum ice crystal in the t-th generation, its initial energy is Then determine a temporary lake center. To avoid external interference during the ice crystal freezing process and prevent it from falling into a local optimum, while accelerating the convergence speed, the position of the lake center is dynamically adjusted by combining the angle and distance. Define the diversity metric equation of the \(i\)-th quantum ice crystal in the \(t\)-th generation as denotes the \(i\)-th quantum ice crystal in the \(t\)-th generation the minimum angle between the \(i\)-th quantum ice crystal and the \(j\)-th crystal, \(i = 1, 2, \cdots, N\) p , \(j = 1, 2, \cdots, N\) p and \(i\neq j\), The larger it is, the more obvious the separation of the crystal from other crystals, and the better its diversity. Among them \(\|\cdot\|_2\) represents the second-order norm of the vector, \((\cdot)\) T represents the transpose. Then define the dynamic adjustment equation of the angle combined with distance of the \(i\)-th quantum ice crystal in the \(t\)-th generation as where \(Z(\rho)\) is defined as the angle influence factor, \(\tau\) is the control coefficient, is the energy value of the \(i\)-th quantum ice crystal in the \(t\)-th generation. Finally, the temporary lake center position is obtained as

[0073] In nature, the freezing of lake water is a complex dynamic energy change process. The freezing of ice crystals is mainly affected by three parts: the lake center, the already frozen ice crystals, and the atmosphere. To simplify the model, this search mechanism summarizes this process into two stages: precipitation and freezing. The main reason is that in the initial stage of lake water freezing, the energy radiated from the lake center dominates and affects the freezing process. In the later stage of freezing, the influence of the already formed ice crystal shell and the atmosphere dominates and absorbs energy. Therefore, in this search mechanism, the difference between the precipitation and freezing stages can be further simplified as the different energy proportions of the influencing factors in different stages. So first judge whether the iteration number \(t\) is less than or equal to \(T\) max / 3. If so, the quantum ice crystal enters the precipitation stage and updates the energy change value of each crystal; otherwise, the quantum ice crystal enters the freezing stage and updates the energy change value of each crystal.

[0074] For the \(i\)-th quantum ice crystal in the \(t\)-th generation, define its energy change equation as which is mainly divided into three parts. The first part is the energy obtained by the \(i\)-th quantum ice crystal in the \(t\)-th generation from the lake center, is the distance from the \(i\)-th quantum ice crystal in the \(t\)-th generation to the lake center, and \(q\) is an energy constant; the second part represents the energy exchange between the \(i\)-th quantum ice crystal in the \(t\)-th generation and the already frozen ice crystals, and \(O\) is an energy constant; the third part \(F\) i t$C_{t,i}$ represents the energy carried away by the $i$-th quantum ice crystal in the $t$-th generation by the wind. $C$ is a fixed energy value, and $\lambda_1$, $\lambda_2$, and $\lambda_3$ are all proportionality coefficients. By adjusting the three proportionality coefficients $\lambda_1$, $\lambda_2$, and $\lambda_3$, it is possible to distinguish whether it is in the precipitation stage or the freezing stage.

[0075] Step Six: Update the energy of the quantum ice crystals. The $N$ q quantum ice crystals with the lowest energy values start to precipitate and freeze, are added to the existing ice crystal shell, and the historical quantum position space is updated

[0076] Update the energy of the quantum ice crystals. The energy of the $i$-th quantum ice crystal in the $(t + 1)$-th generation is After the energy update, the $N$ q quantum ice crystals with the lowest energy values precipitate and freeze. At the same time, record the q quantum positions corresponding to the first mapping positions with the largest fitness values among these $N$ ice crystals to form the historical quantum position space The space represented by its corresponding mapped state position is where As the iteration progresses, continuously update and replace the quantum positions in the historical quantum position space of the quantum ice crystals.

[0077] Step Seven: According to the quantum evolution rule, use the simulated quantum rotation gate to update the quantum positions and the corresponding mapped state positions of the quantum ice crystals that have not precipitated.

[0078] The $j$-th dimensional quantum rotation angle of the $i$-th quantum ice crystal in the $(t + 1)$-th generation is defined as where is the transfer factor, is the energy value of the $i$-th crystal in the $(t + 1)$-th generation, $\mu$ is the proportionality coefficient, is a random number uniformly distributed between $[0, 1]$, $\nu$ t is the inertial weight coefficient in the $t$-th generation, $\nu$ max is the maximum value of $\nu$ t $v$ min is the minimum value of $v$ t . is used to find the $j$-th dimensional quantum position in the $t$-th generation closest to . Finally, update the quantum position of each quantum ice crystal through the simplified simulated quantum rotation gate. Then the $j$-th dimensional quantum position of the $i$-th quantum ice crystal in the $(t + 1)$-th generation is updated to Meanwhile, the dimensional mapped state position of the i-th quantum ice crystal in the (t + 1)-th generation is obtained through the mapping equation as

[0079] Step 8: Calculate the fitness value of each quantum ice crystal, select quantum ice crystals according to the roulette wheel selection strategy to generate the final positions of the quantum ice crystals in the new generation, and then update the global optimal quantum position.

[0080] Through the roulette wheel selection strategy, several individuals are selected from the quantum ice crystals after the simulated quantum rotation gate update with a certain probability to guide the generation of new quantum ice crystals. Therefore, first calculate the probability that the i-th quantum ice crystal individual is selected For the i-th quantum ice crystal in the (t + 1)-th generation, the probability that it is selected is represents the fitness value of the i-th quantum ice crystal in the (t + 1)-th generation, and then determine the cumulative probability of the i-th individual For the i-th ice crystal, its cumulative probability Then generate a random number uniformly distributed between [0, 1] If is greater than or equal to the cumulative probability and less than then the i-th quantum ice crystal is selected to guide the generation of a new quantum ice crystal, and its dimensional quantum position generation formula is needs to be constrained between [0, 1]. If it exceeds the boundary value, take the boundary value. is the weight factor, is a random number uniformly distributed between [0, 1], and then calculate its corresponding mapped state position Repeat this process to generate N q quantum ice crystals, supplement the quantum ice crystals that precipitate and freeze during the update process to ensure that the number of quantum ice crystal individuals remains unchanged. Finally, substitute the mapped state positions of the generated ice crystals into the fitness value formula, calculate the corresponding fitness values, and put them together with the fitness values of the original updated quantum ice crystals. Select the maximum value among them and compare it with the fitness value of the global optimal quantum position in the previous generation. Finally, select the quantum position corresponding to the larger fitness value as the global optimal quantum position guiding the quantum ice crystals in the (t + 1)-th generation, denoted as Its corresponding mapped state position is denoted as

[0081] Step 9: Determine whether the maximum number of iterations T max is reached. If not, let t = t + 1 and return to Step 5 to continue the iteration; if the maximum number of iterations is reached, select the mapped state position of the current global optimal quantum position as the final result and output the mapped state position at this time.

[0082] Step 10: Determine whether the signal source l exists. The quantum ice crystal optimization mechanism belongs to the maximum value search and according to the search objective function, if there is an unknown signal source, then the fitness value of the output mapping state position should be in an absolutely dominant position compared to other positions in the iterative process. Therefore, by setting a threshold Determine whether the fitness value of the output mapping position is much better than other positions, and then determine whether it is an unknown source. It is considered that the output mapping state position has not reached a state far superior to the state, so it is judged that there are no unknown sources in the space, and then the number of unknown sources measured in the space and the corresponding two-dimensional wave arrival direction estimation value are output. The number of sources is l-1. Then select the current mapping state position As the estimated value of the two-dimensional arrival direction of the source, let l=l+1, return to step 3, and bring the obtained arrival direction estimate as the obtained prior information into the simplification of the objective function for the next solution, and then search again until the output fitness value is less than the threshold and there are no unknown sources in the space.

[0083] For the threshold The setting is to determine a multiple by calculating the average fitness δ of all mapping positions during the iteration process. set up At the same time, determine a lower threshold based on the actual problem situation To ensure that the randomness or misjudgment of the search mechanism is eliminated as much as possible, the threshold can be used to determine whether the fitness value of the output position is much better than other positions, and then determine whether the output mapping state position is a two-dimensional wave direction of arrival estimate of an unknown source.

[0084] For ease of description, the weighted signal subspace fitting method based on the quantum ice crystal mechanism and the Laplace kernel function correlation entropy is abbreviated as QCEO-LLOC-WSSF, and the MUSIC array direction of arrival estimation method based on infinite norm normalization is denoted as IN-MUSIC. The simulation experiment parameters are designed as follows: for the Massive MIMO array, the number of array elements in the x-axis direction is M = 16, the number of array elements in the y-axis direction is N = 16, the array element spacing in both the x- and y-directions is d = λ / 2, the Laplace kernel function kernel length η = 129, and the maximum number of iterations is set to T. max =150, The total number of quantum ice crystals N p =80, spatial dimension Control coefficient τ = 2.5, the number of quantum positions in the ice crystal shell The energy constant value q = 1, O = 0.6, C = 0.6. For the proportionality coefficients λ1, λ2, and λ3, at the initial stage of the algorithm, when in the precipitation stage, λ1 = 1, λ2 = 0, λ3 = -0.5; when in the precipitation and freezing stage, λ1 = 0.1, λ2 = -0.85, λ3 = -0.5, and the number of precipitates N q = 20, and the weighting factor The proportionality coefficient μ = 0.005, and the multiple ν max = 2, ν min = 0, and the impact noise is set with α = 1.5. For the IN-MUSIC method, a spectral peak search method with a scanning interval of 1° is adopted.

[0085] Figure 2 In it, the number of signal sources B = 5 is set, and they are incident from the directions of {60°, 20°}, {100°, 50°}, {160°, 30°}, {160°, 50°}, {300°, 75°} respectively. Under the background of impact noise, with a generalized signal-to-noise ratio of 10 dB and the number of snapshots K = 20, it is a schematic diagram of the joint estimation of the number of signal sources and the direction of arrival. From the simulation Figure 2 It can be seen that the present invention can accurately measure the number of signal sources with high precision and full coverage of the angle measurement range.

[0086] Figure 3 In it, the number of signal sources B = 2 is set, and the two signal sources are incident from the directions of {60°, 20°} and {250°, 70°} respectively. The impact noise characteristic index α = 1.5, the generalized signal-to-noise ratio is 10 dB, the number of snapshots K = 20, and the number of Monte Carlo experiments is 30 times. From the simulation Figure 3 It can be seen that the direction finding effect is relatively ideal, the deviation between the estimated value and the true value is small, and there is only a certain deviation at very few points. Therefore, it can be seen that the designed QCEO-LLOC-WSSF can perform effective DOA estimation in an impact noise environment.

[0087] Figure 4 In it, the number of signal sources B = 2 is set, and the two signal sources are incident from the directions of {250°, 70°} and {100.5°, 30°} respectively. The QCEO-LLOC-WSSF in the present invention and the classical IN-MUSIC method are processed simultaneously to compare the influence of the number of snapshots on the direction finding effects of the two methods. The impact noise characteristic index is 1.5, the generalized signal-to-noise ratio is 10 dB, and the number of Monte Carlo tests is 30 times. From the simulation Figure 4From the experimental results, it can be seen that in the background of impulsive noise and a small number of snapshots, the designed QCEO-LLOC-WSSF method has greater superiority and can better complete the direction-of-arrival estimation in a harsh environment with higher accuracy. On the contrary, the IN-MUSIC algorithm cannot adapt to such an environment. Therefore, it also shows that the invented QCEO-LLOC-WSSF method can still perform effective DOA estimation under poor conditions such as small snapshots and impulsive noise.

[0088] Figure 5 In [reference], the number of signal sources B is set to 2, and the two signal sources are incident from directions {250°, 70°} and {100.5°, 30°} respectively. The QCEO-LLOC-WSSF in the present invention and the classical IN-MUSIC method are processed simultaneously to compare the influence of the generalized signal-to-noise ratio on the direction-finding effects of the two methods. The impulsive noise characteristic index α = 1.5, the number of snapshots K = 20, and the number of Monte Carlo trials is 30 times. From the simulation Figure 5 From the experimental results, it can be seen that the designed QCEO-LLOC-WSSF method still has good superiority, and it can be seen that at low signal-to-noise ratio, the joint estimation of the number of signal sources and the direction of arrival can still be completed under poor communication quality.

Claims

1. A method for jointly estimating the number of signal sources and direction of arrival in a large-scale MIMO square array, characterized in that, Including: Step 1: Establish a Massive MIMO square matrix model in an impulsive noise environment and obtain the snapshot data received by the array; Step 2: Use the received data to construct a low-order covariance matrix based on the Laplacian kernel correntropy and obtain a general weighted signal subspace fitting equation based on the Laplacian kernel correntropy; Step 3: Use a piecewise-by-piece solution method to obtain the final objective function of the incoming wave to be searched; Step 4: Initialize and estimate the number of ice crystal individuals of the quantum ice crystal mechanism for the l-th incoming wave, the quantum positions of each quantum ice crystal and the corresponding mapped state positions, and obtain the global optimal quantum position; Step 5: Initialize the quantum ice crystal energy value, determine a temporary lake center position, and exchange energy with the outside world; Step 6: Update the quantum ice crystal energy. The N q quantum ice crystals with the lowest energy values start to precipitate and freeze, are added to the existing ice crystal shell, and the historical quantum position space is updated Step 7: According to the quantum evolution rule, use a simulated quantum rotation gate to update the quantum positions of the un-precipitated quantum ice crystals and the corresponding mapped state positions; Step 8: Calculate the fitness value of each quantum ice crystal, select quantum ice crystals according to the roulette wheel selection strategy to generate the final positions of the new generation of quantum ice crystals, and then update the global optimal quantum position; Step Nine: Determine whether the maximum number of iterations T has been reached max , if not, let t = t + 1, and return to Step Five to continue the iteration; if the maximum number of iterations is reached, select the mapped state position of the current globally optimal quantum position as the final result, and output the mapped state position at this time Step Ten: Determine whether the signal source l exists: Set a threshold If then it is determined that there are no unknown signal sources in the space, and then the number of measured unknown signal sources in the space and the corresponding estimated values of the two-dimensional direction of arrival are output, and the number of signal sources is l-1; if then select the current mapping state position as the estimated value of the two-dimensional direction of arrival of the signal source, then let l = l + 1, return to Step Three, and use the obtained estimated value of the direction of arrival as the prior information obtained and substitute it into the simplification of the objective function for the next solution, and then search again.

2. The joint estimation method of the number of signal sources and direction of arrival for a large-scale MIMO array according to claim 1, wherein: The snapshot data received by the array in Step 1 satisfies: Assume that there are B narrowband far-field signal sources in space, with azimuth angles θ b and elevation angles incident on a massive MIMO horizontal uniform square array composed of M×N array elements. For b = 1, 2,..., B, the element spacing is d, and the incident wavelength is λ. Then the mathematical model of the k-th snapshot sampling data received by the array is as follows: where \(k = 1, 2, \cdots, K\), \(K\) is the maximum number of snapshots, \(z(k)=[z_1(k), z_2(k), \cdots, z MN (k)] T is the \(k\)-th snapshot data vector received by an \(MN\times1\) dimensional array, the superscript \(T\) represents transpose, is an \(MN\times B\) dimensional manifold matrix, where \(\theta = [\theta_1, \theta_2, \cdots, \theta B and are the azimuth angle vector and elevation angle vector of the signal sources respectively, represents the steering vector of the \(b\)-th signal source, where is the steering vector of the \(b\)-th signal source of the manifold matrix in the \(x\)-axis direction, is the steering vector of the \(b\)-th signal source of the manifold matrix in the \(y\)-axis direction, is the Kronecker product, \(s(k)\) is a \(B\times1\) dimensional signal vector, and \(n(k)\) is an \(MN\times1\) dimensional impulse noise vector subject to the \(S_{\alpha}S\) stable distribution.

3. A method for jointly estimating the number of signal sources and direction of arrival in a large-scale MIMO array according to claim 2, characterized in that: The low-order covariance matrix based on the Laplacian kernel correntropy and the general weighted signal subspace fitting equation in Step 2 are specifically: The element in the row and column of the low-order covariance matrix R based on Laplacian kernel correlation entropy is expressed as: wherein, and respectively represent the -th and -th dimensions of the received k-th snapshot signal data vector, (·) * represents the conjugate operation, and η is the kernel length of the kernel function; Given the number B of unknown signal sources in space, then perform eigenvalue decomposition on the low-order covariance matrix R: Among them, U S is the signal subspace spanned by the eigenvectors corresponding to the B larger eigenvalues, and Σ S is the diagonal matrix composed of the B larger eigenvalues, U N is the noise subspace spanned by the eigenvectors corresponding to the remaining smaller eigenvalues, and Σ N is the diagonal matrix composed of the remaining smaller eigenvalues, (·) H represents the conjugate transpose operation, and the orthogonal projection matrix of a uniform square matrix is (·) -1 represents the inverse operation; Then the general formula of the weighted signal subspace fitting equation based on the Laplacian kernel correntropy is: Among them, is the optimal weight matrix of the signal subspace, where σ 2 represents the noise power, and tr(·) is the matrix trace function.

4. A method for jointly estimating the number of signal sources and directions of arrival in a large-scale MIMO array according to claim 3, characterized in that: The final objective function of the incoming wave to be searched in Step 3 includes: Among them, the label of the l-th incoming wave in space that is unknown is l. Initially, the label of the unknown incoming wave l = 1. is the defined normalized vector, satisfying: Among them, and respectively represent the azimuth parameter vector and the elevation angle parameter vector of the (l - 1)-dimensional that have been searched out. The low-order covariance matrix R is decomposed as: Among them is the signal subspace spanned by the eigenvectors corresponding to the first l larger eigenvalues, is a diagonal matrix composed of the first l larger eigenvalues, and σ 2 represents the noise power.

5. A method for jointly estimating the number of signal sources and directions of arrival in a large-scale MIMO array according to claim 4, characterized in that: Initializing and estimating the number of ice crystal individuals of the quantum ice crystal mechanism for the l-th incoming wave, the quantum positions of each quantum ice crystal and the corresponding mapped state positions, and obtaining the global optimal quantum position in Step 4 includes: Set the number of individual quantum ice crystals to N p , and the spatial dimension of each quantum ice crystal is The maximum number of iterations is T max , where t represents the iteration number. For N p quantum ice crystals, the quantum positions are randomly initialized within the quantum domain [0, 1]. Then, the quantum position of the i-th quantum ice crystal in the t-th generation when estimating the l-th source is defined as where represents the -th dimension of the quantum position of the i-th quantum ice crystal in the t-th generation, and the corresponding mapped state position is The mapping relationship is where and represent the lower and upper bounds of the -th dimensional variable respectively, and i = 1, 2,..., N p , Calculate the fitness value of the mapped state position of each quantum ice crystal. When estimating the l-th source, substitute the mapped state position of the i-th quantum ice crystal in the t-th generation into the fitness function expression to obtain Record the quantum position with the maximum fitness value among the quantum ice crystals guiding the t-th generation as the global optimal quantum ice crystal position The corresponding mapped state position is denoted as 6. A method for jointly estimating the number of signal sources and directions of arrival in a large-scale MIMO array according to claim 5, characterized in that: Initializing the quantum ice crystal energy value, determining a temporary lake center position, and exchanging energy with the outside world in Step 5 includes: For the i-th quantum ice crystal of the t-th generation, its initial energy is Then, a temporary lake center is determined. By using the method of combining angle and distance, the position of the lake center is dynamically adjusted. The diversity metric equation of the i-th quantum ice crystal of the t-th generation is defined as Denote the i-th quantum ice crystal of the t-th generation as The minimum angle between the i-th quantum ice crystal of the t-th generation and the j-th crystal, where i = 1, 2,..., N p , j = 1, 2,..., N p and i ≠ j, The larger the value, the more obvious the separation of the crystal from other crystals, and the better its diversity. Among them, ||·||2 represents the second-order norm of the vector, and (·) T represents the transpose; then the dynamic adjustment equation of the angle combined distance of the i-th quantum ice crystal of the t-th generation is defined as where Z(ρ) is defined as the angle influence factor, τ is the control coefficient, is the energy value of the i-th quantum ice crystal of the t-th generation. Finally, the temporary position of the lake center is obtained as The process of lake water freezing is summarized into two stages: precipitation and freezing. First, it is judged whether the iteration number t is less than or equal to T max / 3. If so, the quantum ice crystal enters the precipitation stage and updates the energy change value of each crystal; otherwise, the quantum ice crystal enters the freezing stage and updates the energy change value of each crystal; For the i-th quantum ice crystal in the t-th generation, its energy change equation is defined as Divided into three parts, the first part is the energy obtained by the i-th quantum ice crystal in the t-th generation from the center of the lake. is the distance from the i-th quantum ice crystal in the t-th generation to the center of the lake, and q is a fixed energy value; the second part represents the energy exchange between the i-th quantum ice crystal in the t-th generation and the ice crystals that have completed freezing. O is a fixed energy value; the third part F i t = C represents the energy carried away by the wind from the i-th quantum ice crystal in the t-th generation. C is a fixed energy value, and λ1, λ2, and λ3 are all proportionality coefficients. By adjusting the three proportionality coefficients λ1, λ2, and λ3, it is possible to distinguish whether it is in the precipitation stage or the freezing stage.

7. A method for jointly estimating the number of signal sources and direction of arrival in a large-scale MIMO array according to claim 6, characterized in that: Update the quantum ice crystal energy described in Step 6. The N q quantum ice crystals with the lowest energy values start to precipitate and freeze, are added to the existing ice crystal shell, and update the historical quantum position space including: Update the energy of the quantum ice crystal. The energy of the i-th quantum ice crystal in the (t + 1)-th generation is After the energy is updated, the N q quantum ice crystals with the lowest energy values precipitate and freeze. At the same time, record the first q quantum positions corresponding to the mapping positions with the largest fitness values among these N ice crystals to form the historical quantum position space Among them The space represented by its corresponding mapped state position is Among them As the iteration progresses, continuously update and replace the quantum positions in the historical quantum position space in it.

8. A method for jointly estimating the number of signal sources and direction of arrival in a large-scale MIMO array according to claim 7, characterized in that: According to the quantum evolution rule, using a simulated quantum rotation gate to update the quantum positions of the un-precipitated quantum ice crystals and the corresponding mapped state positions in Step 7 includes: Define the -dimensional quantum rotation angle of the $i$-th quantum ice crystal in the $(t + 1)$-th generation as where is the transfer factor, is the energy value of the $i$-th crystal in the $(t + 1)$-th generation, $\mu$ is the proportionality coefficient, is a random number uniformly distributed in $[0, 1]$, $\nu$ t is the inertial weight coefficient in the $t$-th generation, $\nu$ max is the maximum value of $\nu$ t and $\nu$ min is the minimum value of $\nu$, t is used to find the -th dimension of the quantum position in the $t$-th generation that is closest to Finally, update the quantum position of each quantum ice crystal through the simplified simulated quantum rotation gate. Then, the -dimensional quantum position of the $i$-th quantum ice crystal in the $(t + 1)$-th generation is updated to At the same time, obtain the -dimensional mapped state position of the $i$-th quantum ice crystal in the $(t + 1)$-th generation through the mapping equation as ​ 9. A method for jointly estimating the number of signal sources and directions of arrival in a large-scale MIMO array according to claim 8, characterized in that: Step 8: Calculating the fitness value of each quantum ice crystal, selecting quantum ice crystals according to the roulette wheel selection strategy to generate the final positions of the new generation of quantum ice crystals, and then updating the global optimal quantum position includes: First, calculate the probability that the $i$-th quantum ice crystal individual is selected For the $i$-th quantum ice crystal in the $(t + 1)$-th generation, the probability that it is selected is denotes the fitness value of the $i$-th quantum ice crystal in the $(t + 1)$-th generation, and then determine the cumulative probability of the $i$-th individual For the $i$-th ice crystal, its cumulative probability Next, generate a random number uniformly distributed in the range $[0, 1]$ If is greater than or equal to the cumulative probability and less than then the $i$-th quantum ice crystal is selected to guide the generation of a new quantum ice crystal, and its $j$-th quantum position generation formula is It needs to be constrained within $[0, 1]$. If it exceeds the boundary value, take the boundary value is the weight factor is a random number uniformly distributed in the range $[0, 1]$, and then calculate its corresponding mapped state position Repeat this process to generate $N$ q quantum ice crystals, supplement the quantum ice crystals that precipitate and freeze during the update process to ensure that the number of quantum ice crystal individuals remains unchanged; finally, substitute the mapped state positions of the generated ice crystals into the fitness value formula, calculate the corresponding fitness values, and put them together with the fitness values of the original updated quantum ice crystals. Select the maximum value among them and compare it with the fitness value of the global optimal quantum position of the previous generation. Finally, select the quantum position corresponding to the larger fitness value as the global optimal quantum position guiding the $(t + 1)$-th generation of quantum ice crystals, denoted as Its corresponding mapped state position is denoted as 10. A method for jointly estimating the number of signal sources and directions of arrival in a large-scale MIMO array according to claim 1, characterized in that: The threshold value is set as follows: Set where δ is the average fitness calculated for all mapped positions during the iterative process, is a set multiple, and at the same time a lower threshold is determined Take and the larger value in is used as the threshold.

Images