Instant Determination of Polar Motion with Tri-static Common View Lunar Laser Ranging

1.   INTRODUCTION

Determining Earth’s orientation in inertial space, also known as the Earth Rotation Parameters (ERPs), is important in astrometric studies that relate terrestrial and celestial coordinates. The ERP data consists of three numbers: the polar motion (PM) components ${x}_{p}$ and ${y}_{p}$, and the Universal Time (UT1). Currently, ERPs are usually provided using the Very Long Baseline Interferometry (VLBI) technique. The international VLBI service can determine PM to accuracy of 50 $\mathrm{\mu }\mathrm{a}\mathrm{s}$ and UT1 to accuracy of 3 $\mathrm{\mu }\mathrm{s}$^{[1, 2]}, and delivery of the ERP product with temporal resolution of 1 day normally requires 8–10 days^[2].

Determination of ERPs can also be conducted using the lunar laser ranging (LLR) technique, which involves measurement of the round-trip time between a ground-based laser ranging station (hereafter, station) and a lunar retro-reflector array (hereafter, reflector). The use of LLR data to determine ERPs was first explored soon after the data initially became available^[3–5]. Over the past decade, a number of LLR stations have produced precise data, with range precision enhanced to the centimeter level^{[6, 7]}. Recent LLR analysis studies typically selected observation nights with 10–15 normal points (NPs), with the minimum requirement of 5 NPs per night^[8–10]. The most up-to-date studies have achieved accuracy of several milliarcseconds for PM coordinates, with temporal resolution of 1 day^{[8, 10]}.

The concept of organizing LLR observations in a common view manner, or LLRCV, was proposed by Leick^[11]. It was claimed that the LLRCV method could provide ERP components with measurement accuracy on the Earth’s surface of approximately 30 cm for LLR in the early 1980s, equivalent to 10 mas of PM or 0.6 ms of UT1^[11]. Another form of LLRCV, proposed by Müller et al.^[12], involves placement of an optical transponder on the lunar surface to illuminate multiple laser ranging stations simultaneously. However, follow-up studies on the LLRCV concept have been limited.

Major geodynamic events, such as earthquakes, can cause rapid changes in ERPs^{[13, 14]}. While conventional LLR analysis extracts a daily solution using single-night observations, the common view method provides an instant solution at higher temporal resolution, which is defined by the accumulation time of LLR NPs, i.e., typically 15 min or less. This represents a substantial improvement in monitoring capability for transient geodynamic phenomena by LLR.

In this paper, we explore the instant determination of Earth’s PM using LLR data in the form of tri-static common view (TCV) measurements. Linear equations can be formulated upon three simultaneous range residuals from three cooperating stations to solve for PM in near real time.

Three European stations were selected from the International Laser Ranging Service network to form a station group that we named GMW. The stations are located in Grasse in France (GRSM/7845), Matera in Italy (MATM/7941), and Wettzell in Germany (WETL/8834). The baselines between the stations are 753 km (GW), 877 km (GM), and 990 km (MW), see Fig. 1.

Figure  1.  LLR stations in Europe.

DownLoad: Full-Size Img PowerPoint

Our approach is grounded in contemporary geodetic infrastructures, which include the high-precision global Earth reference frame (ITRF)^[15], International Earth Rotation Service^[16], high-precision planetary ephemerides^[17], fundamental astronomical software^{[18, 19]}, and the International Earth Rotation Service (IERS) conventions 2010, which provide key algorithms^{[20, 21]}. The reference ERP data were extracted from the IERS C04 20 series, which was accessed via the IERS official website^[16].

2.   METHOD AND DATA

This paper explores instant determination of Earth’s PM using LLR data in the form of TCV measurements.

We first established the method for determining Earth’s PM from the LLR model, the range derivatives to ERPs, the PM prediction model, and the least squares solution. Then, we conducted a two-phase study to evaluate the method. The first or simulation phase involved generating LLR TCV simulation data to test the robustness of the PM determination method by introducing UT1 and range errors. In the second or augmentation phase, we combined simulated data with actual LLR NPs to create realistic TCV events and solve for PM. These results were compared with reference values to generate solution error data. Statistical analysis was conducted to assess the accuracy of the solutions.

