Open AccessArticle

Regional GNSS Common Mode Error Correction to Refine the Global Reference Frame

Ruyuan Wang

^1,2,3

Junping Chen

^1,2,4,*

Danan Dong

¹,

Weijie Tan

¹ and

Xinhao Liao

^1,3

Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai 200030, China

School of Astronomy and Space Science, University of Chinese Academy of Sciences, Beijing 100049, China

School of Physical Science and Technology, ShanghaiTech University, Shanghai 200120, China

⁴

Shanghai Key Laboratory of Space Navigation and Positioning Techniques, Shanghai 200030, China

Author to whom correspondence should be addressed.

Remote Sens. 2024, 16(23), 4469; https://doi.org/10.3390/rs16234469

Submission received: 16 October 2024 / Revised: 25 November 2024 / Accepted: 26 November 2024 / Published: 28 November 2024

(This article belongs to the Special Issue Advances in GNSS for Time Series Analysis)

Download

Browse Figures

Graphical abstract
"> Figure 1
Distribution of 180 selected IGS stations. "> Figure 2
Number of solutions of each station. "> Figure 3
Number of effective stations at each epoch. "> Figure 4
Two examples for station position raw time series (green) and residual time series (blue). "> Figure 5
Region division for common mode error estimation. "> Figure 6
Mean and maximum distances between stations in each region. "> Figure 7
Spatial response of the first three principal components (Top: the spatial response for PC1. Middle: the spatial response for PC2. Bottom: the spatial response for PC3). "> Figure 8
First principal component time series for eastern Asia. "> Figure 9
Coordinate time series with CME correction. "> Figure 10
Ratio of the first three eigenvalues of the corresponding components for each region. "> Figure 11
Average spatial response to the first three principal components of each region. "> Figure 12
Two examples of site position residual time series before CME filtering (blue) and after CME filtering (red). "> Figure 13
RMS improvement percentage of residual time series before and after CME filtering. "> Figure 14
Velocity difference and uncertainty difference before and after CME filtering. "> Figure 15
Distribution of 365 reference stations in Scheme 2. "> Figure 16
Helmert transformation parameters between IGS and ITRF2020 solutions before (blue) and after (red) CME filtering (left: Scheme 1, right: Scheme 2). "> Figure 17
Uncertainty of transformation parameters before (blue) and after (red) CME filtering (left: Scheme 1, right: Scheme 2). "> Figure 18
RMS of post-transformation residual time series before (blue) and after (red) CME filtering (left: Scheme 1, right: Scheme 2). "> Figure 19
Region division for common mode error estimation (Division 2). "> Figure 20
RMS improvement percentage of residual time series before and after CME filtering (Division 2). ">

Versions Notes

Abstract

Common mode error (CME) arises from various sources, including unknown regional errors, potential geophysical signals, and other factors present in global navigation satellite system (GNSS) coordinate solutions, undeniably affecting the GNSS precision. This research concentrates on the effects of CME correction in global IGS-based reference frame refinement. We first estimated the regional CME with principal component analysis to obtain CME-corrected GNSS coordinate solutions. Subsequently, effects on the global reference frame with the regional CME correction were analyzed in three aspects: accuracy improvement of the coordinate solutions, variation in the velocity field, and accuracy improvement of the Helmert parameters in the reference frame transformation. The results show that after applying CME correction, the GNSS coordinate accuracy was improved by 28.9%, 22.1%, and 29.5% for the east, north, and vertical components, respectively. Regarding the site velocities, the maximum difference in velocity reached 0.48 mm/yr. In addition, the standard deviation of the Helmert transformation parameters between the International Terrestrial Reference Frame (ITRF) and the IGS-based reference frame—exclusively derived from GNSS technology—was reduced by over 30%, indicating CME correction enhanced the accuracy of the transformation parameters and refined the IGS-based reference frame.

Keywords:

common mode error correction; position coordinates time series; IGS solutions; terrestrial reference frame

Graphical Abstract

1. Introduction

An accurate and stable global terrestrial reference frame is critically important to geodynamics and engineering construction, as it provides the geometric and physical data essential for mapping and engineering [1]. Establishing a millimeter precision coordinate reference frame is an important challenge in the current international geodetic community [2]. Global reference frames are typically realized by the estimated position coordinates and velocities obtained by four space geodetic techniques [3]. The global navigation satellite system, which is one of the four space geodetic techniques, bridges the three other techniques by co-locating all Very Long Baseline Interferometry (VLBI) and Satellite Laser Ranging stations and more than half of the Doppler Orbitography and Radiopositioning Integrated by Satellite stations. IGS provides the coordinate time series of more than 400 reference stations around the world. Therefore, IGS greatly promotes the realization and application of the ITRF, and the accuracy of GNSS station solutions directly impacts the accuracy of ITRF. IGS changes its reference datum regularly and is aligned with the updates of ITRF. The latest version is IGS20, which is aligned with the ITRF2020 datum. Compared with the previous version IGS14, IGS20 adopts a number of new models and conventions, which were used in the third reanalysis campaign (Repro3) [4]. However, although various errors in GNSS processing have been studied by many scholars [5,6,7,8], unmodeled errors, especially spatially correlated errors, remain in the site coordinate time series, which may bias the IGS-based reference frame. One of the major regional spatially correlated errors is the common mode error (CME) in the GNSS time series [9], which affects the accuracy and dependability of GNSS solutions [10], and the error can lead to changes in station coordinate time series at the millimeter scale.

The CME is composed of unknown regional errors and potential geophysical signals [11,12,13]. It can be alleviated using a method known as regional filtering. Wdowinski et al. [9] initially utilized the regional stacking filtering method to define a certain similarity in single-day coordinate residuals as CME. Nikolaidis et al. [14] introduced the single-day solution variance as the weight in regional stacking filtering and proposed the weighted mean method. Stacking methods are sufficient for regional GPS (Global Positioning System) networks with intermediate-length baselines (up to several tens of kilometers), but they cannot necessarily reveal the spatial distribution of the CME for larger regions. To describe CME more objectively, Dong et al. [15] employed principal component analysis (PCA) to identify CMEs. PCA relaxes the restriction of a uniform distribution and enables the network data to uncover the spatial distribution of CME [16]. Therefore, PCA is more rigorous than stacking methods, in theory. Many improved filtering techniques have been put forward and used to estimate CME. To ensure that CME estimates are not affected by local nonsecular deformation, Liu et al. [17] conducted CME estimation over a broader area, revealing a decrease in raw time series dispersion. However, volcanic-related deformation signals remained after the CME was removed. Ming et al. [18] removed CME by PCA to extract the useful geophysical information accurately. Tian et al. [19] identified minor tectonic deformation signals by applying the correlation-weighted spatial filtering method, effectively extracting coherent signals from the comprehensive GPS network. Jiang et al. [20] performed a linear regression analysis of the noise amplitude vectors of 126 GPS stations in various components before and after CME was removed. The findings demonstrated that eliminating CME led to an improvement in the correlation between the horizontal and vertical flicker noise amplitude vectors, transitioning from uncorrelated to moderately correlated. Li et al. [21] found that the influence of CME on the estimation of vertical velocities should not be disregarded. Tan et al. [22] used PCA to identify CME in the Chuandian region of China and showed that approximately 15% of the CME can be attributed to daily surface mass loadings. Li et al. [23] proposed variational Bayesian PCA and selected 44 GNSS stations to estimate and extract the CME. After CME filtering, the interstation correlation coefficients were significantly reduced, and the root-mean-square values of the residual time series and the colored noise amplitudes were greatly suppressed. Zhou et al. [24] employed multichannel singular spectrum analysis and PCA to extract CME on a global scale and estimated the velocity from coordinates with and without CMEs.

CME is analyzed on the regional scale and exhibits variations among different regions. Given that CME exists in all regions of the Earth and primarily reflects the spatially correlated errors (including undesired/unmodeled signals) in GNSS solutions, it is natural to ask whether such errors may affect the accuracy of the reference frame determination. In this study, we first estimated the regional CME with principal component analysis. Then, CME correction was applied to GNSS coordinate solutions. Subsequently, the effects on reference frame determination were analyzed in terms of the coordinate noise, velocity field, and Helmert transformation parameters. The structure of the paper is as follows. Section 2 introduces the data source: we used weekly solutions of 180 IGS stations around the world for the period from 2000.0 to 2021.0 and applied PCA to extract regional related CMEs. Section 3 presents CME results of the regions covering the globe. Section 4 discusses the refinement of the global reference frame by applying CME correction. Section 5 further discusses the contribution of CME correction to the reference frame, using different strategies to define the regional CME. Finally, the conclusions are discussed and summarized in Section 6.

2. Data Source

We chose a 21-year span (2000.0–2021.0) of weekly coordinate solutions from Repro3 of 180 globally distributed IGS stations. Figure 1 depicts the distribution of the selected stations, and the number of weekly solutions is shown in Figure 2. Each station consists of at least 788 weekly solutions (more than 70% of the total number of weeks), namely, more than 14 years of data (Figure 2). The number of effective stations in each epoch is listed in Figure 3. In each epoch, the number of available station solutions is greater than 116 (64% of the total solutions). Thus, we believe that the 180 stations can be used as research objects.

We utilized the QOCA software (Version 1.28) of the Jet Propulsion Laboratory (JPL) to perform time-series analysis on all selected reference stations [15]. The model outlined in Appendix A.1 was applied, in which the offset and trend, as well as the annual and semiannual terms for the coordinate time series, were modeled. In addition to abnormal changes due to antenna changes and coseismic and postseismic effects provided by ITRF2020, several jumps caused by earthquakes or unknown causes were manually identified. Subsequently, we removed these components from the coordinate time series, resulting in the residual time series. Weekly coordinate solutions with formal uncertainties and residuals exceeding 15 mm (east/north), or 30 mm (vertical) were discarded. After preprocessing, the residual time series showed relative stability. Taking station SAMO (Fagalii, Samoa) and station WGTN (Wellington, New Zealand) as examples, the green points in Figure 4 are the original time series, and the blue points are the residual time series after preprocessing. The figure shows that the linear term, seasonal term and jump, logarithmic decay, and exponential decay terms for SAMO and WGTN were eliminated after preprocessing.

3. CME Extraction

After preprocessing the raw coordinate time series, which contained both signals and noise with temporal and spatial correlations, we obtained the residual time series. Dong et al. explored a spatiotemporal filtering approach that enables the network stations to have a nonuniform spatial response to a CME source, called principal component analysis (PCA). In this study, we used PCA embedded in the QOCA for CME estimation. The PCA decomposition method is presented in Appendix A.2.

Solutions with uncertainties over 20 mm for all three components were discarded. Residuals exceeded 150 mm (east/north) or 300 mm (vertical) were considered outliers and were discarded. Full and continuous time series are essential for PCA. Accordingly, for the time gaps shorter than 2 weeks, we filled them using three-point Lagrangian interpolation. For larger gaps, we used the spatial average derived from all valid stations in all epochs as the value for the missing epochs. Covariance matrices were constructed to carry out PCA as well as obtain eigenvectors and principal components. We took only the first three principal modes and used Equation (A4) to perform an iterative PCA process to estimate the missing data. Similar to Dong (2006), we set the convergence criterion to be that the average change in the missing point values compared with the previous iteration is less than

10^{- 6}

[15]. Subsequently, we carried out PCA on continuous residual time series. New covariance matrices

B

in Appendix A.2 were constructed, and then eigenvectors and principal components were generated. After arranging the eigenvectors in descending order, the leading principal components reflected the largest contributions to the variance of the regional network residual time series and were associated with the common mode component of the regional network. The other high-order principal components are typically associated with specific local or individual site characteristics. The principal components regarded as common mode errors were subtracted from the coordinate time series to obtain a new coordinate time series, which is the coordinate time series with CME correction.

Given that CME represents regional errors and PCA allows a nonuniform spatial response across network stations, we randomly divided the 180 reference stations into 22 regions to obtain the CME in each respective region (Figure 5). At least three stations were included in each region to ensure that the CME in each region could be estimated via PCA. The distances between the stations in each region are shown in Figure 6. On average, the distance between stations within each region is approximately 1700 km, which is commonly considered sufficient for maintaining common mode correlations [25,26,27]. The largest distance of 5358 km is consistent with the findings of Marquez Azua [28], which suggested CME may still exist when the distance is up to 6000 km.

