A Study of Rank Defect and Network Effect in Processing the CMONOC Network on Bernese

High-precision GPS data processing on Bernese has been employed to routinely resolve daily position solutions of GPS stations in the Crustal Movement Observation Network of China (CMONOC). The rank-deficient problems of the normal equation (NEQ) system and the network effect on the frame alignment of NEQs in the processing of CMONOC data on Bernese still present difficulties. In this study, we diagnose the rank-deficient problems of the original NEQ, review the efficiency of the controlled datum removal (CDR) method in filtering out the three frame-origin-related datum contents, investigate the reliabilities of the inherited frame orientation and scale information from the fixation of the GPS satellite orbits and the Earth rotation parameters in establishing the NEQ of the CMONOC network on Bernese, and analyze the impact of the network effect on the position time series of GPS stations. Our results confirm the nonsingularity of the original NEQ and the efficiency of the CDR filtering in resolving the rank-deficient problems; show that the frame origin parameters are weakly defined and should be stripped off, while the frame orientation and scale parameters should be retained due to their insufficient redefinition from the minimal constraint (MC) implementation through inhomogeneous and asymmetrical fiducial networks; and reveal the superiority of a globally distributed fiducial network for frame alignment of the reconstructed NEQs via No-Net-Translation (NNT) MC conditions. Finally, we attribute the two apparent discontinuities in the position time series to the terrestrial reference frame (TRF) conversions of the GPS satellite orbits, and identify it as the orbit TRF effect.


Introduction
In 1998, China Earthquake Administration led the creation of a national scientific project, namely the Crustal Movement Observation Network of China (CMONOC), to monitor crustal deformation in the China mainland and surrounding areas [1,2].Based on the solution of GPS station positions and corresponding position time series, researchers can capture the kinematic features of crustal deformation from the CMONOC network extent, especially in the active Tibetan zone, and gain insight into their dynamic mechanics [3][4][5].
GAMIT/GLOBK [6], GIPSY/OASIS [7], and Bernese [8] have been widely used to obtain high-precision daily solutions of GPS station positions.Compared with open source GAMIT/GLOBK and GIPSY/OASIS packages, Bernese is a commercial software which provides excellent support and services and regular updates, and shows rapid growth of user communities, including the Japanese Global Navigation Satellite System (GNSS) Earth Observation Network System (GEONET), the Europe Permanent Network (EPN), and the Chinese CMONOC [9,10].With regards to positioning algorithms and strategies on Bernese, station position parameters are directly estimated from independent double-difference baseline phase observation equations rather than through traditional "network adjustments" on GAMIT/GLOBK or "precise point positioning (PPP)" on GIPSY/OASIS.The EPN center [11] and the Center for Orbit Determination in Europe (CODE) [12] have proposed empirical recommendations for the routine computation of station positions for regional GPS networks (such as the CMONOC network) on Bernese, and the pivotal procedure is to keep all the satellite orbits and the Earth rotation parameters fixed to their "standard" values which are extracted from postprocessing products (such as International GNSS Service (IGS) orbit products) in the final least squares estimation of position parameters.In this way, the procedure may introduce external datum contents into the final outcome of the normal equation (NEQ) matrix, which should essentially contain the necessary geometrical correlations among GPS stations and thus have 3-7 rank deficiencies [13,14].
The estimation of daily positions for the GPS stations is based on the employment of the least squares algorithm on the daily NEQ to align station positions (in the Cartesian coordinate system) to a specific terrestrial reference frame (TRF).The conventional and classical method to handle this procedure is to apply minimal constraints (MCs) to the NEQ [15,16]; this is called "the MC approach".The aim of the MC approach is to impose several indispensable TRF-related conditions on the NEQ, make the NEQ invertible, and obtain the solution upon the condition of "free network" [16].However, due to the peculiar processing strategy on Bernese, the NEQ matrix is already invertible before the MC implementation, and the additional imposing conditions of the MC approach generally yield the solution of station positions, as has been elaborately verified in Kotsakis and Chatzinikos (2017) [13].Kotsakis and Chatzinikos (2017) have signified the existence of the so-called "network distortions" in several GNSS networks which are attributed to the rank-deficient problems in the processed NEQs, and they point out that the distortions can be several millimeters when they apply the MC approach to the original NEQ [13].As GPS data from the CMONOC network are also processed on Bernese, there should be similar rank-deficient problems in the NEQs.
The remedy for such rank-deficient problems is to strip off the external datum contents from the original NEQs in data processing on Bernese, and keep the contained geometrical correlations of the GPS stations intact in the meanwhile [13,14,17].The conventional methods to remove this issue are inspired by Sillard and Boucher (2001) [16]-filtering out the selected TRF-related contents with a similarity transformation matrix which implicitly contains the datum information.Rebischung et al. (2016) have successfully removed the redundant ranks of the position solutions before the combination of solutions from IGS centers establishing the IGS14 reference frame.Their strategy is efficient when the position solution SINEX file lacks adequate constraint information, and guarantees the validity of the procedure of combination on the solution level [14].Kotsakis and Chatzinikos (2017) have proposed an NEQ-level method to dispose of the rank-deficient problems, namely the controlled datum removal (CDR) filtering approach [13].They perform the CDR filtering before the implementation of the MCs, filter out the frame-origin-related datum contents from the original NEQ, and obtain the "ideal" reconstructed NEQ with redundant ranks eliminated and without affecting the intrinsic geometrical correlations for the GPS stations.After the CDR filtering, the rank defect character of the NEQ gets verified, and the obtained position solutions show no network distortions compared with the "free network" solution [13].Their elaborate analysis and discussion have demonstrated that there is an urgent need to properly dispose of the rank-deficient problems of NEQs, especially for the ones that are processed on Bernese.To further evaluate the reliabilities of the other two datum contents, the frame-orientation-and frame-scale-related datum contents, we should increase the extent of the removal of datum contents in the original NEQ and investigate corresponding effects in detail.
Apart from previous rank-deficient problems, the other fatal issue in the frame alignment of NEQs is the network effect [18].The network effect can cause the alias of erratic periodic or even secular signals in the "true" GPS position records, bias the outcome of station position time series, and thus confuse the study of the deformation fields.Strategies have been investigated to mitigate these network effects, including the use of well-distributed subnetworks in the frame alignment [19,20], the exclusion of the scale factor in the establishment of MCs [21][22][23], the modification of the stochastic function model in the least squares estimation [24], and the joint estimation of the low degrees of the deformation field with position parameters [22].The results show that the standard approach, by applying MC through the whole fiducial network, is particularly harmful for the realization of the position solution to its "real" value if the fiducial network is not well distributed, and all the advanced improvements can drastically reduce the magnitude of the aliasing network effects [23].As all the strategic improvements are adapted to the GAMIT/GLOBK solutions, specifically for implementing network adjustments, the impacts of the network effect and the related mitigating efforts should be reviewed on Bernese solutions which are computed by the one-step algorithm.Legrand et al. (2010) performed comparisons between Bernese solutions and TRF realized through regional and global fiducial networks, respectively, and their results reveal that the velocity field estimation through regional fiducial networks presents systematic effects with respect to the one estimated through global fiducial network, and corresponding estimations of the Euler pole differ significantly, too [25].With regards to their study, they have not taken the rank-deficient problems of the EPN weekly NEQs into account, and thus have totally ignored the consequential effects on the velocity field and the Euler pole estimation.We should review their comparisons based on the ideal NEQs, and further extend our research scope to the temporal behaviors of obtained position solutions, especially to the seasonal variabilities which are critical outcomes in view of the need to assess unmodeled loading effects in GNSS height time series [26,27].
The aim of this paper is to evaluate the recommended processing strategy (by the EPN center and CODE) for resolving station positions of the CMONOC network on Bernese, and to propose appropriate suggestions on the estimation of station positions from the daily NEQs.In Section 2, we briefly present the GPS dataset processed in this paper, and review the processing strategy for the computation of station positions on Bernese.In Section 3, we diagnose the rank-deficient problems in the NEQs from the processing of the CMONOC network on Bernese, and review the efficiency of the CDR filtering in filtering out the three frame-origin-related datum contents.In Section 4, we increase the extent of the removal of datum contents for the CDR filtering, investigate the reliabilities of both the frame orientation and frame scale parameters in the original NEQs, and analyze the impacts of the network effect on the position time series.In Section 5, we gain insight into the trajectories of the position time series and identify the satellite orbit TRF effect.In Section 6, we give conclusions of this paper.

GPS Data Set
In Figure 1, we present the spatial distribution of the GPS data set.There are more than 260 continuous and 2000 campaign GPS stations in the CMONOC network (Figure 1a).In the routine processing of CMONOC data, we also include 152 globally distributed International Terrestrial Reference Frame (ITRF) stations (Figure 1b) which are candidates for the ITRF2014 reference stations [28].
Specifically, we define two fiducial networks for the study of the network effect in the following sections, named as the S-set and the A-set.The S-set is formed of 70 globally distributed ITRF2014 reference stations which have steady temporal behaviors, good continuity qualities, and homogeneous and symmetric distribution on the Earth's surface as much as possible.The A-set includes all the ITRF2014 stations that are regionally distributed in the China mainland and surrounding areas, and thus envelops the CMONOC network (Figure 1b).

Data Processing Strategy on Bernese
The CMONOC data were processed on Bernese 5.2 software [8], forming a set of independent double-difference baselines and employing general models recommended by the International Earth Rotation Service (IERS) [29].In addition, Finite Element Solutions 2004 (FES2004) ocean tidal loading corrections [30], atmosphere S1-S2 tidal loading corrections [31], and absolute phase-center corrections for satellites and receivers, as issued by IGS [32], were applied to the formulation of the observation equations.The tropospheric delay was modeled by using the Vienna Mapping Function (VMF) [33], with zenith delay hourly estimated for the wet part and VMF model values assigned for the dry part.First-order ionospheric delays were eliminated by means of the ionosphere-free linear combination, and higher-order terms were modeled as recommended in Petrie et al. (2010) [34].Satellite orbits and Earth rotation parameters are recommended to be fixed in the final least squares estimation of the position parameters, with phase ambiguities resolved as much as possible and geometrical correlations among GPS stations fully considered.The TRF realization of daily position solutions was implemented by applying MCs via NNT (No-Net-Translation) conditions through the S-set fiducial network (Figure 1b) upon the NEQ.For a typical NEQ matrix N from the processing of the CMONOC network on Bernese, it is assumed that any nuisance parameters have been eliminated, and that only position parameters are retained in the NEQ.The NEQ is given in linearized form as where N is the NEQ matrix, X is the estimated station position parameters, X 0 is the a priori station coordinates, and u is the right-hand-side normal vector.

Diagnosing Rank-Deficient Problems and CDR Filtering of the Original NEQ
Kotsakis and Chatzinikos (2017) employed several indices and patterns to infer the datum contents of the NEQ, including the eigenvalues, the reference system effects, the geometrical assessment, the spectral assessment, and the "network distortions" [13].To diagnose the rank-deficient problems of the original NEQ from the CMONOC network, we choose to test three items-the eigenvalue analysis, the geometrical assessment, and the network distortion-for one general NEQ on 2014101.Then, we review the efficiency of the CDR filtering by testing the same three items for the reconstructed NEQ on 2014101.

