The two linearly independent solutions of Legendre equation are called the first and second kind of Legendre functions, respectively. When the differential equation takes the eigenvalue, the first kind of Legendre function is interpreted as a polynomial, so the independent variable can take any value (except infinity). The second kind of Legendre function is still infinite series, diverging when the independent variable is equal to ±1 and converging when the absolute value is greater than 1. Since Legendre equation belongs to the type of hypergeometric equation, we give the expression of arbitrary order derivatives of different special functions of this type equation. Therefore, the hypergeometric expression of the first kind  of Legendre function and its theoretical relationship with other special functions are given directly. In view of the complexity of solving the second kind of Legendre function, the hypergeometric expression of the second kind of Legendre function  is directly given by using the series expansion method.

Data gaps between the GRACE and GRACE-FO missions have prevented continuous monitoring of the terrestrial water storage changes (TWSC). We use the multiple linear regression model, take the GRACE/GRACE-FO solutions as reference values, the precipitation, air temperature, and model-simulated terrestrial water storage data as prediction parameters, and consider three different strategies to reconstruct the continuous TWSC of nine major river basins in China from April 2002 to December 2021. The results show that the reconstructed results based on the de-trended and de-seasonalized signal reconstruction strategy are slightly better than those dependent on the de-trended signal reconstruction strategy. These two reconstructed results are better than those using the overall signal reconstruction strategy, and this advantage is more evident in the basins with intense human activities or glaciers melt (such as Haihe and Luanhe rivers, Yangtze river, southeast rivers, and inland river basins). In addition, the performance of reconstructed result is also affected by the signal-to-noise ratio of GRACE/GRACE-FO data and the correlation of prediction parameters with GRACE/GRACE-FO data.

We quantitatively analyze the drought events in the Yangtze river basin from 2003 to 2021 by using GRACE data to explore the feasibility of satellite gravimetry in monitoring regional drought. Firstly, we use five GRACE data products(CSR_SH, JPL_SH, GFZ_SH, CSR_M and JPL_M) released by three centers to retrieve the terrestrial water storage anomaly(TWSA) in the Yangtze river basin, and calculate the terrestrial water storage deficit(WSD) and water storage deficit index(WSDI). We then compare the results of five data products in combination with meteorological drought data(SPI, SPEI and scPDSI). Lastly, we analyze the drought events in the Yangtze river basin from 2003 to 2021. The results show that GRACE data products released by different centers have a certain impact on quantifying drought severity in the Yangtze river basin. WSDI has the significant correlation with SPEI on 6-month time scale, with the highest correlation coefficient of 0.66 and the lowest correlation coefficient with scPDSI of 0.54. Precipitation is a dominated factor affecting the change of terrestrial water storage in the Yangtze river basin. The most serious drought event in the Yangtze river basin occurred in the summer and autumn of 2019. The drought intensity was 2.31 and lasted for 10 months. The cumulative loss of water storage reached 415 Gt. The spatial distribution of WSDI of this drought event shows it was most serious in September 2019, with extreme drought areas. WSDI can reflect the temporal and spatial variation of drought distribution in the Yangtze river basin and play an important role in monitoring global and large-scale regional drought.

Based on GRACE and GRACE-FO gravity satellite data, we use the moving T-test technique to determine the characteristic time point of the mutation of terrestrial water storage in the Haihe river basin, and analyze its spatiotemporal evolution characteristics. Combining multivariate trend analysis and contribution rate quantification method, we explore the causes of terrestrial water storage changes in the basin, and discuss the contribution of the south-to-north water diversion project. The results show that before (from January 2004 to January 2015) and after (from January 2015 to October 2020) the south-to-north water diversion project, the decline trends of terrestrial water storage in the Haihe river basin are -17.19 mm/a and -13.49 mm/a, respectively, relieved by about 24%, and the mitigation trend gradually increases from south to north. The ratio of contribution rate of human activities and climate change to the terrestrial water storage changes in the basin is 7∶3, anthropogenic water consumption is the main reason for the loss of water storage in the Haihe river basin all year round. The contribution rate and trend of the water transfer volume show a state of increasing and rising in the annual, indicating that the south-to-north water diversion project has great potential in improving regional terrestrial water storage.

