Analysis of GNSS-Derived Tropospheric Zenith Non-Hydrostatic Delay Anomaly during Sandstorms in Northern China on 15th March 2021

: On the 15th of March 2021, the strongest sandstorm in a decade occurred in northern China, and had a great adverse impact on the natural environment and human health in northern China. Real-time monitoring of dust storms is becoming increasingly important. In order to effectively analyze the non-hydrostatic delay (ZNHD) anomaly during a sandstorm, the method based on GNSS-derived tropospheric ZNHD residual to monitor the sandstorm is proposed at the same time. We studied the relationship between ZNHD/PWV and PM10/PM2.5 in Beijing, Changchun, Pingliang and Zhongwei before and after sandstorms. The ZNHD time series was then decomposed by singular spectrum analysis (SSA) and the residuals were obtained. The relationship between the GNSS-derived ZNHD residual and PM10 was analyzed. The results show that the impact of the sandstorm on PM10 is greater than that on PM2.5. Before the sandstorm, the correlation between PM10 and ZNHD was low, less than 0.25. When the sandstorm occurred, the correlation between PM10 and ZNHD increased signiﬁcantly, and the maximum was greater than 0.7. When the sandstorm ended, the correlation between PM10 and ZNHD decreased signiﬁcantly. Through the relationship between the ZNHD residual and PM10, it can be found that when the peak-to-peak values of the ZNHD residual are all above 80 mm, sandstorms may occur. But Rainfall, snowfall, haze and other abnormal weather can also lead to ZNHD anomalies.


