CN112255629A - Sequential ESPRIT two-dimensional incoherent distribution source parameter estimation method based on combined UCA array - Google Patents

Sequential ESPRIT two-dimensional incoherent distribution source parameter estimation method based on combined UCA array Download PDF

Info

Publication number
CN112255629A
CN112255629A CN202011059330.XA CN202011059330A CN112255629A CN 112255629 A CN112255629 A CN 112255629A CN 202011059330 A CN202011059330 A CN 202011059330A CN 112255629 A CN112255629 A CN 112255629A
Authority
CN
China
Prior art keywords
array
subarray
parameter
dimensional
source
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202011059330.XA
Other languages
Chinese (zh)
Other versions
CN112255629B (en
Inventor
古清月
王惠刚
丁键
崔利杰
周中良
荣少巍
邓翔宇
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northwestern Polytechnical University
Air Force Engineering University of PLA
Original Assignee
Northwestern Polytechnical University
Air Force Engineering University of PLA
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northwestern Polytechnical University, Air Force Engineering University of PLA filed Critical Northwestern Polytechnical University
Priority to CN202011059330.XA priority Critical patent/CN112255629B/en
Publication of CN112255629A publication Critical patent/CN112255629A/en
Application granted granted Critical
Publication of CN112255629B publication Critical patent/CN112255629B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S15/00—Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems
    • G01S15/02—Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems using reflection of acoustic waves
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S3/00—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
    • G01S3/80—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using ultrasonic, sonic or infrasonic waves
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S5/00—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations
    • G01S5/18—Position-fixing by co-ordinating two or more direction or position line determinations; Position-fixing by co-ordinating two or more distance determinations using ultrasonic, sonic, or infrasonic waves
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/52—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
    • G01S7/539—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Optimization (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Data Mining & Analysis (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Acoustics & Sound (AREA)
  • Computing Systems (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

The invention discloses a sequential ESPRIT two-dimensional incoherent distribution source parameter estimation method based on a combined UCA array, which solves a central parameter and an extended parameter sequentially and jointly by utilizing a consistency linear transformation relation between a generalized array manifold and an effective subspace; the statistical characteristic parameters of the random process of the two-dimensional incoherent distribution source are used for modeling expansion parameters, so that the robustness of the high-resolution subspace method to different distribution source distribution forms is improved; the method innovatively provides the maximum characteristic value selection principle of the generalized ESPRIT algorithm by analyzing the effective signal subspace of the two-dimensional incoherent distribution source, thereby effectively solving the problem that the estimation performance of the central parameter estimation precision of the incoherent distribution source is inconsistent to the small spread angle and the large spread angle, and realizing the rapid high-precision estimation of a plurality of two-dimensional incoherent distribution sources.

Description

Sequential ESPRIT two-dimensional incoherent distribution source parameter estimation method based on combined UCA array
Technical Field
The invention relates to the technical field of signal processing, in particular to a sequential ESPRIT two-dimensional incoherent distribution source parameter estimation method based on a combined UCA array.
Background
Orientation estimation of signal targets has been one of the key issues to be solved for signal processing. Early targets were often assumed to be point target models, but in practical application environments, the target was not only a central point, but also had a certain spread width, and appeared to be distributed in space. For example, when a high-base-array antenna is used in suburbs or in the field for wireless communication, local scattering around a mobile phone can cause generalized flat rayleigh channel attenuation, and a signal source is an extended source with a certain spatial distribution for a base station. Under the shallow sea environment of far field, when using sonar to estimate the position of target under water, because the influence of channel, the acoustic signal passes through the reflection of sea level, sea surface wave, reaches the hydrophone through a large amount of scattering routes, and the output data of receiving array is the superposition response to all signals of these routes, therefore the signal source presents the spatial distribution characteristic. This multipath phenomenon is also an urgent problem to be solved in SAR radar, indoor communication, and positioning applications.
According to different practical application environments, the correlation degree between echoes at different angles in the same distributed source is different, and the distributed source can be divided into totally uncorrelated (ICD), Partially Correlated (PCD) and totally Correlated (CD)
Most of the past two decades of research on parameter estimation of totally uncorrelated (abbreviated ID) distributed sources has focused on one-dimensional ID distributed sources, and Trump first proposed a maximum likelihood Method (ML) for a single ID one-dimensional distributed source. Although ML estimation has the performance of unbiased estimation, the AML method with reduced computational complexity, which is more popular because of its great computational burden, also achieves approximate estimation accuracy, can be generalized to multiple one-dimensional distributed sources.
After a covariance fitting (COMET) method is provided for further reducing the computational burden Olivier, the ambiguity problem, the optimal initial value selection problem and the calculation simplification problem are solved successively, but the methods are based on an exponential form distribution source signal model and have no universality.
In the characteristic decomposition formed by the covariance matrix of the received data, the signal characteristic vector and the noise characteristic vector have orthogonal characteristics, so that an important subspace fitting method is developed. In the subspace formed by the covariance of the point source data, the number of characteristic values containing signal energy is the same as the number of point sources in the region of interest. However, the distributed source data covariance is a full-rank matrix, and the signal energy in the subspace of the full-rank matrix is spread to each eigenvalue, so that the performance of the traditional point source subspace fitting algorithm is sharply reduced when the parameter estimation is performed on the distributed source. And the DSPE algorithm proposed by Valee populates the classical MUSIC algorithm to the parameter estimation of a plurality of one-dimensional distribution sources, and the WPSF method proposed by mats has estimation performance close to unbiased estimation. Because the feature subspace of the covariance matrix blends the signal and the noise information together, the accurate division of the signal and the noise dimensions directly determines the parameter estimation precision of the ID distribution source under different extension degrees. Y Meng discusses the energy distribution of the distributed source signal subspace in detail, and the DISPARE algorithm, the root-MUSIC algorithm and the Rank-2 algorithm all select low-Rank approximate signal subspace, so that better precision can be obtained only during small-angle expansion. The method for avoiding the effective dimension selection of the signal subspace can obtain better estimation performance on the small-angle expanded ID one-dimensional distribution source under the condition of low signal-to-noise ratio. An array manifold generated by adopting the improved GHQ (Gauss-Hermite Quadrature) can better approximate a large-spread-angle distribution source than a traditional Taylor series. By adopting a Manifold Separation Technology (MST), better estimation precision can be obtained for parameter estimation of a large-spread-angle ID one-dimensional distribution source.
However, most algorithms need to perform spectral peak search or optimization iterative solution regardless of the covariance matrix fitting or subspace fitting method, the calculation load is large, and the optimization result is easily affected by the initial value and cannot obtain the global optimal solution. To achieve faster computation speeds many scholars aim their eyes at an ESPRIT algorithm that can avoid spectral peak searches and directly obtain a closed-form solution. S Shahbazpanahi populates the point source ESPRIT algorithm to the parameter estimation of the one-dimensional ID distribution source, so that the calculation load is greatly reduced, but the spread angle still needs spectral peak search.
In practical application, as the detection distance is shortened, a large detection target is not a one-Dimensional distribution source any more, but gradually presents the appearance characteristics of a volume target, and the target is more accurately described by a Two-Dimensional uncorrelated distribution (Two-Dimensional incorporated Distributed) source model. The parameter estimation problem of a two-dimensional uncorrelated distribution (TDID) source is that an L-shaped plane array is formed by two ULAs or two UCA arrays which are parallel in the vertical direction are provided, two-dimensional parameters are decoupled and reduced in dimension through a special arrangement mode, and four parameters of a single two-dimensional ID distribution source are searched and estimated through a subspace spectrum peak. This method is not only computationally complex, but also does not yield consistent estimates. The parameter estimation method for expanding the COMET method to the two-dimensional incoherent distribution source has better estimation precision on the central parameter, but the estimation error of the expanded parameter is larger and the calculation burden is large. The equidistant array also has the characteristic of generalized rotation invariance to the received two-dimensional ID distribution source data, so that the low calculation burden is very suitable for carrying out rapid parameter estimation on a multi-parameter two-dimensional distribution source target, and the two-dimensional distribution source ESPRIT algorithm utilizing the L-array two-dimensional decoupling characteristic greatly reduces the calculation burden but requires the array to have a wavelength less than 0.1 time of the distance. By utilizing a uniform rectangular array and adopting a two-dimensional distributed source ESPRIT algorithm, the problem of subarray spacing constraint is solved, but the calculation cost is greatly increased by adopting a large number of array elements.
Disclosure of Invention
Aiming at the existing problems, the planar Uniform Circular Array (UCA) can directly realize direction estimation of two dimensions, the effective estimation coverage angle of the horizontal azimuth angle reaches 360 degrees, and in addition, the estimation performance of the pitch angle is independent of the azimuth angle. Under the same aperture size, UCA can obtain better estimation performance than ULA, so the invention aims to provide a method for carrying out joint estimation on four parameters of a plurality of two-dimensional incoherent distribution sources by solving a closed solution based on UCA array configuration and by utilizing the translational invariance of physical space between UCAs.
In order to achieve the purpose, the technical scheme adopted by the invention is as follows:
the method for estimating the parameters of the sequential ESPRIT two-dimensional incoherent distribution source based on the combined UCA array is characterized by comprising the following steps of:
s1: constructing a multi-subarray array structure model based on uniform circular array UCA to obtain a combinable three-uniform circular subarray structure;
s2: combining the three subarrays UCSA pairwise to form three pairs of subarray UCSA combinations with equal intervals, and constructing received data of three combinations of uniform circular subarrays UCSA to obtain observation data;
s3: based on the received observation data, a Taylor approximate covariance matrix model R of a two-dimensional incoherent distribution source is constructed12、R13、R23For decoupling the central and extended parameters of the two-dimensional uncorrelated distributed sources;
s4: calculating covariance matrix of three sample data through joint data received by three groups of subarrays UCSA
Figure BDA0002711866940000021
And
Figure BDA0002711866940000022
obtaining an estimated value of the central parameter through a linear mapping relation between a covariance matrix space and a feature subspace;
s5: using Taylor approximation covariance model R12、R13、R23Linear mapping relation with the signal characteristic valueTo an estimate of the spread parameter.
Further, the specific operation steps of step S1 are as follows:
s11: establishing three uniform circular sub-arrays A with the same radius1、A2、A3The array types and array parameters of the three sub-arrays are the same, and the number of the array elements is M;
the radius r of each subarray is expressed as: r ═ γ · λ, where λ is the wavelength and γ represents a multiple of the wavelength;
the included angle alpha between adjacent array elements is as follows: α is 2 π/M;
array element spacing d1Comprises the following steps: d1＝2rsin(π/M)；
The distance between every two adjacent array elements is half wavelength, and the ratio of the radius of the array elements to the wavelength is
Figure BDA0002711866940000023
S12: setting a subarray A1The reference array element is an array element with the X-axis coordinate of (0, r), and the sub-array A1The 2 nd to M array elements are arranged anticlockwise in sequence, wherein the included angle between the first array element passing through the Y axis and the Y axis is beta degrees;
s13: at edge subarray A1Forward translation distance d of reference array element on X axisxAt, lay and A1Subarrays A of the same structure2；
S14: at edge subarray A1The first array element passing through the Y axis, namely the first array element is translated upwards by a distance d obliquely from the extension line of an included angle beta DEG of the Y axisyAt, lay and A1Subarrays A of the same structure3；
S15: according to the number M of the sub-array elements, the radius R of the sub-array and the distance d between the two pairs of sub-arraysx、dyJudging the number of types of the multiplexing schemes of the subarray elements among the subarrays;
s16: based on the radius and the number of array elements of each UCSA in the UCSA array layout, the flexible adjustment can be carried out according to communication or detection frequency under the condition of meeting the requirement of S11;
s17: and finally obtaining three pairs of UCSA array combinations with equal spacing.
Further, the specific operation steps of step S15 are as follows:
s151: calculating the number k of array elements falling in a first quadrant;
s152: when k is a multiple of 4, there are k sub-array spacing schemes that can form array element multiplexing:
Figure BDA0002711866940000031
s153: when k is a multiple of 2 instead of a multiple of 4, the sub-array spacing scheme that can form array element multiplexing is:
Figure BDA0002711866940000032
further, the structural model of the UCSA array pair composed of three uniform circular sub-arrays described in step S2 is characterized in that:
s21: subarrays A1Neutral subarray A2M groups of array element pairs formed between corresponding array elements
Figure BDA0002711866940000033
And subarrays A1And subarray A3Formed M groups of array element pairs
Figure BDA0002711866940000034
Have the same time delay therebetween, and the subarray A1And subarray A2Subarray A1And subarray A3Subarray A2And subarray A3Satisfies the following transformation relation between array manifold vectors:
Figure BDA0002711866940000035
wherein ,
Figure BDA0002711866940000037
is a subarray A2And subarray A1In between, the conventional array manifold rotation operator,
Figure BDA0002711866940000036
is a subarray A3And subarray A1The traditional array manifold rotation operator in between, and the expression of the traditional array manifold rotation operator is:
Figure BDA0002711866940000041
wherein ,
Figure BDA0002711866940000042
is a random angle of arrival.
Further, constructing the Taylor approximate covariance matrix R of the two-dimensional incoherent distribution source in step S312、R13、R23The method comprises the following operation steps:
s31: three separate uniform circular sub-arrays A1、A2 and A3The received observation data vector is spread by Taylor, in which the sub-array A1The expansion of (a) is:
Figure BDA0002711866940000043
wherein ,X1(t) is subarray A1Taylor expansion of A1(k0)Is a subarray A1Generalized array manifold, noise n1(t) is a Gaussian complex random variable with zero mean value, cyclic symmetry, independence and same distributionk(t) is a 3X 1 matrix, sk(t) is the kth transmitted signal;
likewise, subarray A2 and A3The received observation data vectors are written in a Taylor expanded form as:
Figure BDA0002711866940000044
Figure BDA0002711866940000045
s32: the combined generalized array manifold consisting of three sets of UCSAs formed by three UCSAs is:
Figure BDA0002711866940000046
Figure BDA0002711866940000047
Figure BDA0002711866940000048
wherein ,B12(0)Is a subarray A1And A2Combined generalized array manifold of (B)13(0)Is a subarray A1And A3Combined generalized array manifold of (B)23(0)Is a subarray A2And A3A joint generalized array manifold of (1);
s33: and (3) calculating a covariance matrix of the received data vector of the joint subarray consisting of the three groups of UCSA, wherein the calculation formulas are respectively as follows:
associative subarrays A1 and A2The calculation formula of the covariance matrix of the received data vector is:
Figure BDA0002711866940000051
wherein ,
Figure BDA0002711866940000052
is a space correlation matrix of the circular array noise field, H represents a conjugate transpose operation symbol, ΛΥExpanding a parameter diagonal matrix for a two-dimensional incoherent distributed source:
Figure BDA0002711866940000053
associative subarrays A1 and A3The calculation formula of the covariance matrix of the received data vector is:
Figure BDA0002711866940000054
associative subarrays A2 and A3The calculation formula of the covariance matrix of the received data vector is:
Figure BDA0002711866940000055
wherein the matrix Λ can be used no matter what type of distribution the two-dimensional incoherent distribution source conforms toΥTwo-dimensional extended parameter of
Figure BDA0002711866940000056
The distribution characteristics are expressed, and the robustness of the algorithm to the distribution type of the distributed source is ensured;
s34: the central parameter for the rotation operator of the formula (4)
Figure BDA0002711866940000057
The deterministic function is expressed as:
Figure BDA0002711866940000058
further, the specific operation of step S4 includes:
s41: covariance matrix of three sample data
Figure BDA0002711866940000059
And
Figure BDA00027118669400000510
performing characteristic decomposition, and respectively taking K eigenvectors corresponding to the first K maximum eigenvalues
Figure BDA00027118669400000511
And
Figure BDA00027118669400000512
obtaining a signal subspace after characteristic decomposition;
s42: separately computing signal subspace rotation operators
Figure BDA00027118669400000513
And
Figure BDA00027118669400000514
the evaluation value of (c) is calculated by the formula:
Figure BDA00027118669400000515
Figure BDA00027118669400000516
Figure BDA00027118669400000517
wherein ,
Figure BDA00027118669400000518
and
Figure BDA00027118669400000519
are respectively as
Figure BDA00027118669400000520
And
Figure BDA00027118669400000521
decomposing to obtain a signal subspace;
s43: separate to signal subspace rotation operators
Figure BDA00027118669400000522
And
Figure BDA00027118669400000523
performing characteristic decomposition to obtain a generalized array manifold rotation operator
Figure BDA0002711866940000061
And
Figure BDA0002711866940000062
s44: obtaining matched K groups of generalized array manifold rotation operators by using a rotation operator matching algorithm, and solving estimated values of K incoherent two-dimensional distribution source central parameters by using the K groups of generalized array manifold rotation operators
Figure BDA0002711866940000063
The closed-form solution of (a) obtains an estimated value of the central parameter, and the formula for calculating the closed-form solution is:
Figure BDA0002711866940000064
wherein ,dx and dyRespectively, the distance between the sub-arrays is,
Figure BDA0002711866940000065
Figure BDA0002711866940000066
where λ is the wavelength, c is the speed of sound in the propagation medium, and ω is the dominant frequency.
Further, the rotation operator matching algorithm of step S44 includes the following steps:
s441: subarrays A1And subarray A2Form a set of rotation operators between the array manifolds
Figure BDA0002711866940000067
Wherein takes out phii,12,i＝1；
Subarrays A1And subarray A3Form a set of rotation operators between the array manifolds
Figure BDA0002711866940000068
Each element composition has K combination schemes:
Figure BDA0002711866940000069
and, the formula for calculating the cost function of the K schemes is:
Figure BDA00027118669400000610
the combination mode corresponding to the minimum value of the cost function obtained by calculation is the correct matching mode
Figure BDA00027118669400000611
S442 at
Figure BDA00027118669400000612
Set of rotation operators
Figure BDA00027118669400000613
And
Figure BDA00027118669400000614
set of rotation operators
Figure BDA00027118669400000615
Removing successfully matched combined elements
Figure BDA00027118669400000616
S443, repeating the steps S441-S442 until all the obtained
Figure BDA00027118669400000617
The correct pairing scheme.
Further, step S5: two-dimensional uncorrelated distributed source Taylor approximate covariance model R provided by utilizing method12、R13、R23And obtaining an estimated value of the extension parameter through a linear mapping relation between the estimated value and the signal characteristic value, wherein the specific operation steps comprise:
s51: the estimated value of the center parameter obtained in step S44
Figure BDA00027118669400000618
In the formula (8), the estimated values of the combined generalized array manifold are obtained respectively
Figure BDA0002711866940000071
And
Figure BDA0002711866940000072
taking any one of the joint generalized array manifold estimates, e.g.
Figure BDA0002711866940000073
Then a joint subarray A can be obtained12The linear transformation relationship between the generalized array manifold and the signal subspace is as follows:
Figure BDA0002711866940000074
s52: by joint subarrays A12For example, and A13 and A23Similarly, the covariance of the samples of the received data
Figure BDA0002711866940000075
The characteristic decomposition is carried out, and the same can be analogized to obtain
Figure BDA0002711866940000076
And
Figure BDA0002711866940000077
taking the first 3K maximum eigenvalues formed by each signal source to form a diagonal matrix
Figure BDA0002711866940000078
In relation to linear transformation
Figure BDA0002711866940000079
Solving to obtain an extended parameter matrix:
Figure BDA00027118669400000710
s53: obtaining estimated value of two-dimensional distribution source expansion parameter by the above formula
Figure BDA00027118669400000711
Extended parameter closed-form solution:
Figure BDA00027118669400000712
s54: and obtaining an expansion parameter estimation value according to the obtained expansion parameter closed-form solution.
The invention has the beneficial effects that:
first, the algorithm of the present invention does not have the limitation that the subarray spacing must be much smaller than the wavelength;
secondly, the algorithm in the invention utilizes the rotation invariant relation of the subspace constructed by the physical translation of the uniform circular sub-array to jointly estimate and obtain the closed solution of the central parameter and the extended parameter, thereby greatly reducing the calculation cost of the high-resolution algorithm.
Thirdly, the algorithm in the invention can achieve higher estimation precision with fewer array elements, and particularly under two different conditions of large-spread angle and small-spread angle scenes, the central parameter estimation consistency is better, namely the robustness on the spread degree of the distributed source is achieved;
fourthly, the algorithm of the invention estimates the extended parameters by utilizing the statistical information of the signal subspace, which not only is robust to the distribution form of the two-dimensional distribution source, but also has better stability for the estimation precision of the extended parameters.
Drawings
FIG. 1 is a far field incident diagram of a two-dimensional incoherent distributed source signal;
FIG. 2 is a diagram of a geometric array structure based on multiple uniform circular sub-arrays;
FIG. 3 is a schematic diagram of a relationship between two-dimensional distributed source extension parameters and energy distribution in a signal subspace;
FIG. 4 is a diagram illustrating a relationship between a large eigenvalue of a two-dimensional distributed source signal subspace and a noise subspace and an angle spread parameter;
FIG. 5 is a diagram illustrating the influence of UCSA subarray element number on the estimation performance of the central parameter and the extended parameter;
FIG. 6 is a graph of the effect of sample fast beat number on the center and extended parameters;
FIG. 7 is a graph of the effect of signal-to-noise ratio on algorithm performance;
FIG. 8 is a graph of the results of the impact of two-dimensional incoherent distributed source spread angle on algorithm performance;
FIG. 9 is a result graph of the effect of the spatial angular separation of two-dimensional incoherent sources on the performance of the algorithm;
FIG. 10 is a graph of the results of two-dimensional incoherent distribution sources with different signal powers on algorithm performance;
fig. 11 is a graph of the effect of the parameter estimation performance of distributed source signals and point source signals with different spatial distribution types as a function of signal-to-noise ratio.
Detailed Description
In order to make those skilled in the art better understand the technical solution of the present invention, the following further describes the technical solution of the present invention with reference to the drawings and the embodiments.
First, referring to the incident diagram of the two-dimensional incoherent distribution source signal shown in FIG. 1, assume that there are K narrow-band incoherent two-dimensional distribution sources s in the far-field spacek(t) as shown in FIG. 1. The distributed source is composed of a cluster of signal beams with random angles according with certain distribution characteristics, and uncertain phase differences exist between the signal beams in different directions in the same distributed signal cluster, namely the different signal beams are completely uncorrelated and are incoherent distributed sources. Assuming that the source consists of L incoherent narrow-band signal beams, the data vector received by the array is:
Figure BDA0002711866940000081
wherein s (t) is a complex signal emitted from a sound source, and the distribution range of the incident signal source is thetakl(t)∈(-π,π)、
Figure BDA0002711866940000082
Each sound source bundle in the same distributed source
Figure BDA0002711866940000083
All having random path complex gain gammakl(t) and each sound source beam at any one time
Figure BDA0002711866940000084
Can be expressed as the sum of the central angle and the fluctuation angle as follows:
Figure BDA0002711866940000085
wherein ,
Figure BDA0002711866940000086
is the central angle of the two-dimensionally distributed source, and θk0Is the central azimuth of the two-dimensional distributed source,
Figure BDA0002711866940000087
is the central pitch angle of the two-dimensional distributed source; angle of fluctuation
Figure BDA0002711866940000088
Random variables that are consistent with certain distribution characteristics;
the spread parameter may be defined as the standard deviation of the Gaussian distribution obeyed by the fluctuation angle
Figure BDA0002711866940000089
Therefore, the statistical properties of any single two-dimensional incoherent distributed source can be used
Figure BDA00027118669400000810
Four parameters to measure;
it can further be seen that the noise and signal are uncorrelated between the source's wave angles, between the path gains, and the array is calibrated, the array manifold is a priori information, and the following assumptions are satisfied:
(1) assuming angle of fluctuation
Figure BDA00027118669400000811
Random variables obeying time-varying independent co-distribution, and the variance of the fluctuation angle can be expressed as:
Figure BDA00027118669400000812
(2) for incoherent distribution sources, sound source beams in the same distribution source are not correlated, and path complex gains are circularly, symmetrically, independently and uniformly distributed with zero mean value:
Figure BDA00027118669400000813
(3) the acoustic source emission signal is a constant value that varies with time and is expected to characterize the signal power of the acoustic source:
Figure BDA00027118669400000814
it can be seen that the signal response data vector received by the array is a response to superposition of a plurality of distributed source signals in a sound field, and is a zero-mean circularly symmetric complex gaussian vector;
secondly, as shown in the attached figure 2, the invention constructs three uniform circular sub-arrays A with the same radius1、A2 and A3Yellow, green and red, three sub-array arrays are same in type and array parameters, the radius of each UCA is expressed by the multiple of wavelength r-gamma-lambda, the number of array elements is M, the included angle between adjacent array elements is alpha-2 pi/M, and the distance between array elements is d12rsin (pi/M), in order to avoid phase winding, the distance between two adjacent array elements should be half wavelength, and the ratio of the radius of the array element to the wavelength is
Figure BDA0002711866940000091
For subarray A1The reference array element is an array element with X-axis coordinate of (0, r), and counterclockwise sequentially form a subarray A1The 2 nd to M array elements, wherein the included angle between the first array element passing through the Y axis and the Y axis is beta degrees;
at edge subarray A1The forward translation distance d of the X axis of the reference array elementxIs placed at the position A1Subarrays A of the same structure2；
Similarly, along subarray A1The first array element passing through the Y axis, namely the first array element and the extension line of the included angle beta DEG of the Y axis, are translated by a distance d in an inclined upward manneryAt, lay and A1Subarrays A of the same structure3In this way, the subarrays A1M array elements and subarrays A2M array elements form M pairs of equally spaced array elements, the same sub-array A1And subarray A3M pairs of array elements with equal spacing are also formed;
whether the array elements are multiplexed among the subarrays or not depends on the number M of the subarray elements, the radius R of the subarrays and the distance d between the two pairs of subarraysx、dy. In order to ensure the distribution uniformity of each array element in the circular array, the number of the array elements is even, when the number of the array elements is a multiple of 4, one array element is arranged on four direction axes in a Cartesian coordinate system, but the UCSA array elements are distributed symmetrically about the Y axis no matter whether the number of the array elements is a multiple of 4 or not. When the number of array elements is only a multiple of 2, not a multiple of 4, there are no array elements on the Y-axis in the cartesian coordinate system. The type of the multiplexing scheme of the subarray elements is equal to that of the multiplexing scheme falling in the first quadrant (0 degree < alpha)mThe number k of array elements is related to the number k of the array elements, when the number of the array elements is a multiple of 4, the sub-array spacing scheme capable of forming array element multiplexing comprises k types:
Figure BDA0002711866940000092
when the number of array elements is a multiple of 2, but not a multiple of 4, the sub-array spacing scheme that can form array element multiplexing is:
Figure BDA0002711866940000093
wherein k is the sequence number of the kth array element which falls in the first quadrant and is sorted according to a reverse clock;
based on the radius and the number of the array elements of each UCSA in the UCSA array layout, the method can be flexibly adjusted according to communication or detection frequency under the condition of meeting the requirement of S11, and the increase of the aperture and the number of the array elements of the array layout has an improvement effect on improving the estimation precision of partial parameters;
to describe the combinatorial relationship of different USCA, use BqRepresenting different combinations of sub-arrays, by combining individual sub-arrays apThe received data vectors, p 1,2,3, may be combined in different ways to form different joint array manifolds:
Figure BDA0002711866940000101
wherein ,B12 and B13Parameter value for estimating TDID, B23For pairing a plurality of sources;
three subarrays A as shown in FIG. 11、A2 and A3Is defined as:
Figure BDA0002711866940000102
Figure BDA0002711866940000103
Figure BDA0002711866940000104
wherein ,Δr＝ω/c·r，
Figure BDA0002711866940000108
c is the propagation speed of the wave in the current medium, and M belongs to (1-M) is the sequence number of the array element; ([ a)1]m,[a2]m) The same time delay exists between M E (1-M), and the subarray A is in the same theory1And subarray A3Formed M array element pairs ([ a ]1]m,[a3]m) M is the same as (1-M) in time delay;
thus, subarray A1And subarray A2Subarray A1And subarray A3Has the following transformation relation between the array manifold vectors:
Figure BDA0002711866940000105
wherein ,Φ12Is a subarray A2And subarray A1Conventional array manifold rotation operator in between, phi13Is a subarray A3And subarray A1Conventional array manifold rotation operator in between, phi23Is a subarray A2And subarray A3The traditional array manifold rotation operator in between, namely:
Figure BDA0002711866940000106
thirdly, decoupling and analyzing the central parameters and the extended parameters of the two-dimensional distributed source:
when there are K TDID distribution sources in the region of interest, in subarray A1For example, at some certain time, the array manifold vector is
Figure BDA0002711866940000107
All L considered to have a certain distribution property for the Kth TDID distribution sourcekDiscrete paths of a cluster sound source beam. In the response expression, the central parameter and the fluctuation angle are coupled together, in order to avoid high calculation cost of simultaneous estimation of multiple parameters, the central parameter and the fluctuation angle need to be decoupled, and the decoupling process is as follows:
at any one time, the subarray A1Performing Taylor expansion on an array manifold vector of a certain sound source beam in a distributed source by using a central parameter:
Figure BDA0002711866940000111
wherein the array manifold vector
Figure BDA0002711866940000112
At the central angle
Figure BDA0002711866940000113
The first partial derivative of (a) is:
Figure BDA0002711866940000114
Figure BDA0002711866940000115
Figure BDA0002711866940000116
are column vectors formed by the angles of fluctuation, and
Figure BDA0002711866940000117
at the central angle
Figure BDA0002711866940000118
The values of (a) together form a taylor expansion approximation array manifold matrix:
Figure BDA0002711866940000119
wherein, the formed generalized array manifold matrix is only dependent on the central parameter of the Kth TDID distribution source, which is represented by a subscript K0; and three of the column vectors a1，
Figure BDA00027118669400001110
Is linearly independent, by Taylor approximationUnfolding can decouple any sound source beam in a distributed source with randomness into a central parameter with determinism
Figure BDA00027118669400001111
And a fluctuation angle with randomness
Figure BDA00027118669400001112
Thereby laying a foundation for estimating the central parameters and the extended parameters of the two-dimensional distribution source step by step; in addition, when
Figure BDA00027118669400001113
Only its first derivative is used;
subarrays A2The array manifold vector of (a) can be written with a first order taylor approximation as:
Figure BDA00027118669400001114
wherein ,
Figure BDA00027118669400001115
equation (17) can be written as:
a2＝A1(k0)·Θ12(k0)·Γkl (19)，
wherein ,Θ12(k0)For a generalized rotation operator matrix, likewise, sub-matrix A3The array manifold vector of (a) can be written as: a is3＝A1(k0)·Θ13(k0)·Γkl，Θ12(k0) and Θ13(k0)The generalized rotation operator matrix is:
Figure BDA0002711866940000121
from this, subarray A2Generalized array manifold and subarray A1Has a generalized rotation relation in the form of a matrix between the generalized array manifolds, and similarly, the sub-array A2Generalized array manifold and subArray A2Also has a generalized rotation relationship between the generalized array manifolds of (1):
Figure BDA0002711866940000122
importantly, AP(k0)、Θq(k0)Only with central parameters of distributed sources
Figure BDA0002711866940000124
In relation, the relationship theta can be derived from a generalized rotation relationshipq(k0)To estimate the center parameter.
Thirdly, analyzing the covariance matrix structure of the single-distribution-source single-subarray received data;
three separate uniform circular sub-arrays A according to equation (14)1、A2 and A3The received observation data vector is spread by Taylor, in which the sub-array A1The expansion of (a) is:
Figure BDA0002711866940000126
wherein ,X1(t) is subarray A1Taylor expansion of A1(k0)Is a subarray A1Generalized array manifold, noise n1(t) is a Gaussian complex random variable with zero mean value, cyclic symmetry, independence and same distributionk(t) is a 3 x 1 dimensional column vector, s (t) is an incoherent two-dimensional signal source;
similarly, A2 and A3Taylor expansion was performed to obtain:
Figure BDA0002711866940000127
Figure BDA0002711866940000131
subarrays A1The covariance matrix of the received data vector is:
Figure BDA0002711866940000132
wherein ,A1(0)＝[A1(10)...A1(k0)...A1(0)]Is a generalized array manifold matrix containing K TDID sources of dimension M x 3K, LambdaΥExpanding a parameter diagonal matrix for a two-dimensional incoherent distributed source;
and diagonal matrix ΛΥ＝diag[ΛΥ1ΛΥ2ΛΥk] (26)；
According to the assumption that the signals, the path gain and the fluctuation angle are not related to each other
Figure BDA0002711866940000133
The method comprises the following steps:
Figure BDA0002711866940000134
according to formulas (3), (4), (5)
Figure BDA0002711866940000135
The diagonal elements of (a) are:
Figure BDA0002711866940000136
wherein ,
Figure BDA0002711866940000137
is the power of the signal, it can be clearly seen that the matrix
Figure BDA0002711866940000138
The diagonal elements of (A) are randomly distributed variances, i.e. spreading parameters, to which the fluctuation angle is measured
Figure BDA0002711866940000139
Thus can be selected from
Figure BDA00027118669400001310
Middle estimation extension parameters: matrix can be used no matter what type of distribution the two-dimensional incoherent distribution source conforms to
Figure BDA00027118669400001311
Two-dimensional extended parameter of
Figure BDA00027118669400001312
Thereby ensuring the robustness of the algorithm to the distribution type of the distributed source;
Figure BDA0002711866940000141
likewise, subarray A2The covariance matrix of the received data vector can be abbreviated as:
Figure BDA0002711866940000142
wherein ,Θ12(0)＝diag[Φ12(10)，...Φ12(k0)，...Φ12(K0)]The system is a generalized rotation operator diagonal matrix consisting of generalized rotation operators;
based on the formulas (22), (23) and (24), the observation data received by the two sub-arrays are combined to form a combined observation data set
Figure BDA0002711866940000144
Bq(k0)The first order Taylor expansion of (1) is:
Figure BDA0002711866940000145
then the covariance matrix for the joint observation made up of all the subarray combinations can be written as:
Figure BDA0002711866940000146
it can be seen that the joint generalized array manifold matrix
Figure BDA0002711866940000147
Is formed by combining different sub-arrays of K TDID sources.
In formula (32)
Figure BDA0002711866940000148
Of positive formula, thus RqAlso a positive definite matrix, whose eigendecomposition can be written as follows:
Figure BDA0002711866940000149
wherein ,
Figure BDA00027118669400001410
respectively corresponding to the characteristic value of the signal
Figure BDA00027118669400001411
And noise eigenvalue
Figure BDA00027118669400001412
The signal subspace and the noise subspace of (1);
comparing equations (32) and (33) it can be seen that the generalized array manifold matrix R is combinedqIs similar to the signal subspace and has a linear mapping relationship:
Figure BDA0002711866940000151
Figure BDA0002711866940000152
and
Figure BDA0002711866940000153
are all diagonal arrays, therefore
Figure BDA0002711866940000154
and Bq(0)Using full-rank matrix T approximately in the same subspaceqTo describe
Figure BDA0002711866940000155
and Bq(0)The linear relationship between:
Figure BDA0002711866940000156
wherein ,
Figure BDA0002711866940000157
and is
Figure BDA0002711866940000158
It is brought into formula (34) to obtain:
Figure BDA0002711866940000159
wherein ,
Figure BDA00027118669400001510
as a characteristic value of the signal, TqAs a generalized array manifold matrix Bq(0)And the signal feature vector, so that the expansion parameters of the distributed source can be estimated through the signal subspace by utilizing the mapping relation.
Thirdly, analyzing the relation between the two-dimensional distribution source expansion parameters and the energy distribution in the signal subspace:
95% of signal energy of the one-dimensional distribution source is concentrated on a plurality of large eigenvalues of the characteristic subspace, the signal space dimension of the one-dimensional distribution source is directly related to the size of the expansion parameter, the signal energy in the two-dimensional distribution source characteristic subspace is similar to that of the one-dimensional distribution source, and is concentrated on one maximum eigenvalue but distributed in a plurality of maximum eigenvalues no longer like a point sound source signal model;
referring again to fig. 3, it can be seen that when the extended parameters are within 10 °, 95% of the energy of the two-dimensional distributed source is concentrated on the first three large eigenvalues;
with reference to the spatial diffusion curves of the four maximum eigenvalues of the TDID source in a certain gaussian noise shown in fig. 4(a) and (b), the maximum eigenvalue decreases with the increase of the spread parameter of the two-dimensionally distributed source, and the second and third maximum eigenvalues increase, which means that the energy of the signal subspace gradually diffuses from the first eigenvalue to the second and third eigenvalue with the increase of the spread parameter.
In fig. 4(a) and (b), when the SNR is 15dB, all eigenvalues of the noise subspace remain stable when the spread angle increases, but when the spread angle parameter of the distributed source is within (0 ° -2), the noise eigenvalue affects the second and third signal eigenvalues in a noisy environment or in a more noisy environment;
thirdly, analyzing the circular array noise field:
assuming that the noise n (t) is a zero-mean circularly symmetric, independent and identically distributed Gaussian complex random variable, the M × 1-dimensional noise vector received by the array is:
n(t)＝[n1(t),...,nm(t),...,nM(t)]T (37)，
wherein nm(t) is the background noise received from the m-th array element, then the covariance matrix of the noise data received from a single subarray UCSA is:
Figure BDA0002711866940000161
for a spatially uniform noise field, the noise correlation coefficient between any two points in space is:
[ρM]ij＝sinc(2πdij/λ) (39)，
wherein ,dij2rsin (| i-j | pi/M) is the distance between any two points in any single subarray UCSA. The circular array has super-gain characteristic, and for small-aperture array, super-gain processing can obtain much higher spatial gain than conventional processing, so that a UCA-based multi-subarray array arrangement mode can obtain better estimation accuracy;
finally, analyzing an ESPRIT sequential algorithm to solve a single-distribution-source four-parameter joint estimation method;
the traditional ESPRIT algorithm solves the manifold rotation operator of the point source array by using the linear rotation invariance between different sub-arrays, and how to accurately estimate the TDID distribution source parameters by using the ESPRIT algorithm will be shown below.
Investigating subarrays A1And subarray A2The linear transformation of equation (36) can be expressed as:
Figure BDA0002711866940000162
to form a subarray A1 and A2Corresponding to the corresponding signal subspace, the generalized array manifold of (2) can be obtained:
Figure BDA00027118669400001610
Figure BDA0002711866940000163
due to T12Is non-singular and is obtained according to equation (41)
Figure BDA0002711866940000164
It can be brought into formula (42):
Figure BDA0002711866940000165
define a new matrix:
Figure BDA0002711866940000166
equation (43) is rewritten as:
Us1Ψ12＝Us2 (45)，
wherein ,
Figure BDA0002711866940000167
is a subarray A1、A2Combining to obtain a signal subspace rotation operator, and thus obtaining a signal subspace U by decomposing the covariance matrix characteristics of the joint observation datas1And Us2Also has a rotationally invariant relationship Ψ12；
In order to reduce the disturbance of noise to the signal, the least square idea is adopted to solve psi12The unconstrained cost function of (a):
Figure BDA0002711866940000168
obtaining an analytic solution of a signal subspace rotation operator:
Figure BDA0002711866940000169
it can be seen that the central angle of the two-dimensional incoherent distribution source can be estimated from the estimated signal subspace rotation operator. The effective subspace dimension t of the conventional generalized ESPRIT algorithm when there is only one TDID sources(47) Selecting 3 to obtain 3 x 3 generalized subspace rotation operator matrix
Figure BDA0002711866940000171
In fact, in the case of small angular spread, the second and third largest eigenvalues are very small and more sensitive to noise, as shown in fig. 4, and therefore the effective dimension of the signal subspace may be chosen as t s1. When the spread angle is within 10 degrees, the main energy of the signal is concentrated on the maximum eigenvalue, so that it is reasonable and feasible to discard the second third eigenvalue, which can improve the robustness of the estimation. Therefore, only the vector subspace corresponding to the maximum characteristic value in the sequential ESPRIT algorithm provided by the method
Figure BDA00027118669400001710
Is used to estimate the central angle parameter, so the maximum eigenvalue is selected when the central angle estimates the rotation operator:
Φ12＝maxE12 (48)。
example (b):
1. approximate Clarithrome (CRB) analysis
The method estimates four arrival parameters of a single distributed source, and the four parameters jointly form a matrix
Figure BDA0002711866940000172
Defining the method provides that all parameters of the mathematical model are zeta ═ uT,vT]T∈R(5k+1)×1Approximate (finite sample) FIM information matrix J of model parametersζ,ζ∈R(5k+1)×(5k+1)Comprises the following steps:
Figure BDA0002711866940000173
the approximate CRB is then:
Figure BDA0002711866940000174
the approximate CRB gives the limit of the variance matrix of unbiased estimation of all the distributed source parameter vectors u to be estimated;
defining the variance matrix of the estimation error of the method as follows:
Figure BDA0002711866940000175
measuring the estimation error of the mathematical model estimation distribution source parameter vector u provided by the method by using the approximate CRB, wherein the estimation error satisfies the relation;
the central parameters and the extended parameters of a plurality of TDID sources are decoupled by using a first-order Taylor function, and for calculating CRLB conveniently, the signal energy of the TDID distributed sources is further separated from the diagonal matrix of the extended parameters, and then equation (32) can be written as:
Figure BDA0002711866940000176
wherein ,Bq(0)Including the central parameters of the system, including the central parameters,
Figure BDA0002711866940000177
Figure BDA0002711866940000178
2. performance analysis of sequential ESPRIT algorithms on a single TDID source
The performance of the proposed matrix and algorithm was analyzed from several different aspects by monte carlo simulation (800 times) and the algorithm performance was investigated by comparison with CRB. Setting the random scattering path L of each TDID source as 400, wherein a single two-dimensional distribution signal source follows Gaussian distribution, and the central angle of the two-dimensional distribution signal source is
Figure BDA0002711866940000179
The signal-to-noise ratio is 15 dB;
(1) influence of array element number on algorithm estimation performance
First, three uniform circular arrays as shown in fig. 2 are used, the array element radii are λ, and the number of array elements in a single UCA is M-12, where β is 30 °. The distance between every two array elements in the subarray is d1＝0.5λ；
As can be seen in fig. 5, the abscissa increases the 3 UCSA array elements that are not multiplexed from 24 elements to 72 elements to 12 elements, and the spread angle sets two angles: first, large angle
Figure BDA0002711866940000181
Its second, small angle
Figure BDA0002711866940000182
Both of which are labeled lsprad and Sspread, respectively, in fig. 5. It can be seen from fig. 6 that, in the case of large divergence angles, the central parameter isThe precision is hardly increased along with the increase of the number of array elements in the array, and the phenomenon still exists under the high signal-to-noise ratio of 35 dB; this is because in the ESPRIT algorithm, the estimation of the central parameter is only related to the translation invariant relationship between different UCSAs, in the sequential ESPRIT algorithm, only the first eigenvector corresponding to the largest eigenvalue is selected, and as the number of array elements increases, the energy concentrated on the largest eigenvalue gradually decreases, which means that the eigenvector corresponding to the largest eigenvalue has slight distortion, and it is worth noting that the larger the number of array elements, the larger the aperture of the subarray, the higher the estimation precision of the extension angle parameter is, the higher the precision of the extension angle parameter is than that of the central parameter, however, the estimation performance of the central parameter is mainly affected by the size of the extension angle of the TDID source, and is not related to the number of subarray elements;
(2) impact of fast sample beat number on algorithm estimation performance
Based on the above results, 12 isotropic array elements were arranged in each UCSA, and each UCSA had a radius λ, where α was 30 °, and d wasx＝dyIn the structure, no array element multiplexing exists, 36 isotropic array elements are arranged, and the SNR is 5 dB. Two spread angle size scenarios were compared in conjunction with fig. 6: large spread angle
Figure BDA0002711866940000183
Labeled Lspread in FIG. 6, small spread angle
Figure BDA0002711866940000184
Which is labeled Sspread in fig. 6; as shown in fig. 6, as the number of sample snapshots increases, the improvement of the central parameter estimation accuracy is significantly better than the improvement of the angular spread parameter estimation accuracy; in addition, in the central parameter estimation, the small spread angle is more sensitive to the increase of the sample snapshot, and the parameter estimation precision is improved along with the increase of the samples; and when T is 400, good estimation accuracy is obtained.
3. Performance analysis of sequential ESPRIT algorithms on multiple TDID sources
(1) Effect of signal-to-noise ratio on algorithm estimation performance
Using the same simulation conditions as above, the root mean square error and Cramer-Lo lower bound (CRLB) of the parameter estimation using the generalized ESPRIT algorithm for the UCSA array configuration and URA (uniform rectangular array) were compared only when the SNR was increased from-15 dB to 5dB, as shown with reference to FIG. 7: for the central parameter estimation, the UCSA using only 36 array elements is more stable than the URA array using 100 array elements. The URA array geometry for 100 array elements almost doubles the array aperture compared to the UCSA array configuration. The performance of the central angle estimation by the sequential ESPRIT algorithm cannot be improved with the increase of the signal-to-noise ratio, because the error generated by the first-order taylor approximation of the UCSA array manifold constitutes the main component of the estimation error at high signal-to-noise ratio. The degradation in estimation performance at low signal-to-noise ratios is mainly caused by noise, while at high signal-to-noise ratios and the taylor approximation error of the array manifold does not decrease with increasing signal-to-noise ratio. Even so, under high signal-to-noise ratio, the Root Mean Square Error (RMSE) of estimation reaches below 0.2 °, and good estimation accuracy is maintained. In addition, for the spread angle parameter, the estimation accuracy of the algorithm herein is better than that of the generalized ESPRIT algorithm employing the URA array.
(2) Effect of spread angle size on algorithm estimation performance
The central angle of a two-dimensional distributed signal source is fixed, the spread angles of two distributed sources are changed from 0.1 degrees to 3.1 degrees at the same time, each time, the spread angles are increased by 0.2 degrees, and from the simulation result of the attached figure 8, the sequential ESPRIT algorithm provided by the invention can obtain stable estimation precision under the two conditions of small spread angles and large spread angles, particularly for smaller spread angles, the estimation performance of the central parameter is closer to CRLB, and with the increase of the spread angles, the approximation error is increased, and the estimation performance of the central angle is reduced. When at a small spread angle σθLess than or equal to 1.5 degrees or
Figure BDA0002711866940000191
The RMSE difference of the two algorithms (the present algorithm and the generalized ESPRIT algorithm) is 10 times different, and due to the precise choice of the signal subspace dimension, the algorithm proposed herein can achieve better performance; simulation results show that the sequential ESPRIT algorithm considers a large expansion angle and a small expansion angle at the same time, and improves the robustness of the TDID source expansion value
(3) Analysis of influence of two TDID distribution source space angle separation distance changes on parameter estimation performance
In the present simulation, referring to fig. 9, the second TDID gaussian source starts from an angular distance interval of 4 ° in both azimuth and elevation dimensions, and is gradually increased to an angular distance interval of 18 ° in steps of 2 °. Even if the two TDID spaces are separated by only 4 degrees, the sequential ESPRIT algorithm adopting the UCSA array configuration can show higher estimation precision. With the spatial angle separation of the two TDID sources being more and more distant, the estimation precision is characterized by first descending and then ascending.
(4) Energy difference of two TDID distribution sources has performance influence
Referring to fig. 10, the present simulation will investigate the parameter estimation performance of the algorithm as the difference in power of the two source signals varies for two sources of TDID distribution with the same gaussian distribution. The configuration is the same as the above configuration, the power of the first TDID source is fixed to 1, and the SNR is set to 10 dB; and the power of the second TDID source is gradually increased from 1 to 10 by 1 increments, namely the environmental signal-to-noise ratio of the second TDID source is gradually increased to 30 dB. As can be seen from fig. 11, the sequential ESPRIT algorithm provided by the present invention has better robustness than the conventional ESPRIT algorithm, because the second eigenvalue of the high-power TDID distribution source is large enough to affect the estimation accuracy of the eigenvector of the low-power TDID distribution source; this is why the estimated performance of generalized ESPRIT using URA arrays drops dramatically when the TDID source powers are not equal. The power difference of the two TDID signal sources affects the power distribution of the signal subspace, in particular for obtaining two relatively small eigenvalues of the angular spread parameter. In addition, the generalized ESPRIT algorithm selects three traditional rotation operators to calculate the central parameters at the same time, so that mutual interference of signals is introduced, and poor evaluation performance is obtained.
(5) Robust performance evaluation of distribution morphology of different distributed sources
The two sources of TDID distribution are uniformly and Gaussian spatially distributed, respectively, and the third source is a two-dimensional point source. Spatial range of distributed arrival angles around central angle for spatially uniformly distributed TDID sourcesEnclose as
Figure BDA0002711866940000192
For spatially distributed TDID sources, over 95% of the discrete arrival angles are concentrated
Figure BDA0002711866940000193
Within the range of (1). When there are more TDID signal sources in the region of interest, more array elements and larger aperture can be used to obtain better estimation accuracy, so the radius of each USCA is increased to 2 lambda in the simulation, and the number of array elements of each USCA is 20.
Figure BDA0002711866940000194
The diagonal elements of (2) are only related to the standard deviation of the fluctuation angle distribution without any prior information about the spatial distribution type, so that the extended parameter estimation performance has stronger robustness to different distribution types compared with an early integral model; as can be seen from fig. 11, the estimation performance of the central parameter and the extended parameter can maintain better accuracy for both the uniform distribution and the gaussian distribution. By comparing the estimation performance of the spread angle parameter, it can be seen that the sequential ESPRIT algorithm is degraded when the signal-to-noise ratio is low, and better results can be obtained as the signal-to-noise ratio increases. And for traditional point sources without distribution characteristics, sequential ESPRIT can obtain higher estimation precision, so that the robustness of the algorithm to the distribution source and the point source is proved. Different angular spatial distribution patterns for two-dimensional uncorrelated distributed sources are also robust.
The foregoing shows and describes the general principles, essential features, and advantages of the invention. It will be understood by those skilled in the art that the present invention is not limited to the embodiments described above, which are described in the specification and illustrated only to illustrate the principle of the present invention, but that various changes and modifications may be made therein without departing from the spirit and scope of the present invention, which fall within the scope of the invention as claimed. The scope of the invention is defined by the appended claims and equivalents thereof.

Claims (8)

1. The method for estimating the parameters of the sequential ESPRIT two-dimensional incoherent distribution source based on the combined UCA array is characterized by comprising the following steps of:
s1: constructing a multi-subarray array structure model based on uniform circular array UCA to obtain a combinable three-uniform circular subarray structure;
s2: combining the three subarrays UCSA pairwise to form three pairs of subarray UCSA combinations with equal intervals, and constructing received data of three combinations of uniform circular subarrays UCSA to obtain observation data;
s3: based on the received observation data, a Taylor approximate covariance matrix model R of a two-dimensional incoherent distribution source is constructed12、R13、R23For decoupling the central and extended parameters of the two-dimensional uncorrelated distributed sources;
s4: calculating covariance matrix of three sample data through joint data received by three groups of subarrays UCSA
Figure FDA0002711866930000011
And
Figure FDA0002711866930000012
obtaining an estimated value of the central parameter through a linear mapping relation between a covariance matrix space and a feature subspace;
s5: using Taylor approximation covariance model R12、R13、R23And obtaining an estimated value of the expansion parameter through a linear mapping relation between the estimated value and the signal characteristic value.
2. The method for estimating parameters of a sequential ESPRIT two-dimensional incoherent distributed source based on a joint UCA array according to claim 1, wherein the specific operation steps of step S1 are as follows:
s11: establishing three uniform circular sub-arrays A with the same radius1、A2、A3The array types and array parameters of the three sub-arrays are the same, and the number of the array elements is M;
the radius r of each subarray is expressed as: r ═ γ · λ, where λ is the wavelength and γ represents a multiple of the wavelength;
the included angle alpha between adjacent array elements is as follows: α is 2 π/M;
array element spacing d1Comprises the following steps: d1＝2rsin(π/M)；
The distance between every two adjacent array elements is half wavelength, and the ratio of the radius of the array elements to the wavelength is
Figure FDA0002711866930000021
S12: setting a subarray A1The reference array element is an array element with the X-axis coordinate of (0, r), and the sub-array A1The 2 nd to M array elements are arranged anticlockwise in sequence, wherein the included angle between the first array element passing through the Y axis and the Y axis is beta degrees;
s13: at edge subarray A1Forward translation distance d of reference array element on X axisxAt, lay and A1Subarrays A of the same structure2；
S14: at edge subarray A1The first array element passing through the Y axis, namely the first array element is translated upwards by a distance d obliquely from the extension line of an included angle beta DEG of the Y axisyAt, lay and A1Subarrays A of the same structure3；
S15: according to the number M of the sub-array elements, the radius R of the sub-array and the distance d between the two pairs of sub-arraysx、dyJudging the number of types of the multiplexing schemes of the subarray elements among the subarrays;
s16: based on the radius and the number of array elements of each UCSA in the UCSA array layout, the flexible adjustment can be carried out according to communication or detection frequency under the condition of meeting the requirement of S11;
s17: and finally obtaining three pairs of UCSA array combinations with equal spacing.
3. The method according to claim 2, wherein the step S15 comprises the following steps:
s151: calculating the number k of array elements falling in a first quadrant;
s152: when k is a multiple of 4, there are k sub-array spacing schemes that can form array element multiplexing:
Figure FDA0002711866930000022
s153: when k is a multiple of 2 instead of a multiple of 4, the sub-array spacing scheme that can form array element multiplexing is:
Figure FDA0002711866930000023
4. the method according to claim 1, wherein the step S2 is a structural model of a UCSA array pair composed of three uniform circular sub-arrays, and wherein the method comprises:
s21: subarrays A1Neutral subarray A2M groups of array element pairs formed between corresponding array elements
Figure FDA0002711866930000031
And subarrays A1And subarray A3Formed M groups of array element pairs
Figure FDA0002711866930000032
Have the same time delay therebetween, and the subarray A1And subarray A2Subarray A1And subarray A3Subarray A2And subarray A3Satisfies the following transformation relation between array manifold vectors:
Figure FDA0002711866930000033
wherein ,
Figure FDA0002711866930000034
is a subarray A2And subarray A1Flow-shaping rotation of traditional array in betweenThe number of the operators is converted into a number of operators,
Figure FDA0002711866930000035
is a subarray A3And subarray A1The traditional array manifold rotation operator in between, and the expression of the traditional array manifold rotation operator is:
Figure FDA0002711866930000036
wherein ,
Figure FDA0002711866930000037
is a random angle of arrival.
5. The method of claim 1, wherein the step S3 is performed by constructing taylor-approximated covariance matrix R of two-dimensional incoherent distribution source12、R13、R23The method comprises the following operation steps:
s31: three separate uniform circular sub-arrays A1、A2 and A3The received observation data vector is spread by Taylor, in which the sub-array A1The expansion of (a) is:
Figure FDA0002711866930000038
wherein ,X1(t) is subarray A1Taylor expansion of A1(k0)Is a subarray A1Generalized array manifold, noise n1(t) is a Gaussian complex random variable with zero mean value, cyclic symmetry, independence and same distributionk(t) is a 3X 1 matrix, sk(t) is the kth transmitted signal;
likewise, subarray A2 and A3The received observation data vectors are written in a Taylor expanded form as:
Figure FDA0002711866930000041
Figure FDA0002711866930000042
s32: the combined generalized array manifold consisting of three sets of UCSAs formed by three UCSAs is:
Figure FDA0002711866930000043
Figure FDA0002711866930000044
Figure FDA0002711866930000045
wherein ,B12(0)Is a subarray A1And A2Combined generalized array manifold of (B)13(0)Is a subarray A1And A3Combined generalized array manifold of (B)23(0)Is a subarray A2And A3A joint generalized array manifold of (1);
s33: and (3) calculating a covariance matrix of the received data vector of the joint subarray consisting of the three groups of UCSA, wherein the calculation formulas are respectively as follows:
associative subarrays A1 and A2The calculation formula of the covariance matrix of the received data vector is:
Figure FDA0002711866930000046
wherein ,
Figure FDA0002711866930000047
is null of a circular array noise fieldThe inter-correlation matrix, H represents the sign of the conjugate transpose operation, ΛΥExpanding a parameter diagonal matrix for a two-dimensional incoherent distributed source:
Figure FDA0002711866930000051
associative subarrays A1 and A3The calculation formula of the covariance matrix of the received data vector is:
Figure FDA0002711866930000052
associative subarrays A2 and A3The calculation formula of the covariance matrix of the received data vector is:
Figure FDA0002711866930000053
wherein the matrix Λ can be used no matter what type of distribution the two-dimensional incoherent distribution source conforms toΥTwo-dimensional extended parameter of
Figure FDA0002711866930000054
The distribution characteristics are expressed, and the robustness of the algorithm on the distribution type of the distributed source is ensured;
s34: the central parameter for the rotation operator of the formula (4)
Figure FDA0002711866930000055
The deterministic function is expressed as:
Figure FDA0002711866930000056
6. the method according to claim 1, wherein the step S4 comprises the following steps:
s41: covariance matrix of three sample data
Figure FDA0002711866930000057
And
Figure FDA0002711866930000058
performing characteristic decomposition, and respectively taking K eigenvectors corresponding to the first K maximum eigenvalues
Figure FDA0002711866930000059
And
Figure FDA00027118669300000510
obtaining a signal subspace after characteristic decomposition;
s42: separately computing signal subspace rotation operators
Figure FDA00027118669300000511
And
Figure FDA00027118669300000512
the evaluation value of (c) is calculated by the formula:
Figure FDA00027118669300000513
Figure FDA00027118669300000514
Figure FDA00027118669300000515
wherein ,
Figure FDA0002711866930000061
and
Figure FDA0002711866930000062
are respectively as
Figure FDA0002711866930000063
And
Figure FDA0002711866930000064
decomposing to obtain a signal subspace;
s43: separate to signal subspace rotation operators
Figure FDA0002711866930000065
And
Figure FDA0002711866930000066
performing characteristic decomposition to obtain a generalized array manifold rotation operator
Figure FDA0002711866930000067
And
Figure FDA0002711866930000068
s44: obtaining matched K groups of generalized array manifold rotation operators by using a rotation operator matching algorithm, and solving estimated values of K incoherent two-dimensional distribution source central parameters by using the K groups of generalized array manifold rotation operators
Figure FDA0002711866930000069
The closed-form solution of (a) obtains an estimated value of the central parameter, and the formula for calculating the closed-form solution is:
Figure FDA00027118669300000610
wherein ,dx and dyRespectively, the distance between the sub-arrays is,
Figure FDA00027118669300000611
where λ is the wavelength, c is the speed of sound in the propagation medium, and ω is the dominant frequency.
7. The method according to claim 6, wherein the rotation operator matching algorithm of step S44 comprises the following steps:
s441: subarrays A1And subarray A2Form a set of rotation operators between the array manifolds
Figure FDA00027118669300000612
Wherein takes out phii,12,i＝1；
Subarrays A1And subarray A3Form a set of rotation operators between the array manifolds
Figure FDA00027118669300000613
Each element composition has K combination schemes:
Figure FDA00027118669300000614
and, the formula for calculating the cost function of the K schemes is:
Figure FDA0002711866930000071
the combination mode corresponding to the minimum value of the cost function obtained by calculation is the correct matching mode
Figure FDA0002711866930000072
S442 at
Figure FDA0002711866930000073
Set of rotation operators
Figure FDA0002711866930000074
And
Figure FDA0002711866930000075
set of rotation operators
Figure FDA0002711866930000076
Removing successfully matched combined elements
Figure FDA0002711866930000077
S443, repeating the steps S441-S442 until all the obtained
Figure FDA0002711866930000078
The correct pairing scheme.
8. The method according to claim 6, wherein the step S5 is performed by using a sequential ESPRIT two-dimensional incoherent distribution source parameter estimation method based on a joint UCA array: two-dimensional uncorrelated distributed source Taylor approximate covariance model R provided by utilizing method12、R13、R23And obtaining an estimated value of the extension parameter through a linear mapping relation between the estimated value and the signal characteristic value, wherein the specific operation steps comprise:
s51: the estimated value of the center parameter obtained in step S44
Figure FDA0002711866930000079
In the formula (8), the estimated values of the combined generalized array manifold are obtained respectively
Figure FDA00027118669300000710
And
Figure FDA00027118669300000711
taking any one of the joint generalized array manifold estimates, e.g.
Figure FDA00027118669300000712
Then a joint subarray A can be obtained12The linear transformation relationship between the generalized array manifold and the signal subspace is as follows:
Figure FDA00027118669300000713
s52: by joint subarrays A12For example, and A13 and A23Similarly, the covariance of the samples of the received data
Figure FDA00027118669300000714
Performing characteristic decomposition and analogizing to obtain
Figure FDA00027118669300000715
And
Figure FDA00027118669300000716
taking the first 3K maximum eigenvalues formed by each signal source to form a diagonal matrix
Figure FDA00027118669300000717
In relation to linear transformation
Figure FDA00027118669300000718
Solving to obtain an extended parameter matrix:
Figure FDA00027118669300000719
s53: obtaining estimated value of two-dimensional distribution source expansion parameter by the above formula
Figure FDA00027118669300000720
Extended parameter closed-form solution:
Figure FDA0002711866930000081
s54: and obtaining an expansion parameter estimation value according to the obtained expansion parameter closed-form solution.
CN202011059330.XA 2020-09-30 2020-09-30 Sequential ESPRIT two-dimensional incoherent distribution source parameter estimation method based on united UCA Active CN112255629B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011059330.XA CN112255629B (en) 2020-09-30 2020-09-30 Sequential ESPRIT two-dimensional incoherent distribution source parameter estimation method based on united UCA

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011059330.XA CN112255629B (en) 2020-09-30 2020-09-30 Sequential ESPRIT two-dimensional incoherent distribution source parameter estimation method based on united UCA

Publications (2)

Publication Number Publication Date
CN112255629A true CN112255629A (en) 2021-01-22
CN112255629B CN112255629B (en) 2023-06-02

Family

ID=74233507

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011059330.XA Active CN112255629B (en) 2020-09-30 2020-09-30 Sequential ESPRIT two-dimensional incoherent distribution source parameter estimation method based on united UCA

Country Status (1)

Country Link
CN (1) CN112255629B (en)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113111304A (en) * 2021-04-01 2021-07-13 哈尔滨工程大学 Coherent distribution source direction finding method based on quantum ray mechanism under strong impact noise
CN113468844A (en) * 2021-06-17 2021-10-01 浙江大学 Coupled array beam comprehensive analysis method
CN115407270A (en) * 2022-08-19 2022-11-29 苏州清听声学科技有限公司 Sound source positioning method of distributed array

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4965732A (en) * 1985-11-06 1990-10-23 The Board Of Trustees Of The Leland Stanford Junior University Methods and arrangements for signal reception and parameter estimation
US5262789A (en) * 1992-04-30 1993-11-16 General Electric Company Source identification system for closely separated spatial sources
US20070091723A1 (en) * 2005-07-22 2007-04-26 Institute Of Acoustics, Chinese Academy Of Sciences Method of signal processing for high resolution bathymetric sidescan sonar
CN104931923A (en) * 2015-04-02 2015-09-23 刘松 Grid iterative estimation of signal parameters via rotational invariance techniques (ESPRIT), namely, extensible rapid estimation algorithm capable of being used for uniform circular array 2-dimensional direction of arrival (2D DOA)
CN111123192A (en) * 2019-11-29 2020-05-08 湖北工业大学 Two-dimensional DOA positioning method based on circular array and virtual extension
CN111711500A (en) * 2020-05-06 2020-09-25 中国人民解放军63892部队 Simulation antenna array calibration and radio frequency signal monitoring system

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4965732A (en) * 1985-11-06 1990-10-23 The Board Of Trustees Of The Leland Stanford Junior University Methods and arrangements for signal reception and parameter estimation
US5262789A (en) * 1992-04-30 1993-11-16 General Electric Company Source identification system for closely separated spatial sources
US20070091723A1 (en) * 2005-07-22 2007-04-26 Institute Of Acoustics, Chinese Academy Of Sciences Method of signal processing for high resolution bathymetric sidescan sonar
CN104931923A (en) * 2015-04-02 2015-09-23 刘松 Grid iterative estimation of signal parameters via rotational invariance techniques (ESPRIT), namely, extensible rapid estimation algorithm capable of being used for uniform circular array 2-dimensional direction of arrival (2D DOA)
CN111123192A (en) * 2019-11-29 2020-05-08 湖北工业大学 Two-dimensional DOA positioning method based on circular array and virtual extension
CN111711500A (en) * 2020-05-06 2020-09-25 中国人民解放军63892部队 Simulation antenna array calibration and radio frequency signal monitoring system

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
崔维嘉;代正亮;王大鸣;李祥志;: "一种自动匹配的分布式非圆信号二维DOA快速估计方法" *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113111304A (en) * 2021-04-01 2021-07-13 哈尔滨工程大学 Coherent distribution source direction finding method based on quantum ray mechanism under strong impact noise
CN113111304B (en) * 2021-04-01 2022-09-27 哈尔滨工程大学 Coherent distribution source direction finding method based on quantum ray mechanism under strong impact noise
CN113468844A (en) * 2021-06-17 2021-10-01 浙江大学 Coupled array beam comprehensive analysis method
CN113468844B (en) * 2021-06-17 2023-08-04 浙江大学 Analysis method for coupling array wave beam synthesis
CN115407270A (en) * 2022-08-19 2022-11-29 苏州清听声学科技有限公司 Sound source positioning method of distributed array
CN115407270B (en) * 2022-08-19 2023-11-17 苏州清听声学科技有限公司 Sound source positioning method of distributed array

Also Published As

Publication number Publication date
CN112255629B (en) 2023-06-02

Similar Documents

Publication Publication Date Title
Hu et al. An ESPRIT-based approach for 2-D localization of incoherently distributed sources in massive MIMO systems
CN112255629A (en) Sequential ESPRIT two-dimensional incoherent distribution source parameter estimation method based on combined UCA array
CN106054123A (en) Sparse L-shaped array and two-dimensional DOA estimation method thereof
CN106526530A (en) Propagation operator-based 2-L type array two-dimensional DOA estimation algorithm
CN107919535B (en) three-dimensional array antenna based on directional double circular arrays and construction method thereof
Akbari et al. MUSIC and MVDR DOA estimation algorithms with higher resolution and accuracy
CN107703478B (en) Extended aperture two-dimensional DOA estimation method based on cross-correlation matrix
CN106526531A (en) Improved propagation operator two-dimensional DOA estimation algorithm based on three-dimensional antenna array
CN112379327A (en) Two-dimensional DOA estimation and cross coupling correction method based on rank loss estimation
CN109557504B (en) Method for positioning near-field narrow-band signal source
CN112929962B (en) Positioning method, positioning device, computer equipment and storage medium
CN108872930B (en) Extended aperture two-dimensional joint diagonalization DOA estimation method
CN107290732A (en) A kind of single base MIMO radar direction-finding method of quantum huge explosion
Zhang et al. An RCB-like steering vector estimation method based on interference matrix reduction
CN111352063B (en) Two-dimensional direction finding estimation method based on polynomial root finding in uniform area array
CN113189592A (en) Vehicle-mounted millimeter wave MIMO radar angle measurement method considering amplitude mutual coupling error
CN113835063B (en) Unmanned aerial vehicle array amplitude and phase error and signal DOA joint estimation method
Bai et al. Association of DOA estimation from two ULAs
CN113671439A (en) Unmanned aerial vehicle cluster direction finding system and method based on non-uniform intelligent super-surface array
Mondal Studies of different direction of arrival (DOA) estimation algorithm for smart antenna in wireless communication
CN116224219A (en) Array error self-correction atomic norm minimization DOA estimation method
CN112946615B (en) Phased array system amplitude and phase error correction method
WO2022166477A1 (en) Positioning method and apparatus, base station, computer device, and storage medium
Gu et al. A sequential ESPRIT algorithm based on a novel UCSA configuration for parametric estimation of two-dimensional incoherently distributed source
Ahmed et al. Simulation of Direction of Arrival Using MUSIC Algorithm and Beamforming using Variable Step Size LMS Algorithm

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant