Under-sampled Wideband Array Signal Processing Method Based on Improved Hadamard Matrix

Abstract

The present invention discloses an undersampled wideband array signal processing method based on an improved Hadamard matrix. The method includes establishing an array receiving narrowband signal model and a spatial domain compressive sensing observation model under the array, and proposing a compressive zeroing observation matrix based on sequential partial Hadamard; using the spatial domain sparse adaptive OMP algorithm and the spatial domain sparsity of the array received signal to jointly reconstruct and estimate the arrival direction of the underwater target signal; and extending this method to process underwater wideband signals, and proposing an undersampled wideband array signal processing method based on an improved Hadamard matrix with engineering practical significance. This method can be used as a means for active sonar to reduce the number of array elements in actual engineering when using LFM signals to detect targets, and can be used in fields such as underwater target azimuth capture and tracking.

Description

Technical Field

[0001] The present invention belongs to the technical field of underwater acoustic array signal processing, and relates to a signal processing method for array undersampling, in particular to an undersampling broadband array signal processing method based on an improved Hadamard matrix. Background Art

[0002] Underwater acoustic array signal processing technology involves many fields such as military and national economy, and its applications have been quite mature and extensive. Underwater acoustic array signal processing technology refers to the theory and technology of placing multiple sonar sensors at different spatial positions in the underwater acoustic environment to form a sonar sensor array, inducing and receiving underwater acoustic signals by the sonar sensor array, and then further processing the received signals to recover useful information therefrom.

[0003] For the development of some existing linear arrays with a relatively large total length of underwater sonars, there are other problems such as the huge amount of data stored in the received data and the relatively large impact on the overall engineering test caused by data reception errors of some array elements. For example, a towed linear array sonar is composed of multiple, dozens or even hundreds of hydrophones arranged, and a large amount of data needs to be received and processed through a large number of hardware devices. This not only requires a high hardware cost, but also has many inconveniences in software processing, and the processing process is relatively cumbersome.

[0004] At present, most of the signal processing of compressive sampling arrays (CSA) based on the theory of compressive sensing (CS) uses random measurement matrices, which are difficult to implement in hardware and have no engineering practicality. How to apply the theory of compressive sensing to the signal processing of compressive sampling arrays is a relatively new research direction and a key technology to be broken through in the current underwater acoustic field. Summary of the Invention

[0005] In view of the deficiencies of the prior art, the present invention proposes an undersampling broadband array signal processing method based on an improved Hadamard matrix, mainly extending the theory of compressive sensing to spatial domain array compression. Based on the arrangement of non-uniform linear arrays in engineering, by constructing a deterministic measurement matrix, the compression relationship between the actual engineering linear array and the uniform linear array is established. Then, using the spatial sparsity of the array and through relevant signal reconstruction algorithms, the information of the target is reconstructed and restored from these small amounts of non-uniform linear array data, achieving an effect similar to that of the original uniform linear array. In engineering practice, it has advantages such as greatly reducing the number of linear arrays and the transportation volume.

[0006] The undersampling broadband array signal processing method based on an improved Hadamard matrix specifically includes the following steps:

[0007] Step 1: Establish an array spatial domain signal reception model

[0008] A Uniform Vertical Linear Array (UVLA) consists of N uniformly distributed array elements and is used to receive far-field narrowband target signals from K directions {θ1, θ2, …, θ K}. The vector equation is as follows:

[0009] X(t) = A0·s0(t) + v(t) (1)

[0010] where is the scalar data received by the uniform vertical linear array; t = 1, 2, …, L represents the discrete sampling time; is the amplitude matrix of the underwater echo signal; is the scalar data of the noise received by the uniform vertical linear array; is the K-array vector manifold matrix, and A0 = [a(θ1), a(θ2), …, a(θ K )]. represents a real matrix, represents a complex matrix.

[0011] In actual situations, the underwater target echo signal received by the uniform vertical linear array can be regarded as sparse in the spatial domain. Therefore, an array manifold A can be constructed according to all possible potential signal angle information in the underwater space to establish a spatial domain sparse model.

[0012] After grid division, the K-array vector manifold matrix A0 can be extended to A = [a(θ1), a(θ2), …, a(θ K ), …, a(θ P )], where K << P. The spatial domain sparse model shown in formula (1) is rewritten as:

[0013] X(t) = A·S(t) + v(t) (2)

[0014] is the amplitude information of the K-sparse underwater echo signal; is the array vector manifold matrix under grid division, and the p-th column vector is expressed as:

[0015]

[0016] where c is the propagation speed of sound waves underwater, f0 is the frequency of the underwater echo signal received by the array, d is the element spacing, and e is the natural constant.

[0017] The array manifold matrix A after grid division is equivalent to an overcomplete dictionary, which contains the angle information of all potential signals in the divided spatial domain. Therefore, if there is an actual target echo signal θ at a(θ p )k , then S p = s k , otherwise S p = 0.

[0018] Step 2: Construct a deterministic measurement matrix for the compressible sampling array

[0019] s2.1. Construct the deterministic measurement matrix Φ

[0020] Use the measurement matrix Φ to perform spatial domain compressive sampling on the array. After compressive sampling observation of the N - element linear array data X, obtain the M - element non - uniform linear array data Y:

[0021] Y = ΦX (4)

[0022] Wherein, is a measurement matrix that satisfies the Restricted Isometry Property (RIP), and M << N, represents the received data of the non - uniform linear array after compressive sampling by the measurement matrix Φ. Substitute formula (2) into formula (4) to get:

[0023] Y = Φ(SA + v) = ΦAS + Φv (5)

[0024] Substitute the conversion relationship shown in formula (6) into formula (5), and we can get:

[0025]

[0026] Y = ΦAS + V (7)

[0027] Y M×L = B M×P S P×L + V M×L (8)

[0028] Wherein, B is the vector manifold of the compressive sampling array; V is the noise scalar data received by the underwater array. S is the amplitude information of the underwater target signal that we ultimately want to obtain. Solving for S can be regarded as an NP - hard non - convex optimization problem, but it is difficult to solve in polynomial operations. Therefore, it can be simplified to a minimum l1 - norm problem, and S is reconstructed and recovered through formula (9):

[0029]

[0030] The specific reconstruction and recovery process is carried out in Step 3.

[0031] There are many construction methods for the measurement matrix Φ, but all of them need to satisfy the Restricted Isometry Property (RIP) as a prerequisite. This criterion gives the equivalent conditions between the l1 norm and the l0 norm, that is, assuming that the signal conforms to a certain specific sparsity condition, then the original data can be accurately or approximately recovered with high probability from a small amount of sampled data.

[0032] The deterministic measurement matrix Φ constructed by this method is a compressed zeroing measurement matrix based on the sequential partial Hadamard. The elements of the Hadamard matrix are 1 or -1, and its construction process is as follows:

[0033] H1 = [1] (10)

[0034]

[0035]

[0036]

[0037] It can be seen from the above construction steps that this construction method first generates a Hadamard orthogonal matrix of size U·U, where U = 2 γ , γ = 1, 2, … ∞. Based on the above Hadamard matrix, first intercept the first N column submatrix from it to obtain a partial Hadamard measurement matrix with lower coherence and good partial orthogonality. Then, continuously select the first M row vectors in sequence to form a new submatrix. Finally, use the array compression zeroing method in s2.2 to zero the elements φ n at the positions where there are actually no array elements in the new submatrix, and obtain the final measurement matrix Φ. Therefore, we call this measurement matrix the compressed zeroing measurement matrix based on the sequential partial Hadamard. The simulation results show that this construction method can achieve a better reconstruction effect and has practical significance in engineering.

[0038] s2.2. Array Compression Zeroing Method

[0039] The matrix expansion form of the array compression sampling process shown in formula (4) under single snapshot data is:

[0040]

[0041] For convenience, the single snapshot compressed sampling data y is denoted as [y1, y2, …, y M T , x = [x1, x2, …, x N T .

[0042] ​​Here, we introduce the concept of sparse decomposition, which is to select as few column vectors as possible from the redundant dictionary so that their linear combination is equal to or approximately equal to y. The idea of sparse decomposition is introduced into compressed sampling, that is, the compressed sampling data y under a single snapshot is obtained by selecting M column vectors from the measurement matrix Φ and linearly combining them with the uniform linear array data:

[0043]

[0044] Among them, φ n =[φ 1n ,φ 1n ,…,φ Mn ] T

[0045] Formula (15) can be considered as the conversion relationship between the compressed sampling data y and the echo signal data x. The compressed sampling data y is actually the sampling data of the non-uniform linear array, that is, the data of several array elements arranged relative to the array elements of the uniform linear array. Therefore, the data at the corresponding positions of the non-existent array elements can be set to zero. In other words, the process of obtaining the M-element compressed data Y actually corresponds to discarding the corresponding element array data in the original data, that is, the n-th array data x n The corresponding φ n Set all elements in to zero.

[0046] The measurement matrix Φ is processed by the array compression and zeroing method, and the measurement matrix Φ after processing still satisfies the constrained isometry criterion (RIP) and column non-correlation.

[0047] Step 3: Underwater target detection under narrowband signals

[0048] This method mainly uses the spatial sparse adaptive OMP reconstruction algorithm to solve the reconstruction and recovery problem of S. The specific reconstruction process is as follows:

[0049] a) Initialize the residual r0=y, the index set Λ= , Ι= , the initial value of the number of iterations is z = 1;

[0050] b) Calculate index value u = arg max i | <r,b i >|, and store the index value u in I;

[0051] c) Update support set B Λ =Φ Λ ∪{b u} and index set Λ=Λ∪Ι;

[0052] d) Calculation Update residual

[0053] e) Until the iteration termination condition is met, output as the reconstruction estimate; otherwise, iterate z+1 and jump to step b.

[0054] The iteration termination condition is:

[0055]

[0056] Where: ‖r z ‖2 is the residual value of the iteration number z; Γ z is the sparsity of the z-th iteration.

[0057] When formula (16) is satisfied, it means that the residual decreases slowly after the z-1th iteration and an “inflection point” appears, then Γ z-1 This is the true spatial signal sparsity, or the number of signals from different incoming directions. Simply output the z-1th estimate. After all atoms corresponding to nonzero elements in the sparse signal are included, the remaining atoms are those with low correlation with the array element observation data. If the iteration process continues, only atoms with low correlation will be selected, and the residual energy will change very little. Therefore, this property can be used to determine whether the iteration has achieved the correct sparsity.

[0058] Step 4: Underwater target detection under broadband signals

[0059] The methods described in steps 1 to 3 are suitable for underwater target detection where the incoming wave signal is a far-field narrowband. However, in active sonar detection, facing the increasingly complex signal transmission environment, broadband signals can obtain more information. Therefore, they are extended to underwater target detection using broadband signals. The entire broadband signal received by each element of the array can be divided into several non-overlapping narrowband data according to the frequency band. Then, each of the divided narrowband data is processed using the above method. Finally, the processing results of each frequency band are averaged and integrated to obtain the final azimuth estimation result. The specific method is as follows:

[0060] The broadband signal X(t) received by the uniform linear array is regarded as the superposition of J narrowband signals:

[0061]

[0062] The L snapshot received time domain signal is divided into Q non-overlapping segments according to the frequency band. Then the DFT transform is performed on each segment to obtain Q groups of mutually uncorrelated narrowband frequency domain components. In broadband processing, Q is called the frequency domain snapshot, and L / Q is generally used as an integer value for segmentation. Therefore, the array signal model under a single snapshot in the frequency domain is:

[0063]

[0064] in,

[0065] The deterministic measurement matrix Φ constructed in step 2 is introduced to form a broadband array spatial compression observation model. The compressed data under the qth snapshot is expressed as:

[0066]

[0067] At this point, the wideband signal has been segmented into Q snapshots of narrowband data in the frequency domain. The narrowband reconstruction method in step 3 can be used to reconstruct the Q snapshots and estimate the target's direction of arrival. Finally, the processing results of the Q frequency bands are combined by power averaging as follows to obtain the final direction estimate P(θ):

[0068]

[0069] The present invention has the following beneficial effects:

[0070] The present invention establishes a connection between a non-uniform linear array and a uniform linear array by constructing a compressed zeroing measurement matrix based on sequential partial Hadamard, combines the data of the non-uniform linear array with the sparse matrix through the constructed deterministic measurement matrix, and uses an adaptive reconstruction algorithm to estimate the incoming wave direction of the underwater target. It has good target detection performance, can not only greatly reduce the number of array elements, but also reduce the amount of calculation and hardware processing costs, has good engineering practice significance, and provides an important idea for active sonar in the array field for detecting targets with a small number of array elements. BRIEF DESCRIPTION OF THE DRAWINGS