Introduction
A sandstorm, a collective name for sandstorm and dust storm, arises when a gust front or other strong wind blows loose sand and dirt from a dry surface. During the period of a sandstorm, particles are transported by saltation and suspension, causing turbid air and low visibility of less than 1km. A sandstorm will bring a lot of floating dust, among which particulate matter 10 (PM10), (2.5 to 10 microns particulate matter) and particulate matter 2.5 (PM2.5) (2.5 microns or less particulate matter) are two important indicators of sandstorm pollution. PM2.5 and PM10 consist of a variety of toxic and hazardous substances that can be breathed directly into the lungs and cause respiratory diseases, such as asthma and cardiovascular diseases [1,2]. A sandstorm, as a common natural disaster in northern China, jeopardizes the natural environment, national economy and people's health [2][3][4][5][6]. As can be seen from Figure 1, from 14th of March to the morning of the 15th of March 2021, affected by the cold air, one of the worst sand and dust storms from Mongolia in a decade hit the central and western parts of Gansu, Inner Mongolia, northern Shanxi, northern Hebei, Beijing, Tianjin and other places in China, with a large environmental and economic impact and harmed air quality for millions of people. In many places, levels of the coarse particulate matter, PM10, were hazardous [7][8][9][10][11]. In Inner Mongolia, Gansu, Currently, weather stations and remote sensing satellites equipped with MODIS sensors are mainly used to monitor the occurrence and changes of sandstorms [13][14][15][16][17]. However, due to the limited number of weather stations and the uneven spatial distribution, the transit time of remote sensing satellites equipped with MODIS sensors in China is generally from 10:30 to 12:00 during daytime and from 21:30 to 23:30 at night. Therefore, the above methods may fail to realize a real-time monitoring and forecasting of unforeseen meteorological conditions such as sandstorms. In recent years, the rapid development of the Global Navigation Satellite System (GNSS) has provided large-scale, realtime, and high-precision tropospheric products by improving the temporal and spatial resolution, effectively eliminating the shortcomings of traditional detection technology and promoting the development of GNSS in the field of meteorology [18][19][20][21][22][23]. GNSS technique is widely applied in the monitoring of haze, rainstorm, typhoon and other disasters, as well as forest fires, and has achieved improved results in the field of numerical meteorological forecasting [24][25][26][27][28]. For example, Zhao et al. [29] proposed and applied the method of using PWV data retrieved by GPS to forecast precipitation during typhoons. Choy et al. [30] monitored the storms in Melbourne in 2010 based on the PWV data retrieved from GPS and achieved ideal results. Nykiel et al. [31] successfully monitored the location and development path of the storm by using the GNSS observation network. These advantages of the GNSS technique exactly make up for the insufficiency of using meteorological stations and remote sensing methods to monitor sand and dust storms. During the outbreak of sandstorms, the rapid increase in the amount of particles in the air will inevitably cause changes in the tropospheric delay of GNSS, which makes it possible for GNSS to monitor sandstorms. Therefore, the important practical significance can be accomplished by in-depth analysis of the relationship between sandstorms and tropospheric delay, and research on methods of using the GNSS technique to monitor sandstorms. Currently, weather stations and remote sensing satellites equipped with MODIS sensors are mainly used to monitor the occurrence and changes of sandstorms [13][14][15][16][17]. However, due to the limited number of weather stations and the uneven spatial distribution, the transit time of remote sensing satellites equipped with MODIS sensors in China is generally from 10:30 to 12:00 during daytime and from 21:30 to 23:30 at night. Therefore, the above methods may fail to realize a real-time monitoring and forecasting of unforeseen meteorological conditions such as sandstorms. In recent years, the rapid development of the Global Navigation Satellite System (GNSS) has provided large-scale, real-time, and high-precision tropospheric products by improving the temporal and spatial resolution, effectively eliminating the shortcomings of traditional detection technology and promoting the development of GNSS in the field of meteorology [18][19][20][21][22][23]. GNSS technique is widely applied in the monitoring of haze, rainstorm, typhoon and other disasters, as well as forest fires, and has achieved improved results in the field of numerical meteorological forecasting [24][25][26][27][28]. For example, Zhao et al. [29] proposed and applied the method of using PWV data retrieved by GPS to forecast precipitation during typhoons. Choy et al. [30] monitored the storms in Melbourne in 2010 based on the PWV data retrieved from GPS and achieved ideal results. Nykiel et al. [31] successfully monitored the location and development path of the storm by using the GNSS observation network. These advantages of the GNSS technique exactly make up for the insufficiency of using meteorological stations and remote sensing methods to monitor sand and dust storms. During the outbreak of sandstorms, the rapid increase in the amount of particles in the air will inevitably cause changes in the tropospheric delay of GNSS, which makes it possible for GNSS to monitor sandstorms. Therefore, the important practical significance can be accomplished by indepth analysis of the relationship between sandstorms and tropospheric delay, and research on methods of using the GNSS technique to monitor sandstorms.
The zenith total delay (ZTD) estimated by GNSS is the sum of zenith hydrostatic delay (ZHD) and zenith wet delay (ZWD) (also called zenith non-hydrostatic delay (ZNHD) in surveying). The ZHD can be accurately calculated by modeling. The ZNHD needs to be estimated as an unknown parameter during data processing, which will cause the ZNHD to include the delay caused by both water vapor and particulate matters [27,28]. Therefore, the method of ZNHD is proposed to monitor the occurrence of sandstorms. Firstly, the singular spectrum analysis (SSA) method is used to decompose the time series of the ZNHD in which the principal components are subtracted, and the remaining components are used as residuals. Finally, the method of using the ZNHD is researched to monitor sandstorms. Section 2 describes the different datasets, with their characteristics, and the methodology used in this work. Processing results and analysis are presented in Section 3. Finally, the conclusion is given in Section 4.

Materials
In order to study the feasibility of using the GNSS technique to monitor the sandstorm weather in northern China in March 2021, we conducted research on GPS data of BJFS, GSPL, NXZW and CHAN stations from 20th February to 20th March 2021. The data sampling rate was 30 s. The GPS data of BJFS, CHAN and GSPL, NXZW stations was provided by the International GNSS Service (IGS) data center of Wuhan University (igs.gnsswhu.cn) and The China Earthquake Administration [32]. The BJFS, CHAN, GSPL and NXZW stations were located in Beijing, Changchun in Jilin province, Pingliang in Gansu province and Zhongwei in Ningxia Hui Autonomous Region, respectively, which were among the areas severely affected by sandstorms in March 2021.
The GAMIT/GLOBK 10.7 was used to estimate the precipitable water vapor (PWV) and the ZNHD. The Saastamoinen model [33] was used to accurately estimate the ZHD. The ZNHD and PWV were estimated once every hour. The ZNHD uncertainty of 37-day data at four stations output by GAMIT/GLOBK 10.7 was statistically analyzed, and the results are listed in Table 1. It can be seen from Table 1 that the mean of the ZNHD uncertainty of the four stations was between 4.28 mm and 4.81 mm. We took 4.45 mm (mean of the ZNHD uncertainty of four stations) as the uncertainty of the ZNHD. The uncertainty of the ZTD product provided by the International GNSS Service (IGS) was 5 mm [34]. Thus, according to the law of error propagation, the precision of ZHD is about 2.27 mm. The equation of error propagation law is where m ZHD , m ZTD and m ZNHD are the uncertainty of the ZHD, uncertainty of the ZTD and uncertainty of the ZNHD, respectively. By virtue of all data of PM10 and PM2.5 provided by the air quality online monitoring platform with 1 h time resolution of air quality, we studied PM10 data from Langfang Weather station in Hebei, PM2.5 data from Fangshan Meteorological station in Beijing, and PM10/PM2.5 data from Jingyuetan Meteorological station in Changchun, Pingliang and Zhongwei from the 20th of February to the 20th of March 2021.

PWV Inversion Based on GNSS Data
GAMIT/GLOBK 10.7 was used to estimate ZNHD/PWV. When processing data with GAMIT/GLOBK 10.7, the parameters were set as follows: the cutoff height angle was set to 10; the epoch interval was set to 30 s. The tropospheric delay model and tropospheric mapping function were selected as Saastamoinen model and VMF1 model, respectively. Tropospheric parameters were estimated every hour; the FES2004 model and BERNE model were used for the ocean tidal correction model and the solar photo-pressure perturbation model, respectively. Firstly, GNSS data were processed to get the ZTD as an unknown parameter. The tropospheric model was used to calculate the ZHD which was then subtracted from ZTD to obtain ZNHD [35]. Finally, the empirical equation was used to get PWV inversion based on the GNSS data. The GNSS phase observation equation is [36,37] L where i and f denote the receiver and signal frequency, respectively, j denotes the satellite, ρ j i is the geometry distance from receiver to satellite, c is the speed of light in vacuum, dt i is the receiver clock offset, dt j is the satellite clock offset, and M is the troposphere mapping function.
is the phase ambiguity that is not an integer since it includes the receiver and satellite initial phase biases. ε j i, f contains other unmodelled errors for phase observation. ZHD was calculated using the Saastamoinen model [33]. The equation is where P, L and H are the GNSS station atmospheric pressure, latitude and the GPS station geodetic height in the unit of metre, respectively. Finally, ZNHD can be calculated by The PWV can be obtained from the empirical Equation (5) [18,38]: where Π is the conversion factor, g s (gas constant) = 461 J·kg −1 K −1 ; ρ w is the density of liquid water; k 2 , k 3 (refraction constant), k 2 = 16.48 K/hPa, k 3 = (3.776 ± 0.014) * 105 K 2 /hPa; and T m is the weighted average temperature of the troposphere.
where, e is the water vapor pressure and T is the atmospheric temperature (K). The PWV in this paper is calculated by GAMIT/GLOBK 10.7, in which the water vapor pressure (e) and atmospheric temperature (T) were from the global pressure and temperature 2 wet (GPT2w) model. The T m at 4 GNSS stations is shown in the Figure 2.

Singular Spectrum Analysis
As a digital signal processing technique, the SSA [39] has been applied in fields such as climatology, measurement [40], and oceanography [41][42][43]. SSA is an analysis and prediction method for nonlinear and nonstationary time series. It can extract valid information including slowly changing trends, oscillating components, noise, etc., from time series with unknown prior information. Remote Sens. 2022, 14, x FOR PEER REVIEW 5 of 15 Figure 2. The at four GNSS stations.

