A high-precision time difference estimation method based on cross-spectrum and trilinear decomposition

Abstract

The invention discloses a high-precision time difference estimation method combining cross-spectrum and trilinear decomposition. First, L monitoring nodes are used to receive radiation signals, and the cross-power spectrum of the signals received by all nodes and the reference node is calculated; secondly, non-zero frequency band signals are extracted, and the frequency domain cross-spectrum data between each monitoring node and the reference node are integrated to construct a trilinear model; then, a parallel factor method is used to perform decomposition iteration through the trilinear alternating least squares method to obtain a frequency domain time difference estimation matrix; finally, a normalization method is used to eliminate scale ambiguity, and the time difference matrix is ​​analyzed to obtain a high-precision time difference estimation. The invention can break through the limitation of sampling frequency and obtain a time difference estimation with a higher accuracy than the sampling interval, thereby effectively improving the positioning accuracy.

Description

Technical Field

[0001] The invention relates to the technical field of passive time difference positioning, and in particular to a high-precision time difference estimation method combining cross-spectrum and trilinear decomposition. Background Art

[0002] Passive positioning technology refers to the technology that obtains the target's position by receiving the radiation information emitted by the target itself. Passive positioning technology can be divided into single-node and multi-node passive positioning according to the number of nodes involved in positioning. Multi-node passive positioning uses some information (delay difference, angle, etc.) to determine multiple positioning curves (or surfaces), and the target can be located by solving the intersection of the curves (or surfaces). The structure of the multi-node passive positioning system is relatively complex, and each receiver needs to synchronize data transmission and fusion, but it has very high positioning accuracy and is the most important research direction in passive positioning technology.

[0003] The most widely used passive positioning method is the Time Difference of Arrival (TDOA) positioning method, which measures the absolute time difference between the radiation source signal reaching two nodes to obtain a hyperbola with the two nodes as the focus. The target position can be obtained by solving the intersection of the two hyperbolas. This method only requires synchronization between each node, and does not require clock synchronization between the radiation source and the node. The requirements for clock synchronization are further reduced, and the equipment cost problem is optimized. In addition, the arrival time difference algorithm has certain advantages in positioning accuracy and reduces the complexity of the algorithm. These advantages have made TDOA positioning gain attention and widely used in civil and military fields.

[0004] The focus of TDOA positioning is time difference estimation, that is, calculating the time difference between the radiation source signal arriving at each node and the reference node. As a commonly used solution method for time difference estimation, the cross-correlation method has a wide range of application scenarios due to its simple principle and easy implementation. This method obtains the cross-correlation function by processing the radiation source signals received by different nodes, and finds the location of the cross-correlation function peak to obtain the corresponding delay difference. However, the cross-correlation method not only has high requirements on the sampling frequency, but is also sensitive to the noise in the signal. The time difference estimation capability is limited, and the estimation accuracy cannot exceed the sampling interval. Existing solutions include fitting or interpolation, but the increase in computational complexity puts higher requirements on the equipment hardware, which is not conducive to engineering implementation. Therefore, improving the accuracy of time difference estimation on the basis of low complexity has important practical significance and application value. Summary of the invention

[0005] The technical problem to be solved by the present invention is to provide a high-precision time difference estimation method combining cross-spectrum and trilinear decomposition in view of the defects involved in the background technology.

[0006] The present invention adopts the following technical solutions to solve the above technical problems:

[0007] A high-precision time difference estimation method combining cross-spectrum and trilinear decomposition, the method comprising the following steps:

[0008] Step 1), receiving radiation source signals at L monitoring nodes, setting the mth node as the reference node, 1≤m≤L, calculating the cross-power spectrum of the signals received by all nodes and the reference node, L is the preset base station number threshold;

[0009] Step 2), extracting non-zero frequency band signals, and fusing the frequency domain cross-spectrum data between each monitoring node and the reference node to construct a trilinear model;

