Improving the Specular Point Positioning Accuracy of Ship-Borne GNSS-R Observations in China Seas Based on Comprehensive Geophysical Correction

The accurate modeled GNSS-R reflection delay, which is indispensable for the quantification, modeling, and correction of the GNSS-R altimetry sea-state bias, can be obtained based on the accurate modeled position of the specular point. At present, the reflection surface model of the specular point positioning still has the mean dynamic topography (MDT) error and the deviation of the vertical (DOV) error relative to the instantaneous sea surface. In this study, the following studies have been carried out. Based on the ship-borne GNSS-R observations in China seas, we introduced various elevation parameters including the MDT to correct the elevation error of the reflection surface. We introduced the DOV based on the elevation correction, and the DOV correction positioning method was proposed to correct the slope error of the reflection surface. The specular point was positioned on the instantaneous sea reflection surface. We verified the instantaneous sea reflection surface model and the specular point positioning results, analyzed the relationship between the position correction distance and the reflection incident angle, and discussed the spatial distribution characteristics of the MDT correction distance. The results showed that the reflection surface modeling and the specular point positioning were accurate. The positioning error increased to varying degrees with the increase of the reflection incident angle. The MDT correction improved the positioning by 0.91 m, and the DOV correction further improved the positioning by 0.12 m. Based on the combined application of the two kinds of correction positioning, the positioning was comprehensively improved by 0.99 m. The MDT correction of China seas gradually increased from the north to south. While in the regional sea areas, it gradually decreased from the north to south and showed randomness. The relative position between the antennas and their random changes introduced uncertainty, which can be reduced by integration. The new instantaneous sea reflection surface model and the corresponding specular point positioning method can provide accurate modeled reflection delay for the sea-state bias correction of ship-borne GNSS-R observations, and they can be extended to satellite-borne global observations.


INTRODUCTION
GNSS-R altimetry can provide global high-coverage sea surface height (SSH) observations for research on global ocean mesoscale and sub-mesoscale processes and global climate change as an auxiliary means to traditional altimeters (Martín-Neira, 1993;Wu et al., 1998;Stammer et al., 2000;Hajj and Zuffada, 2003;Saynisch et al., 2015;Zuffada et al., 2015;Li et al., 2016;Xie et al., 2018). Due to the waves, the sea surface is rough, skewed, and rapidly changing, especially for the reflection of the GNSS signal considering the wavelength. The signal is not specularly reflected on the sea surface, which leads to the deviation of the specular point position on the reflected power waveform. This deviation introduces bias to the delay of the reflected signal relative to the direct signal, which cannot be ignored for high-precision altimetry (Hajj and Zuffada, 2003;Rius et al., 2010;Yang and Zhang, 2012). Due to the principle of GNSS-R observation and the complexity of sea surface roughness, the quantitative analysis and correction of reflection delay sea-state bias (SSB) has been one of the technical difficulties and constraints to improve the accuracy of GNSS-R altimetry. This is a key problem to be solved for highly accurate SSH retrieval in GNSS-R altimetry satellite missions (Rodriguez, 1988;Rius et al., 2010;Martín-Neira et al., 2011).
Based on the difference between the observation reflection delay and the modeled reflection delay, the SSB is expected to be quantified. The difference can be used as a prior knowledge to construct the empirical parameter model of the SSB, and thus, it can be predicted and compensated. This requires a large number of observations to suppress random errors and to obtain high coverage of the empirical model parameters. On the other hand, obtaining accurate reflection delay based on model calculation is indispensable in the quantification of the SSB. This requires correcting the specular point geometric positioning error introduced by the difference between the modeled sea surface and the instantaneous sea surface. Ship-borne observations have unique advantages for obtaining both observation delay and modeled delay. For modeled delay, the path of the direct and the reflected signals passing through the atmosphere in shipborne scenario can be considered the same, and no additional delay is caused. The effects of hull's attitude change such as pitch and roll on observation delay and modeled delay can be considered to be the same and offset. In addition, the voyage of large research vessel (RV) is usually long and the route covers different sea areas, and this can support the study of the spatial distribution characteristics of specular point positioning correction. Based on the ship-borne GNSS-R, we conducted sea surface altimetry experiments in China seas (Gao et al., 2020) and carried out modeling and correction of delay SSB. This article focuses on the important basis for obtaining accurate modeled delay-the research on sea reflection surface modeling and specular point positioning.
The modeling of sea reflection surface has experienced the process of refinement of standard sphere, the earth ellipsoid, the geoid, and the ocean tidal surface (Wu et al., 1997;Wagner and Klokocnok, 2003;Kostelecky et al., 2005;Gleason and Gebre-Egziabher, 2009;Rius et al., 2010;Semmling et al., 2014;Jales and Unwin, 2017;Wu et al., 2019a;Wu et al., 2019b). However, the mean dynamic topography (MDT) error has not been corrected. The mean sea surface (MSS) is the average sea surface height after excluding interannual, semiannual, seasonal, and other periodic sea surface height changes over a longer period of time. The MDT is the fluctuation of the MSS relative to the geoid with a global amplitude of −2 ∼2 m. The MDT is the change of sea surface height caused by the external forces of global average flow field, marine hydrological factors, atmospheric pressure, and other nontidal factors. The most important influence on the MDT is caused by the global average flow field, and its highest point is the west Pacific affected by the Kuroshio (Andersen, 2011;Liu, 2014). In this study, the ship-borne experiment's route passes through the influence area of the Taiwan warm current and the Yellow Sea warm current, tributaries of the Kuroshio. Furthermore, due to the difference of the earth's gravity field, the MSS has different slopes relative to the ellipsoid at different locations, that is, the geoid deviation of the vertical (DOV). Martín-Neira's analysis of the effect of the sea surface slope is based on the assumption that the position of the specular point remains unchanged (Martín-Neira, 1993). The resulting slope error needs to be corrected in the modeling of the sea reflection surface and the positioning of the specular point.
There are differences in the MDT of different oceans and seas. Our RV passed through the Yellow Sea, the East China Sea, and the South China Sea. These sea areas have significant sea surface topography differences (Andersen et al., 2016), which will inevitably lead to different specular point positioning corrections. Understanding the spatial distribution characteristics of the position correction distance can help develop targeted strategies of postprocessing and positioning error correction in different sea areas. This study provides a regional approach to acquire this prior knowledge, and it is expected to be extended to satellite-borne global observations.
In this study, we used ship-borne GNSS-R observations in China seas. Based on the geoid and the ocean tidal reflection surface model constructed in our previous research, we sequentially introduced the MDT and the DOV to correct the elevation error and the slope error with the corresponding specular point positioning method. The specular point is finally positioned on the instantaneous sea surface, and the positioning accuracy is improved. This study has laid the foundation for obtaining accurate modeled reflection delay and for the quantification and correction of the SSB. Data and Model introduces the ship-borne data and the geophysical models used, Methodologies introduces the reflection surface modeling and the specular point positioning methods, Results and Discussion discusses the results of the positioning, and Conclusion summarizes and prospects.

Ship-Borne GNSS-R Equipment and Data
We carried Xiang Yang Hong 06 RV and used GNSS-R equipment to carry out a 3-week sea surface altimetry Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 720470 experiment. The route traverses most areas of China seas, including the Yellow Sea, the East China Sea, and the South China Sea. The hardware of the GNSS-R receiver system mainly included two antennas and a GNSS IF raw-data recorder. The up-looking antenna received the direct GPS/ BDS signals, and the down-looking antenna received the signals reflected from the sea surface. The GNSS-R antennas are about 12 m high from the water surface. Figure 1 shows the side-section geometry of the ship-borne GNSS-R antennas. The center line of the up-looking antenna is vertical. The down-looking antenna is installed under the direct antenna, pointing diagonally downward to the sea surface. The angle c between the center lines of the down-looking antenna and the up-looking antenna is 150°. The distance d between the two antennas' phase centers is 0.283 m. We randomly selected 17,000 samples. The positioning of the specular point is based on the GNSS position, the position of the GNSS-R antenna, and the reflection surface. We regarded the phase centers of the direct antenna and the reflection antenna as one same position. The position is calculated from the geodetic coordinates of the ship-borne GNSS navigation antenna combined with the relative position of the navigation antenna and the GNSS-R equipment in the hull coordinate system. The relative position of the two is calculated by their coordinates in the hull coordinate system measured by the total station. The bow direction is the ship's geodetic coordinate at the sampling time pointing to the ship's next geodetic coordinate. The GNSS position is obtained from the ephemeris file.

Geophysical Models
The instantaneous sea reflection surface model is constructed by introducing a series of geophysical parameters into the earth ellipsoid. The used geophysical models included the geoid undulation of the EGM2008 model, the ocean tidal heights of the TPXO model, and the DTU15 MDT elevations. On this basis, we introduced the DOV from the GGMPlus gravitational field to correct the sea surface slope errors.

The EGM2008 Geoid Undulation
The EGM2008 model order is up to 2,159, equivalent to a spatial resolution of about 5′ × 5′. The commission error of the geoid undulation in the ocean area where the latitude is less than 66°is 5.8 cm. The commission error implied by EGM2008 geoid undulation. We used the highest spatial resolution product which is interpolated to a 1′ × 1′ grid, and the difference of interpolated values from those obtained via harmonic synthesis does not exceed ± 1 mm (Pavlis and Saleh, 2005;Pavlis et al., 2012).

The TOPEX/POSEIDON Tidal Model (TPXO)
The TPXO tidal model has performed harmonic analysis along the track on the altimetry data of satellites and incorporates data of tide gauge and satellites in the shallow water areas. The nonlinear 1/4-period day tidal constituent has also been considered to improve the accuracy in the offshore. High-resolution regional assimilation models are developed and added to TPXO global model calculation result with a resolution of 1/6°. These regions are mainly closed and semienclosed oceans, and most of the continental parts shelve coastal areas. The resolution in China seas is 1/30°. TPXO also uses the 1′ bathymetric data in available offshore areas to improve accuracy and spatial resolution (Egbert et al., 1994;Egbert and Ray, 2000;Egbert and Erofeeva, 2002). The average deviation between TPXO and the Global Ocean Tide (GOT) Model is 0.25 cm, and the standard deviation and the RMSE are both approximately 1.5 cm (Liu, 2014). The RMSE of the main tidal constituents of TPXO in China seas is of centimeter level (Wang et al., 2010). The TPXO model is suitable for our high-resolution ship-borne experiments which were mainly carried out in offshore.

The DTU15 MDT
The DTU15 MDT is obtained from the MSS height based on the satellite data from 1993 to 2015 minus the EGM2008 geoid fluctuation, with a spatial resolution of 1′ × 1′. In the study area, the short-wavelength residual MDT signal in the DTU15 MDT associated with EGM2008 ranges within ± 5 cm (Andersen et al., 2019). The difference between the MSS DTU15 and CNES15 models in the study area is basically within the range of ± 2 cm (Andersen et al., 2015).

The GGMplus DOV
The GGMplus model is a synthesis of GRACE and GOCE satellite gravity and EGM2008 and short-wave terrain gravity, with a spatial resolution of 0.002°, approximately 220 m. The DOV data include meridian component and prime component, and the accuracy is about one arc-second (Hirt et al., 2013).

The Ephemeris
GNSS orbital information is obtained from GNSS ephemeris files provided by the International GNSS Service (IGS) (Montenbruck et al., 2017). Unless otherwise specified, the position information used in this study is based on the ECEF WGS-84 coordinate system. Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 720470

METHODOLOGIES
The modeled reflection path of GNSS-R starts from the GNSS transmitter to the specular point on the modeled reflection surface and then to the receiver's antenna. For the modeled reflection delay, the elevation and the slope of the reflection surface determine the position of the specular point and then determine the modeled reflection path. In order to obtain the modeled reflection delay without sea-state error, it is necessary to construct an ideal smooth sea surface model with the elevation and the slope close to the instantaneous sea surface at the moment and the position of the specular point. We introduced the geophysical parameters that are one order of magnitude lower than the delay SSB or more into the reflection surface modeling (see Geophysical Models). In the previous research, we have gradually constructed the geoid reflection surface model and the ocean tidal reflection surface model. On this basis, we further introduced the MDT and the DOV to correct the elevation and slope errors. We constructed the instantaneous sea reflection surface model and positioned the specular point on it.