We conducted PCA on the residual time series to extract the CME for each region. When the normalized response of the first principal component of a certain station was 1, and the responses of all other stations were less than 0.25, we considered the station to have a strong local signal, and therefore, the station should be deleted. Through iterative PCA, we ultimately identified 180 selected stations.

Eastern Asia (Figure 5, Region 2), where IRKJ, ULAB, BJFS, DAEJ, KGNI, and AIRA are located, was taken as an example. The contributions of the first three principal components (PCs) were 44.8%, 16.0%, and 14.0% for the east component; 39.7%, 21.0%, and 15.2% for the north component; and 33.8%, 20.5%, and 15.2% for the vertical component. The contribution of the first principal component (PC1) was much greater than that of the second and third, which means that PC1 accounts for the most energy in the data.

In addition, the spatial response of the first three principal components is shown in Figure 7. For PC1, the east component had mean and minimum spatial responses of 0.76 and 0.54, respectively. The north component’s responses were 0.73 and 0.53, while the vertical component’s responses were 0.65 and 0.43. All stations showed spatial responses above 0.25 for the three components. However, both positive and negative terms were present in the spatial responses of PC2 (the second principal component) and PC3 (the third principal component), with minimum responses of −0.96 and −0.88, respectively. Consequently, we treated PC1 as the CME in Region 2.

Figure 8 shows the time series of PC1 in the region. The root-mean-square (RMS) of PC1 was 0.82 mm, 1.02 mm, and 2.88 mm for the three components. Compared with that for the east and north components, the fluctuation of PC1 for the vertical component was significant, which reveals that the influence of CME on the GNSS solutions for the vertical component is more obvious than that for the horizontal components. By considering PC1 as the CME, and combined with the corresponding spatial response, the coordinate time series with CME correction were obtained and are shown in Figure 9.

For all the stations in the defined regions around the globe, the ratio of the first three eigenvalues of the corresponding components to the summation of all eigenvalues for the E, N, and U coordinate components is shown in Figure 10. The average ratio of PC1 was 46.9%, 48.1%, and 43.5% for the east, north, and vertical components, respectively. PC2 accounted for 20.9%, 20.0%, and 22.0%, respectively. PC3 accounted for 14.7%, 14.3%, and 15.4%, respectively. As shown in the figure, the ratio of the eigenvalue of PC1, which exceeded 40% for all components, was much higher than those of the other principal components.

The average spatial responses to the first three principal components of each region are shown in Figure 11. For all regions, the mean spatial responses corresponding to PC1 had the same sign, and the average responses were 0.70, 0.71, and 0.66 for the east, north, and vertical components, respectively. Clearly, the normalized spatial response of PC1 had a strong consistency, which reflects the common law of the entire region. The spatial responses of the other PCs included positive and negative terms, and the average responses in most regions were less than 0.25; in other words, the standardized spatial responses varied drastically by station and showed a poor consistency, indicating nonuniform spatial responses. The spatial response of PC1 in all regions met the CME definition from Dong et al. [15] that more than 50% of stations have a significant normalized response (>0.25). In the case of the above analysis, we considered PC1 as the CME as applied in Dong et al. [15] to perform regional filtering and removed it from the raw coordinate time series in this study.

4. CME Impact on Coordinate Time Series and Contribution to Reference Frame Refinement

According to the discussion in the preceding section, the coordinate time series with CME correction was derived. We define this approach as CME filtering in later parts of this paper. In the subsequent section, we analyze the impact of CME filtering on the reference frame, with a focus on the following three aspects: (1) accuracy improvement of the coordinate solutions, (2) station velocity differences before and after CME correction, and (3) accuracy improvement of the IGS-based global reference frame.

4.1. Accuracy Improvement of Coordinate Solutions

We further applied CME correction to the residual coordinate time series and derived the new residual time series. Figure 12 shows the time series for the SAMO and WGTN stations. The dispersion of the residual time series after CME filtering was much smaller. For SAMO, the RMS of the residual time series before CME filtering was 1.30 mm, 1.09 mm, and 3.94 mm for the three components, respectively. The values were 0.93 mm, 0.91 mm, and 3.64 mm, respectively, after CME filtering. The RMS decreased by 28.7%, 16.9%, and 7.77% for the three components. For WGTN, the RMS of the residuals before filtering was 1.22 mm, 0.88 mm, and 2.57 mm, respectively, and the values were 0.87 mm, 0.69 mm, and 1.81 mm, respectively, after CME filtering. The RMS decreased by 28.9%, 22.1%, and 29.5% for the three components. The station position precision clearly improved after CME filtering.

The RMS of the residual time series in all epochs before and after CME filtering of the 180 IGS reference stations was calculated. The RMS improvement percentages of the residual position time series for each station are shown in Figure 13. A clear improvement in the RMS among all stations was observed, and the improvement showed no obvious regional distribution law. However, in regions with a large distance between stations, such as Regions 3, 14, and 18, the accuracy improvement of most stations was not significant. By applying CME correction, the mean RMS of the residual time series was reduced by 26.4%, 26.7%, and 23.9% for the east, north, and vertical components, respectively. The value is slightly lower than Li’s result [21], which may be due to the proximity of the stations in Southern California in their study, and is roughly equivalent to the results from Li [23].

4.2. Analysis of Velocity Estimation

We estimated the velocity of the stations along with its uncertainty before and after CME filtering using the model in Appendix A.1. A colored noise model was used to generate the uncertainty in the velocity solutions [29]. Figure 14 displays the velocity difference (

∆ V = V_{unfiltered} - V_{filtered}

) and the corresponding uncertainty difference (

∆ σ = σ_{unfiltered} - σ_{filtered}

) [21]. Table 1 presents their statistics. The mean velocity differences were −0.00 mm/yr, 0.01 mm/yr, and 0.00 mm/yr for the east, north, and vertical components, respectively. Although the mean velocity differences were within 0.01 mm/yr, the maxima reached 0.17 mm/yr, 0.48 mm/yr, and 0.45 mm/yr for the three components, which is much larger than the up-to-date velocity accuracy. The effect of CME on velocity estimation cannot be overlooked. After CME was removed, the velocity uncertainties of almost all stations decreased for all components. The mean reduction in the velocity uncertainty was 0.01 mm/yr, 0.01 mm/yr, and 0.03 mm/yr for the east, north, and vertical components, respectively, and the maximum difference reached 0.09 mm/yr, 0.05 mm/yr, and 0.13 mm/yr, respectively. The effect of CME on velocity estimation therefore cannot be ignored.