Diagnosing Rank-Deficient Problems of the Original NEQ
Mathematically, the matrix condition number (the ratio between the maximum eigenvalue and the minimum one) can indicate the rank-deficient problems of a matrix.In Figure 2a, we present the eigenvalue spectra of the "free" NEQ.The matrix condition number is about 10 4 , which indicates that the NEQ matrix is invertible, and there is no rank deficiency in the original NEQ.There is an apparent gap between the 3rd and 4th eigenvalues, and according to Kotsakis and Chatzinikos (2017) [13], the smallest three eigenvalues largely reflect the weak frame origin information in the NEQ, which can also be revealed in the geometrical assessment.
The geometrical assessment is implemented by multiplying the NEQ matrix N with the so-called similarity transformation matrix G T (see Equation (2) in Kotsakis and Chatzinikos (2017) [13]), which represents the orthogonality level between those two matrices, and further reflects the strength of the contained datum-related contents in the NEQ [35].For more elaborate descriptions and details of the geometrical assessment, refer to Kotsakis and Chatzinikos (2017) [13].In Figure 2b, we present the cosine angles between the NEQ and the seven datum-related parameters, namely, three origin-related, three orientation-related, and one scale-related factors.The closer the cosine angles are to zero, the less information is contained in the NEQ for the corresponding datum parameters.As shown in Figure 2b, the frame-origin-related parts show smaller values at (−0.005, 0.005) than do the orientation-and scale-related parts at (−0.1, 0.1), which indicates that the definition of the datum origin is weaker than the datum orientation and scale in the original NEQ.As regards the influence of the redundant datum contents contained in the original NEQ on the obtained solution of station positions, Kotsakis and Chatzinikos (2017) [13] have proposed the so-called "network distortion" to quantitatively assess the effects.The network distortion is defined by the differences between the "free network" solution (note that here the free network solution is the direct inversion of the free NEQ matrix N in Equation ( 2)) and the "minimally constrained solution" (the solution with MC conditions imposed to the free NEQ matrix N in Equation ( 3)).Conventionally, the differences can be decomposed into the systematic part (the estimated similarity transformation parameters) and the individual part (position residuals of the similarity transformation fitting on the CMONOC stations) through a similarity transformation over all the GPS stations (Equation ( 4)).For our CMONOC network, we present the systematic network distortions in Table 1 and the individual network distortions in Figure 3a-c.Note that the minimally constrained solutions are computed via the NNT conditions with respect to the ITRF2014 frame based on the reference stations shown in Figure 1b (S-set), and the datum-related weight matrix W is set to an accuracy of 10 −5 m for the frame origin fixation.For the systematic part, the estimated rotation parameters are less than 0.01 milliarcseconds (mas) and the scale parameter is −0.24 ppb (Table 1).As there is no constraint imposed on the frame orientation and scale in the MC implementation, their emergences indicate that the redundant internal datum contents in the free NEQ interfere with the imposed MCs and thus corrupt the estimation of the station positions.Similar conclusions can be inferred from the individual position residuals.As shown in Figure 3a-c, the NNT solution is distorted by several millimeters compared with the free network solution, with both horizontal and vertical variation of up to 2 mm.The horizontal residuals show a compression-like behavior toward the center of the CMONOC network (Figure 3a), whereas the vertical residuals show inhomogeneous downward shift which gradually decreases in the northern part (Figure 3b).
In Equations ( 2)-( 4) above, XFREE is the free network solution of the station positions, X 0 is the a priori station positions, N is the original NEQ matrix, u is the right-side vector of the original NEQ system, XMC is the minimally constrained solution of the station positions, H is the coefficients matrix of the MCs, c is the residuals matrix of the MCs, W is a datum-related weight matrix, G FREE is the coefficients matrix of the similarity transformation derived from the free network solution of the station positions, θ is the seven similarity transformation parameters, and V is the residual of the transformation.
Therefore, we can conclude that the original NEQ from the CMONOC network has no rank defect character and contains three weakly defined frame origin datum contents, which can generate conspicuous systematic and individual network distortions in the minimally constrained solution of the station positions.

CDR Filtering of the Original NEQ
According to Kotsakis and Chatzinikos (2017), CDR filtering is used to reconstruct a singular NEQ and conforms to two particular properties: filtering out selected datum parameters and upholding the original data-related information [13].The reconstructed NEQ can be described as

CDR Filtering of the Original NEQ
According to Kotsakis and Chatzinikos (2017), CDR filtering is used to reconstruct a singular NEQ and conforms to two particular properties: filtering out selected datum parameters and upholding the original data-related information [13].The reconstructed NEQ can be described as where N is the original NEQ matrix, N is the reconstructed NEQ matrix, u is the original right-side normal vector, u is the reconstructed right-side normal vector, I is the identity matrix, and E is formed by selected rows of the similarity transformation matrix G and those rows correspond to the ill-defined datum parameters which we intend to redefine via MCs.As has been discussed in Section 3.1, the removal of the frame origin datum contents, which are weakly defined by the fixation of the GPS satellite orbits and the Earth rotation parameters on Bernese, is the primary need.Therefore, the matrix E in Equation ( 5) contains only the first three rows of the matrix G.
The spectral behavior of the reconstructed NEQ matrix after the CDR filtering is depicted in Figure 2c.We can see that the smallest three eigenvalues are in the order of 10 −13 and the apparent gap is in the order of 10 12 , which leads to the NEQ matrix being singular with three rank deficiencies.The cosine angle plots for the reconstructed NEQ are presented in the Figure 2d.We can see that the three rows of the matrix E are fully orthogonal to the reconstructed NEQ matrix N with cosine angles being less than 10 −16 , and the frame orientation and scale parameters seem to retain the same orthogonality level.
The network distortions of the reconstructed NEQ are presented in Table 1 and Figure 3d.We can see that the estimated orientation and scale parameters in the similarity transformation are exactly zero (Table 1), and the position residuals are also zero on all the CMONOC network stations (Figure 3d).In this case, the MC implementation is the same with the one in Section 3.1, and there is no systematic or individual network distortion, which indicates no interference between the data-related information in the reconstructed NEQ and the external datum-related information defined via NNT MCs.
Therefore, we can conclude that the CDR filtering is efficient in stripping off the redundant datum contents in the original NEQ, and the implementation of MCs via NNT conditions on the reconstructed NEQ is free of network distortion.

Frame Alignment Analysis
In Section 3.2, only frame origin parameters are removed by the CDR filtering to realize three rank deficiencies in the reconstructed NEQ, and the frame orientation and scale parameters are still retained as the same within the original free NEQ.In this section, we investigate the reliabilities of the other two sets of frame-related parameters by diagnosing the rank defect property of the reconstructed NEQ, the network distortions on the minimally constrained solutions, and the temporal behaviors of the corresponding position time series.Besides this, the so-called "network effect" will be analyzed to select the optimal fiducial network on the frame alignment of the NEQs.

Frame Orientation
By using a similar CDR filtering algorithm in Section 3.2, we filter out both origin and orientation datum contents from the original NEQ by establishing the matrix E containing the upper six rows of the similarity transformation matrix G.In rank defect analysis, the spectral behavior of the reconstructed NEQ with CDR filtering employed has six rank deficiencies and the six smallest eigenvalues have a gap of 10 12 relative to the rest of the spectrum.The geometrical datum assessment shows orthogonality in both frame origin and orientation parameters, which is similar to the results in Figure 2. To avoid adding to the length of this paper, corresponding figures are waived.
With CDR filtering employed for the original NEQ, the minimally constrained position solutions are computed by applying MCs to the reconstructed NEQ via No-Net-Translation and No-Net-Rotation (NNTR) conditions with respect to the ITRF2014 frame based on the globally distributed reference stations shown in Figure 1b (S-set).A full similarity transformation model is applied for fitting the free network solution and the NNTR solution to reveal network distortions.The estimated transformation parameters and corresponding residuals on 2014101 (as an example) are given in Table 2 and Figure 4. Comparing with the results in Figure 3d (perfect fitting with zero residuals), the similarity fitting is not perfect with post-fit residuals reaching up to 1-2 mm.The horizontal distortions show a rotated pattern, with clockwise rotation on the eastern part and anticlockwise rotation on the western part of the CMONOC network (Figure 4a), while the vertical distortions are symmetrical (Figure 4b).There is conspicuous scale leakage estimated by the similarity transformation (Table 2), and the estimated origin parameters are also different from previous ones shown in Table 1, although the fiducial network for the MC implementation is the same in both cases.similarity transformation (Table 2), and the estimated origin parameters are also different from previous ones shown in Table 1, although the fiducial network for the MC implementation is the same in both cases.
Those results reveal that, with the removal of both origin and orientation datum information from the original NEQ, the MC implementation via the NNTR conditions upon the reconstructed NEQ is not appropriate due to the interference between frame orientation parameters and origin/scale parameters.We conclude that the datum orientation information redefined by applying MCs via NNTR conditions to the reconstructed NEQ through an external fiducial network (S-set) is insufficient and not better than the one inherited from the original NEQ.The interference between estimated parameters is attributed to the inhomogeneous distribution of reference stations in the S-set fiducial network, which leads to imbalance in the matrix  and cross-correlations between datum orientation-related rows (the 4th to 6th rows of ) and other rows (the 1st to 3rd and 7th rows of ).Therefore, unless perfect fiducial networks can be established and employed, we should not filter out the frame orientation datum contents from the original NEQ, but retain those which are inherited from the fixation of the GPS satellite orbits and the Earth rotation parameters in processing the CMONOC network on Bernese.Those results reveal that, with the removal of both origin and orientation datum information from the original NEQ, the MC implementation via the NNTR conditions upon the reconstructed NEQ is not appropriate due to the interference between frame orientation parameters and origin/scale parameters.We conclude that the datum orientation information redefined by applying MCs via NNTR conditions to the reconstructed NEQ through an external fiducial network (S-set) is insufficient and not better than the one inherited from the original NEQ.The interference between estimated parameters is attributed to the inhomogeneous distribution of reference stations in the S-set fiducial network, which leads to imbalance in the matrix E and cross-correlations between datum orientation-related rows (the 4th to 6th rows of G) and other rows (the 1st to 3rd and 7th rows of G).Therefore, unless perfect fiducial networks can be established and employed, we should not filter out the frame orientation datum contents from the original NEQ, but retain those which are inherited from the fixation of the GPS satellite orbits and the Earth rotation parameters in processing the CMONOC network on Bernese.

