Large aperture array direction finding method supported by detection information

Abstract

This invention discloses a direction finding method for a large-aperture array supported by acquired information. The method includes the following steps: Step 1: Acquire dual-channel data to be processed and acquired information; Step 2: Perform signal preprocessing and Fourier transform; Step 3: Determine the signal type; if it is a narrowband pulse signal, proceed to Step 4; if it is a wideband pulse signal, proceed to Step 6; Step 4: Extract the envelopes of the two narrowband pulse signals; Step 5: Estimate the arrival time of the two pulse signals and calculate the time delay difference; Step 6: Calculate the cross-correlation spectrum of the two wideband pulse signals; Step 7: Search for the peak value of the cross-correlation spectrum and calculate the time delay difference; Step 8: Calculate the direction of arrival of the signal and output the direction finding result. This method is applicable to large-aperture dual-element arrays, uses acquired information as support, and realizes a direction finding method that matches the signal type. It has strong engineering practicality and is suitable for real-time signal processing.

Description

Technical Field

[0001] This invention belongs to the field of underwater acoustic signal processing technology, specifically relating to a large-aperture array direction finding method supported by detection information. Background Technology

[0002] Accurate direction finding of underwater acoustic pulse signals is of great significance in underwater target detection, underwater acoustic countermeasures, and marine environmental monitoring. Specifically, with the increasing demands for direction finding accuracy, increasing the array aperture, i.e., increasing the element spacing, has become a key means to improve the system's azimuth resolution. However, while large-aperture arrays bring high accuracy, they also introduce challenges such as phase ambiguity. Furthermore, how to achieve highly robust direction finding of non-cooperative pulse signals (such as CW narrowband pulse signals and LFM broadband pulse signals) using front-end detection information in complex underwater acoustic channels is a problem that urgently needs to be solved.

[0003] Underwater acoustic direction finding methods are often classified into three categories: (1) direction finding methods based on phase difference, which use the carrier phase difference of the received signals of array elements to calculate the azimuth. The accuracy is high but there is ambiguity problem; (2) direction finding methods based on time delay difference (TDOA), which calculate the azimuth by estimating the time difference of the signal arriving at each array element. There is no ambiguity but the accuracy of the time delay estimation is limited; (3) direction finding methods based on beamforming, which use weighted summation of the received signals of each array element to form a spatial directional beam to search for the target azimuth. However, its azimuth resolution is limited by the beam width. Usually, more array elements are needed to achieve high-precision direction finding, resulting in a large system hardware and high computational complexity.

[0004] However, current direction-finding methods mainly focus on conventional aperture arrays or high signal-to-noise ratio environments, and research on direction-finding for large aperture arrays in complex reconnaissance scenarios is relatively limited, with the following main problems:

[0005] (1) Phase ambiguity of large aperture arrays is difficult to solve: For mid-to-high frequency signals, the element spacing of large aperture arrays is much greater than half a wavelength, which leads to the ambiguity of traditional methods based on phase difference. The phase is blurred, making it difficult to directly calculate the unique true orientation.

[0006] (2) Poor performance of time delay estimation under low signal-to-noise ratio and multipath interference: Existing time delay estimation methods use relatively simple features. For example, for narrowband pulse signals, the commonly used Hilbert envelope extraction method is easily affected by out-of-band noise and produces spikes under low signal-to-noise ratio, and cannot effectively remove envelope distortion caused by multipath signals; for wideband pulse signals, the full-band cross-correlation method is prone to introducing out-of-band noise, and the full-time domain search is easily affected by strong multipath peak interference, resulting in insufficient direction finding robustness.

[0007] (3) Low utilization of prior detection information: Most existing direction finding algorithms do not utilize prior information, and there is no method to focus on how to deeply integrate prior information such as signal type, time and frequency range and approximate direction provided by the front-end detection system into the signal processing flow.

[0008] It is evident that current direction finding methods are limited in their ability to achieve high-precision, unambiguous direction finding for various types of underwater acoustic pulse signals under conditions of large aperture, low signal-to-noise ratio, and multiple interference paths. Summary of the Invention

[0009] Technical problem: The purpose of this invention is to provide a large-aperture array direction finding method supported by detection information. This method can solve the phase ambiguity problem under large-aperture array conditions, and can also make full use of prior detection information in low signal-to-noise ratio and multipath interference environments to achieve high-precision and high-robustness direction finding for narrowband and broadband pulse signals.

