Multi-Station Long Baseline Interferometry Passive Localization Method

Abstract

This paper proposes a multi-station long baseline interferometry passive localization method. Multiple stations simultaneously receive target signals, and through inter-station differential interferometry processing, effective signal de-modulation is achieved. Through the relative motion between observation stations and the mission area, two-dimensional sampling of the target radiation signal spatial spectrum is realized, enabling inversion of the target distribution and obtaining two-dimensional ultra-high-precision positioning results. This paper models the multi-station long baseline interferometry scenario, theoretically derives the correspondence between the target distribution map and the bi-static interferometry data, presents the processing flow of the multi-station long baseline interferometry passive localization method, analyzes the positioning accuracy of this method, and verifies through simulation that the method possesses high-resolution positioning capability.

Graphical abstract

Keywords

multi-station / long baseline interferometry / passive localization / high positioning accuracy / two-dimensional moving virtual aperture

1 Introduction

Single-station synthetic aperture passive localization technology is a one-dimensional interferometric positioning method. It achieves ultra-long motion virtual aperture through azimuth interferometric synthesis, possessing high-precision azimuth positioning capability. However, due to the frequency difference between the target transmitter and the receiver platform’s local receiver, residual frequency offset is generated after de-slope processing of the local system. The longer the coherent integration time, the greater the phase error accumulated, which limits the coherent integration duration of the synthetic aperture, suppresses coherent integration gain, and severely degrades the system’s two-dimensional positioning accuracy [1–3]. Subsequently, multi-sub-aperture synthetic aperture passive localization technology was proposed, which assumes the residual frequency offset is approximately a fixed value, and reconstructs the residual frequency offset by solving equations constructed through multi-sub-apertures. However, since the time-varying nature of the residual frequency offset is ignored, the residual error cannot be fully compensated, and the system’s positioning accuracy remains limited [4–6].

Inspired by the very long baseline interferometry (VLBI) imaging technique in radio astronomy [7,8]. Traditional VLBI technology places two radio astronomical telescopes on the ground to receive the brightness temperature of celestial radiation. It exploits the Earth’s rotation to produce relative motion between the observation baseline and the target, and then performs interferometric processing between the observation stations to achieve high-resolution and high-sensitivity observational imaging. This method applies this technology to the field of passive localization. Two moving observation stations are used to form an interferometric baseline to passively receive electromagnetic signals radiated by the target, thereby achieving imaging of electromagnetic energy and thus detecting and locating the target.

This paper proposes a multi-station long baseline interferometry passive localization method. Multiple stations simultaneously receive target signals, and residual frequency offset is canceled through inter-station differential interferometry processing, effectively achieving signal de-modulation. Through the relative motion between observation stations and the mission area, two-dimensional sampling of the target radiation signal spatial spectrum is realized, enabling inversion of the target distribution and obtaining two-dimensional ultra-high-precision positioning results. This paper models the multi-station long baseline interferometry scenario, theoretically derives the correspondence between the target distribution map and the spatial spectrum sampling data, and achieves rapid inversion of the target distribution map through discrete Fourier transform.

2 Multi-Station Localization Model

Multi-station long baseline interferometry passive localization technology utilizes multiple stations moving in the observation scene to achieve two-dimensional sampling of the target radiation spatial spectrum. The greater the relative motion between stations, the longer the constructed virtual aperture, and the higher the positioning accuracy. This paper selects a simple bi-static scenario for principle demonstration. A coordinate system is established with the regional center of the imaging scene as the origin. Two observation platforms flight altitude is the same, H. The slant range from Platform 1 to the scene center is R_1\left({\theta }_1,{\phi }_1\right), with angle {\theta }_1 to the X-axis and angle {\phi }_1 to the Z-axis; the slant range from Platform 2 to the scene center is R_2\left({\theta }_2,{\phi }_2\right), with angle {\theta }_2 to the X-axis and angle {\phi }_2 to the Z-axis, as shown in Fig. 1.

Let the target scene region be X ≥ x ≥ –X, Y ≥ y ≥ –Y, with target point position P\left(x_0,y_0\right). The slant range from station to target can be expressed as \left|R\left(\theta ,\phi \right)-P\right|. The spectrum of the target radiation signal is S\left(f\right), and the carrier frequency of the signal is f_{\mathrm{c}}. The signal spectra received by Station 1 and Station 2 can be expressed as

(1)S\left(f,\theta ,\phi \right)=S\left(f\right)\exp \left\{-2\text{π}{f}_{\mathrm{c}}\frac{\left|R\left(\theta ,\phi \right)-P\right|}{c}\right\}

3 Multi-Station Long Baseline Interferometry Principle

The position of Observation Station 1 P_X\left(X,Y,Z\right) in the scene can be expressed as \left[R_1\sin {\phi }_1\cos {\theta }_1,R_1\sin {\phi }_1\sin {\theta }_1,R_1\cos {\phi }_1\right]. The slant range from Station 1 to the target can be expressed as

(2)\begin{array}{rl}|R_1\left(\theta ,\phi \right)-& P|=\sqrt{{\left(X-x_0\right)}^2+{\left(Y-y_0\right)}^2+{\left(Z\right)}^2}=\\ & \sqrt{R_1{}^2-2x_{0}X-2y_{0}Y+x_0{}^2+y_0{}^2}\approx \\ & R_1-\left(x_0\sin {\phi }_1\cos {\theta }_1+y_0\sin {\phi }_1\sin \theta \right)+\\ & \frac{x_0^2+y_0^2}{2R_1}\end{array}

When the target distribution region ensures \left(x_0{}^2+y_0{}^2\right)/2R_1⩽0, the slant range from Station 1 to the target can be approximated as

(3)\begin{array}{rl}& \left|R_1\left({\theta }_1,{\phi }_1\right)-P\right|\approx \\ & R_1\left({\theta }_1,{\phi }_1\right)-\left(x_0\sin {\phi }_1\cos {\theta }_1+y_0\sin {\phi }_1\sin {\theta }_1\right)\end{array}

Similarly, the slant range from Station 2 to the target point can be expressed as

(4)\begin{array}{rl}& \left|R_2\left({\theta }_2,{\phi }_2\right)-P\right|\approx \\ & R_2\left({\theta }_2,{\phi }_2\right)-\left(x_0\sin {\phi }_2\cos {\theta }_2+y_0\sin {\phi }_2\sin {\theta }_2\right)\end{array}

The bi-static passive reception of the same area signal, with frequency range f_1\sim f_2, the interferometry processing of the bi-static signals can obtain the interferometric signal S_{\mathrm{c}}(f,{\theta }_{1,2},{\phi }_{1,2}).

(5)\begin{array}{rl}& S_{\mathrm{c}}\left(f,{\theta }_{1,2},{\phi }_{1,2}\right)=\\ & \underset{f=f_1}{\overset{f_2}{\int }}{S}_1\left(f,{\theta }_1,{\phi }_1\right)S_2^{\ast }\left(f,{\theta }_2,{\phi }_2\right)\mathrm{d}f=\\ & \underset{f=f_1}{\overset{f_2}{\int }}{\left|S\left(f\right)\right|}^2\mathrm{d}f\exp \left\{-\mathrm{j}2\text{π}f\frac{\left|R_1-P\right|-\left|R_2-P\right|}{c}\right\}=\\ & \underset{f=f_1}{\overset{f_2}{\int }}{\left|S\left(f\right)\right|}^2\mathrm{d}f\exp \left\{-\mathrm{j}2\text{π}f\frac{R_1-R_2}{c}\right\}\times \\ & \exp \left\{-\mathrm{j}2\text{π}{f}_{\mathrm{c}}\frac{x_0\left(\sin {\phi }_2\cos {\theta }_2-\sin {\phi }_1\cos {\theta }_1\right)}{c}\right\}\times \\ & \exp \left\{-\mathrm{j}2\text{π}{f}_{\mathrm{c}}\frac{y_0\left(\sin {\phi }_2\sin {\theta }_2-\sin {\phi }_1\sin {\theta }_1\right)}{c}\right\}\end{array}

Interferometric processing performs signal demodulation, removing the effect of the signal waveform, and simultaneously achieves energy accumulation of the signal in the time-frequency domain, the resulting time-frequency integration gain is \sqrt{BT}, where B is the signal bandwidth, and T is the duration of each slice signal. The wider the signal bandwidth, the higher the SNR of the S_{\mathrm{c}}\left(f,{\theta }_{1,2},{\phi }_{1,2}\right) signal, making it more favorable for target detection.

Redefined variables

(6)\begin{array}{rl}\mathrm{\Delta }\varphi & \left({\theta }_{1,2},{\phi }_{1,2}\right)=\left(f_{\mathrm{c}}/c\right)\left(R_1\left({\theta }_1,{\phi }_1\right)-R_2\left({\theta }_2,{\phi }_2\right)\right)\\ \xi & =\left(f_{\mathrm{c}}/c\right)\left(\sin {\phi }_2\cos {\theta }_2-\sin {\phi }_1\cos {\theta }_1\right)\\ \eta & =\left(f_{\mathrm{c}}/c\right)\left(\sin {\phi }_2\sin {\theta }_2-\sin {\phi }_1\sin {\theta }_1\right)\end{array}

The target radiation intensity atP\left(x_0,y_0\right) is A\left(x_0,y_0\right).

(7)A\left(x_0,y_0\right)=\underset{f=f_1}{\overset{f_2}{\int }}{\left|S\left(f\right)\right|}^2\mathrm{d}f

Substituting into S_{\mathrm{c}}\left(\xi ,\eta \right), we obtain

(8)\begin{array}{rl}& S_{\mathrm{c}}\left(\xi ,\eta \right)=\\ & A\left(x_0,y_0\right)\exp \left\{-\mathrm{j}2\text{π}\mathrm{\Delta }\varphi \left({\theta }_{1,2},{\phi }_{1,2}\right)+x_0\xi +y_0\eta \right\}\end{array}

Since the positions of the observation stations are known, the slant ranges R_1\left({\theta }_1,{\phi }_1\right) and R_2\left({\theta }_2,{\phi }_2\right) from Observation Stations 1 and 2 to the region center can be precisely calculated, and the phase \mathrm{\Delta }\varphi \left({\theta }_{1,2},{\phi }_{1,2}\right) can be precisely compensated. The compensated signal is

(9)\begin{array}{rl}{S}_{cd}\left(\xi ,\eta \right)=& S_{\mathrm{c}}\left(\xi ,\eta \right)\ast \exp \left\{j2\text{π}\mathrm{\Delta }\varphi \left({\theta }_{1,2},{\phi }_{12}\right)\right\}=\\ & A\left(x_0,y_0\right)\exp \left\{-\mathrm{j}2\text{π}\left({\mathrm{x}}_0\xi +y_0\eta \right)\right\}\end{array}

Considering that the received signals at the observation stations consist of all electromagnetic emissions within the mission area, the interferometric compensated signal within the field of view is given by

(10)\begin{array}{rl}& S_{\mathrm{c}\mathrm{d}}\left(\xi ,\eta \right)=\\ & \underset{x=-X}{\overset{X}{\int }}\underset{y=-Y}{\overset{Y}{\int }}A\left(x_0,y_0\right)\exp \left\{-\mathrm{j}2\text{π}\left(x_0\xi +y_0\eta \right)\right\}\mathrm{d}x\mathrm{d}y\end{array}

From the above equation, it can be seen that S_{\mathrm{c}\mathrm{d}}\left(\xi ,\eta \right) is the two-dimensional Fourier transform of the target positions P\left(x_0,y_0\right) in the scene. Thus, the two-dimensional target distribution map A\left(x,y\right) can be obtained by performing an inverse two-dimensional Fourier transform, where the magnitude at each point in the map represents the radiation intensity of the target.

(11)\begin{array}{rl}& A\left(x,y\right)=\\ & \underset{\xi ={\xi }_1}{\overset{{\xi }_2}{\int }}\underset{\eta ={\eta }_1}{\overset{{\eta }_2}{\int }}{S}_{\mathrm{c}\mathrm{d}}\left(\xi ,\eta \right)\exp \left\{\mathrm{j}2\text{π}\left(x\xi +y\eta \right)\right\}\mathrm{d}\xi \mathrm{d}\eta \end{array}

This method achieves two-dimensional sampling of the observation region spatial spectrum through multi-station interferometry processing, and inverts the target distribution from the spatial spectrum. The positioning accuracy depends on the size of the spatial spectrum sampling window. The sampling ranges for \xi and \eta used in the computation are \mathrm{\Delta }\xi ={\xi }_{min}\sim {\xi }_{max}, \mathrm{\Delta }\eta ={\eta }_{min}{\eta }_{max}. The theoretical resolution is

(12)\begin{array}{rl}\delta x& =0.886\times 2/\mathrm{\Delta }\xi \\ \delta y& =0.886\times 2/\mathrm{\Delta }\eta \end{array}

Compared with the traditional method of directly detecting and localizing via single-channel signals, this method achieves energy accumulation of the signal in both the spatial and time-frequency domains, effectively improving the system gain. The gain N includes the virtual array synthesis gain and the coherent signal synthesis gain between the two stations \sqrt{BT}, where N is the number of slices in the synthesis processing. Since the system ultimately performs target detection in the image domain, the system sensitivity is the image noise variance. The noise variance N_0 at pixel \left(x,y\right) is calculated, i.e., the system sensitivity is as follows

(13)N_{\mathrm{s}}=\frac{N_0}{N\sqrt{BT}}

Considering that the image-domain SNR detection threshold is γ, the required SNR of the signal is as follows

(14){\mathrm{S}\mathrm{N}\mathrm{R}}_1=\frac{\gamma }{N\sqrt{BT}}

4 Proposed Method

Considering that the multi-station long baseline interferometry passive localization model must satisfy the following condition, It is usually required that d\varphi be less than \lambda /10.

(15)\mathrm{d}\varphi =\left[\left(x_0^2+y_0^2\right)/2R_1-\left(x_0^2+y_0^2\right)/2R_2\right]>>0

Considering that the observation field of view of the passive localization system is extremely large, if the entire field of view mission area uses one reference center, the above condition may not be satisfied in the far field of the observation area. In this case, the mission area needs to be divided into N sub-mission areas, with one reference center set for each sub-mission area. The size of each sub-mission area must satisfy the above condition, and multi-station long baseline interferometry passive localization processing is performed separately for each sub-mission area.

For positioning processing of each sub-region, with the observation region center as the origin, an XY orthogonal coordinate system is established on the observation region plane, and the observation region is gridded as (x, y), where x=x_{min}:\mathrm{\Delta }x:x_{max}, y=y_{min}:\mathrm{\Delta }y:y_{max}. The bi-static system must simultaneously receive electromagnetic signals within the observation region. The bi-static system moves relative to the observation region center O(0,0), maintains high-precision time synchronization between the two stations, and continuously outputs the position information P_X1 and P_X2 of the observation stations, as well as outputting sampled data s_1\left(t\right) and s_2\left(t\right). The processing flowchart is shown in Fig. 2, the basic steps are given below.

Step 1, data segmentation. Divide the data s_1\left(t\right) and s_2\left(t\right) into N × M segments respectively. Take N consecutive data slices s_{1,2}\left(t,n\right) as the data for one differential processing. t is fast time, n is slow time.

Step 2, Fourier transform. Perform Fourier transform on each of the N data slices s_{1,2}\left(t,n\right) to obtain the spectrum S_{1,2}\left(f,n\right). A total of M positioning processing operations are completed. To increase the number of positioning processing operations, a sliding segmentation method can also be adopted.

Step 3, variable transformation. Extract the sampling start time of each data slice, obtain the corresponding bi-static positions P_X1(n) and P_X2(n), calculate the slant ranges R_{1,2}\left(n\right) from the observation stations to the scene center, as well as the observation angles {\theta }_{1,2}\left(n\right) and {\phi }_{1,2}\left(n\right), perform variable transformation to obtain \xi \left(n\right) and \eta \left(n\right).