2.1   LLR Model and Data

The lunar ranging measurement can be formulated as follows:

where $\tau$ represents the time of flight measured at the ranging station. The distances ${s}_{12}$ and ${s}_{23}$ represent the light path from the station to the reflector and from the reflector to the station, respectively, accounted in barycentric coordinates. The relativistic delay ${\tau }_{\mathrm{G}\mathrm{R}}$ represents the gravitational delay effect on the flying pulse. The tropospheric delay ${\tau }_{\mathrm{t}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{o}}$ represents atmospheric refraction delay during the propagation of light. Tags 12 and 23 represent the outbound flight and the inbound flight, respectively. An empirical term $\Delta {\tau }_{\mathrm{e}\mathrm{m}\mathrm{p}}$ was introduced and fitted to compensate for the remaining periodic effect.

Time-of-flight measurements are performed by ground stations in the clock rate of International Atomic Time (TAI), while planetary motions are calculated using barycentric dynamical time (TDB). The difference in the clock rate between TAI and TDB is considered negligible in this paper.

Light paths ${s}_{12}$ and ${s}_{23}$ are the Euclidean distances in barycentric inertial vacuum space, where ${s}_{12}$ represents the path from the station at transmit time ${t}_{1}$ to the reflector at hit time ${t}_{2}$, and ${s}_{23}$ represents the path from the reflector at ${t}_{2}$ to the station at receive time ${t}_{3}$. Thus, we have the following decompositions:

where ${\overrightarrow{r}}_{ES}$, ${\overrightarrow{r}}_{BE}$, ${\overrightarrow{r}}_{BM}$, and ${\overrightarrow{r}}_{MA}$ are vectors in barycentric inertial space with B/E/M/S/A meaning barycenter, Earth center, moon center, Earth observing station, and lunar reflector array, respectively. The barycentric vectors ${\overrightarrow{r}}_{BE}$ and ${\overrightarrow{r}}_{BM}$ and the selenocentric vector ${\overrightarrow{r}}_{MA}$ are computed with position vectors and lunar libration angles provided by DE430 planetary ephemeris, using TDB. The geocentric vector ${\overrightarrow{r}}_{ES}$ is computed with station coordinates from ITRF2020 and IERS Earth Orientation Parameter data. Minor adjustments were made to the terrestrial and lunar coordinates to minimize the range residuals.

We selected 9747 LLR NP data records, which were retrieved from the EUROLAS Data Center ¹. The data were acquired by European LLR stations GRSM, MATM, and WETL, with the observation time spanning from May 2, 2012 to December 7, 2022.

The data were recorded in consolidated ranging data format. For each NP data record, the consolidated ranging data format includes the following details: the observing station, the observed reflector, the mean epoch of the data, the NP duration (or NP bin size), the number of photoelectrons accumulated within the NP bin, and the root mean square (RMS) value. The NP RMS value is equivalent to 1-sigma standard deviation of the range measurements in the NP.

We discarded 0.56% (55 of 9747) outlier NPs for which the range residuals were >0.4 m. In the remaining 99.44% (9692 of 9747) of the NP data, the overall range measurement residual against our LLR model was 35.9 mm in RMS. The RMS values were 36, 48, and 22 mm for stations GRSM, MATM, and WETL, respectively ². The number of NPs used was 8191, after applying the lunar declination limit, as mentioned in Section 3.1.

2.2   ERPs

The ERPs represent the rotation axis of Earth in the terrestrial reference frame and the rotation angle in the celestial reference frame. The three ERP components, ${y}_{p}$, ${x}_{p}$, and ${\psi }_{ERA}$, are related to frame rotations about the X, Y, and Z axes, respectively.

Earth’s rotation angle, denoted as ${\psi }_{ERA}$, is quantified within the celestial reference frame and represents the angle between Earth’s prime meridian and the equinox. It is conventionally represented by UT1. The formula for conversion from the Earth rotation angle (modulo $2\pi$) to UT1 (modulo 24h) is ${\psi }_{ERA}=15\times 1.002\,737\,811\,911\,354\,48\times \text{UT}1$, with units of seconds and arcseconds for UT1 and the Earth rotation angle, respectively.

The rotation axis is represented by celestial intermediate pole coordinates, or PM, denoted as ${x}_{p}$ and ${y}_{p}$. The celestial intermediate pole is the point of intersection between Earth’s axis of rotation and its ellipsoid, as measured within the ITRF, and PM denotes the angular deviation of the celestial intermediate pole from the north pole of the ITRF. The PM predictions used by modeling and determination are adopted from the IERS Conventions 2010 mean pole model^[20]. After epoch 2010.0, the model is formulated as follows:

where $t$ and ${t}_{0}$ are Besselian epochs of current time and 2000.0, respectively. The results are presented with units are milliarcseconds ³.

In this paper, the reference Earth Orientation Parameter data, including the PM components and the UT1–UTC series, were derived from the IERS C04 20 series using linear interpolation. The data files were retrieved from the IERS website^[16].

During 2012–2022, the PM prediction errors relative to the reference data were both within $\pm 200$ mas.

2.3   Determination of PM by LLR

To investigate the variation $\delta \rho$ of lunar range measurement $\rho$ with respect to variations in ERP components ${y}_{p}$, ${x}_{p}$, and ${\psi }_{ERA}$, or the partial derivatives, we consider the following geometrical relation:

where $H\left(\theta ,\phi \right)=\left(-cos\theta cos\phi ,-sin\theta cos\phi ,-sin\phi \right)$ is the unit vector along the line of sight from the reflector to the station, transposed; angles $\theta$ and $\phi$ are the celestial right ascension and declination of the reflector, respectively, as seen from the station; and the variation vector $\delta {\overrightarrow{r}}_{GCRS}$ is the variation of the station’s geocentric inertial vector, caused by ERP variation.

According to the IERS Conventions 2010^[20], the formula for transformation from an international terrestrial reference system (ITRS) to a geocentric celestial reference system (GCRS) is as follows:

where the bracketed term ${R}_{3}\left(-{\psi }_{ERA}\right){R}_{2}\left({x}_{p}\right){R}_{1}\left({y}_{p}\right)$ represents the effects of the ERPs. The term $Q\left(t\right)$ is the transformation matrix arising from the motion of the celestial pole in the celestial reference system, i.e., precession and nutation, while ${R}_{1}$, ${R}_{2}$, and ${R}_{3}$ are matrices for the coordinate frame rotation about axes X, Y, and Z, respectively. We adapted the equation to its current form using commutativity between the terrestrial intermediate origin locator ${R}_{3}\left(-s{'}\right)$ and the Earth rotation ${R}_{3}\left(-{\psi }_{ERA}\right)$.

The variations in astronomical angles $\theta$ and $\phi$ due to ERP changes are negligible. If we apply partial derivative rules on Equation (7), and let the ERP components be predicted values $\left({y}_{p}^{0},{x}_{p}^{0},{\psi }_{ERA}^{0}\right)$, then we have the following partial derivatives:

where ${M}_{HQR}=H\left(\theta ,\phi \right)Q\left(t\right){R}_{3}\left(-s{'}\right)$, and the matrices $R{{'}}_{k}\left(k=\mathrm{1,2},3\right)$ denote the derivatives of each rotation matrix with respect to the rotation angle.

We consider that the variation in ERPs represented by the vector $\overrightarrow{p}={\left(\delta {y}_{p},\delta {x}_{p},\delta {\psi }_{ERA}\right)}^{T}$ is the exclusive source of the variation in the laser ranging measurements, $\delta \rho$, i.e., $\delta \rho =\dfrac{\partial \rho }{\partial {y}_{p}}\delta {y}_{p}+\dfrac{\partial \rho }{\partial {x}_{p}}\delta {x}_{p}+\dfrac{\partial \rho }{\partial {\psi }_{ERA}}\delta {\psi }_{ERA}$, where $\delta \rho$ is the change in the laser ranging measurement, and the partial derivatives represent the sensitivity of the ranging measurement to changes in each ERP component. To determine the PM components, we assume UT1 is known; hence, $\delta {\psi }_{ERA}=0$. At the epoch of the common view measurement, with the group of stations A, B, and C, we can form a linear equation $M\overrightarrow{p}=\overrightarrow{q}$, where the following holds:

Matrix $M$ is called the design matrix, while $\overrightarrow{q}$ is the range residual vector, i.e., the difference between the measurements acquired and the measurements predicted by the theoretical model. It is an overdetermined problem, with three equations to solve for two unknowns. We apply the least squares technique to the problem, and the normal equation is formed as ${M}^{T}M\overrightarrow{p}={M}^{T}\overrightarrow{q}$, giving the least squares estimate $\overrightarrow{p}={\left({M}^{T}M\right)}^{-1}{M}^{T}\overrightarrow{q}$. If any perturbation $\delta \overrightarrow{q}$ were added to the measurement residual, then the solution would be perturbed correspondingly, i.e., $\delta \overrightarrow{p}={\left({M}^{T}M\right)}^{-1}{M}^{T}\delta \overrightarrow{q}$.

The PM predictions ${\left({\bar {y}}_{p},{\bar {x}}_{p}\right)}^{T}$ were calculated with a linear model as described in subsection 2.4. After solving for $\overrightarrow{p}$, we adjust the prediction with $\overrightarrow{p}$ to form estimates of ${x}_{p}$ and ${y}_{p}$ , i.e., $\left({\widehat{y}}_{p},{\widehat{x}}_{p}\right)=\left({\bar {y}}_{p}+\delta {y}_{p},{\bar {x}}_{p}+\delta {x}_{p}\right)$. The estimated ERP components were then compared with the reference data to form the solution error.

2.4   TCV Simulation

Regarding the simulation phase, we simulated TCV events (hereafter, events) for the GMW station group, according to the existing NP data at the European stations. First, we identified the available NPs at the European stations in the LLR data, taking each NP epoch as the TCV epoch. Second, ranging measurement data were generated at the NP epochs from the GMW stations to the observed reflector. Finally, the same number of PM solutions were generated from the simulated TCV events by solving the underlying equation $M\overrightarrow{p}=\overrightarrow{q}$.

To assess the robustness of the determination method, we introduced perturbations to the range residual vector $\overrightarrow{q}$, and then computed the perturbed solution $\overrightarrow{p}$. The perturbed solution was compared with the original solution to determine the effects of the perturbation.

Table 1 lists the settings of all the simulation cases. We started with the control case denoted as S0. Regarding the perturbation on UT1, an error of $\pm 100\mathrm{\mu }\mathrm{s}$ was introduced to the UT1 prediction to observe its effect, referred to as SU1 and SU2 in Table 1. Regarding the perturbation on the ranging measurement, we have six cases ranging from SG1 to SW2, for which uniform range errors of $\pm 50\mathrm{m}\mathrm{m}$ were introduced to each of the three stations. These perturbations were applied independently, although the linear nature of the model allows them to accumulate.

Table  1.  Settings of the simulation and augmentation cases

Phases/Stages	Case ID	Range Residuals			UT1 bias	Number of Solutions
		GRSM	MATM	WETL
Simulation	S0	0	0	0	0	8191
	SU1	0	0	0	+100 us	8191
	SU2	0	0	0	−100 us	8191
	SG1	+50 mm	0	0	0	8191
	SG2	−50 mm	0	0	0	8191
	SM1	0	+50 mm	0	0	8191
	SM2	0	−50 mm	0	0	8191
	SW1	0	0	+50 mm	0	8191
	SW2	0	0	−50 mm	0	8191
Augmentation Stage 1+2	AG1	real data	+ 48 mm	+ 22 mm	0	7646
	AG2	real data	+ 48 mm	− 22 mm	0	7646
	AG3	real data	− 48 mm	+ 22 mm	0	7646
	AG4	real data	− 48 mm	− 22 mm	0	7646
	AM1	+ 36 mm	real data	+ 22 mm	0	331
	AM2	+ 36 mm	real data	− 22 mm	0	331
	AM3	− 36 mm	real data	+ 22 mm	0	331
	AM4	− 36 mm	real data	− 22 mm	0	331
	AW1	+ 36 mm	+ 48 mm	real data	0	214
	AW2	+ 36 mm	− 48 mm	real data	0	214
	AW3	− 36 mm	+ 48 mm	real data	0	214
	AW4	− 36 mm	− 48 mm	real data	0	214
Augmentation Stage 2+1	AGM1	real data	real data	+ 22 mm	0	126
	AGM2	real data	real data	− 22 mm	0	126
	AGW1	real data	+ 48 mm	real data	0	55
	AGW2	real data	− 48 mm	real data	0	55

 | Show Table

DownLoad: CSV

2.5   TCV Augmentation

In the augmentation phase, we combined the existing NP data with the simulated data. The aim was to assess the accuracy of the solutions with common view events occurring in real-world observation epochs. In this phase, two augmentation strategies were employed in what were termed the '2+1' stage and the '1+2' stage; their settings are listed in Table 1.

In the '1+2' stage, we identified the available NPs at the European stations in the LLR data, which were regarded as monostatic events. Taking each NP epoch as the TCV epoch, we simulated the missing data with range residuals of $\pm$36, $\pm$48, and $\pm$22 mm for the stations GRSM, MATM, and WETL, respectively. The '1+2' stage yielded 12 cases, from AG1 to AW4, which are detailed in Table 1.

In the '2+1' stage, we first identified the bi-static common view events at the European stations in the LLR data. The search for bi-static common view events was based on the criterion that the temporal separation between two NPs at different stations should be <1 h. We then calculated the mid-point of the pair of NP epochs to be the TCV event epoch. Finally, we simulated the range residual at the missing station in the same way as in the '1+2' stage. The '2+1' stage yielded four cases, as listed in Table 1.

The outputs of both stages were evaluated using statistical methods to estimate how the TCV solution would perform under real-world conditions. If a TCV event occurs in reality, the solution should behave similarly to that of our results.

3.   RESULTS AND ANALYSIS

3.1   Preliminary Analysis

Prior to solution of the determination problem, two important attributes can be assessed: the sensitivity and the stiffness. The sensitivity is represented by the partial derivatives, and the stiffness is represented by the condition number.

The sensitivity of a station’s ranging measurement to variations in the ERPs can be quantified by the mean absolute values of the partial derivatives, which are measured in units of meters per arcsecond (m/arcsec). In our dataset, the mean absolute value of $\dfrac{\partial \rho }{\partial {x}_{p}}$ varies from 14 to 17 m/arcsec among the stations, larger than that of $\dfrac{\partial \rho }{\partial {y}_{p}}$, which varies from 8 to 10 m/arcsec. The WETL station exhibited the largest mean absolute values for both partial derivatives. The difference can be attributed to the geographic location of the stations, specifically their distance from the X, Y, and Z axes of the ITRF. The lower sensitivity of the range measurement to ${y}_{p}$ compared with those to ${x}_{p}$ is due to the stations being located closer to the X axis than to the Y axis. ⁴

Another key attribute is the stiffness of the least squares estimate problem ${M}^{T}M\overrightarrow{p}={M}^{T}\overrightarrow{q}$, which can be quantified by the condition number ${c}_{Normal}$ of the normal matrix ${M}^{T}M$, i.e., ${c}_{Normal}=\Vert {M}^{T}M\Vert \Vert {\left({M}^{T}M\right)}^{-1}\Vert$^[22]. Here, the matrix norm $\Vert \cdot \Vert$ is the 2-norm. The condition number ${c}_{Normal}$ is inversely related to ${\delta }_{moon}$, which is the lunar declination. When ${\delta }_{moon}$ approaches zero, the condition number ${c}_{Normal}$ soars to infinity. As an example, the relation between ${c}_{Normal}$ and ${\delta }_{moon}$ is shown in Fig. 2, where the condition number is plotted on a logarithmic scale. In applying the criterion $\left|{\delta }_{moon}\right| > {5}^{\circ }$, 8191 remained out of the original 9692 TCV events. The condition numbers of the remaining TCV events could be confined to an upper limit of $3.2\times {10}^{4}$. Among the remaining 8191 events, the GRSM station accounted for 93% (7646 events), the MATM station accounted for 4% (331 events), and the WETL station accounted for 3% (214 events).

Figure  2.  Relation between the condition number and the lunar declination.

DownLoad: Full-Size Img PowerPoint

3.2   Solutions Overview

A total of 106,845 TCV solutions were obtained for simulation and augmentation cases, as listed in Table 1. In each solution, a set of PM components ${x}_{p}$ and ${y}_{p}$ from a TCV event was evaluated.

In the simulation part, 8191 simulated TCV events were solved one time in the unperturbed control case S0, and 8 more times in the perturbed cases from SU1 to SW2.

In the augmentation part, two stages were investigated. In the '1+2' stage, 8191 TCV events were augmented in four combinations. Augmentation solutions were categorized by station and residual signs into 12 cases, from AG1 to AW4. In the '2+1' approach, there were 181 TCV events, each solved twice.

3.3   Solutions and Assessment of TCV Simulations

In the simulation phase, we simulated TCV events at 8191 epochs selected from the actual LLR data. Then, we compared the solved PM to the reference values to generate the solution error data. In the control case S0, the solution errors were trivially zero. The positively perturbed cases (SU1, SG1, SM1, and SW1) are depicted in Fig. 3. The negatively perturbed cases are symmetrical with the positive cases about the origin. It can be seen from the plots that each perturbation type displays a distinct directional pattern in the error distribution.

Figure  3.  Perturbation effects of UT1 (A) and range at GRSM (B), MATM (C), and WETL (D); the corresponding case IDs are SU1, SG1, SM1, and SW1, respectively.

DownLoad: Full-Size Img PowerPoint

The focus in the simulation phase was to examine the effects of perturbations in UT1 and range on the solution. The effects were quantified by using the mean absolute error (MAE), RMS error (RMSE), and absolute maximum (MaxAbs) metrics for the PM components ${x}_{p}$, ${y}_{p}$, and the pole position error $\sqrt{\delta {x}_{p}^{2}+\delta {y}_{p}^{2}}$. With the introduction of the UT1 error of $\pm$100 μs, the resulting ${x}_{p}$ and ${y}_{p}$ solution errors were only several milliarcseconds, as indicated by both the RMSE and the MAE. The maximum angular separation between the perturbed pole and the reference pole reached 18 mas, with an average of 6 mas. For comparison, range perturbations of $\pm$50 mm resulted in solution errors that were 4–10 times greater. Detailed statistics are provided in Table 2.

Table  2.  Variation in the solutions of the simulation cases

	Solution Variation (mas)
	$\left|\delta {x}_{p}\right|$				$\left|\delta {y}_{p}\right|$				$\sqrt{\delta {x}_{p}^{2}+\delta {y}_{p}^{2}}$
Perturbations	MAE	RMSE	MaxAbs		MAE	RMSE	MaxAbs		MAE	RMSE	MaxAbs
Control Case (S0)	(trivial zero)
UT1 ±100 μs (SU1, SU2)	3	4	15		5	6	18		6	7	18
GRSM ±50 mm (SG1, SG2)	15	22	99		32	45	164		37	50	180
MATM ±50 mm (SM1, SM2)	15	18	69		34	44	150		40	48	150
WETL ±50 mm (SW1, SW2)	17	23	114		19	24	99		27	34	133

 | Show Table

DownLoad: CSV

The results imply that the UT1 prediction error of $\pm$100 μs could, on average, introduce an error of several milliarcseconds into the PM solutions. Furthermore, the range error of $\pm$50 mm could result in RMSEs of around 20 mas in the ${x}_{p}$ solution and nearly 50 mas in the ${y}_{p}$ solution. The range error at the WETL station contributed less to the pole position error because of its higher latitude and therefore the greater distance from the X and Y axes of the ITRF, which is related to the fact that the partial derivatives $\dfrac{\partial \rho }{\partial {x}_{p}}$ and $\dfrac{\partial \rho }{\partial {y}_{p}}$ were larger at the WETL station. The perturbation results are useful for assessing solution errors under a certain amount of ranging error.

3.4   Solutions and Assessment of TCV Augmentations

The augmentation phase focuses on evaluating solution errors under real-world settings. Table 3 displays all the augmentation results by case ID, the settings of which are listed in Table 1.

Table  3.  Solution errors in the augmentation cases

Case ID	Solution Error (mas)
	$\left|\delta {x}_{p}\right|$				$\left|\delta {y}_{p}\right|$				$\sqrt{\delta {x}_{p}^{2}+\delta {y}_{p}^{2}}$
	MAE	RMSE	MaxAbs		MAE	RMSE	MaxAbs		MAE	RMSE	MaxAbs
AG1	15	22	206		37	55	730		42	59	752
AG2	23	30	256		43	59	713		52	66	734
AG3	21	28	184		39	53	554		48	59	570
AG4	14	21	210		34	49	537		38	54	552
AM1	16	24	164		21	45	563		29	51	586
AM2	31	42	210		34	58	596		50	72	620
AM3	30	40	186		36	52	339		51	66	353
AM4	16	24	135		22	39	372		30	46	387
AW1	7	10	39		15	17	36		17	19	53
AW2	12	14	33		37	42	95		40	44	96
AW3	12	14	50		37	41	94		40	44	95
AW4	7	11	43		15	17	50		17	20	62
1+2 Summary	18	26	256		37	53	730		44	59	752
AGM1	23	44	251		40	129	866		50	137	899
AGM2	22	35	198		32	104	748		44	109	760
AGW1	10	13	26		36	41	73		39	43	73
AGW2	10	14	40		18	23	52		23	27	53
2+1 Summary	19	34	251		33	99	866		42	105	899

 | Show Table

DownLoad: CSV

In the ‘1+2’ stage of the augmentation phase, we simulated two range residuals for each selected range residual data. This stage yielded 12 cases from AG1 to AW4. Among the cases, the maximum pole position error was 752 mas. When comparing the pole position data to the linear pole prediction, 88.5% of the events (29,000 out of 32,764) yielded more accurate results than the prediction. The ${y}_{p}$ component errors were markedly larger than those of the ${x}_{p}$ component owing to the geographical locations of the GRSM and MATM stations, i.e., the distance to the X-axis of the ITRF is less than that to the Y-axis. The WETL cases, labeled AW1–AW4, yielded the lowest errors, while the MATM cases ranked in the middle and the GRSM cases exhibited the highest errors. This is partly attributable to the low range of the residual values of the WETL station, and also to the limited size of the WETL dataset.

In the ‘2+1’ stage of the augmentation phase, we identified 181 bi-static common view events from the actual NP data without applying the lunar declination filter. In the selected bi-static events, the range residuals varied from −0.097 m to +0.141 m and the RMS was 0.034 m. Of all the bi-static events, 69.6% (126 of 181) were between GRSM and MATM and 30.4% (55 of 181) were between GRSM and WETL. This stage yielded four cases from AGM1 to AGW2, as listed in Table 1. The statistical values of the ${y}_{p}$ component are obviously larger than those of the ${x}_{p}$ component owing to the same geographical reasons. Comparison of the linear PM predictions with the two perturbed cases combined reveals that 91.2% of the solutions provided values that were more accurate than the predictions. Among all cases, the maximum pole position error was 899 mas, which was attributable primarily to the high value of the maximum ${y}_{p}$. The AGW1&2 cases yielded lower errors owing to the low range of the residual values of the WETL station.

In summarizing the two stages, it is evident that the MAE statistics of ${x}_{p}$, ${y}_{p}$, and the pole position, as well as the percentage of the solutions that outperformed the predictions, were comparable. This indicates that the two augmentation strategies did not introduce substantial bias in simulating reality.

Fig. 4 illustrates the solution errors for the two augmentation strategies, where Fig. 4A displays 99.35% (32,552 of 32,764) of the solutions in the ‘1+2’ augmentation, and Fig. 4B displays 95.6% (346 of 362) of the solutions in the ‘2+1’ augmentation. Note the difference in the scale of the axes.

Figure  4.  Solution errors in the augmentation phase.

DownLoad: Full-Size Img PowerPoint

4.   Discussion and Conclusions

This paper presents a method for determining Earth’s PM using tri-static range measurements in real time. We formulated the linear equations for this purpose, and applied the least squares method to solve them. The solutions were initiated using linear PM predictions. Preliminary analysis showed that the partial derivatives were associated with the geographical locations of the stations, and that the condition numbers were strongly correlated with the lunar declination. The solutions were compared with the IERS Earth Orientation Parameter 20 C04 time series to generate solution error data. Our method was examined through a two-phase study, including a simulation phase and an augmentation phase.

During the simulation phase, we simulated TCV events to the actual range data epochs. Predefined errors were introduced in the UT1 and range to observe the consequent variations in the solution. The results showed that the UT1 prediction error of $\pm$100 μs led to minor variations of several milliarcseconds in the solutions, and that the uniform range errors of $\pm$50 mm could lead to errors of 18–23 mas in the ${x}_{p}$ solution and 24–45 mas in the ${y}_{p}$ solution (RMSE). Perturbation patterns were revealed in the solution variations.

During the augmentation phase, two strategies were used to simulate realistic TCV scenarios. In the first strategy named ’1+2’, we simulated two additional measurements for each NP entry. Realistic range errors were introduced to yield solution errors of 26 and 53 mas in ${x}_{p}$ and ${y}_{p}$ (RMSE), respectively. In the second strategy named as ’2+1’, we identified the actual bi-static common view events and augmented them into tri-static events using simulated data, to yield solution errors of 34 and 99 mas in ${x}_{p}$ and ${y}_{p}$ (RMSE), respectively. Both strategies revealed that the solution error of the ${y}_{p}$ component was markedly larger than that of the ${x}_{p}$ component. The largest pole position error was 899 mas. The realistic scenarios within the '1+2' and '2+1' augmentation strategies yielded solutions that were more accurate than the linear pole prediction for 88.5% and 91.2% of the cases, respectively. Cases involving WETL data had notably smaller error statistics owing to the low range of the residuals. The results presented in this paper are compared with those of other reported methods in Table 4. Our method could yield errors that are 20–60 times larger than LLR long-term solutions, or 300–2000 times larger than VLBI solutions. However, the temporal resolution of our method is the LLR NP bin size, i.e., approximately 15 min, which is much shorter than that of the other methods considered.

Table  4.  Comparison of PM accuracy of the presented method with that of other existing methods

Parameter	LLR TCV (This paper)		LLR long-term Singh et al.[10]	VLBI Schuh & Behrend[2]
	Aug 1+2	Aug 2+1
${x}_{p}$ RMSE (mas)	26	34	1.30	0.050 to 0.080
${y}_{p}$ RMSE (mas)	53	99	1.63
Time Resolution	~ 15 minutes		1 day	1 day

 | Show Table

DownLoad: CSV

Overall, these results indicate that tri-static determination of PM coordinates is possible. Our method requires one set of data from three stations tracking the moon simultaneously, instead of a long period of data accumulation. This makes it a potentially viable tool for detecting transient changes in Earth’s rotation. It could also prove useful for validating ERPs to accuracy of several milliarcseconds.

However, our method cannot replace traditional long-term determination methods, partly because of its dependence on precise geodetic data and partly because it lacks the capability to estimate additional parameters.