[0071] Figure 1 This is a diagram of the uniform vertical linear array receiving signal;

[0072] Figure 2 This is a schematic diagram of airspace grid division;

[0073] Figure 3 This is a schematic diagram of spatial array compression sampling;

[0074] Figure 4 Reconstructing the bearing estimation map for the uncompressed spatial array under wideband signal;

[0075] Figure 5 Reconstructing bearing estimation maps for spatially compressed arrays under wideband signals. DETAILED DESCRIPTION

[0076] The present invention will be further explained below with reference to the accompanying drawings:

[0077] This example is based on the ideal condition that the incoming signal is extended from far-field narrowband to far-field broadband. The N-element uniform vertical linear array is arranged as follows. Figure 1As shown in the figure, the water depth is h = 5m, the speed of sound in water is c = 1500m / s, and far-field narrowband target signals are received from K directions. The array element spacing is d = λ / 2, where:

[0078]

[0079] The kth target signal and its delay can be expressed in complex exponential form as:

[0080]

[0081]

[0082] In the above formula, u k (t), are the amplitude and phase of the signal respectively, and τ is the delay time.

[0083] Since the target is a far-field narrowband signal, it propagates in the water to the array and is received by the array in the form of a plane wave, so And the distance between the signal and the array is much larger than the maximum aperture of the array, so the difference in the source amplitude at each array element can be ignored, that is, u k (t-τ)≈u k (t).

[0084] Therefore, formula (23) can be approximately expressed as:

[0085]

[0086] In order to facilitate theoretical research, it is assumed that the target signals are pairwise uncorrelated echo signals, and that the ocean ambient noise is independent of the target signal source. It is also assumed that the ambient noise received by each array element is a stationary Gaussian white noise with a finite bandwidth and a mean of zero and a variance of σ. 2 Therefore, the noise vector v received by the i-th array element is i (t) and the noise vector v received by the ι-th array element ι (t) Satisfy:

[0087] E[(v i (t)v ι T (t))]=0 (i≠ι) (25)

[0088] W[(v i (t)v l H (t))]=σ 2 I(i=l) (26)

[0089] Based on the characteristics of underwater echo signals in the marine environment, the basic assumption conditions for the array model are summarized as follows:

[0090] a) The array arrangement position is far from the underwater echo source signal, and the echo signals are all far-field narrowband signals and propagate in the form of plane waves.

[0091] b) All array elements are isotropic and can be regarded as point particles, and there is no mutual coupling phenomenon between adjacent array elements.

[0092] c) The total number of array elements is greater than the number of underwater target echo signals to be received, and it is assumed that there are K underwater echo signals.

[0093] The array collects a total of L snapshot data, and the sampling time points are t = 0, 1 / f s , 2 / f s , …, T, where f s represents the signal sampling frequency. Then, the received output signal of the i-th array element under the L snapshot condition can be expressed as the vector superposition of K incoming signals and the underwater environmental noise:

[0094]

[0095] where τ ik represents the reception time difference of the k-th incoming signal at the i-th array element, v i (t) represents the noise of the underwater environment received by the i-th array element, and s k (t - τ ik ) represents that the reception of the k-th incoming signal by the i-th array element can be represented by its relative time delay.

[0096] The time delay τ ik of the echo signal received from direction k between the i-th array element and the reference array element is:

[0097]

[0098] Figure 2 is a schematic diagram of spatial grid division. If the spatial domain is divided by the exhaustive method, an over-complete angle set Θ = {θ1, θ2, …, θ P} can be obtained. To reflect the feasibility of practical applications, in the simulation, the spatial domain is defined as the [0°, 180°] region, and an equally spaced grid division is performed with a step size of 1°. Then, an over-complete angle set Θ = {0°, 1°, …, 180°} can be obtained. The hollow circles represent the grid points for dividing the spatial domain, and the solid circles represent the underwater target echo signals from K directions. When the hollow circle coincides with the solid circle, it means that the direction of arrival of the signal falls exactly on the grid without error. Through such grid division, since it is not an exhaustive division, if the direction of arrival does not fall on the divided grid, a certain error will be introduced.