4.3. Accuracy Improvement of the Reference Frame

In order to assess the accuracy of the reference frame with CME correction, we employed Helmert transformation parameters, which account for quasi-instantaneous translational, rotational, and scale discrepancies between two geometrically similar reference frames, and its estimates may be influenced by the network distributions [30]. When aligning IGS coordinate time series, which contain seasonal signals, with ITRF, significant aliasing errors can arise in the transformation parameters. Research conducted by Chen et al. suggested that the augmented ITRF2020 has a positive effect on frame alignment [31]. However, the focus of this paper is to analyze the impact of CME filtering on the reference frame. Therefore, in this study, we aligned the IGS solutions to the nonlinear ITRF, incorporating postseismic deformation models, and provided seasonal signals in the Center of Figure (CF) frame without considering other methods for augmenting ITRF2020.

Our approach involved treating ITRF2020 as the ground-truth reference frame and deriving corresponding weekly solutions based on the IGS weekly solutions. Subsequently, we obtained the Helmert transformation parameters between ITRF2020 and the IGS reprocessing solutions with and without CME correction, consisting of three translation parameters, three rotation parameters, and scale parameters. The calculation methodology for Helmert transformation parameters is detailed in Appendix A.3. We set outlier rejection thresholds at 10 mm (east/north) and 30 mm (vertical) during the Helmert transformation process [30,32]. To validate the efficacy of the CME filtering, we devised two distinct schemes. The 180 stations participating in estimating CME were used to calculate transformation parameters in Scheme 1, which we call 180 core stations. Other 185 stations that did not participate in estimating CME were added to estimate transformation parameters in Scheme 2, which we call 185 user stations. The distribution of the total 365 reference stations in Scheme 2 is shown in Figure 15. Red spots represent newly added user stations. We calculated the CME of 185 user stations distributed in 22 regions using the average corresponding of each region and obtained the CME-filtered solutions. Results showed that outlier rejection rates are 0.63 and 0.60 with and without CME filtering in Scheme 1 and 1.57 and 1.54 with and without CME filtering in Scheme 2 during the Helmert transformation process.

4.3.1. Helmert Transformation Parameters

The time series of translation (Tx, Ty, Tz), rotation (Rx, Ry, Rz), and scale (Scale) between the ITRF2020 and IGS solutions in Scheme 1 are presented in Figure 16 (left), and those of Scheme 2 are presented in Figure 16 (right). Visible biases and rates can be observed in the seven parameters. Due to differences in the station selection and processing strategies, the values of the seven parameters are inconsistent with the results of Altamimi [33]. Notably, at epoch 2015.0, the fitting scale parameter is −0.68 ppb with a rate of −0.016 ppb/yr, which is consistent with the scale parameter from IGS to ITRF2020 [33]. The near-annual signals present in the rotation parameters and scale parameters could be the aliasing errors in different techniques, periodicity errors in IGS solutions, or unapplied geophysical signals in ITRF [31,33,34].

The mean values and the standard deviation (STD) of the transformation parameters after subtracting the linear fitting terms, both before and after CME filtering, are presented in Table 2. The average of three translation parameters are within 0.05 mm, indicating that the systematic deviation between ITRF2020 and IGS Repro3 was effectively eliminated by the linear fitting. There were noticeable decreases in both the mean and STD of the parameters after CME filtering. The STDs of the seven transformation parameters were reduced by approximately 65%, 58%, 53%, 42%, 33%, 56%, and 46% in Scheme 1 and by 60%, 54%, 51%, 47%, 39%, 62%, and 47% in Scheme 2. The reduction in the STD demonstrates that the transformation parameters exhibited reduced fluctuations and higher stability after CME filtering.

4.3.2. Uncertainty of the Helmert Transformation Parameters

The uncertainty of each parameter of Helmert transformation was also evaluated at each epoch, and the time series are shown in Figure 17. In general, the uncertainty of Scheme 2 was found to be lower than that in Scheme 1, which could be attributed to the increased number of reference stations used for calculating the transformation parameters. Moreover, the values after CME filtering were smaller. The accuracy of the translation parameters improved by approximately 0.04 mm, 0.04 mm, and 0.03 mm in Scheme 1 and 0.02 mm, 0.03 mm, and 0.03 mm in Scheme 2. The rotation parameters were improved by approximately 1.28

μ a s

, 1.22

μ a s

, and 1.32

μ a s

in Scheme 1 and 0.80

μ a s

, 0.79

μ a s

, and 0.77

μ a s

in Scheme 2. The scale parameters were improved by 0.01 ppb and 0.00 ppb in Scheme 1 and Scheme 2, respectively. The accuracy of the Helmert transformation parameters was increased due to the improved accuracy of the IGS coordinate solution after CME filtering. Consequently, the credibility of the transformation parameters between the two frames was enhanced with CME filtering.

4.3.3. Post-Transformation Residuals Time Series

The coordinate difference after Helmert transformation, known as post-transformation residuals, is also an accuracy indicator for evaluating reference frames. We calculated the post-transformation residuals between ITRF2020 and the unfiltered/filtered IGS solutions. The results are depicted in Figure 18, with the blue and red lines representing the post-transformation residual series before and after CME filtering, respectively. It is evident that the post-transformation residuals all decreased after filtering, with the vertical component exhibiting the most notable improvement. It decreased by 0.23 mm, 0.19 mm, and 0.58 mm in Scheme 1 and by 0.21 mm, 0.15 mm, and 0.44 mm in Scheme 2, respectively, in the three directions. The lower post-transformation residuals further affirm the refinement of the reference frame achieved through the utilization of the CME-filtered coordinate time series.

5. Discussion

To explore whether the different divisions of the globe will affect the CME filtering results, we defined a second region division (Division 2) (Figure 19). Most regions are different from those in the first division definition given in Section 3 (Division 1). Similarly, at least three stations were included in each region to ensure that the CME in each region could be estimated. For each region, the average distance between the stations was approximately 1770 km, and the largest distance between two stations reached 5358 km.

