Joint Doppler Center Estimation for Synthetic Aperture Positioning Using Starlink Signals of Opportunity

Abstract

To address the degradation of global navigation satellite system (GNSS) positioning performance in challenging environments such as signal blockage and interference, this paper investigates a high-precision positioning method based on Starlink satellite signals of opportunity. A self-developed receiving platform is used to receive and analyze a Starlink downlink narrowband beacon signal at 11.45 GHz. To overcome the challenges of unknown prior signal parameters, unknown prior frequency offset under low-SNR (signal-to-noise ratio) conditions, and the difficulty of robustly estimating the reference time, a boundary-constrained integrated stationary phase method (BISPM) is proposed to achieve joint and accurate estimation of the reception time and carrier frequency. On this basis, a tomographic synthetic aperture positioning model based on virtual steering vectors is constructed. By coherently integrating short-time observations and performing back-projection imaging, the proposed method improves positioning resolution under single-satellite conditions. Experimental results using the actually received signal from STARLINK-32494 show that, with a synthetic aperture time of 3 s and only a single satellite, positioning accuracies of approximately 0.4 m in azimuth and 8 m in range can be achieved. These results verify the feasibility of using Starlink narrowband beacon signals for opportunistic positioning and provide a useful reference for rapid and high-precision positioning with low-Earth-orbit satellite signals of opportunity.

Graphical abstract

Keywords

Starlink signals of opportunity / Doppler center estimation / boundary-constrained integrated stationary phase method / synthetic aperture positioning / tomographic imaging / virtual steering vector

1 Introduction

Global navigation satellite system (GNSS) is the main enabler of modern positioning, navigation, and timing (PNT) services, but its performance can degrade significantly in challenged environments because of signal blockage, weak received power, and vulnerability to interference and jamming [1]. To overcome these limitations, signals of opportunity (SOP) have been widely studied as an alternative navigation source. Among them, low-Earth-orbit (LEO) satellite signals are especially attractive due to their wide coverage, strong received power, and large constellation size [2]. In particular, compared with Iridium, Starlink offers broader coverage, a larger constellation, and continuous beacon signals, making Starlink-based opportunistic positioning an important research direction [3].

In recent years, many studies have investigated parameter estimation and positioning using LEO satellite SOP signals. Existing approaches mainly include Doppler- or Doppler-rate-based positioning, blind estimation for unknown signal structures, and methods that exploit phase information or coherent integration to improve performance. For Iridium and other LEO systems, previous studies have shown that weighted least-squares positioning based on instantaneous Doppler can achieve a three-dimensional accuracy better than 70 m [4], while instantaneous Doppler positioning with height aiding can provide horizontal accuracy better than 200 m [5]. Joint pseudorange-Doppler processing using iridium satellite time and location (STL) burst signals has achieved a positioning accuracy better than 100 m [6], and a single-search Doppler matching method for a single-satellite scenario has reported an accuracy of about 193.3 m [7]. For Starlink signals, prior work has reported a two-dimensional positioning error of 25.9 m using six satellites over about 800 s [8]. By treating Starlink downlink signals as SOP signals with unknown structures, matched subspace detection and multi-satellite Doppler/range-rate estimation have achieved about 10 m horizontal accuracy [9]. A cognitive beacon detection and estimation framework has further demonstrated 4.15 m navigation accuracy over a 1 km ground-vehicle trajectory [10]. In addition, short-time Fourier transform (STFT) and maximum-likelihood estimation have achieved horizontal accuracy of about 15 m using four satellites over 260 s with height aiding [11], while differential Doppler methods have reported two-dimensional errors of 5.6 m and 2.6 m under different baseline conditions [12]. Frequency-domain sliding-window estimation combined with Doppler-based least-squares solving has achieved 11.8 m horizontal error using six satellites and 100 s of observations [13]. Joint pseudorange-Doppler processing has achieved 2.7665 m two-dimensional accuracy and 6.5015 m three-dimensional accuracy under short-baseline conditions [14]. Iterative blind estimation of full Starlink beacon signals has yielded positioning errors of about 10 m and 2 m within 8 s and 20 s, respectively, using only three satellites [3]. Extended Kalman filter (EKF)-based fusion of synchronization symbols and pilot information has also achieved a two-dimensional error of 6.1 m and a three-dimensional error of 13.5 m under simulation conditions [15]. These studies verify the feasibility of Starlink-based positioning, but most still depend on multiple satellites, long observation times, or external aiding information.

To address these limitations, this paper focuses on Starlink downlink narrowband beacon signals and proposes a boundary-constrained integrated stationary phase method combined with a synthetic-aperture positioning approach based on tomographic imaging and virtual steering vectors. By coherently integrating the received signal, the proposed method improves the signal-to-noise ratio and enhances positioning performance under low input signal-to-noise ratio (SNR) [16]. Experimental results using actual signals from STARLINK-32494 show that, with only 3 s of synthetic aperture time, the proposed method can achieve positioning accuracies of about 0.4 m in azimuth and 8 m in range, demonstrating its potential for single-satellite, short-observation positioning applications.

2 Characteristics of Starlink Downlink Signals

The Starlink system provides broadband satellite Internet services over multiple frequency bands, among which the Ku-band downlink is one of the most accessible signal sources for ground-based reception and signal-of-opportunity processing. This paper focuses on the downlink signals transmitted from Starlink satellites to ground user terminals. In the Ku band, Starlink adopts a multi-channel frequency allocation scheme, where each channel occupies 250 MHz, including 240 MHz of effective bandwidth and two 5 MHz guard bands, as illustrated in Fig. 1.

Previous studies have shown that Starlink Ku-band downlink signals employ an orthogonal frequency-division multiplexing (OFDM) waveform and contain characteristic narrowband tones near the channel center. Compared with the broadband data component, these narrowband features are easier to identify in the frequency domain and exhibit better continuity and stability in the time-frequency domain, making them more suitable for parameter estimation and opportunistic positioning under limited prior signal knowledge. Based on the received data, a stable narrowband characteristic signal was identified near 11.45 GHz, with an adjacent subcarrier spacing of approximately 44 kHz [13, 17].

The use of this narrowband characteristic signal is motivated by two main considerations. 1)Its phase and frequency evolution are relatively clear, facilitating the extraction of Doppler center time and carrier frequency using the proposed boundary-constrained integrated stationary phase method. 2)The continuous Doppler variation caused by satellite motion supports time-domain synthetic aperture construction and coherent integration, thereby providing a basis for subsequent high-resolution positioning[16, 18].

3 System Platform

To enable stable reception and analysis of Starlink downlink narrowband characteristic signals, an observation platform consisting of a horn antenna, a low-noise amplifier (LNA), a down-conversion module, an X310 data acquisition device, and a GPS-disciplined oscillator (GPSDO) clock module was developed. A block diagram of the system platform is shown in Fig. 2. The platform is capable of Ku-band signal reception, frequency down-conversion, and highly stable signal sampling. In particular, the GPSDO provides a unified time-frequency reference for the receiving system, which helps suppress frequency drift and improve timing accuracy. Based on this platform, the Starlink downlink narrowband characteristic signal near 11.45 GHz was successfully received. After signal reception, the source satellite was identified by combining publicly available orbital ephemeris information with satellite matching. Specifically, two-line element (TLE) data of candidate satellites were obtained from CelesTrak, and candidate screening and matching were performed in systems tool kit (STK) based on reception time, antenna pointing direction, and Doppler-curve search results (Fig. 3). This procedure provides the position and velocity of the target satellite during the observation interval, which serve as prior orbital information for subsequent Doppler modeling, the proposed frequency-estimation method, and tomographic synthetic aperture positioning using virtual steering vectors.

4 Boundary-Constrained Integrated Stationary Phase Method (BISPM)

4.1 Signal Model and Problem Formulation

Consider a received downlink narrowband satellite signal. After down-conversion and bandpass filtering, its complex baseband representation can be expressed as

(1)s(t)=A(t)\exp [\mathrm{j}\varphi (t)]+w(t)

where A(t) denotes the slowly varying amplitude envelope, w(t) is additive white Gaussian noise, and the phase term \varphi (t) can be further written as

(2)\varphi (t)=2\text{π}{f}_{\mathrm{p}}t+2\text{π}{\int }_0^t{f}_{\mathrm{D}}(\tau )\mathrm{d}\tau +{\varphi }_0

where f_{\mathrm{p}} denotes the unknown residual carrier-frequency offset, f_{\mathrm{D}}(t) denotes the geometric Doppler frequency induced by the relative motion between the satellite and the receiver, and {\varphi }_0 is the initial phase. Under relative satellite motion, the geometric Doppler frequency can be approximately modeled as

(3)f_{\mathrm{D}}(t)=f_{\mathrm{D}0}+{\dot{f}}_{\mathrm{D}}t+O(t^2)

where f_{\mathrm{D}0} denotes the Doppler centroid frequency, and \dot{f}{}_{\mathrm{D}} denotes the Doppler rate. The time instant satisfying f_{\mathrm{D}}(t)=0 corresponds to the geometrically symmetric moment at which the radial velocity between the satellite and the receiver becomes zero. This instant is of key physical significance for synthetic aperture positioning and orbital parameter estimation. The objective of this study is to accurately estimate the Doppler zero-crossing time and the corresponding carrier frequency offset under unknown initial frequency offset, time mismatch between the predicted satellite-orbit time and the actual signal reception time, finite observation duration, and noisy conditions.

4.2 Principle of the Integrated Stationary Phase Method

According to the stationary phase principle (SPM) [19], when the derivative of the phase function of the integral kernel becomes zero at a certain time instant, that instant makes the dominant contribution to the integral result. In the presence of the residual carrier-frequency offset f_{\mathrm{p}}, the phase derivative contains both f_{\mathrm{p}} and the geometric Doppler component. By differentiating the signal phase \varphi (t), the instantaneous frequency can be obtained as

(4)\hat{f}{}_{\mathrm{D}}(t)=\frac{1}{2\pi }\frac{\mathrm{d}\varphi (t)}{\mathrm{d}t}=f{}_{\mathrm{p}}+f{}_{\mathrm{D}}(t)

Therefore, the Doppler center time t₀, at which the Doppler frequency becomes zero, satisfies

(5)f_{\mathrm{D}}(t_0)=0,\hat{f}{}_{\mathrm{D}}(t_0)=f_{\mathrm{p}}

However, under practical low-SNR and finite-aperture conditions, directly detecting the zero-crossing point of the instantaneous frequency is highly susceptible to noise, endpoint effects, and phase-unwrapping errors. To improve robustness, the proposed method adopts the idea of the integrated stationary phase method. By integrating the phase derivative and analyzing its stationarity, the Doppler zero-crossing point can be identified indirectly. Define the first-order derivative of the phase, i.e., the instantaneous frequency estimate, as

(6)\hat{f}{}_{\mathrm{D}}(t)=\frac{1}{2\pi }\frac{\mathrm{d}\hat{\varphi }(t)}{\mathrm{d}t}

(7)\hat{\varphi }(t)=\text{unwrap}\left(\arg \{s(t)\}\right)

An integral-type stationarity function is then constructed as

(8)I(t)=\frac{1}{t}{\int }_0^t{f}_D(\tau )\mathrm{d}\tau

Here, f{}_{\mathrm{D}}(t) in the integral stationarity function denotes the geometric Doppler component. This function will exhibit stationary behavior provided that the geometric Doppler component stays locally small near the Doppler center time.

It should be noted that the differentiation–integration operation is not used to recover the original phase. Instead, it constructs a robust stationarity measure from the instantaneous frequency variation. The differentiation step removes the unknown initial phase and extracts the Doppler-related frequency variation. The subsequent integration acts as a cumulative averaging operation, which suppresses local noise fluctuations and reduces the sensitivity to phase-unwrapping errors. In this way, the direct pointwise zero-crossing detection of the instantaneous frequency is converted into a more robust local stationarity detection problem.

4.3 Boundary Constraints and Candidate-Point Scoring Mechanism

Within a finite observation interval t\in [0,T], the endpoint regions typically suffer from the following issues: concentrated phase-unwrapping errors, pronounced filter transients and edge effects, and the failure of the approximate Doppler variation model near the boundaries. To avoid falsely identifying nonphysical stationary points as Doppler zero-crossing points, a boundary-constrained mechanism is introduced into the integrated stationary-phase search process, such that the search is performed only within the following valid time interval:

(9)t\in [\alpha T,(1-\alpha )T],\alpha \in (0,0.5)