Frame Scale
Analogously, we remove both origin and scale datum information from the original NEQ.The corresponding matrix E here contains four rows of the similarity transformation matrix G.The spectral behavior of the reconstructed NEQ after CDR has four rank deficiencies and the four smallest eigenvalues have a gap of 10 12 relative to the rest of the spectrum.The geometrical datum assessment shows orthogonality in both frame origin and scale parameters.
With CDR filtering employed, the minimally constrained solutions are computed by applying MCs to the reconstructed NEQ via No-Net-Translation and No-Net-Scale (NNTS) conditions with respect to the ITRF2014 frame based on the globally distributed reference stations shown in Figure 1b (S-set).A full similarity model is applied for fitting the free network solution and the minimally constrained solution to reveal network distortions.The estimated transformation parameters on 2014101 (as an example) are given in Table 2.The estimated rotation factors between the free network solution and the NNTS solution are exactly zero, and the similarity transformation gives a perfect fit with zero position residuals.
Therefore, by removing both frame origin and scale datum contents from the original NEQ, we can also avoid network distortions in the NNTS solution, and there is no interference between frame origin/scale parameters and frame orientation parameters.However, the estimated translation parameters are not exactly equal to the previous ones in the NNT solution (Table 1).The differences of the X-component, the Y-component, and the Z-component translation parameters are 0.1, 0.1, and 0.3 mm, respectively, which reflects the difference in the frame scale parameter between the estimated value and the inherited one in the NEQ.Thus, we can conclude that the frame origin parameters couple with the frame scale parameter in the NEQ system.
To further elaborate the coupling relationship, two sets of position time series were computed by applying MCs via the NNT condition to the NEQ with only origin-related datum contents removed (the NNT solution), and via NNTS conditions to the NEQ with both origin-and scale-related datum contents removed (the NNTS solution), respectively, and temporal behaviors of their differences were investigated.The comparison considers two basic aspects of the differences, namely the estimated full similarity transformation parameters between the NNT solution and the NNTS solution, and their position differences.For this test, the NNT and NNTS solutions were computed in ITRF2014 using the reference stations shown in Figure 1b (S-set) with a prior weight matrix set to an accuracy level of 10 −5 m for the frame origin and scale fixation.
The estimated similarity transformation parameters are given in Figure 5.The amplitudes of their temporal variation are 1-2 mm for the estimated translation factors, and about 1 ppb for the scale factor.Apart from the X-component translation parameter, all the other parameters show significant periodic behaviors.Note that the rotation parameters are exactly equal to zero in the transformation due to the irrelevant correlation between frame origin/scale parameters and frame orientation parameters.
The time series of position differences at station YANChi (YANC) and maps of corresponding derived signals are shown in Figure 6.Taking YANC as an example, we find that the temporal behaviors of the horizontal position differences are random with the magnitude of the variabilities being at the submillimeter level, while the vertical differences are remarkably large and have significant secular and periodic features (Figure 6a).Extending to other stations in the CMONOC network, both the pattern of the vertical trends and the vertical annual signals that are derived from the vertical time series of position differences are homogeneous; the magnitude of the vertical trends can reach up to 0.3-0.5 mm/a; the amplitude of the annual signals are about 1 mm; and the annual harmonic functions reach their peaks in February for almost all CMONOC stations.
The above comparison results tell us that although the S-set fiducial network is globally distributed, there is still interference between frame origin and scale information in NEQs due to the improper redefinition of both frame origin and scale parameters in the NNTS solution based on unevenly and asymmetrically distributed reference stations.This imperfect spatial distribution leads to cross-correlations between the datum scale-related row (the 7th row of G) and origin datum-related rows (the 1st to 3rd rows of G).In former works, Tregoning and van Dam (2005) [21], Lavallée et al. (2006) [22], and Collilieux et al. (2012) [23] have concluded that the scale parameter should be excluded from the conditions in applying MCs to the NEQ if the fiducial network is not evenly distributed, and the height component of the position time series always suffers the most serious aliasing effects; this has been confirmed in our study.
To further verify those conclusions, we altered the fiducial network from the S-set to the A-set, and the NNT and NNTS solutions were aligned to ITRF2014 through reference stations in the A-set fiducial network (Figure 1b).Analogous similarity transformation and position differences were performed between those two solutions, and the corresponding annual behaviors were enhanced, not only in the vertical component but also in the horizontal component where the amplitudes can reach up to 3 mm.
In conclusion, in the absence of perfect fiducial networks to redefine the scale datum information, we should not remove the frame-scale-related datum contents from the original NEQ, but retain its property that is inherited from the fixation of the GPS satellite orbits and the Earth rotation parameters in establishing the NEQ of the CMONOC network on Bernese.To further verify those conclusions, we altered the fiducial network from the S-set to the A-set, and the NNT and NNTS solutions were aligned to ITRF2014 through reference stations in the A-set fiducial network (Figure 1b).Analogous similarity transformation and position differences were performed between those two solutions, and the corresponding annual behaviors were enhanced, not only in the vertical component but also in the horizontal component where the amplitudes can reach up to 3 mm.
In conclusion, in the absence of perfect fiducial networks to redefine the scale datum but retain its property that is inherited from the fixation of the GPS satellite orbits and the Earth rotation parameters in establishing the NEQ of the CMONOC network on Bernese.

Network Effect
Geodetic position time series may contain several millimeter nonlinear variations in amplitude.These are mainly attributed to loading effects, which are not modeled in the coordinate calculation procedure, and to remaining systematic errors.On employment of MCs to realize the TRF (Equation ( 5)), the conventional coordinates of reference stations are always linearly interpolated or extrapolated to the specific epoch, sometimes with discontinuities and postseismic deformation modifications applied [28].Therefore, nonlinear (mostly seasonal) variations in the reference stations can bias the "standard" coordinates of reference stations, thus biasing the frame alignment of the position solution, and the seasonal variations are aliased into the position solution [36].Naturally, the aliasing effect is dependent on the configuration of the fiducial network and on the behaviors of fiducial candidates.Therefore, it is called the "network effect".

Network Effect
Geodetic position time series may contain several millimeter nonlinear variations in amplitude.These are mainly attributed to loading effects, which are not modeled in the coordinate calculation procedure, and to remaining systematic errors.On employment of MCs to realize the TRF (Equation ( 5)), the conventional coordinates of reference stations are always linearly interpolated or extrapolated to the specific epoch, sometimes with discontinuities and postseismic deformation modifications applied [28].Therefore, nonlinear (mostly seasonal) variations in the reference stations can bias the "standard" coordinates of reference stations, thus biasing the frame alignment of the position solution, and the seasonal variations are aliased into the position solution [36].Naturally, the aliasing effect is dependent on the configuration of the fiducial network and on the behaviors of fiducial candidates.Therefore, it is called the "network effect".
To investigate the network effect and to select the optimal fiducial network for frame alignment of the NEQ, we compute two sets of NNT solutions based on the S-set fiducial network and the A-set fiducial network, respectively (Figure 1b), with a prior weight matrix set to an accuracy level of 10 −5 m for the frame origin fixation.Hereafter, we call them the S-Solution and the A-Solution, respectively.The comparisons focus on the position differences (S/A differences); their corresponding temporal behaviors, including the secular and periodic variabilities derived by using least squares adjustment; and the temporal behaviors of the estimated similarity translation parameters between the S-Solution and the A-Solution, which are shown in Figures 7 and 8, respectively.Conventionally, YANC is employed as an example to show the temporal behaviors of position solutions and position differences (Figure 7a-c), and secular trends and nonlinear annual signals are derived from the time series of position differences of all the CMONOC stations (Figure 7d-g).Note that the position time series and position difference in Figure 7a-c are de-trended.We find that the root mean squares (RMS) of the S-Solution are larger than those of the A-Solution in all three components of the topographic coordinate time series (Figure 7a,b) due to the evident presence of periodic signals in the S-Solution and their absence in the A-Solution.The S/A differences are dominated by periodic features especially in the N-S and vertical directions (Figure 7c), and the corresponding amplitudes are about 1 mm and 3 mm, respectively.Extending to the whole CMONOC network, the horizontal secular trends are homogeneous pointing to the northwest direction with a magnitude of about 0.3 mm/a (Figure 7b), while the vertical trends are random with smaller magnitudes; the patterns of horizontal annual signals are respectively rotational and symmetrical with amplitudes reaching up to 1-2 mm on the N-S component and the E-W component (Figure 7c,d), while the vertical annual signals are much larger in amplitude, reaching up to 3 mm, and show a homogeneous pattern reaching its maximum values in January of every year.As regards the similarity transformation parameters for the comparison between the S-Solution and the A-Solution, the estimated values are all zero for the rotation and scale parameters; the remaining three estimated translation parameters, which reflect the differences of frame origin realization, are considerably larger with temporal variations being in the interval (−15 mm, 15 mm) and show evident annual behaviors, especially in the Y-component and the Z-component, with amplitude being about 5 mm (Figure 8).Note that there are two apparent discontinuities in the E-W component of the S-Solution and S/A difference, and the X-component of the estimated similarity translation parameter (Figure 7a,c and Figure 8a).
From previous investigations, we conclude that the network effect aliases into the A-Solution, and thus mitigates the "true" annual signals in the position time series.These mitigation effects are attributed to the fact that the implicit frame origin parameters (on applying MCs) absorb most of the annual differences between the standard position solution and the NNT solution, because periodic behaviors on the reference stations are nearly homogeneous due to the limited extent of the A-set fiducial network, and thus are treated as holistic behaviors of the whole fiducial network.However, for the S-Solution, as the S-set fiducial network is globally distributed, the reference stations have different (or even opponent) periodic behaviors, and thus the absorption by implicit frame origin parameters does not occur.This conclusion can also be verified by the temporal behaviors of estimated similarity translation parameters (Figure 8) and their annual behaviors shown in Figure 8 should be blamed on the network effect on the A-Solution.Furthermore, the network effect not only affects the periodic signals, which are a critical outcome to assess unmodeled loading effects in position time series, but also the secular trends, which are important to the study of crustal deformation.Therefore, we have reason to believe that the globally distributed and stable fiducial network, the S-set, should be superior to the regionally distributed one, the A-set, on the frame alignment of the NEQs.