Step 4, bi-static interferometry processing. Perform interferometry processing on the bi-static spectrum data S_{1,2}\left(f,n\right) to obtain S_{\mathrm{c}}\left(n\right).

(16)S_{\mathrm{c}}\left(n\right)=\underset{f=f_1}{\overset{f_2}{\int }}{S}_1\left(f,n\right)S_2^{\ast }\left(f,n\right)\mathrm{d}f

Step 5, phase compensation. Calculate the phase compensation amount \mathrm{\Delta }\varphi \left(n\right), compensate the phase deviation, and obtain S_{\mathrm{c}\mathrm{d}}\left(n\right).

(17)\begin{array}{rl}{S}_{\mathrm{c}\mathrm{d}}\left(n\right)& =S_{\mathrm{c}}\left(n\right)\exp \left\{\mathrm{j}2\text{π}\mathrm{\Delta }\varphi \left(n\right)\right\}\\ \mathrm{\Delta }\varphi \left(n\right)& =\left(f_{\mathrm{c}}/c\right)\left(R_1\left(n\right)-R_2\left(n\right)\right)\end{array}

Step 6, complete the inverse Fourier transform. calculate the radiation intensity of each grid point in the scene region, and obtain the target distribution map A\left(x,y\right). Construct the inverse Fourier transform factor F\left(n,x,y\right), and calculate the inverse Fourier transform result for each grid point.

(18)\begin{array}{rl}A\left(x,y\right)& =\sum _n{S}_{\mathrm{c}\mathrm{d}}\left(n\right)F\left(n,x,y\right)\\ F\left(n,x,y\right)& =\exp \left\{\mathrm{j}2\text{π}\left[x\xi \left(n\right)+y\eta \left(n\right)\right]\right\}\end{array}

Step 7, target detection. Perform image detection, set detection threshold, and extract target position information through threshold detection.

(19)P\left(x,y\right)=\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{d}\left|\left|A\left(x,y\right)>\gamma \right|\right|

In actual processing, since the task area spans hundreds of kilometers while the system resolution can reach the meter level, the computational load required for inversion is enormous. In this case, a hierarchical resolution processing method can be adopted to reduce the computational load. In the first step, the grid resolution is set to 1–10 km, and the first localization processing is performed to complete coarse localization. In the second step, the target imaging area is narrowed down to the grid points where targets are detected, the grid resolution is set to 1–100 m, and inverse Fourier transform processing and target detection are carried out to complete the second-stage coarse localization, further narrowing the target distribution area. In the third step, the grid resolution is set to 0.1 m, and inverse Fourier transform processing and target detection are performed to complete the final localization.

5 Simulation Verification

Simulation verification of the multi-station long baseline interferometry passive localization method is performed. The observation scene is shown in Fig. 1. The mission area size is 20 km × 20 km, gridded at 0.1 m resolution. Target points are set at (0,0) and (50, –50) respectively. The carrier frequency of the radiation signal is 3 GHz, and the signal type is a single-tone continuous wave s\left(t\right)=\exp \left(2\text{π}{f}_{\mathrm{c}}t\right). The flight altitudes of Station 1 and Station 2 are fixed at 80 km. The initial positions of Station 1 and Station 2 are (0, 0, 80000). The flight speed of Station 1 is 128 m/s, with the flight direction at an angle of 51.3° to the X-axis. The flight speed of Station 2 is 112 m/s, with the flight direction at an angle of 63.4° to the X-axis. Station 1/2 receive signals respectively, continuously collecting 100 s of data. The data is segmented into 1000 data segments at 0.1 s intervals. The observation station positions corresponding to each sampling time are shown in Tab. 1.

The 1000 data segments are processed for differential passive localization. A hierarchical localization processing method is adopted to reduce the computational load. Coarse localization is completed using a grid resolution of 1 km within a 20 km task area, the obtained target distribution map is shown in Fig. 3(a). Coarse localization is performed again using a grid resolution of 1 m within a 1 km suspected area, the obtained target distribution map is shown in Fig. 3(b). Finally, final localization is completed using a grid resolution of 0.1 m within a 20 m suspected area, the obtained target distribution map is shown in Fig. 3(c). Through image detection, the positions of 2 targets can be extracted, located at (0,0) and (50,-50) respectively.