Singular Spectrum Analysis
As a digital signal processing technique, the SSA [39] has been applied in fields such as climatology, measurement [40], and oceanography [41][42][43]. SSA is an analysis and prediction method for nonlinear and nonstationary time series. It can extract valid information including slowly changing trends, oscillating components, noise, etc., from time series with unknown prior information.
Firstly, building the time-delay matrix, the ZNHD time series is = [ 1 , 2 , 3 , … , ]. The window size is generally greater than 1 and less than N/2, and a suitable M is selected to structure the trajectory matrix X of × ( -+ 1). Let = -+ 1. The trajectory matrix is expressed as follows: Secondly, singular values decomposition (SVD), there was a matrix S = XX T , 1 , … , the eigenvalues of arranged in descending order of magnitude as 1 ≥ 2 ≥ ⋯ ≥ ≥ 0, and the corresponding eigenvectors are 1 , 2 , ⋯ , respectively. The matrices have rank 1; such matrices are sometimes called elementary matrices. The SVD of the time-delay matrix can be expressed as: where = ( ), = √ , and = X T /√ ( = 1,2, ⋯ ). √ is a singular value of . The new matrix ( ) ( < ) can be constructed by selecting appropriate eigenvalues and corresponding eigenvectors. The collection (√ , , ) will be called i-th eigentriple (abbreviated as ET) of the SVD.
Thirdly, eigentriple grouping. The grouping procedure partitions the set of indices {1, … , } into disjoint subsets 1, … , . Let = { 3 , … , }. Then the resultant matrix Firstly, building the time-delay matrix, the ZNHD time series is X i = [x 1 , x 2 , x 3 , . . . , x N ]. The window size M is generally greater than 1 and less than N/2, and a suitable M is selected to structure the trajectory matrix X of M × (N -M + 1). Let K = N -M + 1. The trajectory matrix X is expressed as follows: Secondly, singular values decomposition (SVD), there was a matrix S = XX T , λ 1 , . . . , λ M the eigenvalues of S arranged in descending order of magnitude as λ 1 ≥ λ 2 ≥ · · · ≥ λ M ≥ 0, and the corresponding eigenvectors are U 1 , U 2 , · · · U M , respectively. The matrices X i have rank 1; such matrices are sometimes called elementary matrices. The SVD of the time-delay matrix X can be expressed as: √ λ i is a singular value of X i . The new matrix X (r) (r < d) can be constructed by selecting appropriate eigenvalues and corresponding eigenvectors. The collection ( √ λ i , U i , V i ) will be called i-th eigentriple (abbreviated as ET) of the SVD.
Thirdly, eigentriple grouping. The grouping procedure partitions the set of indices {1, . . . , d} into m disjoint subsets I1, . . . , Im. Let I = {i 3 , . . . , i p }. Then the resultant matrix X I corresponding to group I is defined as X I = X I1 + . . . + X I p . The resultant matrices are computed for the groups I = I1, . . . , Im and the expansion (7) leads to the decomposition The procedure of choosing the sets I1, . . . , Im is called eigentriple grouping. If m = d and I j = { j }, j = 1, . . . , d, then the corresponding grouping is called elementary.
Fourthly, diagonal averaging. The purpose of diagonal averaging was to transform the decomposed elementary matrix X i into a new time series of length N, called the reconstruction component (RC), and the sum of all RCs was equal to the raw ZNHD time series x 1 , x 2 , x 3 , . . . , x n . There is a matrix Z = X i , and z 1 , z 2 , · · · , z N is the time series obtained by diagonal averaging of Z. Then the formula for diagonal averaging can be expressed as: where The various signals in the ZNHD time series after preprocessing are separated by SSA. The w-correlation method [44] was used to analyze the correlation between each reconstructed component. The reconstructed component with the same signal features was grouped, and the higher correlation part of the residual sequence was reconstructed. The reconstructed time series of each elementary matrix X k was defined as Y k which consists of elements y k 1 , y k 2 , · · · , y k N . The correlation coefficients between each reconstructed time series can be expressed as: where

Relationship between ZNHD/PWV and PM2.5/PM10
It can be seen from Equations (6) and (7) that during the conversion from ZNHD(ZWD) to PWV, T m is a key coefficient, which changes with the water vapor pressure and atmospheric temperature in each region. Therefore, theoretically, ZNHD(ZWD) and PWV are different. Therefore, we analyzed the relationship between ZNHD/PWV and PM2.5/PM10, respectively.
Firstly, the change trends and correlations of ZNHD, PWV, PM10 and PM2.5 were analyzed before, during, and after the occurrence of sandstorms at the BJFS, GSPL, NXZW and CHAN stations. The results are shown in Figures 3-6 and Table 2.
As can be seen from Figure 3, the concentrations of PM10 in Beijing was greater than 2300 µg/m 3 , about 11.2 times higher than usual. However, the change in PM2.5 was small, and the concentration of PM2.5 in Beijing reached approximately 330 µg/m 3 when the sandstorm occurred. Between the 20th of February and the 14th of March the trends of ZNHD/PWV and PM10/2.5 at the BJFS station were relatively consistent; when the sandstorm occurred from the 15th of March to the 16th of March, the change consistency of ZNHD/PWV and PM10/2.5 reached an extremely high level; after the sandstorm from the 16th of March to the 20th of March the change of ZNHD/PWV and PM10/2.5 remained high. During the period from the 25th of March to the 27th of March ZNHD and PWV at BJFS station showed an obviously increasing trend. We know from the website of the China Meteorological Administration that there was obvious haze in Beijing during this period, which can also be reflected by PM2.5 data.
As can be seen from Figure 4, the concentrations of PM10 in Changchun are greater than 800 µg/m 3 , about seven times higher than usual. However, the change of PM2.5 in Changchun was not obvious, demonstrating that sandstorms mainly transfer large particles, with a greater impact on PM10, but less impact on PM2. 5 As can be seen from Figure 3, the concentrations of PM10 in Beijing was greater than 2300 μg/m 3 , about 11.2 times higher than usual. However, the change in PM2.5 was small, and the concentration of PM2.5 in Beijing reached approximately 330 μg/m 3 when the sandstorm occurred. Between the 20th of February and the 14th of March the trends of ZNHD/PWV and PM10/2.5 at the BJFS station were relatively consistent; when the sandstorm occurred from the 15th of March to the 16th of March, the change consistency of ZNHD/PWV and PM10/2.5 reached an extremely high level; after the sandstorm from the 16th of March to the 20th of March the change of ZNHD/PWV and PM10/2.5 remained high. During the period from the 25th of March to the 27th of March ZNHD and PWV at BJFS station showed an obviously increasing trend. We know from the website of the China Meteorological Administration that there was obvious haze in Beijing during this period, which can also be reflected by PM2.5 data.  As can be seen from Figure 3, the concentrations of PM10 in Beijing was greater than 2300 μg/m 3 , about 11.2 times higher than usual. However, the change in PM2. 5    with PM2.5 failed to come up to a high level. From the 16th of March to the 20th of March after the sandstorms, the trends of ZNHD and PM10/2.5 had a relatively high consistency During the period from the 26th of March to the 28th of March, ZNHD and PWV at CHAN station showed an obviously increasing trend. According to the website of the China Me teorological Administration, there was an obvious rainfall in Changchun during this pe riod.    As can be seen from Figure  In summary, the occurrence of sandstorms can be detected by using the GNSS-derived ZNHD/PWV. However, GNSS-derived ZNHD/PWV is also affected by abnormal weather such as rainfall, snowfall, and haze. Therefore, when using the GNSS-ZNHD to detect sandstorms, it is necessary to consider the impact of rain, snow and haze.
From the Figures 3-6, sandstorms have a greater impact on PM10 than PM2.5, and the changes in ZNHD and PWV are highly consistent. Therefore, in order to discover the relationship more intuitively between ZNHD/PWV and PM10, we examined statistics of the correlation between ZNHD and PM10 before, during and after the occurrence of sandstorms. At the same time, we tested the significance level of α = 0.05 for the correlation coefficient. Except for the correlation coefficient between PM10 and ZNHD of BJFS station before sandstorm and the correlation coefficient between PM10 and ZNHD of Chan station after the sandstorm, all other correlation coefficients have passed the significance test. The statistical results are shown in Table 2.
According to Table 2, before sandstorms, the correlation between ZNHD and PM10 at four GNSS stations were low. This is mainly because the tropospheric delay caused by PM10 has little contribution to ZNHD, which is mainly caused by water vapor. The correlation between ZNHD and PM10 were increasing in significance when sandstorms occurred, for which the proportion of tropospheric delay caused by PM10 in ZNHD was obviously larger than that before the sandstorms. After the sandstorm, the correlation between PM10 and ZNHD decreased significantly, and the reason was similar to that before the sandstorm.