Discussion
As we have inferred from previous rank defect analysis and frame alignment analysis, only the frame origin datum contents should be removed from the original NEQ by employing the CDR filtering, and datum orientation and scale information should be retained in the reconstructed NEQ.By this means, the reconstructed NEQ fully inherits the geometrical correlations of the GPS stations in the CMONOC network, and datum orientation and scale information defined by the fixation of the GPS satellite orbits and the Earth rotation parameters in the processing of CMONOC data on Bernese, and is thus adequate to compute the solution of station positions by applying NNT MCs.Naturally, the inherited datum information should rely on the TRF of the GPS satellite orbits, and even on the antenna phase center calibration (PCC) models for the satellites and the receivers [37].
In Table 3, we present the historical conversions for the TRF of the GPS satellite orbits declared in the header part of the GPS satellite orbits files (sp3), starting with ITR96 in 1999 and gradually turning to IGb08 till the deadline of the included CMONOC data.Due to the irrelevant relationship between the external redefined datum origin and the internal inherited datum orientation/scale information, the obtained NNT solutions, and the corresponding position time series, are affected by the TRF conversions of the satellite orbits, which has been recognized in Figures 7a and 8a with the emergences of two apparent discontinuities.Those two discontinuities, 2001336 and 2006309, are respectively coincident with the TRF conversion epochs when the TRF of the GPS satellite orbits converted from IGS97 to IGS00, and from IGb00 to IGS05, which are highlighted in bold font in Table 3.As this effect is dependent on the TRF of the GPS satellite orbits, hereafter we name it as the "orbit TRF effect".Altamimi et al. (2002) have signified that the ITRF2000 and its successors (IGS00, IGb00) are the first in the history of the International Terrestrial Reference Frame combining unconstrained space geodesy solutions that are free from any tectonic plate motion model [38], thus leading to an absolutely different definition of the TRF parameters (origin, orientation, and scale) compared with their predecessor ITR97 or IGS97.Towards ITR05 and its successor IGS05, Altamimi et al. (2007) have

Discussion
As we have inferred from previous rank defect analysis and frame alignment analysis, only the frame origin datum contents should be removed from the original NEQ by employing the CDR filtering, and datum orientation and scale information should be retained in the reconstructed NEQ.By this means, the reconstructed NEQ fully inherits the geometrical correlations of the GPS stations in the CMONOC network, and datum orientation and scale information defined by the fixation of the GPS satellite orbits and the Earth rotation parameters in the processing of CMONOC data on Bernese, and is thus adequate to compute the solution of station positions by applying NNT MCs.Naturally, the inherited datum information should rely on the TRF of the GPS satellite orbits, and even on the antenna phase center calibration (PCC) models for the satellites and the receivers [37].
In Table 3, we present the historical conversions for the TRF of the GPS satellite orbits declared in the header part of the GPS satellite orbits files (sp3), starting with ITR96 in 1999 and gradually turning to IGb08 till the deadline of the included CMONOC data.Due to the irrelevant relationship between the external redefined datum origin and the internal inherited datum orientation/scale information, the obtained NNT solutions, and the corresponding position time series, are affected by the TRF conversions of the satellite orbits, which has been recognized in Figures 7a and 8a with the emergences of two apparent discontinuities.Those two discontinuities, 2001336 and 2006309, are respectively coincident with the TRF conversion epochs when the TRF of the GPS satellite orbits converted from IGS97 to IGS00, and from IGb00 to IGS05, which are highlighted in bold font in Table 3.As this effect is dependent on the TRF of the GPS satellite orbits, hereafter we name it as the "orbit TRF effect".Altamimi et al. (2002) have signified that the ITRF2000 and its successors (IGS00, IGb00) are the first in the history of the International Terrestrial Reference Frame combining unconstrained space geodesy solutions that are free from any tectonic plate motion model [38], thus leading to an absolutely different definition of the TRF parameters (origin, orientation, and scale) compared with their predecessor ITR97 or IGS97.Towards ITR05 and its successor IGS05, Altamimi et al. (2007) have pointed out that ITR05 uses data time series of station positions and daily Earth Orientation Parameters (EOPs) instead of long-term solutions (used in ITR00) as input data to establish the frame [39].According to those upgrades and the rapid increase of valid data (including GPS, Satellite Laser Ranging (SLR), Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS), and Very Long Baseline Interferometry (VLBI) data) in the establishment of the TRF, there should be apparent differences in the frame parameters among ITR97, IGS00, and IGS05, to which the discontinuities in Figures 7a and 8a must be attributed.Meanwhile, the upgrades for the TRF of the GPS satellite orbits may also lead to secular trend differences in the position time series, as depicted in Figure 7a,b, with the N-S component of the time series presenting piecewise linear behaviors which are segmented on 2006309.Note that the absence of apparent discontinuities in the A-Solution (Figure 7b) and their presence in the X-component of estimated translation parameters (Figure 8a) signify that the orbit TRF effect is homogeneous in the A-set extent, and is absorbed by the implicit similarity translation parameters between the A-Solution and the standard position solution, thus masking the "true" signals in the corresponding position time series.This result further demonstrates the superiority of the S-set fiducial network, which shows the "true" discontinuities, to the A-set fiducial network for frame alignment of the NEQs.
Since CMONOC data are cumulatively and routinely processed and the orbit-related sp3 files are daily updated, the TRF of the GPS satellite orbits has been gradually changed by IGS analysis centers; thus, our daily position solutions inevitably suffer from the orbit TRF effect.Therefore, in the study of crustal deformation based on the CMONOC solutions, time series analysis should be carefully operated, and appropriate discontinuities (sometimes with trend changes) should be set up on the TRF conversion epochs to counteract this effect.Another effective approach is to reprocess CMONOC data with the GPS satellite orbits in a consistent TRF.