Applying the same analysis as in Section 4.1, PCA was performed on the residual time series of all stations in each region, and the residual time series after filtering were derived. The RMS of the residual time series before and after filtering was also estimated, and its improvement percentages for each station are shown in Figure 20. By applying CME correction, the mean RMS of the residual time series was reduced by 24.3%, 23.4%, and 18.7% for the east, north, and vertical components, respectively. Compared with the results in Figure 13, the accuracy improvement of a few stations differed greatly, but the accuracy improvement of most stations was not substantially different from that under Division 1. Similar to Division 1, the accuracy improvement of the stations in a region separated by long distances was not significant.

Table 3 lists the statistics of

∆ V

and

∆ σ

. The mean

∆ V

values were −0.00 mm/yr, 0.01 mm/yr, and 0.01 mm/yr, and the maximum

∆ V

reached 0.28 mm/yr, 0.41 mm/yr, and 0.42 mm/yr for the east, north, and vertical components, respectively. The maxima in

∆ σ

were 0.08 mm/yr, 0.05 mm/yr, and 0.10 mm/yr for the three components, respectively. In general, there is no almost difference in the mean

∆ V

and mean

∆ σ

between Division 1 (shown in Table 1) and Division 2. The maximum

∆ V

difference is 0.42 mm/yr, and the maximum

∆ σ

difference is 0.10 mm/yr. The influence of different divisions of the globe on velocity can be ignored.

The Helmert transformation parameters between the ITRF2020 and IGS solutions with/without CME filtering under the two divisions were compared. A total of 180 core stations were used to calculate the Helmert transformation parameters. The mean and STD of the Helmert parameters were recalculated after deducting the linear terms under Division 2. The results are displayed in Table 4. In the case of filtering by the two divisions, the mean and STD of the Helmert parameters showed only slight differences. The mean of the translation parameter varied within 0.01 mm, and the rotation parameter varied within 0.001 mas. Compared with the accuracy difference of the Helmert parameters before and after CME filtering, these accuracy differences were much smaller. Therefore, applying regional CME filtering can improve the accuracy of the reference frame. In addition, if the globe is divided into a considerable number of regions, and the divisions are conducted based on randomness without considering geophysical factors like tectonic plates, the difference in refinement of the reference frame between different divisions can be considered negligible.

6. Conclusions

In this study, we utilized Repro3 weekly coordinate solutions spanning a period of 21 years (2000.0–2021.0). The GNSS regional-related CMEs were estimated by applying principal component analysis, and the filtered GNSS coordinate solutions were derived. To assess the extent of the refinement on the global reference frame with regional CME correction, the following three aspects were considered: improvement of the coordinate solutions accuracy, variation in the velocity field, and improvement of the reference frame accuracy before and after CME filtering. Based on the analysis, we reached the following conclusions:

Through the implementation of CME correction, the mean RMS of the residual time series was reduced by 28.9%, 22.1%, and 29.5% for the east, north, and vertical components, respectively. So, more accurate station coordinates can be obtained by CME filtering.
The maximum differences in velocity between unfiltered and filtered solutions were found to be 0.48 mm/yr, significantly exceeding the currently available velocity accuracy. After CME correction, the velocity uncertainties of almost all stations decreased for all components. The maximum difference reached 0.13 mm/yr. Thus, it is evident that CME for velocity estimation cannot be disregarded and should be taken into consideration to ensure the highest attainable accuracy.
The reduction in standard deviation values of the transformation parameters between ITRF2020 and the unfiltered/filtered IGS solutions exceeding 30% demonstrates that the transformation parameters exhibit fewer fluctuations and greater stability after CME filtering.
The average uncertainty of the Helmert transformation parameters significantly decreased after CME filtering compared with before. The higher precision of the Helmert parameters indicates that filtered IGS coordinate solutions are more reliable.
After CME filtering, the post-transformation residuals, which represent the coordinate difference between the ITRF and the unfiltered/filtered IGS solutions after Helmert transformation, consistently decreased for all three components. This indicates that CME correction enhances the accuracy of the transformation between the IGS reference frame and ITRF, leading to improved alignment between the two frames.
In comparison with the observed accuracy difference in the Helmert parameters before and after CME filtering, the accuracy difference between the Helmert parameters obtained using different divisions was considerably smaller. This finding suggests that when the globe is divided into a considerable number of regions, and the divisions are conducted based on randomness without considering geophysical factors like tectonic plates, the filtered coordinate time series has little difference in refining the entire reference system.

In conclusion, CME correction proves to be an effective method for obtaining more accurate reference station position coordinates and has effects on velocity estimation. It further enhances the refinement of the IGS reference frame. We believe that the application of common mode error correction can be extended to other technologies like VLBI. If we form per-technique long-term solutions by CME correction to contribute to the ITRF combination, overall accuracy can be further improved. Furthermore, our future research is designed to achieve a deeper insight into the physical origins of the CME and develop more accurate correction methods to effectively refine the reference frame.

Author Contributions

J.C. came up with the concept of the study. R.W. conducted the experiments and drafted the manuscript. D.D. and X.L. offered feedback. J.C. and W.T. corrected the manuscript to be more appropriate. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by National Natural Science Foundation of China (Grant No. 42474034), Natural Science Foundation of Shanghai (Grant No. 22ZR1472800).

Data Availability Statement