[0010] Technical solutions, such as... Figure 1 As shown, this invention proposes a large-aperture array direction finding method supported by detection information, which includes the following steps:

[0011] (1) For a large-aperture dual-element array with an element spacing d greater than half the wavelength of the signal, obtain the element spacing d and the dual-channel data sequences x1(m) and x2(m) to be processed, and obtain the detection information including the pulse coarse time interval [t1,t2], the pulse coarse frequency interval [f1,f2] and the a priori interval of the incoming wave direction [θ1,θ2].

[0012] (2) Based on the pulse coarse time interval [t1, t2], the dual-channel data sequences x1(m) and x2(m) to be processed are truncated and normalized to obtain the dual-channel signal data sequences s1(n) and s2(n). The discrete Fourier transforms S1(k) and S2(k) of s1(n) and s2(n) are calculated. Based on the pulse coarse frequency interval [f1, f2], bandpass filtering is performed on S1(k) and S2(k) to obtain the truncated first and second channel spectra S. 1_cut (k) and S 2_cut (k);

[0013] (3) Calculate the signal bandwidth B based on the pulse coarse frequency range [f1, f2], and compare it with the bandwidth threshold B. th Compare and determine the signal type. If it is a narrowband pulse signal, proceed to step (4); if it is a wideband pulse signal, proceed to step (6).

[0014] (4) Based on the truncated first and second channel spectra S 1_cut (k) and S 2_cut(k), calculate the signal envelope data sequences y1(n) and y2(n) of the dual-channel signal data sequences s1(n) and s2(n);

[0015] (5) Calculate the precise arrival time n of the two narrowband pulse signals. start1 n start2 And calculate the time delay difference between the two pulse signals. , and proceed to step (8);

[0016] (6) Based on the truncated first and second channel spectra S 1_cut (k) and S 2_cut (k) Calculate the time-domain cross-correlation function r of the two-channel signal data sequences s1(n) and s2(n). 12 (n);

[0017] (7) Search for the modulus of the time-domain cross-correlation function |r 12 The maximum value of (n)| and its corresponding peak index n peak And calculate the time delay difference between the two pulse signals. ;

[0018] (8) Based on the time delay difference between the two pulse signals Calculate the direction of arrival θ of the signal and output the direction finding result.

[0019] Furthermore, in step (1), for a large-aperture dual-element array where the element spacing d is much larger than half the signal wavelength, the following method is used to obtain the element spacing d and the dual-channel data sequences x1(m) and x2(m) to be processed. The detection information obtained includes the pulse coarse time interval [t1, t2], the pulse coarse frequency interval [f1, f2], and the a priori interval of the incoming wave direction [θ1, θ2]. Specifically, the following steps are included:

[0020] (1-1) Obtain the element spacing d of the large aperture dual-element array in meters; synchronously receive the real-time acquired digital signal from two underwater acoustic sensors with a spacing of d, or extract the digital signal containing the target signal from the memory to obtain the dual-channel data sequence to be processed x1(m) and x2(m), where m = 0, 1, 2, …M-1 is the original sampling time index, and M is the original data length, requiring M≥1024;

[0021] (1-2) Prior detection information of the target pulse signal received from the front-end detection system, specifically including:

[0022] ① Coarse pulse time interval [t1, t2]: where t1 is the coarse start time of the target pulse signal and t2 is the coarse end time of the target pulse signal, in seconds, and 0 ≤ t1 is required. <t2≤(M-1) / f s And t1 and t2 are both real numbers, fs is the sampling frequency, with the unit of Hertz;

[0023] ② Pulse rough frequency interval [f1, f2]: where f1 is the lower limit of the rough frequency of the target pulse signal, and f2 is the upper limit of the rough frequency of the target pulse signal, with the unit of Hertz. It is required that 0 ≤ f1 < f2 ≤ f s / 2, and both f1 and f2 are real numbers;

[0024] ③ Prior interval of wave arrival direction [θ1, θ2]: where θ1 and θ2 are the estimated starting angle and ending angle of the wave arrival direction of the target signal respectively, with the unit of degree. It is required that 0° ≤ θ1 < θ2 ≤ 180°, and both θ1 and θ2 are real numbers.

[0025] Furthermore, in step (2), based on the time interval [t1, t2], the dual-channel data sequence to be processed x1(m) and x2(m) are intercepted and normalized to obtain the dual-channel signal data sequences s1(n) and s2(n). Calculate the discrete Fourier transforms S1(k) and S2(k) of s1(n) and s2(n), and perform band-pass filtering on S1(k) and S2(k) based on the frequency interval [f1, f2] to obtain the spectra S 1_cut (k) and S 2_cut (k), which specifically includes the following steps:

[0026] (2-1) Based on the pulse rough time interval [t1, t2] and the sampling frequency f s , calculate the interception start index m start and the interception length L:

[0027]

[0028]

[0029] where min(·) is the function to take the minimum value, max(·) is the function to take the maximum value, floor(·) is the function to round down, and M is the length of the original data;

[0030] (2-2) Respectively, from the dual-channel data sequence to be processed x1(m) and x2(m), intercept data segments with the length of L starting from the index m start . Set the length of the data to be processed N = 2 κ , where κ is the smallest positive integer that satisfies 2 κ ≥ L. If the interception length L < N, then perform zero-padding extension at the end of the intercepted data segment to make its length reach N, obtaining the intercepted data sequences x 1_trunc (n) and x 2_trunc(n), where n = 0, 1, …, N-1, and n is the discrete-time index of the data sequence after dual-channel truncation;

[0031] (2-3) The data sequence x after truncation of the first and second channels 1_trunc (n) and x 2_trunc (n) are normalized to obtain the first and second channel signal data sequences s1(n) and s2(n) respectively:

[0032]

[0033] Where |·| is the modulo function;

[0034] (2-4) Calculate the N-point discrete Fourier transform of the first and second channel signal data sequences s1(n) and s2(n) respectively to obtain the first and second channel spectra S1(k) and S2(k):

[0035]

[0036] Where k = 0, 1, …, N-1 are discrete frequency indices, and j is the imaginary unit, i.e. ;

[0037] (2-5) Based on the pulse coarse frequency range [f1,f2] and the sampling frequency f s Calculate the lower limit k of the fundamental frequency index range. lower and upper limit k upper :

[0038]

[0039] Where ceil(·) is the floor function;

[0040] (2-6) Set the number of frequency domain truncation expansion points Δk, requiring 1≤Δk≤10, and Δk is an integer; based on the basic frequency index range, expand to the left and right by the preset number of sampling points Δk, to obtain the final frequency domain truncation interval [k]. start ,k end ]:

[0041] ;