[0099] However, the model in Step 1 is just a traditional simple sparse model constructed by taking advantage of the natural sparsity of the directivity of array signals in the spatial domain. If the array is arranged in the above-mentioned uniform vertical linear array manner, under the requirement of a large aperture, the number of required array elements will be too large, which will not only increase the computational load in hardware arrangement and processing systems, but also bring a relatively large measurement cost. Therefore, the compressive sensing (CS) theory is applied to array signal processing, and the array is compressed and sampled by using a measurement matrix, as shown in Figure 3 is a schematic diagram of spatial domain array compressive sampling. The compression process corresponds to the construction of the deterministic measurement matrix Φ proposed in this method and is closely related to the compression and zeroing method proposed in Step 2. By compressing the array data in this compression manner, a compressible array (CSA) structure can be obtained. This structure can not only reduce the hardware cost of measurement, but also greatly reduce the computational load, and can break the limitation that the element spacing in the spatial domain is less than half of the signal wavelength. Using the data of a smaller number of array elements, a larger sonar aperture can be compressed and reconstructed to obtain a higher azimuth estimation resolution.

[0100] This embodiment mainly demonstrates the simulation experiment under broadband incident signals. The relevant simulation parameter configurations are as follows: the spatial domain grid is divided with a step length of 0.5°, and the azimuth angle [0°, 180°] is divided into 361 grids. The number of original receiving array elements N = 27, the incident direction of the underwater target signal is 40°, the center frequency is 150KHz, and the linear frequency modulation signal with a bandwidth of 100KHz. The sampling frequency f of the array s is 1000KHz, the number of time domain snapshots is 50000, the sampling time T = 0.05 seconds, and the received signal is divided into Q = 100 segments in the time domain. Q is the number of snapshots in the frequency domain, that is, each frequency domain snapshot is obtained by DFT transformation of the time domain data with L / Q = 500 time domain sampling points. The frequency domain resolution under this transformation is In the frequency domain processing, the broadband signal with a bandwidth of 100Khz is divided into J = 50 segments of narrowband data, and finally the non-compressed spatial domain array reconstruction method under the broadband signal as shown in Figure 4 is obtained by estimation. Using the measurement matrix proposed in this method, the 27-element uniform linear array is compressed into non-uniform linear array data with M = 11 elements, as shown in Figure 5 is the azimuth estimation of underwater broadband incident signals based on CS. In the case of a compression ratio of 0.41, the reconstruction effect is significantly restored. Although there is a certain error (the error range is within 0.5°), this does not affect the final target azimuth estimation. This not only greatly reduces the number of array elements, but also reduces the computational load, which has relatively important significance in engineering practice.

Claims