MDT Correction and Positioning
The modeling of the sea reflection surface is realized in the process of specular point positioning. The specular point is initially positioned on the reference ellipsoid to obtain its initial longitude l and latitude b, and the elevation is 0 (Wu et al., 1997). In the process of transforming geodetic coordinates to space coordinates, we introduced the geophysical parameter elevation at the position and the time of the specular point (Wu et al., 2019a). The spatial coordinate of the specular point is shown in where H G , H T , and H MDT are the geoid undulation, the ocean tidal height, and the MDT elevation of the specular point. N a/ [1 − e 2 sin 2 (b)] , where a is the long radius of the WGS-84 reference ellipsoid and e is the first eccentricity of the ellipsoid. Then, we calculated the incident angle, the emergence angle, and the geocentric angle and iterated them with weight to correct the position of the specular point (Wu et al., 1997;Wu et al., 2019a). The MDT elevation correction components ρ X , ρ Y , and ρ Z of the specular reflection point in the X, Y, and Z directions at each iteration are, respectively, H MDT cos(b)cos(l), H MDT cos(b) sin(l), and H MDT sin(b). And, λ X , λ Y , and λ Z are the sum of the geoid and the tidal correction components. Based on the comprehensive consideration of positioning accuracy and iteration times, the iterative cutoff threshold is set to the modulo of the difference between the incident angle and the emergence angle and is less than 10 −8 rad. When the iteration reaches this condition, the elevation correction positioning ends. The threshold is satisfied after n iteration corrections, and the total MDT correction components σ X , σ Y , and σ Z are obtained as follows: The correction distance of the positioning accuracy of the MDT correction D MDT is the distance between the specular points before and after the MDT correction, given as follows:

DOV Correction and Positioning
The essence of the DOV correction is to use the prime component η and the meridian component ζ as the correction to correct the normal direction of the specular point successively. In the space coordinate system, we corrected the ellipsoid normal direction (x 1 , y 1 , z 1 ) to the geoid normal direction (x 2 , y 2 , z 2 ) to correct the reflection geometry and the position of the specular point. The steps are as follows: (1) η correction: solving |x 2 ′| and |y 2 ′| in plane XOY.
We obtained the geoid normal (x 2 , y 2 , z 2 ) based on the above calculation. We used the nonapproximate normal projection correction method in the following correction. By directly solving the spatial geometric relationship between the projection of the normal on the plane and the reflection path, the specular point is corrected to the vertical plane of the geoid. The positioning error caused by the radial normal difference is reduced, and the influence of approximate substitution is reduced (Wu et al., 2019a). The positioning accuracy is further improved towards normal direction.

MDT Correction and Positioning
The average MDT elevation of our samples is 0.66 m, the maximum is 0.70 m, and the minimum is 0.58 m. In the space coordinate system, the correction distance is 0.91 m. The average error of the DTU15 MDT at the specular point is 2.67 cm. We calculated the difference between the specular point positions before and after adding the MDT error, and the mean value is 3.57 cm. The correction distance in the X, Y, and Z directions are −0.36, 0.53, and 0.18 m, respectively, and the corresponding modulus are 0.46, 0.53, and 0.39 m, respectively.

DOV Correction and Positioning
The GGMplus DOV data do not cover all the global oceans. There are 4,246 samples with DOV data. The average values of the prime components and the meridian components of the DOV are −0.0023°and 0.0012°, respectively. In the space coordinate system, the correction distance is 0.12 m. The error of the DOV model is one arc-second (Hirt et al., 2013). We calculated the difference between the positions before and after adding the one arc-second DOV error. The mean value of the difference is 1.96 cm. In the X, Y, and Z directions, the correction distances are −0.04, −0.03, and −0.02 m, respectively, and the moduli of the correction distance are 0.04, 0.03, and 0.04 m, respectively.

Combined Correction
We combined the MDT and DOV correction positioning, and the comprehensive positioning correction distance is the final position compared to the position without these two kinds of corrections. For the samples with the DOV data, the combined correction distance is 0.99 m in the space coordinate system. We calculated the difference between the positions of the specular points before and after adding the MDT error and the DOV error, and the mean value is 5.23 cm. In the X, Y, and Z directions, the  Figure 3 shows the specular point positioning correction distance of the MDT correction D MDT and the corresponding reflection incident angle θ. It can be seen that the θ of most segments changes from large to small and then to large, which is the usual process of GNSS-R equipment from being visible to invisible to a GNSS satellite. And, some of the θ changes from small to large and then to small, which is another relative motion mode. The change of D MDT is very consistent with θ for almost all the segments.