where \alpha is the boundary constraint coefficient. It determines the proportion of the observation interval excluded from both ends. A larger \alpha can more effectively suppress endpoint effects caused by filtering transients and phase-unwrapping errors, but it also reduces the effective search interval. In this work, \alpha =0.1 is used, which means that the first and last 10% of the observation interval are excluded from the candidate search. Within this interval, zero-crossing detection is performed on the derivative of the integral function,{\displaystyle \frac{\mathrm{d}I(t)}{\mathrm{d}t}}, and candidate points are further screened in combination with a local stationarity criterion. Based on the above analysis, a candidate Doppler zero-crossing time t_k is defined to satisfy the following conditions. Integral stationarity condition:

(10){\frac{\mathrm{d}I(t)}{\mathrm{d}t}|}_{t=t_k}\approx 0

Local stationarity condition:

(11)

Boundary constraint condition:

(12)t_k\in [\alpha T,(1-\alpha )T]

where \mathcal{N}(t_k) denotes a local time window in the vicinity of t_k, which is used to evaluate the local stationarity around each candidate point. A larger window provides stronger suppression of local noise fluctuations, but may smooth out the Doppler variation near the center time. A smaller window preserves local variations better, but is more sensitive to noise. ϵ is the stability threshold, which controls the strictness of the local stationarity test. A smaller ϵ rejects unstable candidates more strictly, whereas a larger ϵ may retain noise-induced candidates. To further improve reliability, a center-weighted candidate scoring function is introduced as follows

(13)\begin{array}{rl}\Gamma (t_k)=& \frac{1}{|{\hat{f}}_{\mathrm{D}}(t_k)|+\delta }\cdot \frac{1}{\text{std}\left\{{{\displaystyle \frac{\mathrm{d}I(t)}{\mathrm{d}t}}|}_{t\in \mathcal{N}(t_k)}\right\}+\eta }\cdot \\ & \left(1-\frac{|t_k-T/2|}{T/2}\right)\end{array}

where \delta and \eta are regularization constants introduced to avoid numerical singularities when the denominator terms approach zero. When they are chosen to be much smaller than the corresponding denominator terms, they have little influence on the ranking of valid candidates. The last term in \Gamma (t_k) serves as a center-weighting factor, assigning higher scores to candidates closer to the center of the observation interval. This weighting reduces the probability of selecting edge-induced false stationary points. Finally, the candidate point with the highest score is selected as the Doppler zero-crossing time, i.e., the Doppler center time,

(14){\hat{t}}_0=\mathrm{a}\mathrm{r}\mathrm{g}\underset{t_k}{max}\Gamma (t_k)

4.4 Doppler Center Time and Carrier Frequency Estimation

After obtaining the estimated Doppler zero-crossing time {\hat{t}}_0, the compensated geometric Doppler component satisfies f_{\mathrm{D}}({\hat{t}}_0)\approx 0, and the corresponding raw instantaneous frequency satisfies

(15)\begin{array}{rl}& {\hat{f}}_{\mathrm{D}}({\hat{t}}_0)\approx f_{\mathrm{p}},\hat{f}{}_{\mathrm{p}}={\hat{f}}_{\mathrm{D}}({\hat{t}}_0)\\ & {\hat{f}}_{\mathrm{c}\mathrm{a}\mathrm{r}\mathrm{r}\mathrm{i}\mathrm{e}\mathrm{r}}=f{}_{\mathit{c}}+{\hat{f}}_{\mathrm{p}}=f{}_{\mathit{c}}+\hat{f}{}_{\mathrm{D}}({\hat{t}}_0)\end{array}

where \hat{f}{}_{\mathrm{D}}({\hat{t}}_0) is the instantaneous frequency estimate obtained by BISPM at the Doppler center time. Since this time instant is constrained within a locally stationary and boundary-reliable region, the final carrier-frequency estimation error is mainly affected by phase noise and higher-order Doppler residuals, and thus exhibits good statistical stability. Through the above procedure, the proposed method achieves joint robust estimation of the Doppler center time and carrier frequency for narrowband signals with unknown prior information. These estimates can be used not only for subsequent time-frequency compensation of the received signal, but also to provide a unified phase reference and trajectory modeling basis for tomographic synthetic aperture positioning based on virtual steering vectors.

5 Tomographic Synthetic Aperture Positioning Method Based on Virtual Steering Vectors