Based on the M_W8.0 earthquake in northern Peru on May 26, 2019, using the global ionospheric map(GIM) obtained by the center for orbit determination in Europe(CODE) and the GPS-TEC observation data provided by the international GNSS service, we analyze the changes of ionosphere total electron content(TEC) and equatorial ionospheric anomaly before the earthquake. Using the TEC grid data from May 10 to 26, 2019, we analyze the variation of equatorial ionospheric anomaly and ionospheric TEC with geographic latitude before the earthquake, and statistically analyze the daily variation anomalies of ionospheric TEC before the earthquake using the GPS-TEC data of RIOP and GLPS stations. The results show that three days before the earthquake, the equatorial ionospheric anomaly “double peaks” disappeared, and the ionospheric TEC decreased significantly, reaching about 10 TECu. At 16:00~20:00 UT three days before the earthquake, the ionospheric TEC shows a bimodal curve with latitude, and there is a trough near the epicenter. During this period, the ionospheric TEC at RIOP and GLPS stations also decreased significantly.

We use the original curve analysis and normalized rate method to analyze the changes of three earth resistivity data sets within 300 km from the epicenter of Menyuan earthquake. Judging from the original curve, the EW component of Jinyintan station has a sharp jump and a downward trend of retaining annual change (-3.3%) in the two months before the earthquake. The NS component of Shandan station has a downward trend (-1.9%), and the EW component has an upward trend (3.1%). The NS component of Wuwei Xiaoxigou station has a downward trend (－3.3%), the EW component has an upward trend (2.8%), and the trend of EW component turned flat after the earthquake. The normalized rate curves of the three earth resistivity stations all exceeded the threshold before the Menyuan earthquake, but their normalized rate patterns are different. The sudden change of the hour value, the accelerated change of the daily average and monthly average values, and the turning change, especially the turning change in the stage of rising to the high value, which can be analyzed as an abnormal shape. When the above changes are not obvious, the normalized rate method should be combined to extract the abnormality.

Aiming at the problems of poor timeliness and low efficiency of conventional data processing in large-scale GNSS network, we design a shared memory model data parallel algorithm from two dimensions of multi-period and multi-subnet respectively based on GAMIT/GLOBK software. It further implements a two-layer parallel solution for spatio-temporal integration of large-scale GNSS network in the time and space domains. The results show that the solution breaks through the limitation of traditional serial processing of GNSS data by the software, which has poor timeliness and low utilization of multi-core computing resources. In the test environment, the maximum speedup ratio is up to 19.39, which fully exploits the computing power of the computer and greatly improves the timeliness of large-scale GNSS network data processing.

We use the satellite precise clock offset product of IGS BDS-3 in 2021(GPS week 2 138 to 2 190) to analyze the periodic change and construct a high-precision model function based on the periodic results. The results show that the modeling accuracy of the additional eight periods reaches 0.049 ns, increasing accuracy by about 70% compared with the quadratic polynomial results. Meanwhile, the performance of precise point positioning(PPP) is analyzed based on the modeling results of BDS-3 satellite clock offset. Although the PPP convergence time of the modeling BDS-3 satellite clock offset product is slightly slower than that of the IGS final products, it can achieve centimeter-level positioning. Moreover, the modeling products can greatly save memory space compared with the current IGS discretization products.

To analyze and evaluate the orbit determination performance of the new BDS-3 signals (B1C and B2a), we collect the observation data of 69 MGEX stations that are evenly distributed around the world from May 1 to June 30, 2021 to evaluate the precise orbit determination performance of B1C and B2a signals. The results of the 48 h overlapping arc show that the 3D RMS value of the BDS-3 MEO satellite orbit based on B1C/B2a combined observations is better than 6 cm, the radial value is better than 2 cm, and the laser ranging check residual is better than 4 cm. In addition,under the conditions of single BDS-3 MEO and BDS-3 MEO+Galileo, the orbit determination performance of new and old signals of the former are basically the same, and the difference of orbit determination accuracy is not more than 1 mm; for the latter, B1C/B2a new signal shows better orbit determination performance than B1I/B3I old signal, with its 3D RMS increased by 0.38 cm.

Aiming at the problems of low gait recognition rate, poor universality and stability of traditional peak detection, we propose a pedestrian dead reckoning method based on improved FSM gait detection. The proposed method focuses on the current gait state, devides a gait cycle into five states and introduces the peak threshold, and carries out experiments from different motion states and mobile phone modes.The results show that the algorithm not only has high gait recognition accuracy, but also reduces the time cumulative error. It has certain reliability, stability and can meet the needs of indoor positioning.