[0042] (2-7) Based on frequency domain truncation interval [k start ,k end Bandpass filtering was performed on the dual-channel spectra S1(k) and S2(k) respectively to obtain the truncated first and second channel spectra S. 1_cut (k) and S 2_cut (k):

[0043]

[0044] Where k = 0, 1, …, N-1 is the discrete frequency index, and S 1_cut (k), S 2_cut (k) retains only the spectral components within the expanded signal band, filtering out out-of-band noise.

[0045] Furthermore, in step (3), the signal bandwidth B is calculated based on the frequency range [f1, f2] and compared with the bandwidth threshold B. th Compare and determine the signal type. If it is a narrowband pulse signal, proceed to step (4); if it is a wideband pulse signal, proceed to step (6). Specifically, the steps are as follows:

[0046] (3-1) Based on the signal frequency range [f1, f2] in the detected information described in step (1), calculate the bandwidth B and center frequency f0 of the target pulse signal:

[0047]

[0048]

[0049] Where f0 is the center position of the signal frequency band, and the unit is Hertz;

[0050] (3-2) Set the relative bandwidth ratio coefficient η, and calculate the bandwidth decision threshold B based on the center frequency f0. th :

[0051]

[0052] Where η is a dimensionless constant, 0.01 < η ≤ 0.2, used to define the relative dividing line between narrowband and wideband signals;

[0053] (3-3) Compare the signal bandwidth B calculated in step (3-1) with the bandwidth decision threshold B calculated in step (3-2). th Compare: If B th If the signal type is determined to be a narrowband pulse signal, then proceed to step (4); if B ≥ B th If the signal type is determined to be a broadband pulse signal, then proceed to step (6).

[0054] Furthermore, in step (4), based on the truncated first and second channel spectra S 1_cut (k) and S 2_cut (k) Calculate the signal envelope sequences y1(n) and y2(n) of the dual-channel signal data sequences s1(n) and s2(n), specifically including the following steps:

[0055] ​(4-1) Initialize the parameters for calculating the narrowband pulse signal envelope sequence, specifically including the initialization of the following parameters:

[0056] ① The scaling factor ζ of the iterative smoothing window is initialized to a real number of 10 ≤ ζ ≤ 100;

[0057] ②The initial envelope iteration smoothing window length M2 is initialized as follows:

[0058]

[0059] Where, round(·) is the rounding function, L is the truncation length mentioned in step (2-1), and M2 is a positive odd number;

[0060] ③ Maximum number of iterations for smoothing I max Initialize to: 2≤I max Integers ≤ 10;

[0061] ④ The iterative convergence accuracy threshold δ is initialized to: 10 -5 ≤δ≤10 -3 real numbers;

[0062] ⑤ The number of envelope iterations for smoothing is initialized to I = 0;

[0063] (4-2) The spectrum S after step (2-6) is respectively 1_cut (k), S 2_cut (k) Perform inverse discrete Fourier transform and take the modulus to obtain the initial envelope sequences y of the first and second channels. 1_init (n) and y 2_init (n):

[0064]

[0065] Where n = 0, 1, …, N-1 is the discrete-time index of the dual-channel signal data sequence, N is the data length, and k = 0, 1, …, N-1 is the discrete-frequency index;

[0066] (4-3) Set the initial result y1 for the envelope smoothing of the first and second channels. (0) (n) and y2 (0) (n):

[0067] ;

[0068] (4-4) Update envelope iteration smoothing number I:

[0069] ;

[0070] (4-5) In the I-th iteration, the result y1 calculated in the previous iteration is smoothed using the initial envelope iteration smoothing window length M2. (I-1) (n), y2 (I-1) (n) Smoothing:

[0071]

[0072] When the index n+m<0 or n+m≥N, zero-filling is used to handle the boundary, and y1 is defined. (I-1) (n+m)=0 and y2 (I-1) (n+m)=0;

[0073] (4-6) Calculate the sum of squared residuals ε of the first and second channels in the I-th iteration. y1 (I) and ε y2 (I) :

[0074] ;

[0075] (4-7) Determine whether the iterative loop satisfies the following continuation conditions:

[0076]

[0077] If the condition is met, return to step (4-4); otherwise, proceed to step (4-8).

[0078] (4-8) Output the final envelope sequences y1(n) and y2(n) of the first and second channel narrowband pulse signals:

[0079] .

[0080] Furthermore, in step (5), the precise arrival time n of the two narrowband pulse signals is calculated using the following method. start1 n start2 And calculate the time delay difference between the two pulse signals. Specifically, it includes the following steps:

[0081] (5-1) Pulse edge feature search window length L win Initialize to:

[0082]

[0083] Where, round(·) is the rounding function, max(·) is the maximum value function, L is the truncation length mentioned in step (2-1), ξ is the proportional coefficient of the pulse along the feature search window length, and ξ is a dimensionless positive real number, requiring 10≤ξ≤50;

[0084] (5-2) Calculate the maximum value y1(n) and y2(n) of the envelope sequences of the first and second channel narrowband pulse signals. 1_max and y 2_max and the minimum value y 1_min and y 2_min :

[0085]

[0086] Where n = 0, 1, …, N-1 are the discrete-time indices of the dual-channel signal data sequence;

[0087] (5-3) Set a very small positive real number ε, and calculate the normalized value y of the envelope sequence of the narrowband pulse signals of the first and second channels. 1_norm and y 2_norm :

[0088]

[0089] Where ε is a specified minimal positive real number used to prevent the denominator from being zero, and ε takes values ​​of

[10] . -16 10 -14 Any real number within ];

[0090] (5-4) Calculate the first-order forward difference sequence y of the first and second channels. 1_diff (n) and y 2_diff (n):

[0091]

[0092] Among them, y 1_diff (n), y 2_diff (n) represents the rate of change of the envelope sequence at time n;

[0093] (5-5) Calculate the count values ​​z of the difference sequences within the first and second channel windows that are greater than zero. 1_pos (n) and z 2_pos (n):

[0094]

[0095] Where u(·) is a step function, which takes the value 1 when the input is greater than 0, and takes the value 0 otherwise;

[0096] (5-6) Based on the count value z 1_pos (n), z 2_pos (n) Calculate the consecutive rising weights w of the rising edges of the pulse envelopes of the first and second channels. 1_rise (n) and w 2_rise (n):

[0097]

[0098] Among them, w 1_rise (n), w 2_rise (n) is used to evaluate the monotonically increasing degree of the subsequent signal envelope at time n;

[0099] (5-7) Using the normalized difference sequence y 1_diff (n), y 2_diff (n) as the transient center weights w 1_center (n), w 2_center (n) of the rising edge of the pulse envelope, calculate the comprehensive evaluation factors q 1_rise (n) and q 2_rise (n) of the rising edge of the pulse envelope for the first and second channels:

[0100]

[0101] Among them, q 1_rise (n), q 2_rise (n) comprehensively reflects the instantaneous change rate of the signal at time n and the monotonicity of the subsequent trend, and its peak value corresponds to the direct arrival wave front of the pulse signal;

[0102] (5-8) Within the range of 0 ≤ n < N, search for the discrete sampling times corresponding to the maximum values of the comprehensive evaluation factors q 1_rise (n), q 2_rise (n) of the rising edge of the pulse envelope for the first and second channels, and use them as the precise arrival time n start1 and n start2 :

[0103]

[0104] Among them, represents the index value that makes the function reach the maximum value within the range of 0 ≤ n < N;

[0105] (5-9) Calculate the time delay difference :

[0106]

[0107] Among them, f s is the sampling frequency.

[0108] Further, in step (6), based on the intercepted first and second channel spectra S 1_cut (k) and S 2_cut (k), calculate the time-domain cross-correlation function r 12 (n) of the two-channel signal data sequences s1(n) and s2(n), which specifically includes the following steps:

[0109] (6-1) Based on the truncated first and second channel spectra S obtained in step (2-7) 1_cut (k) and S 2_cut (k) will extract the spectrum of the second channel S. 2_cut (k) and the complex conjugate S of the truncated first channel spectrum * 1_cut (k) Multiply to calculate the frequency domain cross-correlation spectrum R. 12 (k):

[0110]

[0111] Among them, S * 1_cut (k) represents S 1_cut The complex conjugate of (k), where k = 0, 1, …, N-1 is the discrete frequency index;

[0112] (6-2) For the frequency domain cross-correlation spectrum R 12 (k) Perform an N-point discrete Fourier inverse transform to obtain the time-domain cross-correlation function r. 12 (n):

[0113]

[0114] Where n = 0, 1, …, N-1 is the discrete-time index of the dual-channel signal data sequence, and N is the data length.

[0115] Furthermore, in step (7), the modulus of the time-domain cross-correlation function |r is searched using the following method. 12 The maximum value of (n)| and its corresponding peak index n peak And calculate the time delay difference between the two pulse signals. Specifically, it includes the following steps:

[0116] (7-1) Set the underwater acoustic velocity c = 1500 m / s; based on the prior interval [θ1, θ2] of the incoming wave direction in the detected information described in step (1), and combined with the array element spacing d and the underwater acoustic velocity c, calculate the lower bound τ of the theoretical time delay search interval of the two signals. min and the upper bound τ max :

[0117]

[0118] Where θ1 and θ2 are the estimated starting angle and ending angle respectively, and θ1 < θ2, and cos(·) is the cosine function;

[0119] (7-2) Based on the periodic extension property of the discrete Fourier transform, a piecewise mapping function between the discrete time index n and the physical delay τ(n) is established:

[0120]

[0121] Where n = 0, 1, …, N-1 are the discrete-time indices of the dual-channel signal data sequence, N is the data length, and f s The sampling frequency;

[0122] (7-3) Search interval [τ] based on physical time delay min ,τ max ] and mapping function τ(n), construct an efficient search index set Ω search :

[0123]

[0124] Among them, Ω search It contains all discrete-time indices that fall within the prior physical delay constraints;

[0125] (7-4) In the set of valid search indexes Ω search Within, search for the modulus of the time-domain cross-correlation function |r 12 The maximum value of (n)| and its corresponding peak index n peak :

[0126]

[0127] in, This indicates that in n∈Ω search The index value within the range that allows the function to reach its maximum value;

[0128] (7-5) Based on peak index n peak Based on the mapping relationship in step (7-2), calculate the time delay difference between the two pulse signals. :

[0129]

[0130] in, This is a high-precision time delay difference estimate after prior constraints and in-band extraction optimization.

[0131] Furthermore, in step (8), based on the time delay difference between the two pulse signals... Calculate the direction of arrival θ of the signal and output the direction finding result, specifically including the following steps:

[0132] (8-1) Receive the time delay difference of the narrowband pulse signal output in step (5) Or the time delay difference of the broadband pulse signal output in step (7) Combining the array element spacing d and underwater acoustic velocity c obtained in step (1), the intermediate variable V of the direction cosine is calculated. temp :

[0133] ;

[0134] (8-2) The intermediate variable V of the direction cosine obtained in step (8-1) temp Numerical constraints are applied to force the direction cosine value V to be confined to the interval [-1, 1]. clip :

[0135]

[0136] Among them, V clip This is used to eliminate the problem of the inverse cosine function's independent variable going out of bounds due to measurement errors, ensuring the mathematical validity of subsequent angle calculations;

[0137] (8-3) Based on the constrained direction cosine value V clip The direction of arrival θ of the signal is calculated using the inverse cosine function:

[0138]

[0139] Where arccos(·) is the inverse cosine function. θ is the conversion factor for radians to degrees, with the unit being degrees;

[0140] (8-4) The calculated direction of arrival θ is output as the final direction finding result of the large aperture array supported by the detection information.

[0141] Beneficial effects: Compared with the prior art, the technical solution of the present invention has the following beneficial technical effects:

[0142] 1. This invention designs a differentiated direction-finding architecture that matches physical characteristics, improving adaptability to complex environments. As shown in step 3, for narrowband signals, the "envelope rising edge extraction" strategy is adopted to transform phase ambiguity into a time delay estimation problem, thus avoiding physical limitations. For broadband signals, the "frequency domain in-band extraction cross-correlation" strategy is adopted to fully utilize pulse compression gain and achieve high-precision time delay estimation.

[0143] 2. This invention fully utilizes the constraints of the detected frequency range and the envelope iterative smoothing technique, as shown in steps 2 and 4, to achieve high signal-to-noise ratio envelope extraction of narrowband pulse signals. This process abandons the traditional Hilbert transform method and effectively filters out out-of-band noise and eliminates envelope glitches through frequency domain truncation filtering and time domain iterative smoothing, providing a high-quality foundation for time delay estimation.

[0144] 3. This invention fully utilizes the normalized differential and pulse envelope rising edge comprehensive evaluation factors, as shown in step 5, to achieve accurate estimation of the arrival time of narrowband pulse signals. This process comprehensively considers the transient center and continuous rising weight, effectively resisting the envelope tail and distortion caused by the multipath effect of the underwater acoustic channel, and accurately locking the start time of the direct wave.

[0145] 4. This invention fully utilizes the frequency domain sparsity characteristics of the detected frequency information, as shown in steps 2 and 6, to improve the significance of the cross-correlation peak of the broadband pulse signal. By extracting out-of-band noise spectral lines in the frequency domain and setting them to zero, the signal-to-noise ratio under low signal-to-noise ratio conditions is significantly improved, solving the problem that traditional full-band cross-correlation is easily submerged by noise.

[0146] 5. This invention fully utilizes prior angle information and physical aperture constraints, as shown in steps 7 and 8, to solve the phase ambiguity problem of large-aperture arrays. Through constrained time-delay search and numerical anti-boundary constraint mechanisms, it avoids erroneous locking of multipath peaks or grating lobe peaks, ensuring the uniqueness and physical validity of the direction finding results. Attached Figure Description

[0147] Figure 1 This is a schematic flowchart of the method of the present invention;

[0148] Figure 2 This is a schematic diagram of the first channel spectrum S1(k) and the frequency domain truncation interval in Example 1;

[0149] Figure 3 This is a schematic diagram of the envelope data sequence y1(n) after iterative smoothing of the first channel in Example 1;

[0150] Figure 4 The comprehensive evaluation factor q for the rising edge of the first channel pulse envelope in Example 1. 1_rise (n) Schematic diagram;

[0151] Figure 5 This is a schematic diagram of the first channel spectrum S1(k) and the frequency domain truncation interval in Example 2;

[0152] Figure 6 The frequency domain cross-correlation spectrum R of Example 2 12 (k) Schematic diagram;

[0153] Figure 7 The time-domain cross-correlation function magnitude |r in Example 2 12 (n)| and a diagram illustrating the peak search results. Detailed Implementation

[0154] To better understand the purpose, structure, and function of this invention, the invention will be further described below with reference to the accompanying drawings.

[0155] In an embodiment of the present invention, a direction-finding scenario using a large-aperture dual-element array was constructed, with the element spacing d set to 20 meters. For this large-aperture dual-element array, a simulated dual-channel signal model x was used to receive the signal. i (t) is defined as:

[0156]

[0157] Where i is the channel index, i=1,2; A is the signal amplitude; j is the imaginary unit, i.e. τ0 is the signal transmission time; τ is the signal pulse width; τ i ω is the propagation delay of the signal to the i-th array element; i (t) represents a variable with mean 0 and variance σ. 2 Complex Gaussian white noise, variance σ 2 The magnitude depends on the signal-to-noise ratio (SNR), and the relationship is: SNR = 10log 10 (A 2 / σ 2 ).

