Estimation of Surface Water Level in Coal Mining Subsidence Area with GNSS RTK and GNSS-IR

Abstract

Ground subsidence caused by underground coalmining result in the formation of ponding water on the ground surface. Monitoring the surface water level is crucial for studying the hydrologic cycle in mining areas. In this paper, we propose a combined technique using Global Navigation Satellite System Real-Time Kinematic (GNSS RTK) and GNSS Interferometric Reflectometry (GNSS-IR) to estimate the surface water level in areas of ground subsidence caused by underground coal mining. GNSS RTK is used to measure the geodetic height of the GNSS antenna, which is then converted into the normal height using the local height anomaly model. GNSS-IR is employed to estimate the height from the water surface to the GNSS antenna (or, the reflector height). To enhance the accuracy of the reflector height estimation, a weighted average model has been developed. This model is based on the coefficient of determination of the signal fitted by the Lomb-Scargle spectrogram and can be utilized to combine the reflector height estimations derived from multiple GNSS system and band reflection signals. By subtracting the GNSS-IR reflector height from the GNSS RTK-based normal height, the proposed method-based surface water level estimation can be obtained. In an experimental campaign, a low-cost GNSS receiver was utilized for the collection of dual-frequency observations over a period of 60 days. The collected GNSS observations were used to test the method presented in this paper. The experimental campaign demonstrates a good agreement between the surface water level estimations derived from the method presented in this paper and the reference observations.

1. Introduction

Coal has played a significant role in global energy production and industrial development, currently accounting for 40% of global primary energy consumption [1]. Due to the limitation of geological conditions, surface mining is impractical or uneconomical in most circumstances. Hence, the majority of coal production comes from underground mining especially in China [2]. The underground coal-mining can induce ground subsidence over a large area. When the subsided ground falls below the water table level, groundwater can rise to the ground surface, submerging farmland and buildings. Continuous monitoring of surface water levels in areas of ground subsidence caused by underground coal-mining is important for scientific research on the hydrologic cycle in the mining area and for the effective reclamation of submerged land [3,4,5].

Basically, monitoring surface water levels with high temporal resolution based on the conventional leveling surveying method is unfeasible, because the method is labor-intensive and costly. The ground subsidence due to the underground coal-mining can induce considerable spatial position variation for the water level radar fixed on the ground, resulting in a significant error in the water level radar-based measurement. Consequently, the water level radar is basically incapable of measuring the surface water level (and surface water level variation) in coal mining subsidence area. The invalidation of the water level measurement using the water level radar in the coal mining subsidence area is also clearly illustrated in Figure 1. As shown in the figure, at a given time t₁, the vertical distance from the water level radar to the reference datum is H₀, and the radar-measured vertical distance to the water surface is H₁; then the radar-based water level can be easily obtained, which is H_L1 = H₀ − H₁. At the time t₂, the radar-based water level is H_L2 = H₀ − H₂ for the areas of subsidence-free, if the radar-measured distance to the water surface is H₂. However, due to the ground subsidence induced by the underground coal-mining (H_s1 in Figure 1), the in-situ water level at the time of t₂ is H₀ − H_s1 − H₂, instead of H₀ − H₂, in the underground coal-mining induced ground subsidence areas. Hence, accurately and continuously measuring the surface water level in the area of ground subsidence due to underground coal-mining in a cost-effective manner remains a challenging problem.

GNSS RTK (Global Navigation Satellite System Real-Time-Kinematic, GNSS RTK) is a typical positioning technique, which makes use of a combination of the GNSS signal collected by the rover station and the signal collected by the base station to provide accurately position information of the rover station. The technology began to be applied for monitoring ground subsidence due to underground coal mining around 2000 [6]. By fixing the rover to a ground point, GNSS RTK can be used to monitor ground subsidence due to underground coal mining at that point with high temporal resolution. However, because of the high cost of the GNSS receiver and the establishment of the GNSS station, monitoring ground surface subsidence in this manner is still rare to this day. In addition to GNSS RTK, the GNSS Precise Point Positioning (PPP) technique can also be used to monitor the ground subsidence deformation. Basically, the GNSS PPP technique is inferior to GNSS RTK in terms of timeliness, stability, and accuracy; although the technique can reduce the economic cost to some extent [7,8,9]. In recent years, some GNSS devices integrated with solar power modules and wireless network communication modules have been developed. By making use of the integrated GNSS devices and GNSS RTK (or GNSS PPP), automatic and real-time monitoring of ground subsidence deformation in underground coal mining areas has been achieved [10,11,12].

In addition to positioning, GNSS can also be used for remotely sensing a variety of Earth’s surface parameters. The GNSS remote sensing technique is also referred to as GNSS reflectometry [13]. GNSS reflectometry can be divided into two types; that is the reflecting pattern and the interference pattern [14]. The reflecting pattern technique uses a specialized receiver connected with two GNSS antennas. One of GNSS antennas is right-hand circularly polarized and is typically deployed with a zenith-looking orientation to collect the direct GNSS signal. While the other is left-hand circularly polarized and deployed with a nadir-looking orientation to collect the signal reflected from the ground or water surface. Because of the high cost and high energy consumption of the receiving instrument, it is common to install the receiver on a satellite for measuring ocean surface parameters (e.g., ocean surface roughness and wind speed) and ground surface parameters (e.g., vegetation height and soil moisture) [15,16,17,18]. The interference pattern technique uses a geodetic or navigational GNSS receiver equipped with a single antenna to collect the interfered signal. The interfered signal consists of the direct signal and the signal reflected from a medium surface, such as the ground or water surface. By utilizing the frequency variation (or the oscillating period variation) of the interfered signal, the technique can be used to retrieve the height variation of reflection surface relative to the GNSS antenna. The technique is also termed GNSS Interferometric Reflectometry (GNSS-IR) and can be utilized to measure the snow depth, lake and coastal water level variation [19,20,21,22].

Previous studies related to GNSS-based water level estimation mainly focus on circumstances where only the water level varies, such as river or ocean water level variation estimation. In these cases, using the reflected GNSS signal alone (or the GNSS-IR technique) can accurately measure water level variations. For example, in the last few years, Larson et al. proposed using the GNSS-IR method to estimate ocean water levels and developed a model to correct the estimation errors caused by dynamic changes in water levels [23,24,25]. Jin et al. also demonstrated the good performance of BDS signals in ocean water level estimation [26]. Wang et al. utilized multiple GNSS systems and a signal fusion method to achieve millimeter-grade accuracy in estimating river water levels [27,28]. Hitherto, the GNSS-based water level estimation technique has not been utilized to measure water levels in areas affected by underground coal mining-induced ground subsidence. In the context of ground subsidence due to underground coal mining, both the water level and ground surface height vary simultaneously. The alternative method for estimating surface water levels in coal mining subsidence areas is the leveling survey method, which is labor-intensive, time-consuming, and expensive. Therefore, based on GNSS techniques, the development and validation of a new water level estimation method for ground subsidence areas is necessary to improve measuring efficiency and reduce economic costs.

In this paper, a new method for estimating the surface water level in area of ground subsidence due to underground coal-mining is presented. The method combines the GNSS positioning technique (GNSS RTK, specifically) with GNSS-IR, which may also be termed “GNSS+IR”, enabling continuous and accurate estimation of the surface water level in the area of ground subsidence. Additionally, due to the use of a low-cost GNSS chip and antenna, the method presented in this paper is promising for cost-effective monitoring of surface water level variations with high temporal and spatial resolution over a large area.

2. Fundamentals of Ground Subsidence Due to Underground Coal-Mining

When the underground coal seam is not mined out, it carries the pressures from the rocks and soil above it. The pressures will cause the rocks to separate from the bedding planes above the coal seam once it has been fully mined. As the mined-out area (or goaf) increases, the mass of the separated rocks also rises, resulting in increased gravitational force. When this force exceeds the cohesion of the separated rock seam, the rocks will fall into the goaf [29]. As the separated rocks fall into the goaf, the rocks and soil above will experience vertical displacement. After the effects of the displacement reach the surface, ground subsidence due to underground coal mining begins. Basically, the shape of the subsided ground surface due to the underground coal-mining is very similar to a bowl as shown in Figure 2, and is also termed as a subsidence basin. The ground subsidence area due to underground mining in a typical coalfield in China can exceed about 200 km², resulting in serious damage to buildings and land, as well as waterlogging over a larger area, as illustrated in Figure 2.

In general, the subsidence basin area is significantly greater than that of the goaf; subsidence magnitudes are also different for different ground points in the basin. Both the subsidence area and magnitude depend on the mining and geological conditions (e.g., mining speed, coal seam depth and mining thickness); and the ground subsidence is also spatially and temporally related [30,31]. In the spatial dimension, for the ground point within the subsidence basin, the magnitude of subsidence at that point is inversely proportional to the horizontal distance from the goaf center [32]. In other words, as the horizontal distance increases, the subsidence magnitude at the ground point decreases.

In the temporal dimension, ground subsidence due to underground coal-mining develops progressively. At a ground point, a slight downward movement becomes detectable when the mining operation is located approximately 100 to 200 m behind that point [33]. As the mining position passes through the point, the ground subsidence increases significantly. Approximately 80% of the total ground subsidence occurs during this stage (or, termed as “active stage”). Duration of the active stage is typically about 3–10 months, as shown in Figure 3, which shows a typical ground subsidence time series induced by mining of a thick coal seam. It can be observed from Figure 3 that underground coal mining will induce rather larger ground subsidence over a short period (i.e., greater than 5 m within about 6 months). As a result, both the horizontal and vertical positions of the GNSS antenna fixed on the subsided ground will vary considerably.

After the active stage, the subsidence velocity of the ground point gradually closes to 0 mm/day, and the ground subsidence magnitude at the point gradually approaches its maximum. This stage is also referred to as the ‘stable stage’. The duration of the stable stage can last for several years [34].

3. Materials

The experimental campaign was conducted in the Jining coal mining area, which is one of the largest coal-producing areas in China, from 20 August 2023 to 20 October 2023, over about 60 days. The extent of the geographical region of Jining coal mining area is shown in Figure 4. A GNSS station was established at the ground surface above a coal mining working face on 15 August 2023. The relative position between the GNSS station and the coal mining working face is shown in the left part of Figure 5; the right part of Figure 5 shows local environment of the GNSS station. The working face was mined on 17 June 2022. The mining width and length of the working face are about 140 m and 800 m respectively; and the averaged coal seam thickness, mining depth, and mining speed of the working face are about 6.5 m, 550 m, and 2.8 m/d, respectively.

Due to the significant ground subsidence induced by the underground coal mining, the raw ground surface and roads around the GNSS station were submerged by the groundwater after the underground coal-mining, as illustrated in Figure 6. In the experiment, the GNSS station is established on the side of an elevated and repaired road inside the underground coal-mining-induced water pooling area; that is, the GNSS station is deployed inside the water pooling area, as shown in Figure 6.

In addition, a high-frequency water level radar (or, water level gauge, WLG) was also used to measure height from the water surface to WLG (i.e., h_w). The deployment of the WLG is shown in the right part of Figure 6. As shown in the figure, the WLG was fixed on a horizontal crossbeam. The crossbeam is horizontally connected to the GNSS station pole. Besides, a metal tube attached to a metal ruler was vertically fixed in the water pooling area for measuring the surface water level with the leveling surveying method, as shown in Figure 7.

4. Methods

4.1. Geometry of GNSS+IR Based Surface Water Level Estimation

As presented in Section 2, underground coal-mining can induce a considerably subsidence on the ground surface. Once the height of the subsided ground is lower than the water table level, groundwater will be exposed to the air and converted into surface water as illustrated in Figure 8. The geometry of the proposed method for estimating surface water levels in the subsidence area due to underground coal-mining is illustrated in the right figure of Figure 8. A GNSS antenna and receiver is installed on the ground in an area affected by subsidence due to underground coal mining. For a given time, plane coordinates (x_a, y_a) and geodetic height (H_a) of the antenna can be measured by using GNSS RTK. Since the normal height, rather than the geodetic height, is more applicable for realistic engineering applications (such as land reclamation and road repair in coal mining subsidence areas), it is necessary to convert the GNSS RTK-based geodetic height into the normal height for the surface water level estimation in mining subsidence areas.

As studied in [35], the difference between geodetic height and normal height (or, height anomaly) at a local region can be described by a second-order surface function:

$$\xi \left(x,y\right)=b_0+b_{1}x+b_{2}y+b_{3}xy+b_4{x}^2+b_5{y}^2$$

(1)

where {b_i} is the function coefficients, which can be obtained by applying the least-squares method to the field dataset of geodetic and normal height observations at the public points. Substituting the GNSS RTK observed plane coordinate into (1), the height anomaly ξ_a at the station can be obtained; then, the normal height of the GNSS antenna Han can be calculated by:

$$H_{an}=H_a-{\xi }_a=H_a-\left(b_0+b_1{x}_a+b_2{y}_a+b_3{x}_a{y}_a+b_4{x}_a^2+b_5{y}_a^2\right)$$

(2)

Further, if the reflector height (RH) of the GNSS antenna (i.e., the vertical distance from the water surface to the antenna) H_r can be estimated using GNSS-IR, based on the geometry relationship shown in Figure 8, the surface water level in coal mining subsidence area at a given time can be obtained by:

$$H_g=H_{an}-H_r$$

(3)

Note that, in the existing research scenarios (e.g., river, lake, or ocean water level estimation), the normal height of the GNSS antenna can basically be treated as a constant, due to the negligible position variation of the GNSS antenna [21,36]. Thus, the GNSS positioning technique is usually unnecessary for the water level (or, water level variation) estimation in the existing research scenarios. However, in the coal mining subsidence area, the vertical displacement of the GNSS antenna can reach several or even hundreds of centimeters over a short period (e.g., several months). That is, the GNSS antenna normal height and GNSS-IR RH are changing simultaneously and significantly in the underground coal mining subsidence areas. In this scenario, the GNSS RTK must be used to obtain the antenna normal height in near-real-time; and only by combining the antenna’s normal height (H_an) with GNSS-IR-based RH (H_r) can the surface water level (and the surface water level variation) be measured correctly at a given time. Hence, both the GNSS RTK and GNSS-IR are essential for estimating water levels in the underground coal-mining induced ground subsidence scenario.

4.2. Estimation of RH Using Multiple Mode and Band GNSS Signals

When the GNSS signals reach the water surface, part of the signal will be reflected by the surface. If the reflected signal can be captured by a GNSS antenna, the raw SNR observations outputted by the antenna-connected receiver is an interfered signal. This interfered signal includes both the direct GNSS signal and the GNSS signal reflected by the water surface. In general, a direct component of the interfered GNSS signal is a signal with low frequency. The low-frequency direct component can be well-fitted using a second-order polynomial function. Subtracting the direct signal from the raw SNR series, the reflected SNR signal (or, the multipath SNR signal) can be obtained and can be expressed as [37]:

$$SN{R}_m\left(\theta \right)=A_m\cos \left(2\pi \cdot {\displaystyle \frac{2H_r}{\lambda }}\cdot \sin \theta +{\varphi }_m\right)$$

(4)

where A_m and ϕ_m are the amplitude and phase of the reflected SNR signal; θ is the GNSS satellite elevation angle; λ is the GNSS signal wavelength (e.g., 0.1903 m for GPS L1 band signal). Basically, the observed reflection SNR signal is a quasi-sinusoidal signal periodically oscillating with respect to the sine of elevation angle when the elevation angle is lower (e.g., less than 40° for the navigational GNSS receiver), as shown in Figure 9. The figure shows typical examples of the observed multipath GNSS SNR series for three GNSS systems (i.e., GPS, Galileo, and GLONASS).

The multipath series of SNR observations with respect to the sine of satellite elevation angle is an unevenly sampled signal. The Lomb-Scargle periodogram (LSP) method can be exploited to retrieve the amplitude spectrum and phase spectrum of the unevenly sampled signal [38]. Figure 10 shows an example of the observed GPS L1 band multi-path SNR signal, and the corresponding amplitude spectrum and phase spectrum obtained using the LSP method. After determining the peak frequency of the multipath signal from the amplitude spectrum, the RH captured by the multipath signal can be estimated by:

$$H_r=f_{LSP}\cdot \lambda /2$$

(5)

where f_LSP is peak frequency of the multipath SNR signal which can be obtained by exploiting the LSP method to the observed multipath SNR signal. Once the phase, amplitude, and frequency of the observed multipath SNR signal are obtained using the LSP method, the LSP method fitted multipath GNSS SNR signal can also be obtained, which is given by:

$$SN{R}_{LSP}\left(\theta \right)=A_{LSP}\cos \left(2\pi \cdot f_{LSP}\cdot \sin \theta +{\varphi }_{LSP}\right)$$

(6)

where A_LSP and ϕ_LSP are the peak frequency corresponded amplitude and phase

In the field of statistics, the goodness of the fitting result is usually evaluated using the coefficient of determination (CoD) [39]. Similarly, the CoD can also be utilized to indicate the accordance degree between the fitted and the observed multipath SNR signals; and the CoD can be calculated by:

$$R^2=1-{\displaystyle \frac{{\displaystyle {\sum }_{k=1}^n{\left[SN{R}_m\left({\theta }_k\right)-SN{R}_{LSP}\left({\theta }_k\right)\right]}^2}}{{\displaystyle {\sum }_{k=1}^n{\left[SN{R}_m\left({\theta }_k\right)-{\overline{SNR}}_m\right]}^2}}}$$

(7)

where R² is the CoD; SNR_m(θ_k) and SNR_LSP(θ_k) are the observed and fitted multipath SNR at k-th elevation angle; $\overline{{SNR}_m}$ is the average of the observed multipath SNR series. The CoD ranges from zero to one; and the goodness of the LSP fitting result is proportional to the CoD. That is, a larger CoD value indicates a better agreement between the LSP-fitted and the observed multipath signals. Accordingly, the LSP-determined peak frequency with larger CoD, and hence the RH calculated by the peak frequency, should be closer to the ground truth. Therefore, for a RH estimation derived from an observed multipath SNR signal, the accuracy of the RH can be weighted with the CoD of the LSP fitting result. During a given period (e.g., a day), if there are m RH estimations derived from different GNSS systems and band multipath signals, the weighted average of GNSS-IR RH can be calculated by:

$$H_{r\_WA}={\displaystyle \frac{1}{{\displaystyle {\sum }_{i=1}^m{R}_i^2}}}{\displaystyle {\sum }_{i=1}^m\left(R_i^2\cdot H_r^i\right)}$$

(8)

where $H_r^i$ and $R_i^2$ are the LSP peak frequency estimated RH and the RH corresponded CoD for the i-th multipath SNR signal, respectively. By substituting the weighted average of RH given by (8) into (3), the multiple mode and band signals-based GNSS+IR water level estimation in the mining subsidence area can be obtained.

4.3. Data Processing Flow of GNSS+IR Based Surface Water Level Estimation

The data processing flow of the GNSS+IR based surface water level estimation in the area of ground subsidence due to underground coal-mining is shown in Figure 11. The inputs of the estimation include GNSS carrier-phase observation and SNR observation. The former is used to calculate the geodetic height of the GNSS antenna with the GNSS positioning technique (i.e., GNSS RTK). To mitigate affection of the multipath GNSS signal on the positioning accuracy, the GNSS carrier phase observations with a high elevation angle (greater than 40°) are selected for calculation of the GNSS antenna geodetic height. Alternatively, the GNSS SNR observations with an elevation angle of less than 40° are used to calculate RH of the antenna based on GNSS-IR, because the SNR observations with a low elevation angle are easily affected by the reflected GNSS signal and are more suitable for GNSS reflected signal-based RH estimation technique (i.e., GNSS-IR).

As illustrated in Figure 11, the data processing flow of the GNSS+IR-based surface water level estimation in an area of ground subsidence due to underground coal-mining can be described as:

(a)

Select the GNSS carrier phase observations with an elevation angle greater than 40°;

(b)

Calculate the geodetic height (H_a) of the GNSS antenna at a given time with the GNSS RTK technique;

(c)

Convert the geodetic height of the GNSS antenna to the normal height (H_an) using Equation (2);

(d)

Select the GNSS SNR observations with an elevation angle less than 40°;

(e)

Calculate reflector the height (H_r) of the GNSS antenna at a given time with the GNSS-IR technique;

(f)

Based on Equation (8), calculate the weighted average (H_{r_WA}) of the GNSS antenna reflector height estimations derived from different GNSS systems and bands;

(g)

Calculate the surface water level estimation by substituting the normal height and reflector height of the GNSS antenna at a given time into Equation (3).

4.4. GNSS and the Reference Data Processing Methods

Based on the consumer-grade GNSS chip (u-blox F9P) and GNSS antenna (u-blox ANN-MB), a low-cost GNSS receiver was developed. The receiver was utilized to monitor the surface water level with GNSS+IR. The receiver can process dual-frequency band signals transmitted by four GNSS constellations including GPS L1 and L2, Galileo E1 and E5b, GLONASS G1 and G2, and BDS B1 and B2 signals. Details related to the chip and the structure of the receiver can be referred to [40,41]. Because the strength of the reflected GNSS signal is inversely proportional to the elevation angle, it is beneficial to mitigate the impact of multipath error on GNSS positioning accuracy by using observations with high elevation angles. However, the received GNSS signal is more sensitive to geophysical variations of the reflecting surface, such as changes in the water surface, when the elevation angle is lower. In this respect, GNSS SNR observations with low elevation angles would be more suitable for GNSS-IR-based reflector height estimation. Thereby, the carrier-phase observations with elevation angle greater than 40° were used to determine the plane coordinate and geodetic height of the GNSS antenna based on GNSS RTK. The SNR observations within an elevation angle range of 20° to 40° were used to calculate GNSS-IR RH. The GNSS observation types and corresponding functions in the water level estimation can also be seen clearly in

Table 1. The GNSS observation types and corresponding functions in the experimental campaign.

Observation Types	Function	Utilized Technique	Utilized GNSS Information
			System	Observation Code	Elevation Angle Range
Carrier phase	Measure Ha	GNSS RTK	GPS	LlC, L2L	40° to 90°
			Galileo	L1C, L7Q
			GLONASS	L1C, L2C
			BDS	L2I, L7I
SNR	Measure Hr	GNSS-IR	GPS	SlC, S2L	20° to 40°
			Galileo	S1C, S7Q
			GLONASS	S1C, S2C

.

The observing interval of the GNSS instrument is 1s. By averaging all the GNSS RTK measured geodetic height observations in a day (~25,000 observations with fixed solution), the daily averaged geodetic height of the antenna can be calculated. Then, the height is converted into normal height by using (2), obtaining the daily normal height of the GNSS antenna.

The number of satellites utilized to calculate GNSS-IR RH is 9, 22, and 22 for GPS, Galileo, and GLONASS respectively; these satellite sky trajectories during the experiment are shown in Figure 12. Assuming that the GNSS-IR RH estimations derived from different GNSS system and band signals in a day have the same weight (e.g., 1), and by substituting these weights and estimations into (8), the daily normal average of GNSS-IR RH for combining multiple systems and (or) band signals can be obtained. By calculating the weight of the RH estimation with CoD and substituting the calculated weights and estimations into (8), the daily weighted average of GNSS-IR RH for combining multiple system and (or) band signals is obtained. Substituting the daily normal height of the GNSS antenna and the daily normal average of GNSS-IR RH into (3), the normal average (NA) of the surface water level in the coal mining subsidence area can be obtained. Similarly, the weighted average (WA) of the water level can also be obtained by substituting the daily normal height of the antenna and the daily weighted average of GNSS-IR RH into (3). Formulas related to the calculation of WA and NA of the water level are given by:

$$\left\{\begin{array}{l}{H}_g^{WA}=H_{an}-{\displaystyle \frac{1}{{\displaystyle {\sum }_{i=1}^m{R}_i^2}}}{\displaystyle {\sum }_{i=1}^m\left(R_i^2\cdot H_r^i\right)}\\ H_g^{NA}=H_{an}-{\displaystyle \frac{1}{m}}{\displaystyle {\sum }_{i=1}^m{H}_r^i}\end{array}\right.$$

(9)

Notice that, based on the location of the GNSS station in relation to the water surface shown in Figure 6, the GNSS reflection tracks with an azimuth angle range of 270° to 360° were utilized to estimate the GNSS-IR RH. In the azimuth angle range, the number of BDS reflection tracks is rather limited (less than two tracks in a day), which will result in an unreliable performance assessment of the BDS signals in the water level estimation. The maximum number of channels that can be utilized to receive BDS signals is four for the u-blox F9P chip, which should be the main reason for the limited availability of BDS tracks. Therefore, the BDS multipath signals were not utilized in the experiment.

The height from the GNSS antenna to the WLG (i.e., Δh) was manually measured by using a ruler, which is 0.742 m. Then, the GNSS RTK and WLG (RTK+WLG) measured surface water level can be calculated by:

$$H_{g\_RTK+WLG}=H_{an}-\Delta h-h_w$$

(10)

The measuring interval and accuracy of the WLG are 10 s and ±5 mm, respectively.

A DSZ1 level was used to measure the normal height of the top of the metal tube H_lst. The base point of the leveling survey is located in an area unaffected by ground subsidence. The vertical distance from the metal tube top to the water surface (H_rws) was manually measured using the metal ruler attached to the tube. By subtracting this vertical distance from the normal height of the tube top, the leveling surveying method (LSM) based surface water level (H_{g_ls}) can be obtained, as illustrated in (11).

$$H_{g\_ls}=H_{lst}-H_{rws}$$

(11)

There are seven LSM-based surface water level observations during the experimental campaign. Both the RTK+WLG-based and the LSM-based surface water level observations were compared with the GNSS+IR-derived ones in the experimental campaign.

5. Results

The time series of daily normal averages (NA) of the GNSS+IR water level for GPS, Galileo, and GLONASS dual-frequency signals, and the water level measured by the RTK+WLG method, are shown in Figure 13. Additionally, seven LSM-based water level observations are also displayed in the figure. Using the RTK+WLG based water level as the reference data, the Mean, STD, and RMS errors of the GNSS+IR water level estimations can be obtained, as shown in

Table 2. Mean error, error STD and RMSE of GNSS+IR-based surface water level estimations with reference data of RTK+WLG.

System	Band	Mean (m)	STD (m)	RMSE (m)
GPS	L1	−0.006	0.013	0.014
	L2	0.012	0.016	0.020
Galileo	E1	−0.015	0.023	0.027
	E5b	0.018	0.055	0.058
GLONASS	G1	0.009	0.011	0.014
	G2	0.005	0.019	0.019

. The ground subsidence of the GNSS station due to underground coal-mining was 38 mm during the experiment, and the GNSS RTK-based ground subsidence time series is also presented in Figure 13. The RTK+WLG measured surface water level ranges from 41.725 m to 42.092 m during the experiment.

The scatterplots of the normal average (NA) and weighted average (WA) of the multiple-system-and-band-signals combined water level versus the RTK+WLG-based water level are shown in Figure 14. The scatterplots of the normal average (NA) and weighted average (WA) of the multiple-system-and-band-signals combined water level versus the LSM-based water level are shown in Figure 15. Taking the RTK+WLG and LSM water level as the reference data respectively, the Mean, STD, and RMS of errors of the NA and WA water level estimations can be obtained, which is given by

Table 3. The mean error, error STD, and RMSE of GNSS+IR based surface water level estimation in coal mining subsidence area.

Reference	Method	Mean (m)	STD (m)	RMSE (m)
RTK+WLG	NA	0.004	0.013	0.013
	WA	0.004	0.007	0.008
LSM	NA	0.004	0.011	0.011
	WA	0.002	0.007	0.007

.

6. Discussion

It can be observed from Figure 13 and Table 2 that, the accuracy of GPS based water level estimations is equivalent to that of the GLONASS based ones basically. The performance of the second band signal (e.g., GPS L2) derived estimation is slightly inferior to that of the first band signal (e.g., GPS L1) derived one for these two systems. The inferior performance of the second band signal should be mainly attributed to the quality of the LSP-fitted multipath signal. Because the CoD of the second band multipath signal fitted by the LSP is lower than that of the first band signal (e.g., about 0.50 and 0.55 for typical GPS L2 and L1 band signals respectively), the RH, and thus the water level, estimating error of the former is greater than that of the latter. This also infers the possibility of weighing the GNSS-IR RH with CoD to improve the water level estimating accuracy.

Due to the same reason, the Galileo E5b band signal derived water level estimating error is greater than E1 band signal derived one. Furthermore, Galileo based water level estimating error is obviously larger than other two systems; e.g., RMS of the water level estimating error is 0.027 m and 0.014 m for Galileo E1 and GPS L1 band signals, respectively. The signal transmitted by Galileo satellite makes use of the improved modulation technique (i.e., CBOC and AltBOC for E1 and E5b band signal), which performs better in suppression of multipath signal especially when the satellite elevation angle is larger [42,43]. The improved modulation technique could decrease quality, and thus CoD, of the multipath signal significantly, resulting in a larger RH and water level estimating error. Besides, the less available number of reflection tracks (about four, eight, and six tracks for Galileo, GPS, and GLONASS respectively) in a day should be another reason causing the poor performance of Galileo system in the experiment.

In addition, the mean errors provided in Table 2 show different signs and inter-differences in water level estimations obtained using different GNSS frequencies. The different signs and inter-differences might be induced by the inter-frequency bias of GNSS-IR RH. The inter-frequency bias might be related to the combined interaction of the Fresnel reflection coefficients or the antenna gain pattern [44,45]. However, since the main reason causing the bias has not been well identified and proved, it has been difficult to develop a reliable correction model to accurately compensate for the bias until now.

It can be observed from Figure 14 and Figure 15 that the NA based water level estimation deviates to some extent from the RTK+WLG (and LSM) based one; whereas the deviation can be compensated by making use of the proposed weighted averaged method. This can also be seen clearly in Table 3. It can be seen from the table that, the RMS errors of the normal average results are 0.013 m and 0.011 m for the RTK+WLG and LSM-based reference datasets, respectively. By utilizing the weighted average method, the RMSE is reduced by 38.5% and 36.3% for these two datasets, respectively. Overall, the weighted average method demonstrates good effectiveness in enhancing the accuracy of the surface water level estimations.

7. Conclusions

In this paper, the combined technique of GNSS RTK and GNSS-IR is proposed for estimating surface water level of ground subsidence due to under-ground-coal-mining with the low-cost GNSS receiver. The method can make use of a single GNSS station to monitor ground deformation and surface ponding water level caused by underground coal mining simultaneously. The weighted average model has also been developed to combine multiple GNSS system and band signals-based RH estimations; the model can improve the water level estimating accuracy significantly.

Due to the limited availability of BDS multipath signals, they were not utilized in the research. As studied in [46], there is a noticeable difference in the RH estimating performance between the BDS2 and BDS3 signals. The next work will focus on evaluating the performance of surface water level estimation in the area of ground subsidence due to underground coal-mining using different BDS signals. Besides, comprehensive research on the generation mechanism of the inter-frequency bias and the development of a reliable correction model will also be studied in the future.