The IGS solutions from IGN are available via FTP (ftp://igs.ign.fr/pub/igs/products/repro3/, accessed on 1 October 2024), and the ITRF solutions are available via FTP (ftp://itrf-ftp.ign.fr/pub/itrf/itrf2020, accessed on 1 October 2024).

Acknowledgments

We are grateful to the French National Institute of Geographic and Forest Information (IGN) for providing the GNSS solutions and to the International Earth Rotation and Reference Systems Service for providing the ITRF solutions. We thank the Jet Propulsion Laboratory (JPL) for providing the QOCA software. Ideas and technical support from Xu, X. and Liu, J. significantly improved the quality of the manuscript.

Conflicts of Interest

The authors declare no conflicts of interest.

Appendix A

Appendix A.1. Fitting Model of Coordinate Time Series

The coordinate time series contain both signals and noise in the temporal and spatial domains. For each coordinate time series, the model proposed by Nikolaidis (2002) [14], which has the following form, was used:

x (t) = x_{0} + v \cdot (t - t_{0}) + \sum_{i} [S_{i} \sin (ω_{i} t) + C_{i} \cos (w_{i} t)] + \sum_{k} D_{k} B (t_{k}^{1}, t_{k}^{2}) + \sum_{j} F_{j} (t)

(A1)

where

x (t)

is the site coordinate in epoch

t

;

x_{0}

is the offset;

v

is the velocity;

S_{i}

and

C_{i}

are the amplitudes of the annual and semiannual seasonal changes, respectively;

D_{k}

is the amplitude of the box function;

B

is a box function, which has two time dummy parameters, the start time and the end time; and

F_{j}

is the amplitude of the remaining nonlinear terms.

Appendix A.2. Mathematical Model of Principal Component Analysis

Define the

m \times n

real-valued matrix

X (t_{i}, x_{j})

(

i = 1, 2, \dots, m

and

j = 1, 2, \dots, n

, assuming

m \geq n

) as the daily station coordinate time series of a regional network with n stations and spanning m days. In the

X

matrix, each column represents the residual value for the east, north, and vertical components from a single station in the network, and the rows represent the geodetic component values for all stations in a given epoch. The elements of the covariance matrix

B

of the matrix

X (t_{i}, x_{j})

are defined as

b_{i j} = \frac{1}{m - 1} \sum_{k = 1}^{m} X (t_{k}, x_{i}) X (t_{k}, x_{j})

(A2)

The (

n \times n

) symmetric matrix

B

can be decomposed as

B = V Λ V^{T}

(A3)

where the

n \times n

eigenvector matrix

V^{T}

is a row orthonormal matrix, and the matrix

Λ

has

k

nonzero diagonal eigenvalues

\{λ_{k}\} (n \geq k)

. In real data, the rank of matrix

B

is usually full

(k = n)

, and

B

can be expanded by

n

orthonormal vectors. Therefore, matrix

X (t_{i}, x_{j})

can be expanded by the orthogonal function basis

V

X (t_{i}, x_{j}) = \sum_{k = 1}^{n} a_{k} (t_{i}) v_{k} (x_{j})

(A4)

where

a_{k} (t)

is derived as

a_{k} (t_{i}) = \sum_{j = 1}^{n} X (t_{i}, x_{j}) v_{k} (x_{j})

(A5)

a_{k} (t)

is called the kth principal component (PC) of matrix

X

, and

v_{k} (x)

is its corresponding eigenvector. The principal components represent temporal variations, and the eigenvectors represent the corresponding spatial responses to each principal component.

Appendix A.3. Calculation Formula of Helmert Transformation Parameters

Helmert transformation parameters are calculated with the formula given by [33]:

{(\begin{matrix} x \\ y \\ z \end{matrix})}_{I G S} = {(\begin{matrix} x \\ y \\ z \end{matrix})}_{I T R F 2020} + T + D {(\begin{matrix} x \\ y \\ z \end{matrix})}_{I T R F 2020} + R {(\begin{matrix} x \\ y \\ z \end{matrix})}_{I T R F 2020}

(A6)

where

T

is the translation vector

T = {(T_{x}, T_{y}, T_{z})}^{T}

D

is the scale factor, and

R

is the matrix encompassing the rotation angles, given by

R = (\begin{matrix} \begin{matrix} 0 {- R}_{z} R_{y} \\ R_{z} 0 {- R}_{x} \end{matrix} \\ {- R}_{y} R_{x} 0 \end{matrix})

(A7)

References

Chen, J. Progress in theory and practice for geodetic reference coordinate frame. J. Geod. Geodyn. 2007, 27, 1–6. [Google Scholar]
Dang, Y.M.; Jiang, T.; Yang, Y.X.; Sun, H.; Jiang, W.; Zhu, J.; Xue, S.; Zhang, X.; Yu, B.; Luo, Z.; et al. Research progress of geodesy in China (2019–2023). Acta Geod. Cartogr. Sin. 2023, 52, 1419–1436. [Google Scholar] [CrossRef]
Han, C.H.; Li, L.; Cai, Z.W.; Lin, Y.T. The space-time references of BeiDou Navigation Satellite System. Satell. Navig. 2021, 2, 18. [Google Scholar] [CrossRef]
Rebischung, P.; IGS Reference Frame Working Group Coordinator. IGS Contribution to ITRF2020. 2022. Available online: https://itrf.ign.fr/docs/solutions/itrf2020/IGS-contribution-to-ITRF2020.pdf (accessed on 10 March 2024).
Su, K.; Jiao, G.Q. Estimation of BDS pseudorange biases with high temporal resolution: Feasibility, affecting factors, and necessity. Satell. Navig. 2023, 4, 17. [Google Scholar] [CrossRef]
Ding, J.S. Research on GNSS tropospheric delay modeling and spatial-temporal characteristics analysis of bias. Acta Geod. Cartogr. Sin. 2024, 53, 1659. [Google Scholar] [CrossRef]
Jiang, G.W.; Wang, P.L.; Guo, C.X.; Wang, B.; Yang, Y.X. Short-term GNSS network solution and performance in large height difference region with tropospheric delay constraint. Acta Geod. Cartogr. Sin. 2022, 51, 2255–2264. [Google Scholar] [CrossRef]
Song, Z.Y.; Chen, J.P.; Zhang, Y.Z.; Tang, W.J. Methods of URA evaluation refinement and re-broadcast based on BDS-3 B2b-PPP service. J. Navig. Position. 2024, 12, 58–67. [Google Scholar] [CrossRef]
Wdowinski, S.; Bock, Y.; Zhang, J.; Fang, P.; Genrich, J. Southern California Permanent GPS Geodetic Array: Spatial Filtering of Daily Positions for Estimating Coseismic and Postseismic Displacements Induced by the 1992 Landers Earthquake. J. Geophys. Res. Solid Earth 1997, 102, 18057–18070. [Google Scholar] [CrossRef]
Tian, Y. Anomalous high frequency seasonal noises in GPS positions time series. Sci. Surv. Mapp. 2011, 36, 26–28. [Google Scholar]
Dong, D.; Fang, P.; Bock, Y.; Cheng, M.K.; Miyazaki, S.I. Anatomy of apparent seasonal variations from GPS-derived site position time series. J. Geophys. Res. Solod Earth 2002, 107, ETG 9-1–ETG 9-16. [Google Scholar] [CrossRef]
Gruszczynski, M.; Klos, A.; Bogusz, J. Orthogonal Transformation in Extracting of Common Mode Error from Continuous GPS Networks. Acta Geodyn. Geomater. 2016, 13, 291–298. [Google Scholar] [CrossRef]
Dong, D.N.; Chen, J.Y.; Wang, J.X. The GNSS High Precision Positioning Principle; Science Press: Beijing, China, 2018; pp. 254–265. (In Chinese) [Google Scholar]
Nikolaidis, R. Observation of Geodetic and Seismic Deformation with the Global Positioning System. Ph.D. Thesis, University of California, San Diego, CA, USA, 2002; p. 249. [Google Scholar]
Dong, D.; Fang, P.; Bock, Y.; Webb, F.; Prawirodirdjo, L.; Kedar, S.; Jamason, P. Spatiotemporal Filtering Using Principal Component Analysis and Karhunen-Loeve expansion approaches for regional GPS network analysis. J. Geophys. Res. 2006, 111, B03405. [Google Scholar] [CrossRef]
Williams, S.D.P.; Bock, Y.; Fang, P.; Jamason, P.; Nikolaidis, R.M.; Prawirodirdjo, L.; Miller, M.; Johnson, D.J. Error analysis of continuous GPS position time series. J. Geophys. Res. 2004, 109, B03412. [Google Scholar] [CrossRef]
Liu, Z.; Dong, D.; Lundgren, P. Constraints on time-dependent volcanic source models at Long Valley Caldera from 1996 to 2009 using InSAR and geodetic measurements. Geophys. J. Int. 2011, 187, 1283–1300. [Google Scholar] [CrossRef]
Ming, F.; Yang, Y.X.; Zeng, A.M.; Jing, Y.F. Analysis of seasonal signals and long-term trends in the height time series of IGS sites in China. Sci. China Earth Sci. 2016, 59, 1283–1291. [Google Scholar] [CrossRef]
Tian, Y.; Shen, Z. Extracting the regional common-mode component of GPS station position time series from dense continuous network. J. Geophys. Res. 2016, 121, 1080–1096. [Google Scholar] [CrossRef]
Jiang, W.; Ma, J.; Li, Z.; Zhou, X.; Zhou, B. Effect of removing the common mode errors on linear regression analysis of noise amplitudes in position time series of a regional GPS network & a case of GPS stations in Southern California. Adv. Space Res. 2018, 61, 2521–2530. [Google Scholar] [CrossRef]
Li, Z.; Yue, J.; Li, W.; Lu, D.; Hu, J. Comprehensive Analysis of the Effects of Common Mode Error on the Position Time Series of a Regional GPS Network. Pure Appl. Geophys. 2019, 176, 2565–2579. [Google Scholar] [CrossRef]
Tan, W.; Chen, J.; Dong, D.; Qu, W.; Xu, X. Analysis of the Potential Contributors to Common Mode Error in Chuandian Region of China. Remote Sens. 2020, 12, 751. [Google Scholar] [CrossRef]
Li, W.; Jiang, W.; Li Zhao Chen, H.; Chen, Q.; Wang, J.; Zhu, G. Extracting Common Mode Errors of Regional GNSS Position Time Series in the Presence of Missing Data by Variational Bayesian Principal Component Analysis. Sensors 2020, 20, 2298. [Google Scholar] [CrossRef]
Zhou, M.; Guo, J.; Liu, X.; Shen, Y.; Zhao, C. Crustal movement derived by GNSS technique considering common mode error with MSSA. Adv. Space Res. 2020, 66, 1819–1828. [Google Scholar] [CrossRef]
Yuan, L.; Ding, X.; Chen, W.; Simon, K.; Chan, S.; Hung, P.; Chau, K. Characteristics of Daily Position Time Series from the Hong Kong GPS Fiducial Network. Chin. J. Geophys. 2008, 51, 1372–1384. [Google Scholar]
Tian, Y.; Shen, Z. Progress on Reduction of Non-tectonic Noise in GPS Position Time Series. Acta Seismol. Sin. 2009, 31, 68–81. [Google Scholar]
Tian, Y.; Shen, Z.; Li, P. Analysis on Correlated Nosie in Continuous GPS Observations. Acta Seismol. Sin. 2010, 32, 696–704. [Google Scholar] [CrossRef]
Marquez-Azua, B.; DeMets, C. Crustal velocity field of Mexico from continuous GPS measurements, 1993 to June 2001: Implications for the neotectonics of Mexico. J. Geophys. Res. 2003, 108, 2450. [Google Scholar] [CrossRef]
Williams, S.D.P. The Effect of Coloured Noise on the Uncertainties of Rates Estimated from Geodetic Time Series. J. Geod. 2003, 76, 483–494. [Google Scholar] [CrossRef]
Liu, J.; Chen, J.P.; Liu, P.Z.; Tan, W.J.; Dong, D.N.; Qu, W.J. Comparison and Assessment of Three ITRS Realizations. Remote Sens. 2021, 13, 2304. [Google Scholar] [CrossRef]
Chen, G.; Wei, N.; Wang, H.; Li, T.; Li, M.; Zhao, Q.L. Impacts of non-linear ITRF2020 on reference frame alignment. Adv. Space Res. 2023, 74, 2569–2582. [Google Scholar] [CrossRef]
Fritsche, M.; Sośnica, K.; Rodríguez-Solano, C.J.; Steigenberger, P.; Wang, K.; Dietrich, R.; Rothacher, M. Homogeneous reprocessing of GPS, GLONASS and SLR observations. J. Geod. 2014, 88, 625–642. [Google Scholar] [CrossRef]
Altamimi, Z.; Rebischung, P.; Collilieux, X.; Métivier, L. Chanard KITRF2020: An augmented reference frame reefing the modeling of nonlinear station motions. J. Geod. 2023, 97, 47. [Google Scholar] [CrossRef]
Rebischung, P.; Altamimi, Z.; Ray, J.; Garayt, B. The IGS contribution to ITRF2014. J. Geod. 2016, 90, 611–630. [Google Scholar] [CrossRef]

Figure 1. Distribution of 180 selected IGS stations.

Figure 2. Number of solutions of each station.

Figure 3. Number of effective stations at each epoch.

Figure 4. Two examples for station position raw time series (green) and residual time series (blue).

Figure 5. Region division for common mode error estimation.

Figure 6. Mean and maximum distances between stations in each region.

Figure 7. Spatial response of the first three principal components (Top: the spatial response for PC1. Middle: the spatial response for PC2. Bottom: the spatial response for PC3).

Figure 8. First principal component time series for eastern Asia.

Figure 9. Coordinate time series with CME correction.

Figure 10. Ratio of the first three eigenvalues of the corresponding components for each region.

Figure 11. Average spatial response to the first three principal components of each region.

Figure 12. Two examples of site position residual time series before CME filtering (blue) and after CME filtering (red).

Figure 13. RMS improvement percentage of residual time series before and after CME filtering.

Figure 14. Velocity difference and uncertainty difference before and after CME filtering.

Figure 15. Distribution of 365 reference stations in Scheme 2.

Figure 16. Helmert transformation parameters between IGS and ITRF2020 solutions before (blue) and after (red) CME filtering (left: Scheme 1, right: Scheme 2).

Figure 17. Uncertainty of transformation parameters before (blue) and after (red) CME filtering (left: Scheme 1, right: Scheme 2).

Figure 18. RMS of post-transformation residual time series before (blue) and after (red) CME filtering (left: Scheme 1, right: Scheme 2).

Figure 19. Region division for common mode error estimation (Division 2).

Figure 20. RMS improvement percentage of residual time series before and after CME filtering (Division 2).

Table 1. Velocity difference and uncertainty difference before and after CME filtering (unit: mm/yr).

Component	Velocity Difference			Velocity Uncertainty Difference
Component	E	N	U	E	N	U
Mean	−0.00	0.01	0.00	0.01	0.01	0.03
Max	0.17	0.48	0.45	0.09	0.05	0.13

Table 2. Means and standard deviations of the transformation parameters before and after CME filtering.

Parameters	After Deducting Linear Fitting Terms (Scheme 1)		After Deducting Linear Fitting Terms (Scheme 2)
Parameters	Unfiltered	Filtered	Unfiltered	Filtered
Tx (mm)	0.04	0.00	0.05	0.01
Tx (mm)	±0.60	±0.21	±0.78	±0.31
Ty (mm)	0.02	0.01	0.03	0.01
Ty (mm)	±0.43	±0.18	±0.65	±0.30
Tz (mm)	0.01	0.00	0.01	0.00
Tz (mm)	±0.45	±0.21	±0.47	±0.23
Rx (mas)	0.000	0.000	0.000	0.000
Rx (mas)	±0.012	±0.007	±0.017	±0.009
Ry (mas)	0.001	0.000	0.001	0.001
Ry (mas)	±0.015	±0.010	±0.018	±0.011
Rz (mas)	0.000	0.000	0.000	0.000
Rz (mas)	±0.009	±0.004	±0.013	±0.005
Scale (ppb)	0.01	0.00	0.01	0.00
Scale (ppb)	±0.13	±0.07	±0.15	±0.08

Table 3. Velocity difference and uncertainty difference before and after CME filtering (unit: mm/yr, Division 2).

Component	Velocity Difference			Velocity Uncertainty Difference
Component	E	N	U	E	N	U
Mean	−0.00	0.01	0.00	0.01	0.01	0.03
Max	0.28	0.41	0.42	0.08	0.05	0.10

Table 4. Means and STDs of the transformation parameters before and after CME filtering.

Parameters	With Deducting Linear Fitting Terms
Parameters	Unfiltered	Filtered, Division 1	Filtered, Division 2
Tx (mm)	0.04	0.00	0.00
Tx (mm)	±0.60	±0.21	±0.21
Ty (mm)	0.02	0.01	0.01
Ty (mm)	±0.43	±0.18	±0.17
Tz (mm)	0.01	0.00	0.00
Tz (mm)	±0.45	±0.21	±0.20
Rx (mas)	0.000	0.000	0.000
Rx (mas)	±0.012	±0.007	±0.004
Ry (mas)	0.001	0.000	0.000
Ry (mas)	±0.015	±0.010	±0.005
Rz (mas)	0.000	0.000	0.000
Rz (mas)	±0.009	±0.004	±0.004
Scale (ppb)	0.01	0.00	0.00
Scale (ppb)	±0.13	±0.07	±0.05

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Wang, R.; Chen, J.; Dong, D.; Tan, W.; Liao, X. Regional GNSS Common Mode Error Correction to Refine the Global Reference Frame. Remote Sens. 2024, 16, 4469. https://doi.org/10.3390/rs16234469

AMA Style

Wang R, Chen J, Dong D, Tan W, Liao X. Regional GNSS Common Mode Error Correction to Refine the Global Reference Frame. Remote Sensing. 2024; 16(23):4469. https://doi.org/10.3390/rs16234469

Chicago/Turabian Style

Wang, Ruyuan, Junping Chen, Danan Dong, Weijie Tan, and Xinhao Liao. 2024. "Regional GNSS Common Mode Error Correction to Refine the Global Reference Frame" Remote Sensing 16, no. 23: 4469. https://doi.org/10.3390/rs16234469

APA Style

Wang, R., Chen, J., Dong, D., Tan, W., & Liao, X. (2024). Regional GNSS Common Mode Error Correction to Refine the Global Reference Frame. Remote Sensing, 16(23), 4469. https://doi.org/10.3390/rs16234469

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Regional GNSS Common Mode Error Correction to Refine the Global Reference Frame

Abstract

1. Introduction

2. Data Source

3. CME Extraction

4. CME Impact on Coordinate Time Series and Contribution to Reference Frame Refinement

4.1. Accuracy Improvement of Coordinate Solutions

4.2. Analysis of Velocity Estimation

4.3. Accuracy Improvement of the Reference Frame

4.3.1. Helmert Transformation Parameters

4.3.2. Uncertainty of the Helmert Transformation Parameters

4.3.3. Post-Transformation Residuals Time Series

5. Discussion

6. Conclusions

Author Contributions

Funding

Data Availability Statement

Acknowledgments

Conflicts of Interest

Appendix A

Appendix A.1. Fitting Model of Coordinate Time Series

Appendix A.2. Mathematical Model of Principal Component Analysis

Appendix A.3. Calculation Formula of Helmert Transformation Parameters

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI