A high-precision time delay estimation method based on parallel factor

Abstract

This invention discloses a high-precision time delay estimation method based on parallel factors. First, L monitoring nodes receive radiated signals, and the received signals undergo Fourier transform. Second, non-zero frequency band signals are extracted, and frequency domain data from each monitoring node are fused to construct a trilinear model. Then, the parallel factor method is used, and the trilinear alternating least squares method is employed for decomposition and iteration to obtain the frequency domain time delay estimation matrix. Finally, a normalization method is used to eliminate scale ambiguity, and the Vandermonde feature of the time delay matrix is ​​used for time delay estimation. Compared with traditional methods, this invention achieves higher accuracy in time delay estimation, effectively improving positioning accuracy; it eliminates the need for spectral peak search, resulting in relatively low computational complexity; and by fusing frequency domain data from multiple nodes, multiple time delay estimates can be obtained simultaneously without the need for node-specific processing, making it of significant application value.

Description

Technical Field

[0001] This invention relates to the field of passive positioning technology, and in particular to a high-precision time delay estimation method based on parallel factors. Background Technology

[0002] The time difference of arrival (TDOA) method is an important tool in passive positioning technology. As an important component of signal parameter estimation and extraction in modern signal processing, it is widely used in radar, sonar, wireless communication system positioning, radio monitoring and other fields. Therefore, further improving the performance of passive time difference positioning technology is of great practical significance.

[0003] Time delay estimation is a crucial aspect of time delay difference measurement and positioning. Specifically, it involves using parameter estimation and signal processing theories and methods to estimate the time delay of signals from the same source due to differences in transmission distance, thereby obtaining relevant parameters such as the target's location. Current research on time delay estimation still faces significant challenges. For example, time delay estimation methods based on cross-correlation algorithms obtain the time delay estimate by analyzing the correlation peaks of the received and reference signals. Since the maximum time delay estimation accuracy achievable by cross-correlation algorithms depends on the system's sampling period, the positioning accuracy of the system is limited by the system's signal sampling rate. Furthermore, the statistical characteristics of the signal and noise are unknown, leading to significant errors in this method. While time delay estimation methods based on subspace decomposition provide high-resolution results, they are computationally intensive, their performance is limited to signals with flat or approximately flat spectra, and they are sensitive to noise. Cost function-based time delay estimation algorithms involve extensive parameter optimization, and their high computational complexity prevents them from achieving dynamic estimation requirements. Time delay estimation errors affect the subsequent accuracy of radiation source positioning, making it difficult to achieve asymptotically optimal estimation performance. Therefore, exploring a high-precision time delay estimation method for complex scenarios has significant practical implications and research value. Summary of the Invention

[0004] The technical problem to be solved by the present invention is to provide a high-precision time delay estimation method based on parallel factors, which addresses the deficiencies mentioned in the background art.

[0005] To solve the above-mentioned technical problems, the present invention adopts the following technical solution:

[0006] A high-precision time delay estimation method based on parallel factors, the method includes the following steps:

[0007] Step 1) Receive radiation source signals at L known monitoring nodes and perform Fourier transform on the received signals, where L is a preset threshold for the number of base stations.

[0008] Step 2) Extract non-zero frequency band signals and fuse the frequency domain data of each monitoring node to construct a trilinear model;

[0009] Step 3) Using the parallel factor method, the frequency domain time delay estimation matrix is ​​obtained by decomposing and iterating through the trilinear alternating least squares method.

[0010] Step 4) Eliminate scale ambiguity by normalization and use the Vandermonde feature of the time delay matrix for time delay estimation.

[0011] As a further optimization of the high-precision time delay estimation method based on parallel factors of the present invention, step 1) includes the following steps:

[0012] Step 1.1): Receive radiation source signals at L known monitoring nodes. The received signal x at the l-th monitoring node... l (t)=α l s k (t-τ l )+n l (t), where 1≤l≤L, s k (t) represents the radiation source signal received in the k-th segment, 1≤k≤K, where K is a preset threshold for the number of radiation source signal segments, α l and τ l Let n represent the path attenuation factor and time delay of the radiation source signal reaching the l-th monitoring node, respectively. l (t) represents mutually independent zero-mean additive white Gaussian noise, and the signal and noise are uncorrelated.

[0013] Step 1.2): Perform a Fourier transform on the received signal to obtain the time-domain output matrices of L monitoring nodes, where the time-domain output matrix of the l-th monitoring node is in the form of...

[0014]

[0015] In the formula, x l (t) represents the received signal of the l-th monitoring node, F H Let x be the inverse Fourier transform matrix. l (ω) is x l The Fourier transform form of (t), α l and τ l Represent the path attenuation factor and time delay of the radiation source signal reaching the l-th monitoring node, respectively. k (ω1)…S k (ω N )] T and [N(ω1)…N(ω N )] T They represent s respectively k (t) and n lThe Fourier transform form of (t), where ω n Let N be the nth frequency point, where 1 ≤ n ≤ N, and N is the number of sampling points.

[0016] As a further optimization of the high-precision time delay estimation method based on parallel factors of the present invention, the specific steps of step 2) are as follows:

[0017] Step 2.1), the time-domain output matrices of the L monitoring nodes are processed according to the frequency band w = [ω c1 ,...,ω cP ] T Extract:

[0018]

[0019] in, It is the inverse Fourier transform matrix F H According to the frequency band w = [ω c1 ,...,ω cP ] T The column extraction performed, [S k (ω c1 )…S k (ω cP )] T and [N(ω c1 )…N(ω cP )] T They represent [S] k (ω1)...S k (ω N )] T and [N(ω1)...N(ω N )] T According to the frequency band w = [ω c1 ,...,ω cP ] T The row extraction is performed, where ω cp Let P be the p-th frequency point, 1≤p≤P, where P is the number of frequency points extracted.

[0020] Step 2.2), define the time delay vector of the l-th node as... Then the time delay matrix of L nodes is

[0021]

[0022] Ignoring noise, the received signal of monitoring node l is X. l =F c H diag(h l )S T =F c H Dl (H)S T

[0023] in,

[0024] To collect the source matrix of K-band radiation source signals, For the k-th segment of the radiated signal, according to the frequency band w = [ω c1 ,...,ω cP ] T Extracted non-zero vector; X l Let l be the received signal matrix of the l-th monitoring node. It is the inverse Fourier transform matrix F H According to the frequency band w = [ω c1 ,...,ω cP ] T Column extraction was performed; by combining the received signal matrices from L nodes, a three-dimensional matrix was constructed:

[0025]

[0026] Where X is a three-dimensional matrix of received signals from L monitoring nodes, a parallel factor trilinear model is constructed:

[0027]

[0028] Where, x n,k,l f is the (n,k,l)th element in the three-dimensional matrix X. n,p It is matrix F c H The (n,p)th element, s k,p h is the (k,p)th element in matrix S. l,p It is the (l,p)th element in matrix H.

[0029] As a further optimization of the high-precision time delay estimation method based on the parallel factor of the present invention, the frequency domain time delay estimation matrix obtained in step 3) is:

[0030]

[0031] in, This is the time delay estimation matrix. for Column vectors.

[0032] As a further optimization of the high-precision time delay estimation method based on parallel factors of the present invention, the specific steps of step 4) are as follows:

[0033] Step 4.1): Solve the scale ambiguity problem using a normalization method, and set a...p All elements in p = 1, 2, ..., P are divided by the first element of the column. At the same time, the time delay estimate is also changed to the time delay difference estimate relative to the first monitoring node.

[0034]

[0035] in This is the time delay estimation matrix after row vector normalization. for Column vectors;

[0036] Step 4.2), will Column vector a in the matrix pp Divide p = 2, ..., P by the dot notation between the previous column to eliminate the attenuation factor and prevent phase blurring.

[0037]

[0038] in, This is the time delay estimation matrix after column vector normalization. for The row vectors are obtained by taking their modulus:

[0039]

[0040] Where q is a vector containing frequency difference and time delay difference;