This paper adopts a back-projection positioning method based on tomographic synthetic aperture imaging with virtual steering vectors. In a single-station satellite signal-of-opportunity positioning scenario, the satellite motion along its trajectory forms an equivalent synthetic aperture during the observation interval, while the received signal exhibits a continuously varying Doppler history in the slow-time dimension. By exploiting the continuous satellite trajectory information, the theoretical Doppler history corresponding to each candidate search grid point can be reconstructed and then correlated with the actual received signal to obtain the response output of each grid point. The position with the maximum correlation response is taken as the estimated location of the radiation source. Furthermore, by equivalently modeling the slow-time observations as a one-dimensional virtual array distributed along the satellite trajectory, the positioning process can be formulated as a near-field beamforming and spatial tomographic focusing process based on virtual steering vectors.

5.1 Geometry and Signal Model

Consider a ground receiver positioning scenario based on satellite signals of opportunity, as shown in Fig. 4. The satellite acts as a non-cooperative signal source and moves along a known orbit, with its instantaneous spatial position denoted by {\mathit{r}}_s(t). The ground receiver is located at an unknown position {\mathit{r}}_r(t). Over the slow-time interval t\in \left[-{\displaystyle \frac{T}{2}},{\displaystyle \frac{T}{2}}\right],the received baseband signal can be expressed as

(16)s(t)=A(t)\exp (-\mathrm{j}\frac{2\pi }{\lambda }R(t))+n(t)

where R(t)=\Vert {\mathit{r}}_s(t)-{\mathit{r}}_r\Vert is the instantaneous propagation distance from the satellite to the receiver, \lambda is the carrier wavelength, and n(t) denotes noise. Due to the high-speed motion of the satellite, R(t) varies continuously with time, causing the received signal to exhibit a pronounced Doppler history along the slow-time dimension, which provides the physical basis for synthetic aperture positioning.

5.2 Principle of Near-Field Beamforming

Conventional far-field beamforming is based on the plane-wave assumption, which requires the signal source distance to satisfy the far-field condition R⩾{\displaystyle \frac{2D^2}{\lambda }}. In the single-satellite synthetic aperture positioning scenario, the one-dimensional virtual array formed by the satellite motion along its trajectory exhibits evident range-dependent characteristics. The propagation distance differences from the radiation source to different virtual array elements cannot be neglected, making the conventional far-field plane-wave approximation no longer applicable. Therefore, this paper adopts a near-field spherical-wave model to describe the phase evolution of slow-time observations, and formulates the positioning process as a near-field beamforming and tomographic focusing procedure based on virtual steering vectors. Specifically, near-field beamforming is used to obtain the matched response for each candidate position, while tomographic focusing is employed to form the spatial response distribution over the entire search area. Based on this model, the virtual steering vector of each candidate grid point is constructed and correlated with the actual received signal. The position corresponding to the peak response is then taken as the estimated location of the radiation source.

Consider the virtual array model, as shown in Fig. 5. The satellite is continuously sampled during the synthetic aperture interval, forming a set of virtual array elements. For an arbitrary ground point P, the distances from P to different virtual array elements are not identical, resulting in different received signal phases. Under far-field conditions, the phase difference can be approximated as

(17){\phi }_n=-\frac{2\pi }{\lambda }nd\sin \theta

Under near-field conditions, however, the curvature of the spherical wavefront must be taken into account, yielding

(18){\phi }_n=-\frac{2\pi }{\lambda }\sqrt{n^2{d}^2+R^2+2ndR\sin \theta }

where d is the array-element spacing, R is the distance from the reference point to the target, and \theta is the azimuth angle. Compared with the linear phase term nd\sin \theta in the far-field model, the near-field model introduces an additional quadratic phase term related to d^2/R, enabling more accurate characterization of the nonlinear Doppler history caused by trajectory curvature.

5.3 Construction of Virtual Steering Vectors

The area of interest is discretized into an M\times N grid, as shown in Fig. 6, and the coordinates of each grid point are given by P_{m,n}=(x_m,y_n). For the k-th virtual array element, corresponding to the satellite position at time t_k, the steering-vector component with respect to the grid point P_{m,n} is given by

(19)s_d(k;P_{m,n})=\exp \left[-\mathrm{j}\frac{2\pi }{\lambda }\sqrt{x_m^2+(kd-y_n{)}^2}\right]

For a spatial point P_{m,n}, the complete steering vector formed by all K virtual array elements is

(20){\mathit{a}}_{m,n}=\left[\begin{array}{c}{s}_d(1;P_{m,n})\\ s_d(2;P_{m,n})\\ ⋮\\ s_d(K;P_{m,n})\end{array}\right]

This process essentially transforms the Doppler history in the time domain into phase encoding in the spatial domain, thereby providing prior information for subsequent tomographic imaging.

5.4 Beamforming and Tomographic Imaging-Based Positioning

To improve the computational efficiency of the near-field beamforming and tomographic imaging-based positioning process, the virtual steering vectors corresponding to all candidate grid points in the entire search area are arranged column by column to form a steering matrix. In this way, a matrix-form positioning model can be established as follows

(21)\mathit{r}=\mathit{A}\mathit{\sigma }+\mathit{n}

where \mathit{r} denotes the received signal vector, \mathit{\sigma } denotes the target spatial response vector arranged in column-stacked form, \mathit{A} denotes the steering matrix, and \mathit{n} denotes the additive noise vector.

The steering matrix \mathit{A} is composed of the virtual steering vectors corresponding to all candidate grid points, and its structure can be written as

(22)\mathit{A}=\left[\begin{array}{cccc}{a}_{1,11}& a_{1,12}& \cdots & a_{1,MN}\\ a_{2,11}& a_{2,12}& \cdots & a_{2,MN}\\ ⋮& ⋮& \ddots & ⋮\\ a_{K,11}& a_{K,12}& \cdots & a_{K,MN}\end{array}\right]

On this basis, the grid-by-grid matched-correlation process can be uniformly expressed in matrix form, yielding the beamforming output, i.e., the matched-filter output, for all candidate grid points as

(23)\hat{\mathit{\sigma }}={\mathit{A}}^{\mathrm{H}}\mathit{r}

where each element of \hat{\mathit{\sigma }} corresponds to the response output of a candidate grid point. Furthermore, by reshaping the squared magnitude of \hat{\mathit{\sigma }}, the tomographic imaging response map can be constructed as

(24)\mathit{F}=\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{h}\mathrm{a}\mathrm{p}\mathrm{e}(|\hat{\mathit{\sigma }}{|}^2,[M,N])

Through peak detection on the tomographic response map, the radiation source location estimation can be obtained:

(25){\hat{P}}_0=\arg \underset{m,n}{max}\mathit{F}(m,n)

Unlike the traditional point-by-point matched filter search, this paper constructs a steering matrix that includes the virtual steering vectors of all candidate grid points, thereby unifying the pointwise inner product operations into matrix–vector multiplication. This matrix-based processing does not change the algorithm itself, but merely rewrites the original pointwise computation as standard linear algebraic operations. Consequently, it enables parallel acquisition of the response outputs for all grid points, reduces computational overhead from loops and repeated memory access, and improves the computational efficiency of tomographic localization.

6 Experimental Results

Experiments on tomographic synthetic aperture positioning based on virtual steering vectors were conducted using the receiving platform and frequency estimation method described above. The experimental scenario is shown in Fig. 7. The signal received by the X310 is shown in Fig. 8.

6.1 Satellite Identification

The signal segment with the strongest energy was first selected for satellite identification.

According to the reception time and antenna beam pointing, five candidate satellites were screened, as shown in Fig. 3. A Doppler matching method was then employed to determine the actual source satellite of the received signal, and the matching results are shown in Fig. 9. From the results, it can be concluded that the received signal originated from satellite STARLINK-32494. The corresponding time-frequency representation of the received signal is shown in Fig. 10. As can be seen from the time-frequency plot, the two subcarriers of the Starlink beacon signal are separated by approximately 44 kHz, which is consistent with previously reported results.

6.2 Frequency Estimation

The carrier frequency of the received signal was first coarsely estimated from the spectrum acquired by the X310. The down-converter translated the signal centered at 11.45 GHz to an intermediate frequency centered at 1.8 GHz, and the X310 performed zero-intermediate frequency (IF) sampling. Based on this configuration, the carrier frequency of the received signal was coarsely estimated to be approximately 11.449749340\mathrm{G}\mathrm{H}\mathrm{z}, corresponding to an offset of about -250.660\mathrm{k}\mathrm{H}\mathrm{z} from the nominal 11.45\mathrm{G}\mathrm{H}\mathrm{z} center frequency. The time interval from the start of signal acquisition to the Doppler center time was approximately 537\mathrm{s}. The GPSDO module of the X310 provided a stable and accurate timing reference. A 3 s data segment centered at the Doppler center time was then selected for subsequent processing. For reproducibility, the main implementation parameters of the BISPM are listed in Tab. 1. These parameters include the boundary constraint coefficient, local stationarity window, stability threshold, and regularization constants used in the candidate-point scoring process.

After coarse time-frequency compensation, the received signal was processed using the proposed BISPM to refine the Doppler center time and carrier frequency, as shown in Fig. 11.

Fig. 11(a) shows the signal amplitude and unwrapped phase within the selected 3 s observation interval. The phase exhibits a continuous nonlinear variation, which reflects the Doppler variation induced by the satellite motion. Fig. 11(b) gives the instantaneous frequency estimated from the phase derivative, where the Doppler center is located around the locally stationary region after removing the residual frequency offset. Fig. 11(c) shows the cumulative phase variation rate used in the integrated stationarity analysis. Compared with directly detecting the zero-crossing point of the noisy instantaneous frequency, this integrated quantity provides a smoother criterion for identifying the Doppler center. Fig. 11(d) presents the candidate-point scoring result. A dominant peak appears within the boundary-constrained search interval, indicating that the selected candidate satisfies both the integral stationarity condition and the local stationarity criterion.

The experimental results indicate that the Doppler center time is estimated as 537.395\mathrm{s}, where the time origin is defined as the beginning of the received signal, and the corresponding carrier frequency is estimated as 11449749331.675\mathrm{H}\mathrm{z}. The estimated Doppler center time provides the time reference for synthetic aperture construction, while the corresponding frequency estimate provides the carrier-frequency reference for subsequent coherent integration and positioning.

6.3 Positioning Results

After compensation, the received signal was processed using the tomographic synthetic aperture positioning method based on virtual steering vectors. The positioning results are shown in Fig. 12.

Fig. 12(a) shows the three-dimensional tomographic imaging response over the search area. A distinct focused peak can be observed, indicating that the received signal is coherently accumulated at the correct candidate position after time-frequency compensation. Fig. 12(b) gives the corresponding two-dimensional search result, where the estimated position is close to the reference position. This confirms that the virtual-steering-vector-based synthetic aperture model can effectively convert the slow-time Doppler history into spatial focusing information. The estimated receiver coordinates are X=-2170625.099\mathrm{m}, Y=4388199.865\mathrm{m}, and Z=4074412.870\mathrm{m}. The azimuth-angle error is 4.8648985\times {10}^{-5}, corresponding to an azimuth-direction error of 0.3917\mathrm{m}, while the signed range-direction error is -8.106\mathrm{m} (an absolute error of 8.106\mathrm{m}); the reported root-mean-square error (RMSE) is 32.263\mathrm{m}.

The result shows that the response is more concentrated in the azimuth direction than in the range direction. This is consistent with the single-satellite synthetic aperture geometry, because the satellite motion provides stronger Doppler diversity and coherent accumulation gain along the azimuth direction, whereas the range direction is less constrained under the short-time single-satellite observation condition. Using only 3 s of received signal from a single satellite, the proposed method achieves a positioning accuracy of approximately 0.4\mathrm{m} in the azimuth direction and 8\mathrm{m} in the range direction. These results demonstrate the effectiveness of the proposed method for improving positioning resolution under a short-observation, single-satellite scenario.

7 Conclusion

This paper investigated opportunistic positioning using Starlink downlink narrowband beacon signals under a single-satellite and short-observation scenario. Based on a self-developed receiving platform, a Starlink downlink narrowband characteristic signal near 11.45\mathrm{G}\mathrm{H}\mathrm{z} was successfully received and analyzed, verifying its observability and potential for opportunistic positioning. To estimate the Doppler center time and carrier frequency under unknown prior signal parameters, a BISPM was proposed. Based on the estimated time-frequency reference, a tomographic synthetic aperture positioning method using virtual steering vectors was further developed to achieve spatial focusing from short-time satellite observations.

Experiments using the actually received signal from STARLINK-32494 show that, with a synthetic aperture time of only 3\mathrm{s} and a single satellite, positioning accuracies of approximately 0.4\mathrm{m} in azimuth and 8\mathrm{m} in range can be achieved. These results verify the feasibility of using Starlink narrowband characteristic signals for rapid opportunistic positioning. Future work will focus on multi-satellite cooperative positioning, robust parameter estimation in complex noise environments, and three-dimensional high-precision positioning.