1. An undersampled wideband array signal processing method based on an improved Hadamard matrix, characterized in that: Step 1: Establish an array spatial domain signal reception model For a uniform vertical linear array that receives far-field narrowband target signals from K directions {θ1, θ2, …, θ K} through N uniformly distributed array elements, its vector equation is: X(t) = A0·s0(t) + v(t) (1) Among them, is the scalar data received by the uniform vertical linear array; t = 1, 2, … L represents the discrete sampling time; is the underwater echo signal amplitude matrix; is the noise scalar data received by the uniform vertical linear array; is the K - array vector manifold matrix, A0 = [a(θ1), a(θ2), …, a(θ K )]; represents a real - number matrix, represents a complex - number matrix; After grid division, the K-array vector flow pattern matrix A0 is extended to A = [a(θ1), a(θ2)…, a(θ K ), …, a(θ P ), where K << P, and the spatial domain sparse model shown in formula (1) is rewritten as: X(t) = A·S(t) + v(t) (2) The amplitude information of the underwater echo signal that is K-sparse; It is the array vector manifold matrix under grid division, where the p-th column vector is expressed as: Among them, c is the propagation speed of sound waves underwater, f0 is the frequency of the underwater echo signal received by the array, d is the element spacing, and e is the natural constant; if there is an actual target echo signal at a(θ p ) where θ k , then S p = s k , otherwise S p = 0; Step 2: Construct a deterministic measurement matrix for a compressible sampling array Perform spatial domain compressive sampling on the array using the measurement matrix Φ, and obtain the data of the N-element linear array After compressive sampling observation, obtain the data of the M-element non-uniform linear array Y = ΦX (4) Among them, is a measurement matrix that satisfies the finite isometry criterion, M << N; expanded for single-snapshot data as: s2.1 Array compression and zeroing method Introduce the idea of sparse decomposition into compressive sampling, and regard the compressive sampling data y = [y1, y2, …, y M T obtained by the linear combination of M column vectors in the measurement matrix Φ and the uniform linear array data X:​ where φ n = [φ 1n , φ 1n , …, φ Mn T ;​ Since the compressed sampling data y is the sampling data of a non-uniform linear array, that is, the data of several array elements under the array element arrangement of a uniform linear array, the data at the corresponding positions of the non-existent array elements is set to zero. Such a zeroing method is called the array compression and zeroing method; s2.2 Construct a deterministic measurement matrix Φ Generate a Hadamard orthogonal matrix H of size U·U U : where U = 2 u , u = 1, 2, … ∞; U > N, intercept the first N columns of the sub-matrix from it to obtain a partial Hadamard measurement matrix, then continuously and sequentially select the first M row vectors to form a new sub-matrix, and then through the array compression and zeroing method described in s2.1, zero the elements at the positions where there are actually no array elements in the new sub-matrix to obtain the spatial domain compression sampling measurement matrix Φ; Step 3: Underwater target detection under narrowband signals Substitute formula (2) into formula (4) to get: Y = Φ(AS + v) = ΦAS + Φv (8) Substitute B = ΦA and V = Φv into formula (8) to get: Y = ΦAS + V (9) Y M×L = B M×P S P×L + V M×L (10) Among them, B is the vector manifold of the compressed sampling array; V is the noise scalar data received by the underwater array; S is the amplitude information of the underwater target signal to be reconstructed and recovered; the problem of solving S is simplified to a minimum l1 norm problem: Use the spatial domain sparse adaptive OMP reconstruction algorithm to solve and reconstruct S; Step 4: Underwater target detection under broadband signals The broadband signals received by each array element are divided into several non-overlapping narrowband data according to frequency bands, and then the methods described in Steps 1 to 3 are used to process each divided narrowband data. Finally, the processing results of each frequency band are averaged and synthesized to obtain the final azimuth estimation result.

2. The undersampled wideband array signal processing method based on the improved Hadamard matrix as claimed in claim 1, wherein: The elements of the Hadamard matrix are 1 or -1, and its construction process successively generates a Hadamard orthogonal matrix of size U·U: H1＝[1] (12) 3. The undersampled wideband array signal processing method based on the improved Hadamard matrix according to claim 1, characterized in that: In Step 3, the reconstruction process of the underwater target signal amplitude information S is: a) Initialize the residual \(r_0 = y\), and the index set Initialize the initial value of the iteration count as \(z = 1\); b) Calculate the index value and store the index value u in Ι; c) Update the support set B Λ = Φ Λ ∪ {b u} and the index set Λ = Λ ∪ Ι; d) Calculate Update the residual e) Until the iteration termination condition is met, output as the reconstruction estimate; otherwise increment the iteration count z by 1 and jump to step b); The iteration termination condition is: Among them, ‖r z ‖2 is the residual value of the iteration number z; Γ z is the sparsity at the z-th iteration.

4. The undersampled wideband array signal processing method based on the improved Hadamard matrix as claimed in claim 1, wherein: The specific process of underwater target detection under broadband signals is: Regard the broadband signal X(t) received by the uniform linear array as the superposition of J narrowband signals: Divide the time-domain received signal of L snapshots into Q non-overlapping segments according to frequency bands, and L / Q is an integer value; then perform DFT transformation on each segment to obtain Q groups of uncorrelated narrowband frequency domain components, and obtain the array signal model in the frequency domain with a single snapshot as: Among them, Introduce the measurement matrix Φ constructed in Step 2 to form a broadband array spatial domain compression observation model. The compressed data at the qth snapshot is: Then, for the broadband signal divided into Q snapshot narrowband data in the frequency domain, use the reconstruction method in Step 3 to reconstruct Q segments of snapshots in turn, and finally perform power averaging and synthesis on the processing results of Q frequency bands to obtain the final azimuth estimation result P(θ):

Images