[0010] Step 3), using the parallel factor method, performing decomposition iteration through the trilinear alternating least squares method to obtain the frequency domain time difference estimation matrix;

[0011] Step 4), use the normalization method to eliminate scale ambiguity and perform time difference estimation based on the time difference matrix.

[0012] As a further optimization scheme of the high-precision time difference estimation method combining cross-spectrum and trilinear decomposition of the present invention, the step 1) includes the following steps:

[0013] Step 1.1), receive the radiation source signal at L monitoring nodes with known positions, and the received signal x of the lth 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 kth segment, 1≤k≤K, K is the preset radiation source signal segment number threshold, α l and τ l They represent the path attenuation factor and delay of the radiation source signal reaching the lth monitoring node, n l (t) is an independent zero-mean additive white Gaussian noise, and the radiation source signal and the white Gaussian noise are uncorrelated with each other;

[0014] Step 1.2), calculate the cross power spectrum g of the received signals of all nodes and the reference node l (ω), 1≤l≤L, and perform Fourier transform to obtain the cross-correlation function, where the cross-correlation function between the lth monitoring node and the reference node is:

[0015]

[0016] In the formula, r l(t) is the cross-correlation function between the lth monitoring node and the reference node, F H is the inverse Fourier transform matrix, g l (ω) is the mutual power spectrum between the lth monitoring node and the reference node, α l and τ ml They represent the path attenuation factor of the radiation source signal reaching the lth monitoring node and the delay difference with the reference node, s k (ω)=[G ss (ω1),...,G ss (ω N )] T Represents the radiation source signal s k (t) is the power spectrum function, where ω n is the nth frequency point, 1≤n≤N, N is the number of sampling points.

[0017] As a further optimization scheme of the high-precision time difference estimation method combining cross-spectrum and trilinear decomposition of the present invention, the specific steps of step 2) are:

[0018] Step 2.1), the cross-correlation function between the lth monitoring node and the reference node is calculated according to the frequency band w = [ω c1 ,...,ω cP ] T To extract:

[0019]

[0020] in, is the inverse Fourier transform matrix F H According to the frequency band w = [ω c1 ,...,ω cP ] T Column extraction performed, s k (ω c ) means s k (ω) According to the frequency band w = [ω c1 ,...,ω cP ] T The element extraction is performed, where ω cp is the pth frequency point, 1≤p≤P, P is the number of extracted frequency points;

[0021] Step 2.2), define the time difference vector between the lth node and the reference node as Then the time difference matrix of L nodes is

[0022]

[0023] When the noise effect is not considered, we have

[0024] R l =Fc H diag(h cl )S T =F c H D l (H)S T ,l=1,2,...,L

[0025] Among them, R l is the cross-correlation function matrix between the lth monitoring node and the reference node; To collect the power spectrum matrix of K-segment radiation source signals; combine the cross-correlation function matrix of L nodes to construct a three-dimensional model:

[0026]

[0027] Among them, R is a three-dimensional matrix of the cross-correlation function containing L monitoring nodes and reference nodes, and a parallel factor trilinear model is constructed:

[0028]

[0029] Among them, r n,k,l is the (n,k,l)th element in the three-dimensional matrix R, f n,p is the matrix F c H The (n,p)th element in s k,p is the (k,p)th element in the matrix S, h l,p is the (l,p)th element in the matrix H.

[0030] As a further optimization scheme of the high-precision time difference estimation method combining cross-spectrum and trilinear decomposition of the present invention, the frequency domain time difference estimation matrix obtained in step 3) is:

[0031]

[0032] in, is the time difference estimation matrix, Δ is The existing scale fuzzy matrix can be eliminated by normalization later. l for The row vector of .

[0033] As a further optimization scheme of the high-precision time difference estimation method combining cross-spectrum and trilinear decomposition of the present invention, the specific steps of step 4) are:

[0034] Step 4.1), Each row of a l Divide the first row a1 by 1≤l≤L, and use the normalization method to solve the scale ambiguity problem;

