# Beacon Aided Positioning Method for Lunar Probes Using Temporal Interferometry

Authors: Yang Liu
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

With a series of complex lunar exploration missions carried out, reliable autonomous navigation approaches for lunar probes are in demand. In this paper, aiming to obtain the position of the probes with high accuracy and computational efficiency, a novel passive positioning scheme for lunar probes is presented. We first designed a reference system including four beacons with known positions. Then, we proposed arctangent temporal interferometry (ATI) method for ambiguity resolution in displacement measurement. Moreover, the quad-beacon displacement search (QBDS) algorithm is proposed for 2-dimension (2D) position search, which enables accurate and computational-efficient 2D positioning for lunar probes. The simulation validations show the effectiveness of the proposed positioning scheme and algorithms.


---

## Body

## 1. Introduction

In recent years, with the development of aerospace science and technology, there are increasing space probes launched for lunar, planetary, and interplanetary research, including Chang’e-5, InSight, and Tianwen-1 [1–4]. The moon is the closest celestial body to the earth, and it is also the first stop for humans to carry out deep-space exploration activities. In the past 60 years, more than a hundred spacecrafts have been launched for lunar exploration [5]. The probe navigation system has always been an important part of a lunar exploration mission [6, 7]. At present, lunar probes mainly rely on ground stations to provide telemetry and navigation services. However, as the moon exploration missions expand from the near to the far side of the moon, the ground telemetry is no longer available. Under such circumstance, reliable autonomous navigation approaches must be applied to ensure the probe has the ability to survive and explore during the mission cycle [8].A number of autonomous navigation schemes have been proposed for lunar exploration missions. The inertial navigation is a commonly used method of autonomous navigation [9, 10]. It does not require any external information. However, its cumulative error leads to poor long-term accuracy. The stellar navigation method is another autonomous navigation method which can be used for lunar probes [11–13]. It obtains the absolute position of the probe by measuring the angular differences between the earth and the navigation stars from the probe. Whereas the positioning accuracy using this method is limited by the accuracy of the sensors, which is not sufficiently high. A third method is to use relay satellites for radio navigation, which measures the distance and Doppler frequency to obtain the relative position and velocity information of the probe [14]. But unlike global navigation satellite system (GNSS) which is available all day, the visible time of lunar orbiting satellites is usually short as the number of satellites is very limited. Therefore, stable and continuous positioning services for lunar probes using relay satellites are infeasible at present. GNSS is also applied for lunar probe positioning, but the geometry of the visible satellites deteriorates greatly as the receiver leaves the earth [15]. There are attempts to build a GNSS-like system orbiting moon; however, the orbit accuracy and cost limits the application of such systems [16].As a commonly used ranging method, carrier phase measurement has been utilized in various fields. For GNSS applications, the carrier phase can be measured to obtain the pseudo-range from the satellite to the receiver [17, 18]. The pseudorange has a very high precision but with an ambiguity due to the whole number of cycles between satellite and receiver is unknown. By resolving the integer ambiguity, the pseudorange with high precision can be utilized in high-precision positioning. Similar method is utilized for micromotion measurements using Doppler radar [19–21]. By measuring the phase of the signal reflected and modulated by the moving target, the displacement of the target referring to the radar sensor is obtained. As the range difference exceeds the wavelength of the carrier, the phase difference is greater than 2π, resulting the phase ambiguity. To solve this problem, a number of phase unwrapping approaches have been proposed, including conventional arctangent demodulation (CAD), extended differentiate and cross-multiplying (DACM), and arcsine demodulation [20–23]. However, the disadvantages of these approaches include (1) high computational complexity and (2) significant accumulation error for fast-moving targets after approximation. These disadvantages limit the application of real-time high-precision ranging using carrier phase.To solve the aforementioned problems, in this paper, a novel passive positioning scheme for lunar probes is proposed and investigated. A reference system including four beacons with known positions is designed which transmits continuous waves in a series of subbands from each beacon. Then, an accurate and computational-efficient algorithm for ambiguity resolution, which is called arctangent temporal interferometry (ATI), is presented for displacement measurement. After that, a target function for 2-dimension (2D) positioning based on the displacements referring to the beacons is constructed. A parallel 2D position search algorithm for lunar probes called quad-beacon displacement search (QBDS) was proposed, which enables accurate and computational-efficient 2D positioning results for lunar probes. The simulation results validate the effectiveness of the proposed scheme.The remainder of this paper is organized as follows. The principle of displacement measurements using carrier phase and ATI algorithm is introduced in Section2. The system configuration and proposed QBDS positioning algorithm of lunar probes is explained in Section 3. Section 4 gives the simulation results and analyzes the positioning performance. The conclusion is drawn in Section 5.


## 2. Displacement Measurement Using Carrier Phase

### 2.1. Signal Model

Consider a beacon transmitting a continuous wave (CW), neglecting the amplitude vibration, we can normalize the amplitude of the transmit signal to 1. Hence, the beacon transmits the signal with the following form,(1)Tt=cos2πft+φt,where f and t denote the carrier frequency of the signal and time, respectively. And ϕ denotes the arbitrary initial phase of the transmitter.After transmission in the free space, the signal is received by the receiver on the lunar probe. The transmission induces a delayτt with respect to the transmit signal, that is,
(2)Rt=Arcos2πft−τt+φt+εt,where Ar is the amplitude of the received signal, εt is the noise term, and delay τt is determined by
(3)τt=rtc,where rt denotes the instantaneous range from the beacon to the probe, and c is the speed of light.Then, the received signal is quadrature demodulated to baseband. Combining Equation (2) and Equation (3), the in-phase (I) and quadrature (Q) outputs are represented by
(4)It=Arcosφt−φr−2πfrtc+εIt,=AIcosφ−2πrtλ+εIt,Qt=Arsinφt−φr−2πfrtc++εQt,=AIsinφ−2πrtλ+εQt,where φr is the arbitrary phase of the receiver, λ=c/f is the carrier wavelength, and +εIt,+εQt denote the noise terms in I/Q channels. We can observe from Equation (4) that the range rt can be resolved by extracting the phase of the received signal Φt, that is,
(5)Φt=φ−2πrtλ.Theoretically, we can obtain the range from the beacon to the probe with high accuracy via carrier phase measurement. And with these range information, the precise position of the receiver can be resolved as long as there are sufficient beacons. In this case, only two beacons are required at least for 2D positioning, and three beacons are required at least for 3D positioning (regardless of finite multiple solutions). Therefore, by placing several reference beacons with known exact coordinates on the moon face, the position of the lunar probes can be resolved by measuring the ranges from the beacons to the probe.However, there are two practical problems for direct range measurement. First, the arbitrary initial phase of the transmitter and receiver is difficult to obtain. Meanwhile, the moving area of the probe is usually far larger than the scale of the microwave wavelength, which leads to2π phase ambiguity. In this scenario, the unambiguous range cannot be obtained directly.In order to solve this problem, we apply range difference measurement instead of direct range measurement. Although we do not have the actual range, we can still resolve the variation of the range by continuously calculating the instantaneous phase of the received signal. The wrapped phaseΦwt can be calculated by taking the arctangent of the received signal, that is,
(6)Φwt=mod2πφ−2πrtλ.To obtain unambiguous phase without2π jumping, unwrapping procedure is required, which is illustrated in Figure 1. A conventional unwrapping method is to perform direct phase shift, which is called conventional arctangent demodulation. Whenever the difference between consecutive phases is greater than or equal to 2π, the following phases are shifted by adding multiples of ±2π until the difference is smaller than π. Taking the first sampled phase as a reference, we have the unwrapped phase
(7)Φut=unwrapΦwt.Figure 1
Illustration of original wrapped phase and unwrapped phase.Therefore, we can calculate the pseudorangerpt from the beacon to the probe.
(8)rpt=−Φutλ2π.It is called pseudorange as the measurement has an unknown offset. However, this offset is not important as it is subtracted when calculating the range difference. Therefore, the range differenceδrt1,t2 between the beacon and the probe from t1 to t2 can be determined by the following equation
(9)δrt1,t2=rpt2−rpt1.The position of the probe at timet1 and t2 can then be resolved using the range differences from the probe to different beacons.


### 2.2. ATI Demodulation

As is shown above, the key to resolve range difference is to perform phase unwrapping. However, the aforementioned conventional arctangent demodulation algorithm has a series of drawbacks. First, the computational complexity is high as the entire data block from the current time to the end needs to be frequently shifted with period2π, which might be computational costly for the embedded system on the probe, as the computing resources are very limited. Also, this algorithm leads to serious accumulative error for a single wrong unwrapping point, where the latter data will have an entire phase offset.In order to reduce the impact of wrong unwrapping, extended DACM algorithm is proposed. This algorithm calculates the phase increment between two adjacent samples and then accumulates the phase increments. Extended DACM algorithm greatly improves the performance facing the wrong unwrapping problem. The basic procedure of the extended DACM algorithm is as follows.To obtain the phase increment, the continuous form is considered first. The angular frequency is the time derivative of the instantaneous phase, that is,(10)ωt=ddtarctanQtIt=ItQ̇t−İtQtIt2+Qt2,where İt and Q̇t denote the time derivative of It and Qt. The angular frequency can then be approximated in a discrete form, which is,
(11)ωk=IkQk−Qk−1/Δt−QkIk−Ik−1/ΔtIk2+Qk2,=Ik−1Qk−IkQk−1Ik2+Qk2Δt,where Δt is the sampling interval. Here, we approximate the time derivative by the slope of adjacent samples. Then, we can derive the unambiguous phase increment ΔΦk=Φk−Φk−1, which has the expression
(12)ΔΦk=ωkΔt=Ik−1Qk−IkQk−1Ik2+Qk2.TakingΦun=0, the unambiguous phase is represented as follows
(13)Φun=∑k=2nΔΦk,=∑k=2nωkΔt,=∑k=2nIk−1Qk−IkQk−1Ik2+Qk2,wheren≥2.In the extended DACM algorithm, the time derivative ofIt and Qt is approximated by the difference between two samples. Then, the phase increment is calculated by multiplying the instantaneous angular frequency ωk and the time interval Δt. These approximations may introduce additional error to the measurement, especially in low signal-to-noise ratio (SNR) scenarios.In order to solve this problem, herein we propose an algorithm which calculates the phase increment between two adjacent samples regardless of the aforementioned assumptions, namely, ATI algorithm.To simplify the phase calculation, we rewrite the baseband signal in complex form,(14)Xt=It+jQt,wherej is the imaginary unit. In discrete domain, it has the form
(15)Xk=Ik+jQk.The phase incrementΔΦk can be resolved by taking the difference of Φk and Φk−1. Assuming that the sampling rate is high enough that the absolute value of phase increment is smaller than π, the phase increment ΔΦk can be derived as the following equation:
(16)ΔΦk=Φk−Φk−1=argXk−argXk−1=argXkX∗k−1=argIk+jQkIk−1+jQk−1∗=argIk+jQkIk−1−jQk−1=argIk−1Ik+Qk−1Qk+jIk−1Qk−IkQk−1=arctanIk−1Qk−IkQk−1Ik−1Ik+Qk−1Qk.Comparing Equation (16) with Equation (11), we can observe that the extended DACM does not consider the variation of signal amplitude over samples. If we assume Ik≈Ik−1,Qk≈Qk−1, and the phase increment satisfies small angle assumption, i.e., ΔΦk≈ΔΦk when ΔΦk is sufficiently small, Equation (16) and Equation (11) can have a uniform form.Then, similar to Equation (13), taking Φun=0, the unambiguous phase can be represented by
(17)Φun=∑k=2nΔΦk,=∑k=2narctanIk−1Qk−IkQk−1Ik−1Ik+Qk−1Qk,wheren≥2.As we already have the unambiguous phase, we can calculate the range difference according to Equation (9), that is,
(18)Δrm,n=rpn−rpm,=−λ2πΦun−Φum.


## 2.1. Signal Model

Consider a beacon transmitting a continuous wave (CW), neglecting the amplitude vibration, we can normalize the amplitude of the transmit signal to 1. Hence, the beacon transmits the signal with the following form,(1)Tt=cos2πft+φt,where f and t denote the carrier frequency of the signal and time, respectively. And ϕ denotes the arbitrary initial phase of the transmitter.After transmission in the free space, the signal is received by the receiver on the lunar probe. The transmission induces a delayτt with respect to the transmit signal, that is,
(2)Rt=Arcos2πft−τt+φt+εt,where Ar is the amplitude of the received signal, εt is the noise term, and delay τt is determined by
(3)τt=rtc,where rt denotes the instantaneous range from the beacon to the probe, and c is the speed of light.Then, the received signal is quadrature demodulated to baseband. Combining Equation (2) and Equation (3), the in-phase (I) and quadrature (Q) outputs are represented by
(4)It=Arcosφt−φr−2πfrtc+εIt,=AIcosφ−2πrtλ+εIt,Qt=Arsinφt−φr−2πfrtc++εQt,=AIsinφ−2πrtλ+εQt,where φr is the arbitrary phase of the receiver, λ=c/f is the carrier wavelength, and +εIt,+εQt denote the noise terms in I/Q channels. We can observe from Equation (4) that the range rt can be resolved by extracting the phase of the received signal Φt, that is,
(5)Φt=φ−2πrtλ.Theoretically, we can obtain the range from the beacon to the probe with high accuracy via carrier phase measurement. And with these range information, the precise position of the receiver can be resolved as long as there are sufficient beacons. In this case, only two beacons are required at least for 2D positioning, and three beacons are required at least for 3D positioning (regardless of finite multiple solutions). Therefore, by placing several reference beacons with known exact coordinates on the moon face, the position of the lunar probes can be resolved by measuring the ranges from the beacons to the probe.However, there are two practical problems for direct range measurement. First, the arbitrary initial phase of the transmitter and receiver is difficult to obtain. Meanwhile, the moving area of the probe is usually far larger than the scale of the microwave wavelength, which leads to2π phase ambiguity. In this scenario, the unambiguous range cannot be obtained directly.In order to solve this problem, we apply range difference measurement instead of direct range measurement. Although we do not have the actual range, we can still resolve the variation of the range by continuously calculating the instantaneous phase of the received signal. The wrapped phaseΦwt can be calculated by taking the arctangent of the received signal, that is,
(6)Φwt=mod2πφ−2πrtλ.To obtain unambiguous phase without2π jumping, unwrapping procedure is required, which is illustrated in Figure 1. A conventional unwrapping method is to perform direct phase shift, which is called conventional arctangent demodulation. Whenever the difference between consecutive phases is greater than or equal to 2π, the following phases are shifted by adding multiples of ±2π until the difference is smaller than π. Taking the first sampled phase as a reference, we have the unwrapped phase
(7)Φut=unwrapΦwt.Figure 1
Illustration of original wrapped phase and unwrapped phase.Therefore, we can calculate the pseudorangerpt from the beacon to the probe.
(8)rpt=−Φutλ2π.It is called pseudorange as the measurement has an unknown offset. However, this offset is not important as it is subtracted when calculating the range difference. Therefore, the range differenceδrt1,t2 between the beacon and the probe from t1 to t2 can be determined by the following equation
(9)δrt1,t2=rpt2−rpt1.The position of the probe at timet1 and t2 can then be resolved using the range differences from the probe to different beacons.


## 2.2. ATI Demodulation

As is shown above, the key to resolve range difference is to perform phase unwrapping. However, the aforementioned conventional arctangent demodulation algorithm has a series of drawbacks. First, the computational complexity is high as the entire data block from the current time to the end needs to be frequently shifted with period2π, which might be computational costly for the embedded system on the probe, as the computing resources are very limited. Also, this algorithm leads to serious accumulative error for a single wrong unwrapping point, where the latter data will have an entire phase offset.In order to reduce the impact of wrong unwrapping, extended DACM algorithm is proposed. This algorithm calculates the phase increment between two adjacent samples and then accumulates the phase increments. Extended DACM algorithm greatly improves the performance facing the wrong unwrapping problem. The basic procedure of the extended DACM algorithm is as follows.To obtain the phase increment, the continuous form is considered first. The angular frequency is the time derivative of the instantaneous phase, that is,(10)ωt=ddtarctanQtIt=ItQ̇t−İtQtIt2+Qt2,where İt and Q̇t denote the time derivative of It and Qt. The angular frequency can then be approximated in a discrete form, which is,
(11)ωk=IkQk−Qk−1/Δt−QkIk−Ik−1/ΔtIk2+Qk2,=Ik−1Qk−IkQk−1Ik2+Qk2Δt,where Δt is the sampling interval. Here, we approximate the time derivative by the slope of adjacent samples. Then, we can derive the unambiguous phase increment ΔΦk=Φk−Φk−1, which has the expression
(12)ΔΦk=ωkΔt=Ik−1Qk−IkQk−1Ik2+Qk2.TakingΦun=0, the unambiguous phase is represented as follows
(13)Φun=∑k=2nΔΦk,=∑k=2nωkΔt,=∑k=2nIk−1Qk−IkQk−1Ik2+Qk2,wheren≥2.In the extended DACM algorithm, the time derivative ofIt and Qt is approximated by the difference between two samples. Then, the phase increment is calculated by multiplying the instantaneous angular frequency ωk and the time interval Δt. These approximations may introduce additional error to the measurement, especially in low signal-to-noise ratio (SNR) scenarios.In order to solve this problem, herein we propose an algorithm which calculates the phase increment between two adjacent samples regardless of the aforementioned assumptions, namely, ATI algorithm.To simplify the phase calculation, we rewrite the baseband signal in complex form,(14)Xt=It+jQt,wherej is the imaginary unit. In discrete domain, it has the form
(15)Xk=Ik+jQk.The phase incrementΔΦk can be resolved by taking the difference of Φk and Φk−1. Assuming that the sampling rate is high enough that the absolute value of phase increment is smaller than π, the phase increment ΔΦk can be derived as the following equation:
(16)ΔΦk=Φk−Φk−1=argXk−argXk−1=argXkX∗k−1=argIk+jQkIk−1+jQk−1∗=argIk+jQkIk−1−jQk−1=argIk−1Ik+Qk−1Qk+jIk−1Qk−IkQk−1=arctanIk−1Qk−IkQk−1Ik−1Ik+Qk−1Qk.Comparing Equation (16) with Equation (11), we can observe that the extended DACM does not consider the variation of signal amplitude over samples. If we assume Ik≈Ik−1,Qk≈Qk−1, and the phase increment satisfies small angle assumption, i.e., ΔΦk≈ΔΦk when ΔΦk is sufficiently small, Equation (16) and Equation (11) can have a uniform form.Then, similar to Equation (13), taking Φun=0, the unambiguous phase can be represented by
(17)Φun=∑k=2nΔΦk,=∑k=2narctanIk−1Qk−IkQk−1Ik−1Ik+Qk−1Qk,wheren≥2.As we already have the unambiguous phase, we can calculate the range difference according to Equation (9), that is,
(18)Δrm,n=rpn−rpm,=−λ2πΦun−Φum.


## 3. 2D Positioning Scheme for Lunar Probes

### 3.1. System Configuration

As mentioned in Section II, the position of the probe at timet1 and t2 can be resolved using the range differences from the probe to different beacons as long as there are sufficient number of beacons.Suppose the probe moves on a 2D plane, the trajectory starts fromrt1=xt1,yt1 and ends at rt2=xt2,yt2. As we can see, different from direct range positioning where only one coordinate needs to be resolved, there are two unknown coordinates, comprising four parameters xt1,yt1,xt2,yt2. Therefore, if we can obtain at least four range differences, it is possible to determine the position of the probe.In order to obtain a fine geometric dilution of precision (GDOP) and simplify the configuration, we set four beacons on the four vertices of a square, and the moving range of the probe is limited inside that square area. The positioning scene for lunar probe positioning is illustrated in Figure2.Figure 2
Positioning scene for lunar probe positioning.


### 3.2. The Description of the Positioning Problem

Suppose the four beacons are distributed on the four vertices of a square with side-lengthD, we can establish a Cartesian coordinate system referring to the beacons. The locations of the beacons rbk=xbk,ybk are as follows: rb1=0,0,rb2=D,0,rb3=D,D,rb4=0,D. The probe moves from the initial position ri=xi,yi to the terminal position rt=xt,yt, and the range difference from the initial position to the terminal position with respect to the four beacons is δrk,k=1,2,3,4. We can locate the probe by solving the following equations
(19)rpt−rbk2−rpi−rbk2=Δrk,wherek=1,2,3,4.Equation (19) can be further written in a scalar form.
(20)xt−02+yt−02−xi−02+yi−02=δr1,xt−D2+yt−02−xi−D2+yi−02=δr2,xt−D2+yt−D2−xi−D2+yi−D2=δr3,xt−02+yt−D2−xi−02+yi−D2=δr4.As we can see, by solving the four equations in Equation (20), the initial position rpi=xi,yi and the terminal position rpt=xt,yt can be estimated. However, Equation (20) is a set of nonlinear equations, which are difficult to solve the initial and terminal positions of the probe directly. In order to solve the initial and terminal positions of the probe, the following algorithm can be adapted.Consider the probe is originally located at the initial positionri. For the four beacons k=1,2,3,4, we can obtain the initial range rk from the probe to the beacons.
(21)rk=xt−xbk2+yt−ybk2.Then, we can then obtain the terminal range from the probe to the beacons, which can be calculated by(22)rk′=rk+δrk,where the range difference δrk can be measured by applying ATI algorithm. Given two fixed reference points on the plane and the range from the third point to the two points, we can easily solve for the coordinates of the third point using law of cosines and we can divide the beacons into two separate groups to solve for the terminal position rt, respectively. For instance, Group A consists of beacons 1 and 2, and Group B consists of beacons 3 and 4.Here, we derive the expression of terminal position calculated by beacons Group A. The geometry of the beacons and the probe is illustrated in Figure3. The terminal angle of the probe with respect to Beacon 1 can be calculated by
(23)cosθ=D2+r1′2−r2′22Dr1′.Figure 3
The geometry of the beacons and the probe of Group A.Then, we derive the coordinates of terminal positionrptA=xtA,ytA. Theoretically, there are two solutions of the terminal position. The solution inside the feasible region is as follows.
(24)xtA=r1′cosθ,=D2+r1′2−r2′22D,ytA=r1′2−xtA2,=r1′2−D2+r1′2−r2′224Dr1′2,=2D2r1′2+2D2r2′2+2r1′2r2′2−D4−r1′4−r2′42D.Similarly, the coordinates of terminal position solved using beacons Group BrptB=xtB,ytB are
(25)xtB=D2−r3′2+r4′22D,ytB=D−2D2r3′2+2D2r4′2+2r3′2r4′2−D4−r3′4−r4′42D.If the initial positionrpi is correct, the terminal position calculated by two groups should be identical. Therefore, the problem is transformed into finding an initial position rpi in the feasible region, so that the terminal positions rptA and rptB overlap. In other words, the distance between rptA and rptB can be used to measure the accuracy of the given initial position rpi. Here, we define a target function δp, which is defined by the following equation
(26)δp=rptA−rptB22,=xtA−xtB2+ytA−ytB2.Combining Equations (21), (22), (24), (25), and (26), the cost function δp can be represented by a function of rpi. As the specific form of the equation is very complicated, here, the cost function is simply denoted by
(27)δp=δrpi.And the problem becomes an optimization problem, that is,(28)minrpiδprpisubjecttorpi,rptA,rptB∈Swhere S is the feasible region, which is the area enclosed by four beacons.


### 3.3. Newton-Raphson Method Based Algorithm

We can apply Newton-Raphson method to solve this optimization problem. Finding the minimum of the functionδrpi requires finding the critical point which satisfies ∇δrpi=0. After setting an initial solution rpi0, we can use the following formula for iteration. Each step generates a new iterate rpik+1 from the previous iterate rpik according to
(29)rpik+1=rpik−Hδrpik−1∇δrpik,where Hδrpik is the Hessian matrix defined by
(30)Hδrpi=∂2δrpi∂x2∂2δrpi∂x∂y∂2δrpi∂y∂x∂2δrpi∂y2,and ∇⋅ is the gradient of δrpi defined by,
(31)∇δrpi=∂δrpi∂x∂δrpi∂y.After each iteration, judge whether the difference between the two iterates is smaller than the threshold, i.e.,(32)∣δrpik+1−δrpik∣<ε0.If Equation (32) stands, output δrpik+1 as the final solution of the initial position. Otherwise, continue the iteration. The whole process of initial position solution process is illustrated in Table 1.Table 1
Initial position solution process using Newton-Raphson method.
Step 1Generate an initial solutionrpi0.Step 2Calculate the target functionδrpik according to Equation (27).Step 3Verify the solution has converged according to Equation (32) (k≥1).Abort iteration and outputrpik if the difference between thetwo iterates is smaller than the threshold.Step 4Calculaterpik+1 according to Equation (29).Step 5Return to step 2 and start next iteration.


### 3.4. QBDS Algorithm

Despite the initial position of the probe can be solved quickly by applying Newton-Raphson method, the calculation of the Hessian function is quite difficult. Also, Newton-Raphson method may not be able to obtain the global optimal solution ofrpi. Within the allowable range of accuracy, we can use a method based on grid search to solve the optimization problem in Equation (28), namely, QBDS.As a simple but effective method, we can perform grid search to the initial positionrpi, and calculate the target function δp corresponding to each on-grid rpi. By finding the minimum δp, we can solve the true initial position rpi. Figure 4 gives an intuitive explanation of the QBDS process.Figure 4
Positioning scene for lunar probe positioning.For convenience in presenting the figures, here, we define another target functionτrpi=1/δprpi. The problem becomes
(33)maxrpiτrpisubjecttorpi,rptA,rptB∈SThen, we can locate the probe in 2D plane whenτ reaches its maximum.We can observe from the QBDS algorithm that after calculating the range differencesδrk, the searching procedure can be run in a parallel form. This feature enables fast position resolution for real-time positioning, which helps the probe to survive in complex lunar environment.


## 3.1. System Configuration

As mentioned in Section II, the position of the probe at timet1 and t2 can be resolved using the range differences from the probe to different beacons as long as there are sufficient number of beacons.Suppose the probe moves on a 2D plane, the trajectory starts fromrt1=xt1,yt1 and ends at rt2=xt2,yt2. As we can see, different from direct range positioning where only one coordinate needs to be resolved, there are two unknown coordinates, comprising four parameters xt1,yt1,xt2,yt2. Therefore, if we can obtain at least four range differences, it is possible to determine the position of the probe.In order to obtain a fine geometric dilution of precision (GDOP) and simplify the configuration, we set four beacons on the four vertices of a square, and the moving range of the probe is limited inside that square area. The positioning scene for lunar probe positioning is illustrated in Figure2.Figure 2
Positioning scene for lunar probe positioning.


## 3.2. The Description of the Positioning Problem

Suppose the four beacons are distributed on the four vertices of a square with side-lengthD, we can establish a Cartesian coordinate system referring to the beacons. The locations of the beacons rbk=xbk,ybk are as follows: rb1=0,0,rb2=D,0,rb3=D,D,rb4=0,D. The probe moves from the initial position ri=xi,yi to the terminal position rt=xt,yt, and the range difference from the initial position to the terminal position with respect to the four beacons is δrk,k=1,2,3,4. We can locate the probe by solving the following equations
(19)rpt−rbk2−rpi−rbk2=Δrk,wherek=1,2,3,4.Equation (19) can be further written in a scalar form.
(20)xt−02+yt−02−xi−02+yi−02=δr1,xt−D2+yt−02−xi−D2+yi−02=δr2,xt−D2+yt−D2−xi−D2+yi−D2=δr3,xt−02+yt−D2−xi−02+yi−D2=δr4.As we can see, by solving the four equations in Equation (20), the initial position rpi=xi,yi and the terminal position rpt=xt,yt can be estimated. However, Equation (20) is a set of nonlinear equations, which are difficult to solve the initial and terminal positions of the probe directly. In order to solve the initial and terminal positions of the probe, the following algorithm can be adapted.Consider the probe is originally located at the initial positionri. For the four beacons k=1,2,3,4, we can obtain the initial range rk from the probe to the beacons.
(21)rk=xt−xbk2+yt−ybk2.Then, we can then obtain the terminal range from the probe to the beacons, which can be calculated by(22)rk′=rk+δrk,where the range difference δrk can be measured by applying ATI algorithm. Given two fixed reference points on the plane and the range from the third point to the two points, we can easily solve for the coordinates of the third point using law of cosines and we can divide the beacons into two separate groups to solve for the terminal position rt, respectively. For instance, Group A consists of beacons 1 and 2, and Group B consists of beacons 3 and 4.Here, we derive the expression of terminal position calculated by beacons Group A. The geometry of the beacons and the probe is illustrated in Figure3. The terminal angle of the probe with respect to Beacon 1 can be calculated by
(23)cosθ=D2+r1′2−r2′22Dr1′.Figure 3
The geometry of the beacons and the probe of Group A.Then, we derive the coordinates of terminal positionrptA=xtA,ytA. Theoretically, there are two solutions of the terminal position. The solution inside the feasible region is as follows.
(24)xtA=r1′cosθ,=D2+r1′2−r2′22D,ytA=r1′2−xtA2,=r1′2−D2+r1′2−r2′224Dr1′2,=2D2r1′2+2D2r2′2+2r1′2r2′2−D4−r1′4−r2′42D.Similarly, the coordinates of terminal position solved using beacons Group BrptB=xtB,ytB are
(25)xtB=D2−r3′2+r4′22D,ytB=D−2D2r3′2+2D2r4′2+2r3′2r4′2−D4−r3′4−r4′42D.If the initial positionrpi is correct, the terminal position calculated by two groups should be identical. Therefore, the problem is transformed into finding an initial position rpi in the feasible region, so that the terminal positions rptA and rptB overlap. In other words, the distance between rptA and rptB can be used to measure the accuracy of the given initial position rpi. Here, we define a target function δp, which is defined by the following equation
(26)δp=rptA−rptB22,=xtA−xtB2+ytA−ytB2.Combining Equations (21), (22), (24), (25), and (26), the cost function δp can be represented by a function of rpi. As the specific form of the equation is very complicated, here, the cost function is simply denoted by
(27)δp=δrpi.And the problem becomes an optimization problem, that is,(28)minrpiδprpisubjecttorpi,rptA,rptB∈Swhere S is the feasible region, which is the area enclosed by four beacons.


## 3.3. Newton-Raphson Method Based Algorithm

We can apply Newton-Raphson method to solve this optimization problem. Finding the minimum of the functionδrpi requires finding the critical point which satisfies ∇δrpi=0. After setting an initial solution rpi0, we can use the following formula for iteration. Each step generates a new iterate rpik+1 from the previous iterate rpik according to
(29)rpik+1=rpik−Hδrpik−1∇δrpik,where Hδrpik is the Hessian matrix defined by
(30)Hδrpi=∂2δrpi∂x2∂2δrpi∂x∂y∂2δrpi∂y∂x∂2δrpi∂y2,and ∇⋅ is the gradient of δrpi defined by,
(31)∇δrpi=∂δrpi∂x∂δrpi∂y.After each iteration, judge whether the difference between the two iterates is smaller than the threshold, i.e.,(32)∣δrpik+1−δrpik∣<ε0.If Equation (32) stands, output δrpik+1 as the final solution of the initial position. Otherwise, continue the iteration. The whole process of initial position solution process is illustrated in Table 1.Table 1
Initial position solution process using Newton-Raphson method.
Step 1Generate an initial solutionrpi0.Step 2Calculate the target functionδrpik according to Equation (27).Step 3Verify the solution has converged according to Equation (32) (k≥1).Abort iteration and outputrpik if the difference between thetwo iterates is smaller than the threshold.Step 4Calculaterpik+1 according to Equation (29).Step 5Return to step 2 and start next iteration.


## 3.4. QBDS Algorithm

Despite the initial position of the probe can be solved quickly by applying Newton-Raphson method, the calculation of the Hessian function is quite difficult. Also, Newton-Raphson method may not be able to obtain the global optimal solution ofrpi. Within the allowable range of accuracy, we can use a method based on grid search to solve the optimization problem in Equation (28), namely, QBDS.As a simple but effective method, we can perform grid search to the initial positionrpi, and calculate the target function δp corresponding to each on-grid rpi. By finding the minimum δp, we can solve the true initial position rpi. Figure 4 gives an intuitive explanation of the QBDS process.Figure 4
Positioning scene for lunar probe positioning.For convenience in presenting the figures, here, we define another target functionτrpi=1/δprpi. The problem becomes
(33)maxrpiτrpisubjecttorpi,rptA,rptB∈SThen, we can locate the probe in 2D plane whenτ reaches its maximum.We can observe from the QBDS algorithm that after calculating the range differencesδrk, the searching procedure can be run in a parallel form. This feature enables fast position resolution for real-time positioning, which helps the probe to survive in complex lunar environment.


## 4. Simulation Results and Analysis

In order to verify the feasibility and effectiveness of the theoretical part in this paper, a series of simulation is presented in this section.


### 4.1. Phase Unwrapping

To evaluate the phase unwrapping performance of ATI method, a scenario with linear range increment is considered. Suppose the lunar probe is originally located at the position of a beacon and moves away from the beacon with constant velocityv=1 m/s. The carrier frequency is set f=1000 MHz. Three phase unwrapping algorithms including extended DACM, arcsine demodulation, and ATI are compared in the simulation. The curves of the mean square error (MSE) of the final position versus SNR are plotted in Figure 5, where the SNR is changing from 10 dB to 40 dB and each point is obtained utilizing 104 Monte-Carlo simulations.Figure 5
The MSE versus SNR after 104 Monte-Carlo simulations.As shown in Figure5, the performance of ATI deteriorates for SNR lower than 12 dB. For SNR higher than 14 dB, ATI performs best among the three algorithms. All three algorithms approach similar SNR in phase unwrapping as the SNR approaches 40 dB. In this scenario, the algorithm for phase unwrapping has little effect on the result. The simulation validates the performance of the presented ATI algorithm for SNR higher than 14 dB.


### 4.2. Positioning

Suppose the lunar probe is restricted in a square area with lengthD=100 m. We set four beacons on the four vertices of that area. Therefore, the locations of the beacons are, rb1=0,0,rb2=100,0,rb3=100,100,rb4=0,100 (m). In order to prevent signals from different beacons from interfering with each other, we adopt a frequency division multiplexing mode to manage the signals. We set the carrier frequency near 1000 MHz and a frequency interval of 1 MHz, so we can separate signals from four beacons with different band-pass filters. The carrier frequency of the beacons is set as follows: f1=1000 MHz, f2=1001 MHz, f3=1002 MHz, and f4=1003 MHz. As the Doppler frequency of the lunar probe is far smaller than the frequency interval, we can recover the received signal from each beacon without overlapping. The simulation parameters are shown in Table 2.Table 2
Simulation parameters for lunar probe positioning.
ParameterValueParameterValueBeacon spacingD100 mFrequency intervalΔf1 MHzLocation of Beacon 1rb1(0 m, 0 m)Carrier frequency of Beacon 1f11000 MHzLocation of Beacon 2rb2(0 m, 100 m)Carrier frequency of Beacon 2f21001 MHzLocation of Beacon 3rb3(100 m, 100 m)Carrier frequency of Beacon 3f31002 MHzLocation of Beacon 4rb4(100 m, 0 m)Carrier frequency of Beacon 4f41003 MHzThree typical trajectories of the probe are investigated, those are (a) straight line; (b) semicircle; and (c) polyline, which are illustrated in Figure6. The simulation time is 10 s, and the detailed parameters of the three trajectories are listed in Table 3.Figure 6
Three typical trajectories investigated in the simulation.
(a)
Straight line(b)
Semicircle(c)
PolylineTable 3
Detailed parameters of the three trajectories.
TrajectoryInitial positionrpi (m)Terminal positionrpt (m)Velocity(a)Straight line(50.50)(70.50)2 m/s(b)Semicircle(50.50)(70.50)π m/s(c)Polyline(50.50)(70.50)4 m/s (1st segment) 2 m/s (2nd segment)To illustrate the influence of the SNR on range difference estimation and positioning, the SNR is set to 10 dB, 15 dB, and 20 dB, respectively. The range differenceδrk resolved by phase unwrapping and the distribution of the target function τ are presented in Figures 7–9. The positioning result of the terminal position rpt indicated by the maximum point of the distribution of τ is shown in Table 4.Figure 7
Range differenceδrk and distribution of target function τ of straight line.
(a)
δrk at 10 dB(b)
τ at 10 dB(c)
δrk at 15 dB(d)
τ at 15 dB](e)
δrk at 20 dB(f)
τ at 20 dBFigure 8
Range differenceδrk and distribution of target function τ of semicircle.
(a)
δrk at 10 dB(b)
τ at 10 dB(c)
δrk at 15 dB(d)
τ at 15 dB(e)
δrk at 20 dB(f)
τ at 20 dBFigure 9
Range differenceδrk and distribution of target function τ of polyline.
(a)
δrk at 10 dB(b)
τ at 10 dB(c)
δrk at 15 dB(d)
τ at 15 dB(e)
δrk at 20 dB(f)
τ at 20 dBTable 4
The positioning result of the terminal positionrpt.
TrajectorySNR (dB)Terminal position (m)Position error (m)(a)Straight line10(69.47)3.2(a)Straight line15(70.50)0(a)Straight line20(70.50)0(b)Semicircle10(69.48)1.6(b)Semicircle15(70.50)0(b)Semicircle20(70.50)0(c)Polyline10(69.64)4.1(c)Polyline15(70.60)0(c)Polyline20(70.60)0It can be observed that QBDS method works for all three trajectories. ForSNR=10dB, the large noise induces wrong unwrapping results, hence reduce positioning accuracy significantly. For SNR=15dB and SNR=20dB, no wrong unwrapping results, and the error of range difference is only induced by the phase noise within the 2π period. As the SNR increases, the distribution map of τ becomes more concentrated, which may provide higher positioning accuracy for smaller search grids.


## 4.1. Phase Unwrapping

To evaluate the phase unwrapping performance of ATI method, a scenario with linear range increment is considered. Suppose the lunar probe is originally located at the position of a beacon and moves away from the beacon with constant velocityv=1 m/s. The carrier frequency is set f=1000 MHz. Three phase unwrapping algorithms including extended DACM, arcsine demodulation, and ATI are compared in the simulation. The curves of the mean square error (MSE) of the final position versus SNR are plotted in Figure 5, where the SNR is changing from 10 dB to 40 dB and each point is obtained utilizing 104 Monte-Carlo simulations.Figure 5
The MSE versus SNR after 104 Monte-Carlo simulations.As shown in Figure5, the performance of ATI deteriorates for SNR lower than 12 dB. For SNR higher than 14 dB, ATI performs best among the three algorithms. All three algorithms approach similar SNR in phase unwrapping as the SNR approaches 40 dB. In this scenario, the algorithm for phase unwrapping has little effect on the result. The simulation validates the performance of the presented ATI algorithm for SNR higher than 14 dB.


## 4.2. Positioning

Suppose the lunar probe is restricted in a square area with lengthD=100 m. We set four beacons on the four vertices of that area. Therefore, the locations of the beacons are, rb1=0,0,rb2=100,0,rb3=100,100,rb4=0,100 (m). In order to prevent signals from different beacons from interfering with each other, we adopt a frequency division multiplexing mode to manage the signals. We set the carrier frequency near 1000 MHz and a frequency interval of 1 MHz, so we can separate signals from four beacons with different band-pass filters. The carrier frequency of the beacons is set as follows: f1=1000 MHz, f2=1001 MHz, f3=1002 MHz, and f4=1003 MHz. As the Doppler frequency of the lunar probe is far smaller than the frequency interval, we can recover the received signal from each beacon without overlapping. The simulation parameters are shown in Table 2.Table 2
Simulation parameters for lunar probe positioning.
ParameterValueParameterValueBeacon spacingD100 mFrequency intervalΔf1 MHzLocation of Beacon 1rb1(0 m, 0 m)Carrier frequency of Beacon 1f11000 MHzLocation of Beacon 2rb2(0 m, 100 m)Carrier frequency of Beacon 2f21001 MHzLocation of Beacon 3rb3(100 m, 100 m)Carrier frequency of Beacon 3f31002 MHzLocation of Beacon 4rb4(100 m, 0 m)Carrier frequency of Beacon 4f41003 MHzThree typical trajectories of the probe are investigated, those are (a) straight line; (b) semicircle; and (c) polyline, which are illustrated in Figure6. The simulation time is 10 s, and the detailed parameters of the three trajectories are listed in Table 3.Figure 6
Three typical trajectories investigated in the simulation.
(a)
Straight line(b)
Semicircle(c)
PolylineTable 3
Detailed parameters of the three trajectories.
TrajectoryInitial positionrpi (m)Terminal positionrpt (m)Velocity(a)Straight line(50.50)(70.50)2 m/s(b)Semicircle(50.50)(70.50)π m/s(c)Polyline(50.50)(70.50)4 m/s (1st segment) 2 m/s (2nd segment)To illustrate the influence of the SNR on range difference estimation and positioning, the SNR is set to 10 dB, 15 dB, and 20 dB, respectively. The range differenceδrk resolved by phase unwrapping and the distribution of the target function τ are presented in Figures 7–9. The positioning result of the terminal position rpt indicated by the maximum point of the distribution of τ is shown in Table 4.Figure 7
Range differenceδrk and distribution of target function τ of straight line.
(a)
δrk at 10 dB(b)
τ at 10 dB(c)
δrk at 15 dB(d)
τ at 15 dB](e)
δrk at 20 dB(f)
τ at 20 dBFigure 8
Range differenceδrk and distribution of target function τ of semicircle.
(a)
δrk at 10 dB(b)
τ at 10 dB(c)
δrk at 15 dB(d)
τ at 15 dB(e)
δrk at 20 dB(f)
τ at 20 dBFigure 9
Range differenceδrk and distribution of target function τ of polyline.
(a)
δrk at 10 dB(b)
τ at 10 dB(c)
δrk at 15 dB(d)
τ at 15 dB(e)
δrk at 20 dB(f)
τ at 20 dBTable 4
The positioning result of the terminal positionrpt.
TrajectorySNR (dB)Terminal position (m)Position error (m)(a)Straight line10(69.47)3.2(a)Straight line15(70.50)0(a)Straight line20(70.50)0(b)Semicircle10(69.48)1.6(b)Semicircle15(70.50)0(b)Semicircle20(70.50)0(c)Polyline10(69.64)4.1(c)Polyline15(70.60)0(c)Polyline20(70.60)0It can be observed that QBDS method works for all three trajectories. ForSNR=10dB, the large noise induces wrong unwrapping results, hence reduce positioning accuracy significantly. For SNR=15dB and SNR=20dB, no wrong unwrapping results, and the error of range difference is only induced by the phase noise within the 2π period. As the SNR increases, the distribution map of τ becomes more concentrated, which may provide higher positioning accuracy for smaller search grids.


## 5. Conclusion

In conclusion, we proposed a positioning method and system for lunar probes. A novel phase unwrapping algorithm named ATI algorithm was presented, which calculates the unambiguous phase with higher precision and lower computational complexity. Then, a 2D positioning scheme for lunar probes is proposed, where we configure a set of square distributed beacons and apply QBDS algorithm for lunar probe positioning. Finally, the positioning method was evaluated by simulation. The feasibility and effectiveness of the method was validated by the simulation results.



---

*Source: 6201703-2022-03-30.xml*