We use the smartphone Xiaomi 11 Lite, which could track the B1I/B1C/B2a frequency of BDS-3 and the L1/L5 frequency of GPS, as the research object, and do static observation experiment and dynamic experiments with different sites respectively to analyze the quality of raw observation data of triple-frequency BDS and dual-frequency GPS output by smartphone and the performance of single-frequency SPP. The experimental results show that the low cost GNSS antenna of smartphone leads to poor signal quality, and the anti-interference of GPS L1 frequency and BDS B1C frequency of smartphone is strong. In static and dynamic experiments, the positioning result of BDS B1C frequency of smartphone is relatively stable. In static experiment, the positioning result of GPS L5 frequency is better than that of GPS L1 frequency, and in dynamic experiment, the positioning result of GPS L1 frequency is better than that of GPS L5 frequency.

The coordinated universal time offset(UTCO) parameter broadcast by the global navigation satellite system(GNSS) is updated with the modernization upgrading navigation systems. In particular, the UTCO parameters for GPS and BDS have been significantly changed. From the perspective of model, the first-order traceability model has been replaced with a second-order traceability model, compatible with the first-order model. From the perspective of resolution, the parameter resolution has generally decreased, improving sensitivity to small changes in UTCO. By investigating the UTCO parameters in the GPS LNAV messages, GPS CNAV messages, BDS D1/D2 messages, BDS B-CNAV messages, Galileo messages and GLONASS messages of the four major navigation systems GPS, BDS, Galileo and GLONASS, and based on the UTCr-UTC(k) data published by international bureau of measurement(BIPM), we integrate the GNSS UTCO parameters into the fast coordinated universal time(UTCr) to analyze the performance of UTCO parameters broadcast by various systems. Meanwhile, we adopt the time offset of GNSS system monitored by the national time service center(NTSC) and the circular T bulletin published by BIPM as references to evaluate the coordinated universal time offset error(UTCOE) of GNSS system. The results show that the UTCO parameters of GPS CNAV messages perform better than GPS LNAV messages, and the UTCO parameters of BDS B-CNAV messages perform better than BDS D1/D2 messages within the experimental period. The UTCO parameters of GPS and Galileo perform better than BDS and GLONASS, and BDS and GLONASS had some systematic errors in the UTC traceability.

We adopt the function model with damping factor and use the genetic algorithm(GA) assisted nonlinear least squares(NLS) to solve the phase parameters. The results show that:1) Compared with the standard cosine function model, the correlation coefficient between the inversion phase and soil moisture is significantly improved, and the inversion results are more stable. In the three elevation angle ranges of 5°-15°, 5°-20°, and 5°-25°, the correlation coefficients can all be greater than 0.68. The difference between the correlation coefficients of different elevation angles is less than 0.07. 2) The inversion accuracy has been improved to varying degrees, with R² increased by 5.72%-76.06%, RMSE decreased by 6.12%-24.24%, MAE decreased by 2.7%-28.3%. By applying the inversion results of this method in multi-satellite linear regression model, the average RMSE decreased by 10%.

To analyze the accuracy of GNSS water vapor inversion performed by surface temperature and pressure products in China based on the latest MERRA-2 reanalysis data released by NASA, we combine the measured temperature and pressure data of 609 surface weather stations in China, 48 GNSS stations and co-location radiosonde stations. We evaluate MERRA-2 temperature and pressure products and their accuracy in GNSS PWV inversion. The results show that: 1) The annual average bias of MERRA-2 pressure and temperature data are -0.01 hPa and 0.38 K, respectively, and the annual average RMSE are 1.08 hPa and 2.66 K, respectively. The MERRA-2 reanalysis temperature and pressure products have high accuracy in China. 2) The temperature and pressure products of MERRA-2 reanalysis data show a negative deviation in most regions of China, and the order of accuracy from high to low is the southern region, the northern region, the northwest region and the Qinghai-Tibet Plateau region. 3) Compared  with the PWV measured at the co-location sites, the average RMSE of GNSS PWV derived from MERRA-2 temperature and pressure products is 2.16 mm, and can better reflect the daily variation of PWV. Therefore, MERRA-2 surface temperature and pressure products have important applications in meteorological research and GNSS water vapor monitoring in China.

We select multiple seismic events from seismic-prone regions in the world. Taking Jinghe earthquake as an example, we study the selection method of fitting region for atmospheric error correction based on DEM. Then, to analyze the accuracy and applicability of the correction, we carry out the atmospheric phase delay correction of the co-seismic deformation field based on DEM and GACOS data. The results show that when the phase-elevation correlation coefficient after deducting the main co-seismic deformation region is greater than 0.4, the atmospheric error correction must be made to the obtained deformation field, and the noise can be reduced by 1-4 cm after correction. When the correlation coefficient is 0.4-0.8, the correction accuracy can be verified by the two methods. When the correlation coefficient is greater than 0.8, only DEM fitting atmospheric error distribution function can be used for correction. The atmospheric error correction based on GACOS is universal and can achieve better results in plain areas where atmospheric noise is more serious and is not correlated with terrain.

GAMIT/GLOBK has been used to process the data of 43 GNSS continuous stations in Yunnan since 2018-2021 to obtain the time series. We use Kriging interpolation method to calculate and obtain the strain rate field and analyze the variation of regional strain field before and after the Yangbi M_S6.4 earthquake. The results show that the seismic risk area should be focused on the compression area of surface strain and the conversion area of high and low values of maximum shear strain. The Yangbi M_S6.4 earthquake occurred in the compression area of surface strain in the early stage, and quickly changed to extension state after the earthquake. Yangbi area has been in the high value area of maximum shear strain rate for a long time, and the shear activity is strong. The regional strain accumulation rate of maximum shear strain near the epicenter decreases rapidly in the short-term before the earthquake.

To improve the prediction accuracy of surface subsidence of the GNSS CORS automatic monitoring system in mining areas, we propose a combined prediction method based on the Kalman filter model optimized by genetic algorithm (GA-KF) and the BP neural network strong prediction model (BP-Adaboost) after phase space reconstruction based on the wavelet analysis. The trend sequence and random sequence of the original monitoring data are obtained by wavelet analysis, and are predicted by the GA-KF model and the BP-Adaboost model respectively. The sum of the two data is the final prediction result. Taking the data of Bozhou Banji monitoring station as an example, the experimental results show that: 1) Compared with the prediction values of using the single GA-KF and BP-Adaboost model after phase space reconstruction, the combined model has higher prediction accuracy. 2) It is found that the combined model is less affected by the length of modeling sequence, and the average relative error is within 0.003%, which is much smaller than the two single models and has a certain anti-interference ability.

Using the original seismic waveform data of Tonga submarine volcanic eruption in the South Pacific on January 15, 2022, recorded by 70 broadband stations of the Shandong seismic network, we analyze the recording characteristics of volcanic eruption seismic wave from time and frequency domains. By studying the surface wave with 20 s dominant period, we calculate the surface wave magnitude and seismic wave radiation energy of the volcanic eruption. The results show that the surface wave magnitude of the Tonga volcanic eruption is M_S5.674 and the energy magnitude is M_e5.704, which is basically consistent with the magnitude measured by GFZ and 0.126 smaller than that measured by USGS. The Tonga volcanic explosivity index (VEI) is grade 5, but the seismic wave radiation energy measured in this paper is 8.155×10¹² J. The volcanic eruption releases less energy in the form of seismic wave radiation, and most of the energy is rapidly released into the air through the atmosphere in the form of shock wave with the speed more than 1 000 km/h, inhibiting the compression and extension of deep stress.

Based on the data from Ningxia regional seismic network, we determine the focal depth of the Zhongning M_L3.6 earthquake using the first arrival P phase locating method, and the focal mechanism and moment tensor solution of the earthquake using gCAP inversion method. Then, using the simulation method of focal mechanisms and stress field, we calculate the relative shear stress and normal stress on the two nodal planes of focal mechanism of the Zhongning M_L3.6 earthquake under current stress systems. The results show that the focal depth of the Zhongning M_L3.6 earthquake measured by the first arrival P phase locating method and gCAP method are both 11 km. The focal mechanism of the earthquake is strike 242°, dip 63°, rake 8° for nodal plane I, and strike 148°, dip 83°, rake 153° for nodal plane II. The moment tensor result reveals that double couple component of the earthquake accounts for 97.07% of total moment tensor solution, indicating that the earthquake is a natural earthquake. Combined with spatial distribution of relocation of M_L≥1.0 earthquakes from 2009 near the epicentre, regional geological structure, and relative shear stress, we conclude that the seismogenic fault of the Zhongning M_L3.6 earthquake may be related to the NNW trending weak fault zone near the southeast segment of Tianjingshan fault. Under the NE direction compression of the northeastern Tibetan plateau, the earthquake is a right-lateral strike-slip seismic event caused by shear slip along NNW direction fault.