References

[1] X. Zhang , F. Ma . Review of leo navigation augmentation for GNSS. Acta Geodaetica et Cartographica Sinica, 2019, 48(9): 1073–1087

[2] C. Yang and A. Soloviev, “Starlink doppler and doppler rate estimation via coherent combining of multiple tones for opportunistic positioning,” in Proceedings of the IEEE/ION Position, Location and Navigation Symposium (PLANS), pp. 1143–1153, 2023.

[3] S. Kozhaya , J. Saroufim , Z. M. Kassas . Unveiling Starlink for PNT. NAVIGATION: Journal of the Institute of Navigation, 2025, 72(1): 1–35

[4] H. Qin , Z. Li , C. Zhao . Integrated positioning technology using iridium/orbcomm signals of opportunity. Journal of Beijing University of Aeronautics and Astronautics, 2022, 48(10): 1845–1853

[5] H. Qin , Z. Tan , L. Cong , C. Zhao . Positioning technology based on iridium signals of opportunity. Journal of Beijing University of Aeronautics and Astronautics, 2019, 45(9): 1691–1699

[6] Z. Cui , F. Yue , R. Tian , S. Zhang . Research on navigation and positioning technology based on iridium burst signals. Global Positioning System, 2021, 46(2): 77–85

[7] Y. Xu, J. Ge, Z. Fang, C. Yuan, X. Chen, M. Gao, and C. Jiang, “A single-search doppler matching positioning method for leo satellite signals of opportunity,” Patent CN 117 148 397 B, Dec. 29, 2023.

[8] J. Khalife , M. Neinavaie , Z. M. Kassas . The first carrier phase tracking and positioning results with starlink leo satellite signals. IEEE Transactions on Aerospace and Electronic Systems, 2022, 58(2): 1487–1491

[9] M. Neinavaie, J. Khalife, and Z. M. Kassas, “Exploiting starlink signals for navigation: First results,” in Proceedings of the 34th International Technical Meeting of the Satellite Division of the Institute of Navigation, pp. 2766–2773, 2021.

[10] S. Kozhaya , S. Hayek , Z. M. Kassas . Cognitive beacon estimation of unknown leo satellite signals of opportunity for PNT. IEEE Journal on Selected Areas in Communications, 2025, 43(2): 456–470

[11] H. Qin , Y. Zhang . Positioning method using starlink signals of opportunity. Journal of Navigation and Positioning, 2023, 11(1): 67–73

[12] M. Neinavaie, Z. Shadram, S. Kozhaya, and Z. M. Kassas, “First results of differential doppler positioning with unknown starlink satellite signals,” in 2022 IEEE Aerospace Conference (AERO), Big Sky, MT, USA, pp. 1–14, 2022.

[13] Y. Yuan , J. Yi , X. Wan , C. Hao . Doppler positioning method and experiment based on starlink beacon signals. Systems Engineering and Electronics, 2024, 46(8): 2535–2545

[14] Y. Xu , Y. Liu , M. Lei , M. Gao , Z. Fang , C. Jiang . Joint pseudo-range and doppler positioning method with leo satellite signals of opportunity. Satellite Navigation, 2025, 6(10): 1–30

[15] Z. Li , D. Wang . Doppler positioning method based on starlink downlink signals. Journal of Ordnance Equipment Engineering, 2025, 46(10): 78–84

[16] Y. Wang , G. Sun , M. Xing , Z. Zhang . Performance analysis and parameter design of passive synthetic aperture localization. Journal of Electronics and Information Technology, 2022, 44(9): 3155–3162

[17] T. E. Humphreys , P. A. Iannucci , Z. M. Komodromos , A. M. Graff . Signal structure of the starlink ku-band downlink. IEEE Transactions on Aerospace and Electronic Systems, 2023, 59(5): 6016–6030

[18] L. Zhang , H. Xun , R. Tao . Single-satellite passive high-precision localization method based on passive synthetic aperture. Aerospace Electronic Warfare, 2020, 36(6): 43–48

[19] N. Bleistein and R. A. Handelsman, Asymptotic Expansions of Integrals. New York, NY, USA: Dover, 1986.

Rights & permissions

The associate editor coordinating the review of this manuscript was Dr. Junhua Yang