Shown in Fig. 4, the equivalent observation variable ranges in the X and Y directions are calculated as \xi =-3.3223\sim 0, \eta =0\sim 1.302. The obtained observation window lengths are \mathrm{\Delta }\xi =3.3223 and \mathrm{\Delta }\eta =1.302 respectively. The theoretical resolution is calculated as

(20)\begin{array}{rlr}\delta x& =0.886\times 2/\mathrm{\Delta }\xi & =0.53\mathrm{m}\\ \delta y& =0.886\times 2/\mathrm{\Delta }\eta & =1.36\mathrm{m}\end{array}

The X-direction and Y-direction slices of the target (0,0) image are obtained as shown in Fig. 5 and Fig. 6 respectively. The extracted 3 dB resolutions are X = 0.5 m and Y = 1.3 m, which agree with the theoretical values, proving that this method possesses two-dimensional high-precision positioning capability.

Through simulation, the multi-station ultra-long baseline interferometry passive localization method is verified. Two-dimensional high-precision positioning of electromagnetic targets can be achieved, and the positioning accuracy agrees with the theoretical analysis.

6 Conclusion

This paper proposes a multi-station ultra-long baseline interferometry passive localization method. The observation stations only need to passively receive target signals through a single antenna. Through interferometry processing between the bi-static stations, two-dimensional sampling of the spatial spectrum region is achieved, the target distribution map is inverted, and two-dimensional high-precision positioning of the target is obtained. This paper establishes a digital model of the interferometry scene, theoretically derives the correspondence between the target distribution map and the spatial spectrum sampling data, analyzes the positioning accuracy of this method, presents the processing flow of this method, and verifies the effectiveness of the method through simulation.

References

[1] S. Li, H. Huan, and J. Yang, “Passive synthetic aperture equivalent model and localization for spaceborne squint scenario, ” in Proceedings of the IEEE Wireless Communications and Networking Conference (WCNC), Kuala Lumpur, Malaysia, pp. 1-6, 2026.

[2] J. Sun, H. Huan, R. Tao, and Y. Wang, “A passive synthetic aperture localization method based on sparse sampling reconstruction, ” IEEE Geoscience and Remote Sensing Letters, vol. 22, pp. 1-5, 2025.

[3] L. Zhang, H. Huan, R. Tao, and X. Tang, “A residual frequency offset estimation method based on range migration fitting, ” Journal of Beijing Institute of Technology, vol. 33, no. 2, pp. 103-110, 2025.

[4] J. Yang, H. Huan, and R. Tao, “Multisubaperture backward projection positioning algorithm for radiation sources by single satellite, ” IEEE Transactions on Aerospace and Electronic Systems, vol. 61, no. 2, pp. 1267-1282, 2025.

[5] X. Liu, H. Huan, and R. Tao, “Synthetic aperture positioning using subaperture interference, ” IEEE Sensors Journal, vol. 25, no. 7, pp. 11378-11391, 2025.

[6] J. Yang, H. Huan, Z. Xu, and R. Tao, “Research on distributed synthetic aperture passive positioning and optimal geometric configuration, ” IEEE Signal Processing Letters, vol. 32, pp. 3685-3689, 2025.

[7] J. Y. Kim , T. P. Krichbaum , A. E. Broderick . et al. Event horizon telescope imaging of the archetypal blazar 3C 279 at an extreme 20 microarcsecond resolution. Astron. Astrophys., 2020, 640: A69

[8] Z. X. Huang , O. Titov , M. K. Schmidt , B. Pope , G. K. Brennen , D. Oi , P. Kok . Quantum-enabled optical large-baseline interferometry: applications, protocols and feasibility. Adv. Phys. X, 2026, 11(1): 2597311

Rights & permissions

The associate editor coordinating the review of this manuscript was Dr. Hao Huan