[0041] Step 4.3), using the known inverse Fourier transform matrix F c H Solve the column swap ambiguity problem to obtain the corresponding frequency band. Then calculate the corresponding frequency difference w. c =[0,ω c2 -ω c1 ,...,ω cP -ω cP-1 ] T The time delay difference estimate can then be obtained:

[0042]

[0043] in, This represents the time delay difference between the l-th monitoring node and the 1st monitoring node (reference node).

[0044] Compared with the prior art, the present invention, employing the above technical solution, has the following technical effects:

[0045] This invention provides a high-precision time delay estimation method based on parallel factors. By fusing frequency domain data from various monitoring nodes, a trilinear model is constructed, and time delay is estimated using trilinear decomposition. This solves the problems of limited accuracy and large computational load in traditional time delay estimation methods, improves the accuracy of time delay estimation, and provides a guarantee for the accuracy of passive positioning. Attached Figure Description

[0046] Figure 1 A flowchart of the high-precision time delay estimation method based on parallel factors provided by the present invention;

[0047] Figure 2 This is a scene diagram of source localization based on TDOA as described in this invention;

[0048] Figure 3 This is a time delay estimation diagram under the simulation scenario of the method described in this invention;

[0049] Figure 4 This is a comparison chart of the time delay estimation performance of the method described in this invention and the generalized cross-correlation method under different signal-to-noise ratios;

[0050] Figure 5 This is a diagram illustrating the distribution of information sources and nodes in the actual test environment described in this invention.

[0051] Figure 6 This is a time delay estimation diagram under the actual test scenario of the method described in this invention. Detailed Implementation

[0052] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings:

[0053] This invention can be implemented in many different forms and should not be considered limited to the embodiments described herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully express the scope of the invention to those skilled in the art. In the drawings, components are enlarged for clarity.

[0054] It should be understood that although the terms first, second, third, etc., may be used herein to describe various elements, components, and / or parts, these elements, components, and / or parts are not limited by these terms. These terms are merely used to distinguish elements, components, and / or parts from one another. Therefore, the first element, component, and / or part discussed below may be a second element, component, or part without departing from the teachings of this invention.

[0055] The detailed flowchart of a high-precision time delay estimation method based on parallel factors provided by this invention is as follows: Figure 1As shown, L monitoring nodes with known locations receive radiation source signals, and Fourier transforms are performed on the received signals. Non-zero frequency band signals are extracted, and frequency domain data from each monitoring node are fused to construct a trilinear model. A parallel factor method is used, and the trilinear alternating least squares method is employed for decomposition and iteration to obtain the frequency domain time delay estimation matrix. Scale ambiguity is eliminated through normalization, and the time delay estimation result is obtained using the Vandermonde eigenvalue of the time delay matrix. The specific implementation is as follows:

[0056] Step 1), use L monitoring nodes with known locations to receive radiation source signals and perform Fourier transform on the received signals:

[0057] Consider as Figure 2 The scenario shown illustrates a positioning problem. Assume there is an unknown location in space, u = [x, y, z]. T The radiation source signal s(t) exists, and there are L monitoring nodes with known locations that can receive the signal, located at p respectively. l =[x l ,y l ,z l ] T If l = 1, 2, ..., L, then the received signal of the l-th (1 ≤ l ≤ L) monitoring node is x. l (t)=α l s k (t-τ l )+n l The expression for the received signals from L base stations is (t).

[0058] x1(t)=α1s k (t-τ1)+n1(t)

[0059]

[0060] x L (t)=α L s k (t-τ L )+n L (t)

[0061] Among them, s k (t)(1≤k≤K) represents the radiation source signal received in the k-th segment, where K is a preset threshold for the number of radiation source receiving segments, and α l and τ l Let n represent the path attenuation factor and time delay of the radiation source signal reaching the l-th monitoring node, respectively. l (t) represents mutually independent zero-mean additive white Gaussian noise, with no correlation between the signal and noise. The received signal x at the l-th monitoring node... l(t) Perform a Discrete Fourier Transform (DFT). Let N be the number of sampling points, i.e., the number of Fourier transform points. Then, the frequency domain output matrix of the l-th monitoring node is in the form of:

[0062]

[0063] Where, x l (ω) is x l The Fourier transform form of (t), α l and τ l Represent the path attenuation factor and time delay of the radiation source signal reaching the l-th monitoring node, respectively. k (ω1)...S k (ω N )] T and [N(ω1)...N(ω N )] T They represent s respectively k (t) and n l The Fourier transform form of (t), where ω n (1≤n≤N) represents the nth frequency point, and N represents the number of sampling points.

[0064] make If the DFT matrix is ​​F, then the IDFT matrix is ​​F. H Therefore, the time-domain output matrix of the l-th monitoring node is in the form of

[0065]

[0066] Where, x l (t) represents the vector form of the signal received by the l-th monitoring node.

[0067] Step 2) Extract non-zero frequency band signals and fuse the frequency domain data from each monitoring node to construct a trilinear model:

[0068] Because the sampling bandwidth of the signal is greater than the signal bandwidth, therefore [S k (ω1),...,S k (ω N )] T Many elements in the matrix are zero or close to zero, assuming they are in w = [ω]. c1 ,...,ω cP ] T If the frequency band is non-zero, the above formula is modified to:

[0069]

[0070] in, It is the IDFT matrix F H According to the frequency band w = [ωc1 ,...,ω cP ] T The column extraction performed, [S k (ω c1 ...S k (ω cP )] T and [N(ω c1 ...N(ω) cP )] T They represent [S] k (ω1)…S k (ω N )] T and [N(ω1)…N(ω N )] T According to the frequency band w = [ω c1 ,...,ω cP ] T The row extraction is performed, where ω cp (1≤p≤P) represents the p-th frequency point, and P is the number of frequency points extracted, so P<N.

[0071] Define the time delay vector of the l-th node as follows: Then the time delay matrix of L nodes is

[0072]

[0073] Neglecting noise, the time-domain vector of the received signal at monitoring node l is:

[0074] x l (t)=F c H diag(h l )s(k) T =F c H D l (H)s(k) T

[0075] in, That is, D l (·) indicates that the l-th row is extracted from matrix H and constructed as a diagonal matrix; For the k-th segment of the radiated signal in discrete frequency domain form, according to the frequency band w = [ω c1 ,...,ω cP ] T Extracted non-zero vectors.

[0076] Define the source matrix as After receiving signals from K radiation sources, the received signal matrix of the l-th monitoring node is:

[0077] Xl =F c H D l (H)S T l = 1, 2, ..., L

[0078] By combining the received signal matrices from L nodes, a three-dimensional matrix can be constructed:

[0079]

[0080] Where X is a three-dimensional matrix of received signals from L monitoring nodes, a parallel factor trilinear model can be constructed:

[0081]

[0082] Where, x n,k,l f is the (n,k,l)th element in the three-dimensional matrix X. n,p It is matrix F c H The (n,p)th element, s k,p h is the (k,p)th element in matrix S. l,p It is the (l, p)th element in matrix H. N is the number of sampling points, L is the number of nodes, and K is the number of signal receiving segments from the radiation source. Multiple diversity spaces exist in the trilinear space of the received signal X, reflecting different diversity methods: spatial diversity, time diversity, and frequency diversity. On the other hand, X... l =F c H D l (H)S T Let l = 1, 2, ..., L be two-dimensional slices of three-dimensional data in spatial directions. Due to the symmetry of the trilinear model, the above equation has the following two reconstructions:

[0083] Y k =HD k (S)(F c H ) T k = 1, 2, ..., K

[0084] Among them, Y k It is the k-th slice in the time direction.

[0085] Z n =SD n (F c H )H T n = 1, 2, ..., N

[0086] Among them, Z n It is the nth slice in the frequency direction.

[0087] Step 3) Using the parallel factor method, the frequency domain time delay estimation matrix is ​​obtained by decomposition and iteration through the trilinear alternating least squares (TALS) method:

[0088] The received signal model in step 2) is re-expressed as

[0089] According to the above formula, the least squares fit is:

[0090] Among them ||·|| F It is the Frobenius norm. Let l = 1, 2, ..., L be a slice of signal containing noise. The S matrix is ​​updated using least squares:

[0091]

[0092] in[·] + This is a generalized inverse. l = 1, 2, ..., L represents the previously obtained H. l The estimate.

[0093] Similarly, the second slice form using the trilinear model: Z n =SD n (F c H )H T n = 1, 2, ..., N, can be rewritten as:

[0094]

[0095] According to the above equation, the least squares fit is:

[0096] in n = 1, 2, ..., N, representing a slice of signal containing noise. The H matrix is ​​updated using least squares:

[0097]

[0098] in This indicates a previously obtained estimate of S.

[0099] Matrix S and H are updated using the least squares method until the algorithm converges, yielding the source matrix S and the frequency domain time delay matrix H. The TALS algorithm used in this patent utilizes a known matrix F. c H This will accelerate the convergence speed.

[0100] Given Z n = SD n (F c H )H T (n = 1, 2,..., N), where F c H and H are Vandermonde matrices composed of non - zero sequences. If N + L + k S ≥ 2P + 2, then F c H , H and S are distinguishable (up to column - modulus exchanges and scaling transformations).

[0101] From the previous analysis, k S is full - rank and has full k - rank, that is, k S = r S = min(K, P). Therefore, the above formula becomes N + L + min(K, P) ≥ 2P + 2. In practical situations, since K ≥ P, the identifiability condition is N + L ≥ P + 2. So when N + L ≥ P + 2 is satisfied, the uniqueness of the parallel factor model decomposition can be guaranteed.

[0102] Any that constitutes X l , l = 1, 2,..., L has the following relationship with H and S: where П is a column - ambiguity matrix, and the scale - ambiguity matrices Δ1 and Δ2 are both diagonal matrices and satisfy Δ1Δ2 = I<着 P . Since the IDFT matrix F c H is known, the column - exchange ambiguity (permutation ambiguity) can be solved by the known matrix F c <00002!8>. In the obtained result[[ID=;5]] the column ambiguity does not affect the estimation accuracy, and the scale ambiguity inherent in the trilinear decomposition can be eliminated by normalization methods.

[0105] After estimating the time - delay matrix using the trilinear decomposition method, the time - delay difference is estimated using the Vandermonde characteristics of this matrix. Let a certain column of be a p , p = 1, 2,..., P. First, normalize the time - delay vector a p , p = 1, 2,..., P so that its first term is 1; then the column vectors a in the obtained pp ​​​​​​p = 2, ..., P divided by the previous column; let the result be... a certain behavior b l l = 2, 3, ..., L, take -angle(b) l ), obtain the product of frequency difference and time delay difference; process the matrix according to the above method. After adjustment, the time delay difference between each monitoring node is estimated using the least squares method. The specific process is as follows:

[0106]

[0107] in Delay matrix Column vectors.

[0108] Step 4) Eliminate scale ambiguity using a normalization method and estimate time delay using the Vandermonde eigenvalue of the time delay matrix:

[0109] First, the scale ambiguity problem is solved using a normalization method. p All elements in p = 1, 2, ..., P are divided by the first element of that column. At the same time, the time delay estimate is changed to the time delay difference estimate relative to the first monitoring node.

[0110]

[0111] in This is the time delay estimation matrix after row vector normalization. for The column vector. Then... Column vector a in the matrix pp The dot division between p = 2, ..., P and the corresponding previous column vector eliminates the attenuation factor and prevents phase ambiguity.

[0112]

[0113] in This is the time delay estimation matrix after column vector normalization. for The row vectors are obtained by taking their modulus:

[0114]

[0115] Where q is a vector containing the frequency difference and time delay difference. Using the known inverse Fourier transform matrix F... c H Solve the column exchange ambiguity (permutation ambiguity) problem to obtain the corresponding frequency band. Then calculate the corresponding frequency difference w. c =[0,ω c2 -ω c1,...,ω cP -ω cP-1 ] T The time delay difference estimate can then be obtained:

[0116]

[0117] in, This represents the time delay difference between the l-th monitoring node and the 1st monitoring node (reference node).

[0118] Figure 3 This is a comparison chart of the estimated latency and the actual latency in a simulated scenario using the method described in this invention. The simulation parameters are set as follows: the monitored base station locations are node 1: p1 = [100, 200, 600]. T Node 2: p2 = [500, 1000, 100] T Node 3: p3 = [1000, 0, 1500] T Node 4: p4 = [1500, 1500, 2000] T Here, node 1 is chosen as the reference node, and the location of the radiation source is u = [35, 170, 85]. T All units used are in meters (m). The signal-to-noise ratio (SNR) of the radiation source signal is 10 dB, and the sampling frequency is f. s =10MHz, signal center frequency f c =5MHz, signal duration T = 0.5ms, number of sampling points N = 5000, number of extracted frequency points P = 4, and K = 100 radiation source signals are collected. As can be seen from the figure, the difference between the time delay difference obtained using the method of this invention and the actual time delay difference is very small, indicating that this method can obtain an accurate time delay difference.

[0119] Figure 4 This paper compares the time delay estimation performance of the method described in this invention with that of a generalized cross-correlation method based on phase transformation under different signal-to-noise ratios. Simulation parameters are set as follows: the monitoring base station locations are node 1: p1 = [100, 200, 600]. T Node 2: p2 = [500, 1000, 100] T Node 3: p3 = [1000, 0, 1500] T Node 4: p4 = [1500, 1500, 2000] T Here, node 1 is chosen as the reference node, and the location of the radiation source is u = [35, 170, 85]. T All units used are in meters (m). The sampling frequency f of the radiation source signal. s =30MHz, signal center frequency f c=10MHz, signal bandwidth B =0.9MHz, signal duration T =0.1ms, number of sampling points N =3000, number of extracted frequency points P =30, collecting K =100 radiation source signals. Figure 4 It can be seen that, under the same conditions, as the signal-to-noise ratio increases, the time delay estimation error of the method of the present invention decreases and remains small, and it has significantly better estimation performance.

[0120] Figure 5 This is a scene diagram of the distribution of radiation sources and monitoring nodes in the actual test environment described in this invention. There are four monitoring nodes distributed in different locations, and the radiation source is placed on the playground.

[0121] Figure 6 This is a comparison chart of the estimated and actual time delay values ​​of the method described in this invention under actual test scenarios. A Cartesian coordinate system is established based on the known latitude and longitude positions of four monitoring nodes, with node three as the origin. Therefore, the node positions are: Node one: q1 = [60.1263, -32.0600], Node two: q2 = [83.7238, -127.4608], Node three: q3 = [0,0], Node four: q4 = [65.5065, -220.4126]. Node three is selected as the reference node, and the position of the radiation source is u = [-37.0008, -75.4746]. The sampling frequency of the radiation source signal is f. s =125MHz, center frequency f c =716MHz, signal bandwidth B=15MHz, number of sampling points N=1024, number of extracted frequency points P=14, collecting K=27 radiation source signals. Figure 6 As can be seen, by processing the measured data, the time delay value estimated by the method of the present invention is very close to the actual time delay value, further verifying the effectiveness and practicality of the method of the present invention.

[0122] It will be understood by those skilled in the art that, unless otherwise defined, all terms used herein (including technical and scientific terms) have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. It should also be understood that terms such as those defined in general dictionaries should be understood to have the same meaning as in the context of the prior art, and should not be interpreted in an idealized or overly formal sense unless defined as herein.

[0123] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A high-precision time delay estimation method based on parallel factors, characterized in that, The method includes the following steps: Step 1) Receive radiation source signals at L known monitoring nodes and perform Fourier transform on the received signals, where L is a preset threshold for the number of base stations. Step 1.1): Receive radiation source signals at L known monitoring nodes. The received signal x at the l-th monitoring node... l (t)=α l s k (t-τ l )+n l (t), where 1≤l≤L, s k (t) represents the radiation source signal received in the k-th segment, 1≤k≤K, where K is a preset threshold for the number of radiation source signal segments, α l and τ l Let n represent the path attenuation factor and time delay of the radiation source signal reaching the l-th monitoring node, respectively. l (t) represents mutually independent zero-mean additive white Gaussian noise, and the signal and noise are uncorrelated. Step 1.2): Perform a Fourier transform on the received signal to obtain the time-domain output matrices of L monitoring nodes, where the time-domain output matrix of the l-th monitoring node is in the form of... In the formula, x l (t) represents the received signal of the l-th monitoring node, F H Let x be the inverse Fourier transform matrix. l (ω) is x l The Fourier transform form of (t), α l and τ l Represent the path attenuation factor and time delay of the radiation source signal reaching the l-th monitoring node, respectively. k (ω1)...S k (ω N )] T and [N(ω1)...N(ω N )] T They represent s respectively k (t) and n l The Fourier transform form of (t), where ω n For the nth frequency point, 1≤n≤N, where N is the number of sampling points; Step 2) Extract non-zero frequency band signals and fuse the frequency domain data of each monitoring node to construct a trilinear model; Step 2.1), the time-domain output matrices of the L monitoring nodes are processed according to the frequency band w = [ω c1 ,...,ω cP ] T Extract: in, It is the inverse Fourier transform matrix F H According to the frequency band w = [ω c1 ,...,ω cP ] T The column extraction performed, [S k (ω c1 ...S k (ω cP )] T and [N(ω c1 ...N(ω) cP )] T They represent [S] k (ω1)...S k (ω N )] T and [N(ω1)...N(ω N )] T According to the frequency band w = [ω c1 ,...,ω cP ] T The row extraction is performed, where ω cp Let P be the p-th frequency point, 1≤p≤P, where P is the number of frequency points extracted. Step 2.2), define the time delay vector of the l-th node as... Then the time delay matrix of L nodes is Ignoring noise, the received signal of monitoring node l is X. l =F c H diag(h l )S T =F c H D l (H)S T in, To collect the source matrix of K-band radiation source signals, For the k-th segment of the radiated signal, according to the frequency band w = [ω c1 ,...,ω cP ] T Extracted non-zero vector; X l Let l be the received signal matrix of the l-th monitoring node. It is the inverse Fourier transform matrix F H According to the frequency band w = [ω c1 ,...,ω cP ] T Column extraction was performed; by combining the received signal matrices from L nodes, a three-dimensional matrix was constructed: Where X is a three-dimensional matrix of received signals containing L monitoring nodes, a parallel factor trilinear model is constructed as follows: Where, x n,k,l f is the (n,k,l)th element in the three-dimensional matrix X. n,p It is matrix F c H The (n,p)th element, s k,p h is the (k,p)th element in matrix S. l,p It is the (l, p)th element in matrix H; Step 3) Using the parallel factor method, the frequency domain time delay estimation matrix is ​​obtained by decomposing and iterating through the trilinear alternating least squares method. Step 4) Eliminate scale ambiguity by normalization and estimate time delay using the Vandermonde feature of the time delay matrix.

2. The high-precision time delay estimation method based on parallel factor according to claim 1, characterized in that, The frequency domain time delay estimation matrix obtained in step 3) is: in, This is the time delay estimation matrix. for Column vectors.

3. The high-precision time delay estimation method based on parallel factor according to claim 2, characterized in that, The specific steps of step 4) are as follows: Step 4.1): Solve the scale ambiguity problem using a normalization method, and set a... p All elements in p = 1, 2, ..., P are divided by the first element of the column. At the same time, the time delay estimate is also changed to the time delay difference estimate relative to the first monitoring node. in This is the time delay estimation matrix after row vector normalization. for Column vectors; Step 4.2), will Column vector a in the matrix pp Divide p = 2, ..., P by the dot notation between the previous column to eliminate the attenuation factor and prevent phase blurring. in, This is the time delay estimation matrix after column vector normalization. for The row vectors are obtained by taking their modulus: Where q is a vector containing frequency difference and time delay difference; Step 4.3), using the known inverse Fourier transform matrix F c H Solve the column swap ambiguity problem to obtain the corresponding frequency band. Then calculate the corresponding frequency difference w. c =[0,ω c2 -ω c1 ,...,ω cP -ω cP-1 ] T The time delay difference estimate can then be obtained: in, This represents the time delay difference between the l-th monitoring node and the first monitoring node.

Images