Model Verification
In order to further verify the corrected positioning result based on the instantaneous sea reflection surface model, D MDT is compared with the simulated positioning correction distance D MDT '. Figure 4 shows the elevation correction geometry of the reflection surface and we can have D MDT ′ H MDT /cosθ. SP is the specular point before correction, and SP' is the corrected specular point. θ ranges from approximately 15°to approximately 70°in this study.
We calculated |ΔD MDT | |D MDT '-D MDT |, and the average value is 1.09 m × 10 -4 m, the standard deviation is 2.29 m × 10 -4 m, and D MDT is very close to D MDT '. The correlation coefficient between D MDT and D MDT ' is 97.66%. Figure 5 shows the fitted straight line of D MDT and D MDT ', the slope is 1.004 ± 0.003, SSE of the fitted straight line is 19.98, and the RMSE is 3.43 × 10 -2 , which are small. The correlation does not decrease significantly with the increase of D MDT or D MDT '. We arranged D MDT and D MDT ' in the ascending order of θ, and the correlation between D MDT and D MDT ' is 99.95%. Since the two are very similar, Figure 6 shows D MDT (the blue dots) and mean (H MDT )/cosθ (the red curve) in order to distinguish their trends of change. It can be seen that the changing trends of D MDT and mean (H MDT )/cosθ are very consistent and θ is the main influencing factor of D MDT . On the far right side of Figure 6, a small number of samples have large θ, resulting in large D MDT , which is consistent with the mean (H MDT )/cosθ of the corresponding θ. The above results have verified the high accuracy of the MDT correction positioning.

Relationship Between the Position Correction Distance and the Reflection Incident Angle
We arranged ΔD MDT in the ascending order of θ, as shown in Figure 7. ΔD MDT is centered at 0 and has the characteristics of positive and negative symmetrical distributions. After approximately the 10,000th sample, as θ increases, ΔD MDT increases. This is because as θ increases, the reflection path lengthens and hence the uncertainty introduced by the relative position between the antennas and their change increases. The symmetry feature gradually disappears after approximately the 10,000th sample; this is because large θ introduces additional increase in ΔD MDT . The symmetry feature still exists after approximately the 15,000th sample, and it can be considered that this symmetrical distribution feature exists in the entire reflection incident angle range covered by the sampling. We arranged |ΔD MDT | in the ascending order of θ and performed linear fitting, as shown in Figure 8. The trend of |ΔD MDT | increases with θ. A small number of samples near the minimum and the maximum values of θ correspond to larger deviations.   Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 720470 7 The DOV correction distance D DOV is arranged in the ascending order of θ (see Figure 9). D DOV shows an increasing trend with the increase of θ, but it is not so consistent with the change of θ like D MDT . For the samples of a single track, in the incident plane, when θ increases, the increasing direction of θ is not necessarily the same as the changing direction of the normal of the incident plane. When the two directions are consistent, D DOV will increase with θ; otherwise, D DOV will decrease. Therefore, the change of D DOV with θ has randomness. The comprehensive position correction distance D MDT+DOV is the distance between the specular points with and without the MDT and the DOV correction. D MDT+DOV is arranged in the ascending order of θ as shown in Figure 10. D MDT+DOV and θ generally have the same trend but are not completely consistent. As introduced above, the inconsistency is introduced by the randomness of D DOV relative to θ.