Conclusions
We have processed 17 years of GPS data from the CMONOC network, including more than 260 continuous stations and 2000 campaign stations, since 1999.Together with another 152 ITRF stations, the underlying GPS network is processed on Bernese 5.2 software, and the daily solution of GPS station positions is computed via the MC approach.The recommended data processing strategy was employed in our work, with the GPS satellite orbits and the Earth orientation parameters kept fixed in the least squares estimation of position parameters, to obtain daily NEQs on Bernese.
Through rank defect analysis of the original NEQ from the CMONOC network, we found that it does not have proper singularities, and contains weakly defined frame-origin-related datum contents.The nonsingularity property of the NEQ disobeys the basic discipline of the MC implementation; it thus yields the estimation of station positions, and generates conspicuous network distortions in the advanced minimally constrained solution.The systematic part of the network distortion reveals that the external additions of the frame origin parameters via NNT MCs interfere with the internal redundant datum contents of the NEQ, and the individual part shows millimeter-level (up to 2 mm in our example of the NEQ on 2014101) position differences which corrupt the geometrical correlations between GPS stations in the NEQ.To dispose of these rank-deficient problems, the CDR filtering method was employed to remove the three weak frame origin parameters.CDR filtering is efficient for realizing the necessary three rank deficiencies in the NEQ, without interfering with internal geometrical and other datum contents or generating any network distortions in the minimally constrained solution.Furthermore, we increased the extent of the removal of datum contents and investigated the reliabilities of the other two datum contents (orientation and scale) in the NEQ.We found that the external redefined frame orientation parameters interfere with both frame origin and scale parameters, and lead to several-millimeter network distortions in the minimally constrained NNTR solution upon the reconstructed NEQ, while the external redefined frame scale parameter only couples with the frame origin parameters, and generates no network distortion in the NNTS solution.To further investigate the coupling relationship, we analyzed the temporal behaviors of the position differences between the NNT solution and the NNTS solution (minimally constrained upon the reconstructed NEQ), and found that the height component of the position time series suffers from conspicuous secular and periodic aliasing effects for the NNTS solution.Hence, we conclude that the weakly defined frame-origin-related datum contents must be filtered out from the original NEQ, and be redefined to guarantee the rightful MC implementation; on the other hand, we should not redefine the frame orientation and scale parameters through an inhomogeneous and asymmetrical fiducial network, but retain their values that are inherited from the fixation of the GPS satellite orbits and the Earth rotation parameters in establishing the NEQ of the CMONOC network on Bernese.On the basis of the minimally constrained NNT solutions, we analyzed the impacts of the network effect on the frame alignment of the NEQs.Through comparisons between the S-Solution and the A-Solution, we found that the network effect aliases into the A-Solution due to the limited extent of the underlying fiducial network, and biases the derived annual signals from the position time series.Therefore, we conclude that the globally distributed S-set fiducial network is superior to the regionally distributed A-set fiducial network on the frame alignment of the NEQs.
Lastly, we identified the orbit TRF effect, which leads to apparent discontinuities or even secular trend changes to the position time series, and attributed it to the TRF conversions of the GPS satellite orbits.The proper method to counteract this effect is either to manually set up discontinuities (and sometimes velocity changes) in the position time series analysis, or to reprocess CMONOC data with the GPS satellite orbits in a consistent TRF.

Figure 1 .
Figure 1.Spatial distribution of the GPS data set.(a) The spatial distribution of GPS stations in the Crustal Movement Observation Network of China (CMONOC) network; (b) the spatial distribution of the included International Terrestrial Reference Frame 2014 (ITRF2014) stations.

Figure 1 .
Figure 1.Spatial distribution of the GPS data set.(a) The spatial distribution of GPS stations in the Crustal Movement Observation Network of China (CMONOC) network; (b) the spatial distribution of the included International Terrestrial Reference Frame 2014 (ITRF2014) stations.

Figure 2 .
Figure 2. Spectral behaviors of the normal equation (NEQ) matrix on 2014101.Subfigures (a) and (c) depict eigenvalues decomposed from the NEQs with and without controlled datum removal (CDR) filtering employed, respectively; (b) and (d) depict the cosine angles between the seven rows of the similarity transformation matrix and the rows of the NEQ matrix with and without CDR filtering employed, respectively.Rx, Ry and Rz are X-, Y-and Z-components of the orientation-related parameters, respectively.Tx, Ty, Tz are X-, Y-and Z-components of the origin-related parameters, respectively.Note that the unit of the cosine angle is 10 -16 in the upper figure of (d).