Monitoring Sandstorms Based on the ZNHD Residual
In order to use the GNSS to monitor sandstorms more effectively and accurately, we have proposed a method of monitoring sandstorms based on the ZNHD residual.
Firstly, the ZNHD time series was decomposed by SSA and the order of principal components was judged by the W-correlation method and eigenvalue ratio. The original series subtracts principal components to obtain the ZNHD residual. The results of Wcorrelation analysis of different orders are shown in Figure 7. It can be seen from Figure 7 that the first nine orders of ZNHD at BJFS station and CHAN station contain five signals with good frequency resolution, which can be effectively separated by SSA. The components after the 9th order have a poor separation effect. The eigenvalue contribution rates of the first nine orders of ZNHD decomposition results at BJFS station, CHAN station, GSPL station and NXZW station reach 98.2%, 98.1%, 98.9% and 99.1%, respectively, which shows that the first nine orders already contain most ZNHD signals. Therefore, we take the first nine orders as the principal components and the remaining components as the residuals. tween PM10 and ZNHD decreased significantly, and the reason was similar to that before the sandstorm.

Monitoring Sandstorms Based on the ZNHD Residual
In order to use the GNSS to monitor sandstorms more effectively and accurately, we have proposed a method of monitoring sandstorms based on the ZNHD residual.
Firstly, the ZNHD time series was decomposed by SSA and the order of principal components was judged by the W-correlation method and eigenvalue ratio. The original series subtracts principal components to obtain the ZNHD residual. The results of W-correlation analysis of different orders are shown in Figure 7. It can be seen from Figure 7 that the first nine orders of ZNHD at BJFS station and CHAN station contain five signals with good frequency resolution, which can be effectively separated by SSA. The components after the 9th order have a poor separation effect. The eigenvalue contribution rates of the first nine orders of ZNHD decomposition results at BJFS station, CHAN station, GSPL station and NXZW station reach 98.2%, 98.1%, 98.9% and 99.1%, respectively, which shows that the first nine orders already contain most ZNHD signals. Therefore, we take the first nine orders as the principal components and the remaining components as the residuals.   from about 9 am on the 15th of March and the concentration of PM10 returned to normal at about 5 am on the 16th of March. About 22:00 on the 14th of March, the ZNHD residuals began to show abnormalities. During the period from the 25th of March to the 27th of March, the ZNHD residuals at BJFS station showed an obvious increasing trend. We know from the website of the China Meteorological Administration that there was obvious haze in Beijing during this period. According to Figure 9, the PM10 concentration in Changchun increased significantly from 12 pm on the 15th of March and the PM10 concentration returned to normal around 2 am on the 20th of March. The ZNHD residual showed obvious abnormalities from the early morning of the 15th of March to the noon of the 19th of March. Judging from the ZNHD residuals in Beijing and Changchun, the abnormal occurrence of ZNHD residuals was ahead of schedule. From the 23rd of March to the 24th of March, PM10 showed minor abnormalities, with PM10 content below 200, but there was no obvious abnormality in the ZNHD residual. This may be because the PM10 content was too low, and the impact on the GNSS signal is submerged in noise.  According to Figures 10 and 11, the variation of ZNHD residuals in GSPL and NXZW stations is very similar to that in BJFS and CHAN stations. When the sandstorm occurs, the fluctuation of ZNHD residuals increases significantly, which can reach about ten times normal. The difference is that the ZNHD residuals anomaly of GSPL station appears before the PM10 anomaly, while the ZNHD residuals anomaly of NXZW station appears after the PM10 anomaly. The PM10 content in Pingliang City and Zhongwei city was abnormal on the 28th of February and the 27th of February, respectively, but there was no obvious abnormality in ZNHD residuals. This may be because the GNSS technology detects the particle anomaly on the signal path of the zenith direction of the observation station, and the meteorological observation station detects the particle anomaly on the ground. The wind direction may also be a factor that affects the difference in the time of anomaly detection by different technologies. The specific impact mechanism needs further research in the future.
In summary, the ZNHD and the PWV estimated by GNSS have obvious correlation with PM10. When sandstorms occur, the correlation increases in an obvious manner. Through the relationship between the ZNHD residual and PM10, it can be found that the peak-to-peak value of ZNHD residuals at the four GNSS stations exceeded 80 mm during  According to Figures 10 and 11, the variation of ZNHD residuals in GSPL and NXZW stations is very similar to that in BJFS and CHAN stations. When the sandstorm occurs, the fluctuation of ZNHD residuals increases significantly, which can reach about ten times normal. The difference is that the ZNHD residuals anomaly of GSPL station appears before the PM10 anomaly, while the ZNHD residuals anomaly of NXZW station appears after the PM10 anomaly. The PM10 content in Pingliang City and Zhongwei city was abnormal on the 28th of February and the 27th of February, respectively, but there was no obvious abnormality in ZNHD residuals. This may be because the GNSS technology detects the particle anomaly on the signal path of the zenith direction of the observation station, and the meteorological observation station detects the particle anomaly on the ground. The wind direction may also be a factor that affects the difference in the time of anomaly detection by different technologies. The specific impact mechanism needs further research in the future.
In summary, the ZNHD and the PWV estimated by GNSS have obvious correlation with PM10. When sandstorms occur, the correlation increases in an obvious manner. Through the relationship between the ZNHD residual and PM10, it can be found that the peak-to-peak value of ZNHD residuals at the four GNSS stations exceeded 80 mm during the sandstorm. During the sandstorm, ZNHD and PWV will be obviously abnormal. This According to Figure 8, the concentration of PM10 in Beijing increased significantly from about 9 am on the 15th of March and the concentration of PM10 returned to normal at about 5 am on the 16th of March. About 22:00 on the 14th of March, the ZNHD residuals began to show abnormalities. During the period from the 25th of March to the 27th of March, the ZNHD residuals at BJFS station showed an obvious increasing trend. We know from the website of the China Meteorological Administration that there was obvious haze in Beijing during this period.
According to Figure 9, the PM10 concentration in Changchun increased significantly from 12 pm on the 15th of March and the PM10 concentration returned to normal around 2 am on the 20th of March. The ZNHD residual showed obvious abnormalities from the early morning of the 15th of March to the noon of the 19th of March. Judging from the ZNHD residuals in Beijing and Changchun, the abnormal occurrence of ZNHD residuals was ahead of schedule. From the 23rd of March to the 24th of March, PM10 showed minor abnormalities, with PM10 content below 200, but there was no obvious abnormality in the ZNHD residual. This may be because the PM10 content was too low, and the impact on the GNSS signal is submerged in noise.
According to Figures 10 and 11, the variation of ZNHD residuals in GSPL and NXZW stations is very similar to that in BJFS and CHAN stations. When the sandstorm occurs, the fluctuation of ZNHD residuals increases significantly, which can reach about ten times normal. The difference is that the ZNHD residuals anomaly of GSPL station appears before the PM10 anomaly, while the ZNHD residuals anomaly of NXZW station appears after the PM10 anomaly. The PM10 content in Pingliang City and Zhongwei city was abnormal on the 28th of February and the 27th of February, respectively, but there was no obvious abnormality in ZNHD residuals. This may be because the GNSS technology detects the particle anomaly on the signal path of the zenith direction of the observation station, and the meteorological observation station detects the particle anomaly on the ground. The wind direction may also be a factor that affects the difference in the time of anomaly detection by different technologies. The specific impact mechanism needs further research in the future.
In summary, the ZNHD and the PWV estimated by GNSS have obvious correlation with PM10. When sandstorms occur, the correlation increases in an obvious manner. Through the relationship between the ZNHD residual and PM10, it can be found that the peak-to-peak value of ZNHD residuals at the four GNSS stations exceeded 80 mm during the sandstorm. During the sandstorm, ZNHD and PWV will be obviously abnormal. This paper puts forward a feasible and effective new method for monitoring sandstorm weather.

Conclusions
In order to analyze the abnormal response of the ZNHD during a sandstorm, the relationship between ZNHD/PWV and PM10/PM2.5 before and after a sandstorm was analyzed. The method with the GNSS-derived ZNHD residual was proposed. The results show that the variation of PM10 is more dramatic than that of PM2.5 when a sandstorm occurs. Before sandstorms, the correlation between ZNHD and PM10 is less than 0.25. The correlation between ZNHD and PM10 increased significantly when sandstorms occurred, and the correlation coefficient can reach greater than 0.7. After sandstorms, the correlation between ZNHD and PM10 decreased significantly. Through the study of this sandstorm, it was found that the peak-to-peak value of the ZNHD residuals at the four GNSS stations exceeded 80 mm during the sandstorm. However, GNSS-derived ZNHD/PWV was also affected by abnormal weather such as rainfall, snowfall, and haze.