[0158] φ(t) is the instantaneous phase of the signal, and according to the signal types covered by this invention, it is expressed as follows:

[0159]

[0160] Where φ0 is the initial phase of the signal; when the signal is a narrowband pulse signal, f c f1 is the carrier frequency; when the signal is a broadband pulse signal, f1 is the signal start frequency, and μ is the modulation frequency of the linear frequency modulated signal, defined as μ=(f2-f1) / τ, where f2 is the signal end frequency.

[0161] With sampling frequency f s For the signal x received in the above simulation i Discrete sampling is performed on (t) to obtain the dual-channel signal sampling data sequence x. i (n) is:

[0162]

[0163] Wherein, the discrete starting point n of the signal in the i-th channel start,i =round{(τ0+τ i f s}; Number of signal sampling points N = round{τ f s}; round{·} is the rounding function.

[0164] The discrete form of the instantaneous phase of the corresponding signal φ(n / f) s ) is represented as:

[0165]

[0166] Based on the above signal model, the present invention has conducted simulation verification for narrowband pulse signals and wideband pulse signals, as detailed in Embodiments 1 and 2 below.

[0167] Example 1

[0168] The simulation signal is a narrowband pulse signal, and its parameters are set as follows: element spacing d = 20m, sound velocity c = 1500m / s, and sampling frequency f. s =75kHz; Signal center frequency f c =3400Hz, signal pulse width τ=0.1s, signal period T=0.3s, target direction of arrival θ tar =135°, target distance R=5000m, signal-to-noise ratio SNR=-10dB.

[0169] The simulation signal is then processed using a large-aperture array for direction finding:

[0170] Based on steps (1) and (2), the detection information is obtained as follows: the pulse coarse time interval [t1,t2]=[0,0.15]s, the pulse coarse frequency interval [f1,f2]=[3350,3450]Hz, and the a priori interval of the incoming wave direction [θ1,θ2]=[0°,180°]. Based on the detection time interval, data segments are extracted from the dual-channel original data and normalized to obtain dual-channel data sequences s1(n) and s2(n) of length N; their discrete Fourier transforms are calculated to obtain the frequency domain spectra S1(k) and S2(k); the number of frequency domain truncation expansion points Δk=5 is set, and the frequency domain truncation interval [K] is constructed based on the detection frequency interval [3350,3450]Hz. start ,K end = [1458, 1513], bandpass filtering is performed on S1(k) and S2(k) to obtain the truncated first and second channel spectra S. 1_cut (k) and S 2_cut (k); the spectrum of the first channel S1(k) and the frequency domain truncation interval are as follows Figure 2 As shown;

[0171] Based on step (3), the signal bandwidth B = f2 - f1 = 100Hz is calculated based on the detected frequency range, and the center frequency f0 = 3400Hz is calculated; the relative bandwidth ratio coefficient η = 0.1 is set, and the bandwidth decision threshold B is calculated. th =ηf0=340Hz; Since B=100Hz th =340Hz, signal type is narrowband pulse signal, proceed to step (4).

[0172] According to step (4), the truncated first and second channel spectra S are respectively...​1_cut (k) and S 2_cut (k) Perform inverse discrete Fourier transform and take the modulus to obtain the initial envelope sequences y1(n) and y2(n); set the scaling factor ζ = 28 for the iterative smoothing window length, the initial envelope iterative smoothing window length M2 = 403, and the maximum number of iterations I. max =5, iteration convergence accuracy threshold δ=10 -4 The initial envelope sequence is iteratively smoothed to obtain the first and second channel narrowband pulse signal envelope sequences y1(n) and y2(n); the smoothed envelope data sequence y1(n) of the first channel is as follows: Figure 3 As shown;

[0173] According to step (5), the envelope sequences y1(n) and y2(n) of the first and second channel narrowband pulse signals are normalized, and the first-order forward difference sequences y1(n) of the first and second channels are calculated. 1_diff (n) and y 2_diff (n); Set the scaling factor ξ = 50 for the pulse-edge feature search window length, then the pulse-edge feature search window length L win =225; Calculate the consecutive rising weights w of the rising edges of the pulse envelopes of the first and second channels. 1_rise (n) and w 2_rise (n), and calculate the comprehensive evaluation factor q of the rising edge of the pulse envelope of the first and second channels. 1_rise (n) and q 2_rise (n), searching for q within the characteristic search interval along the pulse. 1_rise (n) and q 2_rise Find the maximum value of (n) to obtain the precise arrival time n of the two signals. start1 =1583、n start2 =876; Calculate the time delay difference between the two pulse signals. The comprehensive evaluation factor q of the rising edge of the first channel pulse envelope. 1_rise (n) such as Figure 4 As shown;

[0174] Based on step (8), the time delay difference between the two pulse signals is... Using the formula The direction of arrival was calculated; the calculated direction finding result was θ=134.99°, which has a very small error with the actual angle, thus realizing accurate direction finding of narrowband pulse signals by a large aperture array under low signal-to-noise ratio.

[0175] Example 2

[0176] The simulated signal is a broadband pulse signal, with the following parameters set: element spacing d = 20m, sound velocity c = 1500m / s, and sampling frequency f. s =75kHz; Signal start frequency f start=3200Hz, termination frequency f stop =3600Hz, signal pulse width τ=0.1s, signal period T=0.3s, target direction of arrival θ tar =150°, target distance R=5000m, signal-to-noise ratio SNR=-10dB.

[0177] The simulation signal is then processed using a large-aperture array for direction finding:

[0178] Based on steps (1) and (2), the detection information is obtained: the pulse coarse time interval [t1,t2]=[0,0.15]s, the pulse coarse frequency interval [f1,f2]=[3150,3650]Hz, and the a priori interval of the incoming wave direction [0°,180°]. Data is extracted and normalized based on the detection time interval to obtain dual-channel data sequences s1(n) and s2(n); their discrete Fourier transforms are calculated to obtain the spectra S1(k) and S2(k); the number of frequency domain extraction expansion points Δk=5 is set, and the frequency domain extraction interval [K] is constructed based on the detection frequency interval [3150,3650]Hz. start ,K end = [1371, 1600], perform bandpass filtering on S1(k) and S2(k) to obtain the truncated first and second channel spectra S. 1_cut (k) and S 2_cut (k); the spectrum of the first channel S1(k) and the frequency domain truncation interval are as follows Figure 5 As shown;

[0179] Based on step (3), the signal bandwidth B = f2 - f1 = 500 Hz is calculated based on the detected frequency range, and the center frequency f0 = 3400 Hz is calculated; the relative bandwidth ratio coefficient η = 0.1 is set, and the bandwidth decision threshold B is calculated. th =ηf0=340Hz; Since B=500Hz>B th =340Hz, signal type is broadband pulse signal, proceed to step (6).

[0180] According to step (6), the truncated second channel spectrum S 2_cut (k) and the complex conjugate S of the truncated first channel spectrum * 1_cut (k) Multiply to calculate the frequency domain cross-correlation spectrum R. 12 (k); for R 12 (k) Perform the inverse discrete Fourier transform to obtain the time-domain cross-correlation function r. 12 (n); Frequency domain cross-correlation spectrum R 12 (k) such as Figure 6 As shown;

[0181] Based on step (7), and using the detected angle prior interval [0°, 180°], combined with the element spacing and sound speed, the lower bound τ of the physical time delay search interval is calculated. min =-0.0133s and upper bound τ max =0.0133s; Search for the magnitude of the time-domain cross-correlation function |r| within the physical time delay search interval [-0.0133, 0.0133]s. 12 Find the maximum value of (n)|, lock the peak index, and calculate the time delay difference between the two pulse signals. ; Modulus of the time-domain cross-correlation function |r 12 (n)| and peak search results such as Figure 7 As shown;

[0182] Based on step (8), the time delay difference between the two pulse signals is... Using the formula Calculate the direction of arrival; the calculated direction finding result is θ. tar =150°, consistent with the actual angle, verifying that under an extremely low signal-to-noise ratio of -10dB, the present invention effectively solves the problems of noise interference and phase ambiguity through frequency domain extraction and limited space search.

[0183] It is understood that the present invention has been described through some embodiments, and those skilled in the art will recognize that various changes or equivalent substitutions can be made to these features and embodiments without departing from the spirit and scope of the invention. Furthermore, under the teachings of the present invention, these features and embodiments can be modified to adapt to specific situations and materials without departing from the spirit and scope of the invention. Therefore, the present invention is not limited to the specific embodiments disclosed herein, and all embodiments falling within the scope of the claims of this application are within the protection scope of the present invention.

Images