[0035]

[0036] in is the normalized time difference estimation matrix, b p is a matrix Column vector of , 1≤p≤P;

[0037] Step 4.2), Column vectors in a matrix With the corresponding previous column Point division, 2≤p1≤P, 1≤p2≤P-1, that is, eliminating the attenuation factor and preventing phase ambiguity:

[0038]

[0039] in, To eliminate the time difference estimation matrix after the attenuation factor is eliminated, is a matrix The row vector of , 2≤l2≤L, is obtained by taking its modulus:

[0040]

[0041] Among them, q l is a vector containing frequency difference and delay difference;

[0042] Step 4.3), using the known inverse Fourier transform matrix F c H Solve the column exchange ambiguity problem and get the corresponding frequency band Then find the corresponding frequency difference You can get

[0043]

[0044] in, is the estimated value of the delay difference from the lth monitoring node to the reference node (the first monitoring node), 1≤l≤L, 1≤m≤L.

[0045] Compared with the prior art, the present invention adopts the above technical solution and has the following technical effects:

[0046] The present invention provides a high-precision time difference estimation method combining cross-spectrum and trilinear decomposition. By fusing the frequency domain cross-spectrum data between each monitoring node and a reference node, a trilinear model is constructed, and the time difference is estimated using trilinear decomposition. This solves the problems of limited accuracy and large amount of calculation in traditional time difference estimation methods, improves the accuracy of time difference estimation, and provides a guarantee for the accuracy of passive positioning. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1A flow chart of a high-precision time difference estimation method combining cross-spectrum and parallel factor provided by the present invention;

[0048] Figure 2 It is a source positioning scene diagram based on TDOA described in the present invention;

[0049] Figure 3 A time difference estimation diagram in a simulation scenario of the method of the present invention;

[0050] Figure 4 This is a comparison diagram of the time difference estimation performance between the method of the present invention and the generalized cross-correlation method under different signal-to-noise ratios;

[0051] Figure 5 A diagram of information source and node distribution scenarios in an actual test environment described in the present invention;

[0052] Figure 6 This is a time difference estimation diagram under the actual measurement scenario of the method described in the present invention. DETAILED DESCRIPTION

[0053] The technical solution of the present invention is further described in detail below in conjunction with the accompanying drawings:

[0054] The present invention can be implemented in many different forms and should not be considered to be limited to the embodiments described herein. On the contrary, these embodiments are provided to make the present invention thorough and complete, and will fully express the scope of the present invention to those skilled in the art. In the accompanying drawings, components are enlarged for clarity.

[0055] The detailed process of the high-precision time difference estimation method of combining cross-spectrum and trilinear decomposition provided by the present invention is as follows: Figure 1 As shown in the figure, L monitoring nodes are used to receive the radiation signal, and the cross-power spectrum of the received signals of all nodes and the reference node is calculated; the non-zero frequency band signal is extracted, and the frequency domain cross-spectrum data between each monitoring node and the reference node are fused to construct a trilinear model; the parallel factor method is used to perform decomposition and iteration through the trilinear alternating least squares method to obtain the frequency domain time difference estimation matrix; the scale ambiguity is eliminated through the normalization method, and the time difference estimation result is obtained based on the time difference matrix. The specific implementation is as follows:

[0056] Step 1), receive the radiation source signal at L monitoring nodes with known positions, let the mth node be the reference node, 1≤m≤L, and calculate the cross power spectrum of the received signals of all nodes and the reference node:

[0057] Consider Figure 2 The TDOA positioning scenario shown. Assume that there is an unknown position u=[x,y] in space T The radiation source signal s(t), and the L monitoring nodes are located at p l =[x l ,yl ] T ,1≤l≤L, select the mth monitoring node as the reference node, then the signal received by the monitoring node from the radiation source can be expressed as

[0058]