Spatial Distribution Characteristics of the MDT Positioning Correction Distance
The H MDT and D MDT in China seas covered by the sampling (17°N  ∼ 35°N) show a gradual increase from the north to south (see Figure 11 and Figure 12). For both H MDT and D MDT , the Yellow Sea is the lowest, the East China Sea is higher, and the South China Sea is the highest (see Table 1). In the two sea areas from the southern part of the Yellow Sea to the northern part of the East China Sea (32°N ∼ 35°N) and the northern part of the South China Sea (17°N ∼ 21°N), both H MDT and D MDT show a gradual decrease from the north to south. They increase near Xiamen and reach maximum in the sea area from 20°N to 21°N of the South China Sea. In the South China Sea, where the sampling coverage of latitude and longitude are both high, H MDT and D MDT have a tendency to gradually decrease from the northwest to southeast (away from the coast in the northwest). In the entire China seas and some regional seas, D MDT is consistent with the trend of H MDT .
H MDT changes monotonously and smoothly in the entire China seas and the regional sea areas. Different from H MDT , the partial spatial variation of D MDT presents randomness. This feature is more obvious in the northern part of the South China Sea (17°N ∼  21°N) where the sampling coverage of latitudes and longitudes are both high. The randomness of the spatial distribution of D MDT is due to the fact that in addition to H MDT , the determinants also  Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 720470 8 include the reflection incident angle, the distance between the transmitter and the specular point, and the distance between the specular point and the down-looking antenna. For a continuous sample sequence of the same GNSS satellite, the abovementioned parameters in the reflection geometry are also continuously changing, and this change will cause the D MDT of the sequence to have monotonicity. The sequences with different monotonicities are intersected and connected to form the entire measurement segments along the track as shown in Figure 12. This makes the overall spatial variation of the segments of D MDT present randomness. We believe that this feature is also present in the satellite measurement segments. In addition, there is a distance between the up-looking antenna and the down-looking antenna. The relative position of the two antennas changes with the incident plane, the incident angle, and the hull attitude. The uncertainty introduced to the model reflection geometry is estimated to be in the order of decimeters, and it introduces random errors to the correction positioning result. The impact on the MDT model near the coast may introduce errors in the correction results. In the DTU15 MDT error map, the error near Xiamen does not increase significantly compared with other parts of the study area, and the number of samples in Xiamen is small, so this will not be the main source of error.

CONCLUSION
Accurate modeled reflection delay is indispensable for quantifying and correcting GNSS-R sea-state bias. Constructing an accurate instantaneous sea reflection surface model is the key to improve the accuracy of the modeled reflection delay. The MDT and the DOV are nonnegligible errors in the reflection surface modeling and need to be corrected. Based on the geoid and the ocean-tidal reflection surface model, this study has introduced the MDT to further correct the elevation error of the reflection surface model. Then, the DOV is introduced to correct the slope error of the reflection surface model, and the corresponding specular point positioning method is proposed. The specular point is finally positioned on the instantaneous sea surface.
The MDT positioning correction distance is very consistent with the geometric simulation result, the model, and the positioning accuracy. The MDT correction improves the positioning accuracy by 0.91 m, and the DOV correction further improves the positioning accuracy by 0.12 m. Based on the combined application of the two kinds of correction, the positioning accuracy is improved by 0.99 m. The MDT correction positioning error increases with the reflection incident angle. It is presumed that the relative position between the antennas and their change introduces greater uncertainty as the reflection incident angle increases. The correlation between the DOV correction positioning error and the reflection incident angle decreased compared to the MDT correction. This is because the difference between the direction of the DOV and that of the increase in the reflection incident angle introduces randomness. The positioning correction distance by the combined application of the two kinds of correction is consistent with the overall trend of the reflection incident angle. The MDT correction distance of China seas gradually increases from the north to south, the Yellow Sea being the lowest, the East China Sea being higher, and the South China Sea being the highest. In some partial sea areas, the MDT correction distance gradually decreases from the north to south. The MDT correction distance presents randomness locally. The randomness is introduced by the intersection between measurement segments of FIGURE11 | Track of Xiang Yang Hong 06 (the pink curve) and the H MDT along the track.  different GNSS satellites and the random changes of the relative position of the GNSS-R antennas. The distance between the antennas can be effectively reduced by hardware integration. The instantaneous sea reflection surface model constructed in this study is expected to be applied to satellite GNSS-R sea surface altimetry to provide accurate modeled reflection delay for the separation, quantification, and modeling correction of sea-state bias. We hope to build a gridded global two-dimensional reflection surface model based on satellite observations, which can be directly used to position the specular points. The issues of interpolation and griding will be discussed in depth.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary materials; further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
FW was responsible for the study design, data processing, scientific analysis, and manuscript writing. WZ took part in project management and manuscript modification. ZL and XS were involved in data processing. All authors contributed to the article and approved the submitted version.