Figure 2 .
Figure 2. Spectral behaviors of the normal equation (NEQ) matrix on 2014101.Subfigures (a,c) depict eigenvalues decomposed from the NEQs with and without controlled datum removal (CDR) filtering employed, respectively; (b,d) depict the cosine angles between the seven rows of the similarity transformation matrix and the rows of the NEQ matrix with and without CDR filtering employed, respectively.Rx, Ry and Rz are X-, Y-and Z-components of the orientation-related parameters, respectively.Tx, Ty, Tz are X-, Y-and Z-components of the origin-related parameters, respectively.Note that the unit of the cosine angle is 10 −16 in the upper figure of (d).

Figure 3 .
Figure 3. Post-fit residuals between the free network solution and the No-Net-Translation (NNT) solution.Subfigures (a)-(c) depict residuals between the free network solution and the NNT solution upon the NEQ without CDR filtering employed, and (d) depicts similar results but for the NNT solution upon the NEQ with CDR filtering employed; (a) and (b) are horizontal and vertical maps of post-fit residuals in the CMONOC network, respectively; (c) and (d) depict station-by-station residual variations in the topographic coordinates.Note that the units in (d) are nm.

Figure 3 .
Figure 3. Post-fit residuals between the free network solution and the No-Net-Translation (NNT) solution.Subfigures (a)-(c) depict residuals between the free network solution and the NNT solution upon the NEQ without CDR filtering employed, and (d) depicts similar results but for the NNT solution upon the NEQ with CDR filtering employed; (a,b) are horizontal and vertical maps of post-fit residuals in the CMONOC network, respectively; (c,d) depict station-by-station residual variations in the topographic coordinates.Note that the units in (d) are nm.

Table 2 .
Similarity transformation parameters between the free network solution and the NNTR/ NNTS solution.

Figure 4 .
Figure 4. Maps of post-fit residuals between the free network solution and the No-Net-Translation and No-Net-Rotation (NNTR) solution.Subfigures (a) and (b) depict horizontal and vertical residual maps in the CMONOC network, respectively.

Figure 4 .
Figure 4. Maps of post-fit residuals between the free network solution and the No-Net-Translation and No-Net-Rotation (NNTR) solution.Subfigures (a,b) depict horizontal and vertical residual maps in the CMONOC network, respectively.
RemoteSens. 2018, 10, 357   12 of 20 excluded from the conditions in applying MCs to the NEQ if the fiducial network is not evenly distributed, and the height component of the position time series always suffers the most serious aliasing effects; this has been confirmed in our study.

Figure 5 .
Figure 5. Temporal behaviors of estimated similarity transformation parameters between the NNT solution and the No-Net-Translation and No-Net-Scale (NNTS) solution.Subfigures (a)-(c) depict the time series of the estimated translation parameters in three components, respectively, and (d) depicts the time series of the estimated scale parameter.

Figure 5 .
Figure 5. Temporal behaviors of estimated similarity transformation parameters between the NNT solution and the No-Net-Translation and No-Net-Scale (NNTS) solution.Subfigures (a)-(c) depict the time series of the estimated translation parameters in three components, respectively, and (d) depicts the time series of the estimated scale parameter.

Figure 6 .
Figure 6.Time series of position differences and maps of derived vertical signals.(a) The time series of position differences between the NNT solution and the NNTS solution at station YANChi (YANC), which is highlighted as a blue triangle in (b) and (c); (b) the distribution of the derived vertical secular trends; and (c) phase (arrow direction) and amplitude (arrow length) of the derived vertical annual signals.The length of the vector in (c) represents the amplitude, and the phasors point directly in the direction of the maximum ('up' = max in January, 'right' = max in April, 'down' = max in July, 'left' = max in October [27]).

Figure 6 .
Figure 6.Time series of position differences and maps of derived vertical signals.(a) The time series of position differences between the NNT solution and the NNTS solution at station YANChi (YANC), which is highlighted as a blue triangle in (b) and (c); (b) the distribution of the derived vertical secular trends; and (c) phase (arrow direction) and amplitude (arrow length) of the derived vertical annual signals.The length of the vector in (c) represents the amplitude, and the phasors point directly in the direction of the maximum ('up' = max in January, 'right' = max in April, 'down' = max in July, 'left' = max in October [27]).

Figure 7 .
Figure 7. De-trended position time series (and their differences) and maps of horizontal secular trends and annual signals derived from the S/A differences.Subfigures (a,b) depict the de-trended S-Solution and A-Solution at station YANC, which is highlighted as a blue triangle in (d)-(g), respectively; (c) depicts corresponding S/A differences; (d) depicts the horizontal secular trends derived from the S/A differences; (e)-(g) depict phase (arrow direction) and amplitude (arrow length) of periodic signals on N-S, E-W, and vertical components, respectively.The explanation of the arrow in (e)-(g) refers to Figure 6.The brown lines in (a,c) indicate apparent discontinuities in the S-Solution and S/A differences.

Figure 8 .
Figure 8. Temporal behaviors of estimated similarity transformation parameters between the S-Solution and the A-Solution.Subfigures (a)-(c) depict time series of the X-component, the Ycomponent, and the Z-component translation parameters, respectively.The brown lines in (a) indicate apparent discontinuities in the time series of the X-component translation parameter.

Figure 8 .
Figure 8. Temporal behaviors of estimated similarity transformation parameters between the S-Solution and the A-Solution.Subfigures (a)-(c) depict time series of the X-component, the Y-component, and the Z-component translation parameters, respectively.The brown lines in (a) indicate apparent discontinuities in the time series of the X-component translation parameter.

Table 1 .
Similarity transformation parameters between the free network solution and the NNT solutions upon the NEQ with and without CDR filtering employed, respectively.

Table 1 .
Similarity transformation parameters between the free network solution and the NNT solutions upon the NEQ with and without CDR filtering employed, respectively.

Table 2 .
Similarity transformation parameters between the free network solution and the NNTR/NNTS solution.

Table 3 .
The conversion history for the terrestrial reference frames (TRFs) of the GPS satellite orbits.