[0059] Among them, x1(t) and x l (t) are the received signals of the reference node m (1≤m≤L) and the monitoring node l (1≤l≤L), τ m With τ l denote the delay of the received signal from the radiation source to the reference node and the monitoring node l, respectively, and α m With α l represents the corresponding random path attenuation factor, s k (t)(1≤k≤K) represents the unknown radiation source signal received in the kth segment, K is the preset threshold of the number of radiation source receiving segments, n m (t) and n l (t) are independent zero-mean additive white Gaussian noises, and the signal and the noise are uncorrelated. To simplify the analysis, we assume that α m =1. The cross-correlation function of the two received signals is calculated as

[0060]

[0061] in represents the cross-correlation function between the reference node and the monitoring node l, τ ml =τ l -τ m , which represents the delay difference between the received signal from the radiation source to the monitoring node and the reference node. According to the Wiener-Khinchin theorem, the cross-correlation function of two signals and their cross-power spectrum density function are a pair of Fourier transform pairs. Therefore, the cross-correlation function is subjected to discrete Fourier transform (DFT) to obtain the cross-power spectrum of the two received signals.

[0062]

[0063] in is the cross power spectrum between the reference node and the monitoring node l, G ss (ω n ) represents the radiation source signal s k (t) power spectrum, where ω n is the nth frequency point, 1≤n≤N, N is the number of sampling points. s k (ω)=[G ss (ω1),...,Gss (ω N )] T , then the above formula becomes

[0064] g l (ω)=diag(h l )s k (ω)

[0065] make is the DFT matrix, then the IDFT matrix is ​​F H , so the cross-correlation function between the monitoring node l and the reference node is

[0066] r l (t) = F H g l (ω)=F H diag(h l )s k (ω)

[0067] Among them, r l (t) is the vector form of the cross-correlation function between the lth monitoring node and the reference node.

[0068] Step 2), extract the non-zero frequency band signal, and fuse the frequency domain cross-spectrum data between each monitoring node and the reference node to build a trilinear model:

[0069] Since the sampling bandwidth of the signal is greater than the signal bandwidth, s k Many elements in (ω) are zero or close to zero. Assuming that they are in the pre-set w = [ω c1 ,...,ω cP ] T If it is non-zero in the frequency band, the above formula is modified to

[0070]

[0071] in, is the inverse Fourier transform matrix F H According to the frequency band w = [ω c1 ,...,ω cP ] T The column extraction is performed. s k (ω c )=[G ss (ω c1 ),...,G ss (ω cP )] T Yes l ,s k (ω) respectively according to the frequency band w=[ω c1 ,...,ω cP ]T The element extraction is performed, where ω cp (1≤p≤P) is the pth frequency point, P is the number of extracted frequency points, so P<N.

[0072] Then the time difference matrix of L monitoring nodes is

[0073]

[0074] The cross-correlation function matrix between the monitoring node l and the reference node can be expressed as

[0075]

[0076] Where R l is the lth slice in the spatial direction; D l (·) means selecting the lth row from the matrix H and constructing a diagonal matrix; Represents the power spectrum matrix of the K-band radiation source signal. When considering the influence of noise, Where N l represents the noise of the lth slice.

[0077] By combining the cross-correlation function matrix of L nodes and the reference node, a three-dimensional matrix can be constructed:

[0078]

[0079] Among them, R is a three-dimensional matrix of the cross-correlation function containing L monitoring nodes and reference nodes, and a parallel factor trilinear model is constructed:

[0080]

[0081] Among them, r n,k,l is the (n,k,l)th element in the three-dimensional matrix R, f n,p is the matrix F c H The (n,p)th element in s k,p is the (k,p)th element in the matrix S, h l,p is the (l,p)th element in the matrix H. N is the number of sampling points, L is the number of nodes, and K is the number of receiving segments of the radiation source signal. There are multiple diversity spaces in the trilinear space of the received signal R, which reflect different diversity modes: spatial diversity, time diversity, and frequency diversity. On the other hand, R l =F c H D l (H)S T ,l=1,2,...,L is a two-dimensional slice of three-dimensional data in the spatial direction. Due to the symmetry of the trilinear model, the above formula has the following two reconstructions:

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

[0083] Among them, Y k is the kth slice in the time direction.

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

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

[0086] Step 3), using the parallel factor method, the trilinear alternating least squares (TALS) method is used to perform decomposition iteration to obtain the frequency domain time difference estimation matrix:

[0087] The essential characteristic of trilinear decomposition is its uniqueness. Let’s first discuss the identifiability of the algorithm. Among them, F c H , H and S are matrices with full k-rank, then F c H , the sufficient condition for H and S to be distinguishable (column module exchange and scale transformation) is min(N,P)+min(L,P)+min(K,P)≥2P+2. In actual situations, since N≥P,L≤P,K≥P, the distinguishability condition becomes L≥2. Therefore, when L≥2 is satisfied, the trilinear model decomposition is guaranteed to be unique.

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

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

[0090] where ||·|| F is the Frobenius norm. is a noisy signal slice. Use least squares to update the S matrix:

[0091]

[0092] in[·] + is the generalized inverse. Indicates the previously obtained H l The estimated form.

[0093] The second slice form of the trilinear model is also used: Z n =SD n (F c H )H T ,n=1,2,...,N,re-expressed as:

[0094]

[0095] According to the above least squares fit,

[0096] in is a noisy signal slice. Use least squares to update the H matrix:

[0097]

[0098] in It represents the estimated form of S that has been obtained previously. From this, the frequency domain time difference estimation matrix can be obtained

[0099] The matrices S and H are updated by the least squares method until the algorithm converges. This method is used to obtain the source matrix and frequency domain delay matrix The TALS algorithm is optimal when the noise is Gaussian independent white noise. The TALS algorithm has the following advantages: it is easy to implement, can guarantee convergence and can be easily extended to high-order data. The TALS algorithm used in this patent uses the known DFT matrix F c H , which will speed up the convergence.

[0100] Step 4) Use the normalization method to eliminate scale ambiguity and perform time difference estimation based on the time difference matrix:

[0101] Any composition R l ,l=1,2,...,L The relationship between H and S is as follows:

[0102]

[0103] Where Π is the column fuzzy matrix, and the scale fuzzy matrices Δ1 and Δ2 are both diagonal matrices, and satisfy Δ1Δ2=I P Since the IDFT matrix F c H It is known that column exchange fuzziness (permutation fuzziness) can be obtained by using the known matrix F cH can be solved. In the trilinear decomposition, we can further find

[0104]

[0105] Π can be obtained through F c H and After eliminating the permutation ambiguity, the frequency domain time difference estimation matrix becomes

[0106]

[0107] Where Δ is the scale fuzzy matrix, a1(1≤l≤L) is the matrix The row vector of . Using the normalization method to solve the scale ambiguity problem, Each row of the matrix a l Divide (1≤l≤L) by the first row a1 to get the matrix

[0108]

[0109] in is the normalized time difference estimation matrix, b p (1≤p≤P) is a matrix By Column vector of a matrix With the corresponding previous column By dividing 2≤p1≤P, 1≤p2≤P-1, the influence of fading factor and phase ambiguity can be eliminated, and the matrix

[0110]

[0111] in To eliminate the time difference estimation matrix after the attenuation factor is eliminated, is a matrix The row vector of , 2≤l2≤L, can be obtained by taking its phase

[0112]

[0113] Among them, q l is a vector containing frequency difference and delay difference. c H And trilinear decomposition The corresponding frequency point position can be found, and then the permutation fuzzy matrix π is used to obtain By The element ω in the vector cp (2≤p≤P) and the corresponding previous element ω cp (1≤p≤P-1) Subtract and calculate the frequency difference vector The delay difference can be obtained:

[0114]

[0115] in is the estimated value of the delay difference from the lth monitoring node to the reference node, 1≤l≤L, 1≤m≤L.

[0116] Figure 3 This is a comparison diagram of the estimated time difference and the actual time difference in the simulation scenario of the method of the present invention. The simulation parameters are set as follows: the locations of the monitoring base stations are node 1: p1 = [0, 0] T , Node 2: p2 = [200, 600] T , Node 3: p3 = [700,0] T , here we select node 1 as the reference node, and the position of the radiation source is u=[100,100] T , all units used are m (meter). The signal-to-noise ratio SNR of the radiation source signal is 10dB, and the sampling frequency is f s =10MHz, signal center frequency f c =1MHz, signal time width T=0.1ms, number of sampling points N=1000, number of extracted frequency points P=6, and K=200 radiation source signals are collected. It can be seen from the figure that the difference between the delay difference calculated by the method of the present invention and the real delay difference is very small, which shows that the method can obtain accurate delay difference.

[0117] Figure 4 The performance of the time difference estimation of the method of the present invention and the generalized cross-correlation method based on phase conversion under different signal-to-noise ratio conditions is compared. The simulation parameters are set as follows: The locations of the monitoring base stations are node 1: p1 = [0, 0] T , Node 2: p2 = [200, 600] T , Node 3: p3 = [700,0] T , here we select node 1 as the reference node, and the position of the radiation source is u=[100,100] T , all units used are m (meter). Radiation source signal sampling frequency f s =10MHz, signal center frequency f c =1MHz, signal bandwidth B = 0.3MHz, signal time width T = 0.1ms, number of sampling points N = 1000, number of extracted frequency points P = 6, and K = 200 radiation source signals are collected. Figure 4 It can be seen that, under the same other conditions, as the signal-to-noise ratio increases, the time difference estimation error of the method of the present invention gradually decreases and remains small, and obviously has better estimation performance.

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

[0119] Figure 6 This is a comparison chart of the estimated time difference and the real time difference of the method of the present invention in the actual test scenario. A Cartesian coordinate system is established based on the known longitude and latitude position information of the four monitoring nodes, with node three as the coordinate origin. Therefore, the positions of the nodes are node one: q1 = [0, 0], node two: q2 = [60.1263, -32.0600], node three: q3 = [83.7238, -127.4608], node four: q4 = [65.5065, -220.4126]. Node one 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 = 9, and K = 27 radiation source signals are collected. Figure 6 It can be seen that by processing the measured data, the time difference value estimated by the method of the present invention is very close to the real time difference value, which further verifies the effectiveness and practicality of the method of the present invention.

[0120] It will be understood by those skilled in the art that, unless otherwise defined, all terms (including technical and scientific terms) used herein have the same meaning as those generally understood by those skilled in the art in the art to which the present invention belongs. It should also be understood that terms such as those defined in common dictionaries should be understood to have meanings consistent with the meanings in the context of the prior art, and will not be interpreted with idealized or overly formal meanings unless defined as herein.

[0121] The specific implementation methods described above further illustrate the objectives, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above description is only a specific implementation method 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 in the scope of protection of the present invention.

Claims

1. A high-precision time difference estimation method combining cross-spectrum and trilinear decomposition, characterized in that: The following steps are involved: Step 1), receiving the radiation source signal at L monitoring nodes with known positions, setting the mth node as the reference node, 1≤m≤L, and calculating the cross-power spectrum of the signal received by each monitoring node and the reference node; Step 1.1), receive the radiation source signal at L monitoring nodes with known positions, and the received signal x of the lth monitoring node l (t) = α l s k (t-τ l )+n l (t), where 1≤l≤L, s k (t-τ l ) represents the kth segment of the radiation source signal received at the lth monitoring node, 1≤k≤K, K is the preset threshold of the number of segments of the received radiation source signal, α l and τ l They represent the path attenuation factor and delay of the radiation source signal reaching the lth monitoring node, n l (t) is an independent zero-mean additive white Gaussian noise, and the radiation source signal and the white Gaussian noise are uncorrelated with each other; Step 1.2), calculate the cross power spectrum g of the received signals of all monitoring nodes and the reference node l (ω), 1≤l≤L, and perform Fourier transform to obtain the cross-correlation function, where the cross-correlation function between the lth monitoring node and the reference node is: In the formula, r l (t) is the cross-correlation function between the lth monitoring node and the reference node, F H is the inverse Fourier transform matrix, g l (ω) is the mutual power spectrum between the lth monitoring node and the reference node, α l and τ ml They represent the path attenuation factor of the radiation source signal reaching the lth monitoring node and the delay difference with the reference node, s k (ω)=[G ss (ω1),...,G ss (ω N )] T Represents the radiation source signal s k (t) power spectrum function, where ω n is the nth frequency point, G ss (ω n ) is s k The autopower spectrum of (t), 1≤n≤N, N is the number of sampling points; Step 2), extracting non-zero frequency band signals, and fusing the frequency domain cross-spectrum data between each monitoring node and the reference node to construct a trilinear model; Step 2.1), the cross-correlation function between the lth monitoring node and the reference node is calculated according to the preset frequency band w = [ω c1 ,...,ω cP ] T To extract: in, is the inverse Fourier transform matrix F H According to the frequency band w = [ω c1 ,...,ω cP ] T Column extraction performed, s k (ω c ) means s k (ω) According to the frequency band w = [ω c1 ,...,ω cP ] T The element extraction is performed, where ω cp is the pth frequency point, 1≤p≤P, P is the number of extracted frequency points; Step 2.2), define the time difference vector between the lth node and the reference node as Then the time difference matrix of L nodes is When the influence of noise is not considered, we have R l =F c H diag(h cl )S T =F c H D l (H)S T ,l=1,2,...,L Among them, R l is the cross-correlation function matrix between the lth monitoring node and the reference node; To collect the power spectrum matrix of K-segment radiation source signals; combine the cross-correlation function matrix of L nodes to construct a three-dimensional model: Among them, R is a three-dimensional matrix of the cross-correlation function containing L monitoring nodes and reference nodes, and a parallel factor trilinear model is constructed: Among them, r n,k,l is the (n,k,l)th element in the three-dimensional matrix R, f n,p is the matrix F c H The (n,p)th element in s k,p is the (k,p)th element in the matrix S, h l,p is the (l,p)th element in the matrix H; Step 3), using the parallel factor method, performing decomposition iteration through the trilinear alternating least squares method to obtain the frequency domain time difference estimation matrix; Step 4), use the normalization method to eliminate scale ambiguity and perform time difference estimation based on the time difference matrix.

2. The high-precision time difference estimation method of joint cross-spectrum and trilinear decomposition according to claim 1 is characterized in that: The frequency domain time difference estimation matrix obtained in step 3) is: in, is the time difference estimation matrix, Δ is The scale fuzzy matrix that comes with it, a l for The row vector of .

3. The high-precision time difference estimation method of joint cross-spectrum and trilinear decomposition according to claim 2 is characterized in that: The specific steps of step 4) are: Step 4.1) Each row of a l Divide the first row a1 by 1≤l≤L, and use the normalization method to solve the scale ambiguity problem; in is the normalized time difference estimation matrix, b p is a matrix Column vector of , 1≤p≤P; Step 4.2), Column vectors in a matrix With the corresponding previous column Point division, 2≤p1≤P, 1≤p2≤P-1, that is, eliminating the attenuation factor and preventing phase ambiguity: in, is the time difference estimation matrix after eliminating the attenuation factor; is a matrix The row vector of , 2≤l2≤L, is obtained by taking its modulus: Among them, q l is a vector containing frequency difference and delay difference; Step 4.3), using the known inverse Fourier transform matrix F c H Solve the column exchange ambiguity problem and get the corresponding frequency band Then find the corresponding frequency difference Obtain in, is the estimated value of the delay difference from the lth monitoring node to the reference node, 1≤l≤L, 1≤m≤L.

Images