1. Introduction
Lakes, as an important part of the global hydrological system, not only affect local ecosystems, but also affect many aspects of human life [
1,
2]. Lake ice covers are present in the northern hemisphere and alpine regions during winter [
3]. Ice covers play a critical role in physical, biogeochemical, and ecological processes in lakes [
4]. The overall amount of lake ice has shrunk in the past few decades, as evidenced by delayed dates of freezing or earlier thawing, as well as shorter freeze-up periods, although the change has been spatially heterogeneous [
3,
5,
6]. The most significant reductions have been observed in the Qinghai–Tibet Plateau, Eastern Europe, and Alaska [
7]. Changes in regional lake ice act as a sensitive indicator of response to climate change. Ice covers are also important for ice transport, fishing activities, and tourism. In winter, the ice is strong enough to support roads, and the narrow width of the lake means that transport can reach the opposite coast quickly [
8]. As an important physical parameter of ice-covered lakes, ice thickness reflects the intensity of energy exchange and material migration at the water–air interface, which is of great ecological value [
9]. In the field of ice engineering, ice thickness is an important parameter for ice control and ice utilization, and is a difficult physical indicator to monitor.
Ice thickness is usually measured using contact and non-contact methods. Contact measurement methods include manual drilling, and the use of resistance heating wires and pressure sensors [
10]. Although the above methods are accurate, they have shortcomings such as low efficiency, poor operational safety, poor mobility, etc., making it difficult to achieve continuous observation in the region with these approaches. Methods of non-contact ice thickness measurement have been increasingly applied as equipment has improved, including electromagnetic sensors [
11], shallow water ice profilers (SWIP) [
12,
13], upward-looking sonars [
14], ground-penetrating radar (GPR) [
10], and satellite-borne altimeters [
15].
With the continuous improvement in GPR technology, it has been widely used in various research fields. GPR has shown prospects in the observation of ice thickness because of its advantages of being non-destructive, continuous, rapid, efficient, and operability in real time, and it has been successfully employed in relation to sea ice [
16,
17], river ice [
10,
18,
19,
20], glacial ice [
21,
22,
23,
24], lake ice [
25,
26,
27], and reservoir ice [
28]. Stevens et al. [
17] described a multi-frequency GPR system with 50 MHz, 100 MHz, and 250 MHz antennas that identified shallow subsurface features of floating and grounded ice, water depth, and sedimentary formation in the nearshore Mackenzie Delta in northwest Canada. Fu et al. [
10] performed a study using dual-frequency GPR with 100 MHz and 1500 MHz antennas to measure ice thickness and water depth in the Mohe section of the Heilongjiang River and the Togtok section of the Yellow River in Inner Mongolia. Galley et al. [
18] measured snow and ice thickness at the Churchill Estuary using 250 MHz and 1 GHz dual-frequency GPR. Using the positions of the reflected media in GPR images, combined with the radar waveform amplitude and polarity change information, the ice thickness and the changing medium position at the bottom of a temperate glacier were identified by Liu et al. [
24]. Gusmeroli et al. [
26] demonstrated that a commercially available, lightweight, 1 GHz ground-penetrating radar system can detect and map the extent and intensity of overflow. Li et al. [
28] tested the detection of ice thickness using GPR in the Hongqipao reservoir; ice crystals, gas bubbles, ice density, and ice thickness were also determined and validated by concurrent drilling. Researchers also continue to explore ways to make non-contact ice thickness measurements more effective. A UAV is a fast and efficient means of carrying a GPR; this can solve the problem of unstable ice thickness detection, and is fast and safe for ice thickness detection at the beginning of lake ice formation and during the stable freezing period. While still at infancy, remote-sensing UAVs have some advantages over their traditional counterparts.
The Qinghai–Tibet Plateau (QTP) is located at the crossroads of the Indian monsoon, the East Asian monsoon and the westerly circulation, and its special geographical location makes it exceptionally sensitive to climate change. The QTP is home to a large number of freshwater, saltwater, and saline lakes, and the regional differences in lake ice formation and thickness are important indicators of climate change research, as well as being of great value to the regional ecosystem and human production and life. Previous studies have focused on lake ice phenology, but few have used GPR to detect lake ice thickness on the QTP [
6]. Qinghai Lake is the largest saltwater lake in China and is an important water body that maintains the ecological security of the northeastern QTP. To assess the proposed method, observation fields were set up at Qinghai Lake and Gahai Lake, and several flight experiments were conducted to validate the accuracy of the UAV-borne ice-penetrating radar used for measuring ice thickness.
2. Study Area
Qinghai Lake is located on the northeastern edge of the Tibetan Plateau, with a geographic range of latitude of 36.32°N–37.15°N and longitude of 99.36°E–100.45°E. As the largest inland saltwater lake in China, it covers an area of 4500 km
2; it is 109 km long and 40 km wide, it has an average water depth of 18.3 m and a maximum water depth of 26.6 m, and an average altitude of 3200 m. It has a typical plateau semi-arid alpine climate, with an annual average temperature of −0.53 °C. The annual average precipitation and evaporation are 400 mm and 800~1100 mm, respectively. The lake is mainly fed by rivers, followed by precipitation and groundwater. The salinity of the lake water is about 7100 mg/L, and due to the high content of inorganic salts in the lake, its freezing temperature is lower than 0 °C. The local temperature usually drops below 0 °C in mid-November, and the lake begins to freeze in mid-December. The temperature is the lowest in January, forming a stable ice sheet and freezing completely. In late March, the ice sheet tends to break, and the floating ice begins to melt. By early April, the ice sheet will have completely melted, with an annual freeze-up period of about 88 days [
29]. Gahai Lake is a sub-lake formed by the retreat of Qinghai Lake as its water level has declined, located in the northeastern part of the main body of Qinghai Lake at a distance of about 3.5 km, with an area of 47.2 km
2 and a water depth of 8.0–9.5 m.
In this paper, we set up two observation fields on the western shore of Qinghai Lake in January 2019. Field 1 was located on the west side of Bird Island, with specifications of a 200 × 200 m area and sample point intervals of 50 m, for a total of 25 sample sites (
Figure 1b). Field 2 was also located in the Bird Island area, similarly shaped to a rectangle, and had the same 50 m sample intervals with 10 sample sites (
Figure 1c). In total, five sample sites were established during the flight experiment in Gahai Lake in January 2021, and the locations are shown in
Figure 1d.
3. Materials and Methods
3.1. Radar Principle of Ice Thickness Measurement
According to the principle of radar wave propagation, the propagation velocities of radar waves in different media are related to the dielectric constant. The dielectric constants of air, ice, water, and sandstone (slit) are 1, 3–4, 80, and 5–40, respectively. Therefore, there are differences in the propagation velocities of radar waves in different media, and the received signals can reflect the positions of these different media [
30]. A radar wave enters two different types of media at different speeds, which causes not only the amplitude of the reflection coefficient R to change, but also the polarity (positive and negative) to change in relation to the characteristics of the reflection source waveform [
24]. The reflection coefficient R of the radar wave can be calculated using Equation (1):
where
R is the reflection coefficient, and
ε1 and
ε2 are the dielectric constant of two different media. According to Equation (1), the greater the difference between the dielectric constants of the two media, the larger the absolute value of the reflection coefficient. The positive- and negative-phase polarities change in response to a significant change in dielectric constant, which occurs when the radar wave propagates through the interfaces of different media. Therefore, the locations of the different interfaces can be determined by the change in phase polarity of the radar wave image.
As shown in
Figure 2, we have only used a drone with a radar. The radar emits high-frequency radar waves, and when it reaches the air–ice interface and the ice–water interface, areflection occurs. The drone was flown at a constant speed, and the reflected radar wave signals were received by the antenna. We assumed that the radar waves took the same amount of time in transmitting to the interfaces and returning to and being received by the antenna, so the two-way travel time through the media can be calculated using Equation (2):
where
t is the two-way travel time of radar waves,
H is the ice thickness,
v is the propagation velocity of the radar wave in ice, and
d is the distance between the transmitting antenna and the receiving antenna, i.e., it is the distance moved by the drone in time
t.
According to the radar theory, the propagation velocity of radar waves in a specific media can be expressed as:
where
v represents the propagation velocity of the radar waves in the medium (m/ns),
c denotes the propagation velocity of radar waves in a vacuum (0.3 m/ns), and
ε denotes the ice layer’s real complex dielectric constant.
Then, the thickness of the measured media (
H) can be obtained using Equation (4):
when performing in situ measurements, the value of
d is small relative to
H, and can be ignored, so the equation can be simplified as follows:
Since the velocity of radar wave propagation in ice is influenced by the dielectric constant, the key to determining the thickness of ice is to determine the equivalent dielectric constant. When the in situ thickness of the ice and the radar wave propagation time are known, the equivalent dielectric constant can be deduced.
3.2. Radar Equipment
We used an IGPR (ice- and ground-penetrating radar) system (IGPR-30, Dalian Zoroy Scientific Instrument Co., Ltd., Dalian, China) to measure ice thickness. The ice detection accuracy reached the mm level, and the maximum detection depth was 6 m. Four modes of sampling points were available: 526, 512, 1024 and 2048, and 1024 was selected for sampling. The minimum time window was 2 ps, and the radar dynamic range was greater than 120 dB. The radar was connected to a GPS antenna, while the latitude and longitude coordinates were positioned with cm-level accuracy, and the UAV was equipped with a 200 W pixel camera to record video. The radar was powered with lithium batteries, and its working ambient temperature was between −30 °C and +60 °C. A GPR with a center frequency of 400 MHz was used for the January 2019 experiment, and one with a frequency of 900 MHz was used in January 2021. Both frequencies of GPR are able to meet the depth range requirements of lake ice detection.
The radar was carried on a multi-rotor UAV platform (Matrice 600 Pro, Shenzhen DJI Innovation Technology Co., Ltd., Shenzhen, China), as shown in
Figure 3. The UAV had a maximum payload of 6 kg, meaning it could carry the weight of the radar module (5 kg). During the airborne measurements, the maximum flight altitude was 20 m, the maximum range of each flight was 5 km, and the maximum remote controlling distance was 2 km.
3.3. Data Acquisition
Two observation fields were set up on the western shore of Qinghai Lake in January 2019, as shown in
Figure 4a (an overhead image of Bird Island). We used manual drilling and airborne radar measurements to determine the ice thickness. Sample sites were placed in the observation fields at equal intervals of 50 m. The locations of the sample sites were marked with yellow tape stickers. We then used the DJI GS Pro 2.0 software for route planning based on the locations of the sample sites. To reduce radar signal attenuation along the flight path and ensure the attitude stability of the UAV, we set the altitude to 2 m and the flight speed to 2 m/s, and a wind speed of force 5 or less was assumed when flying (
Figure 4b). Markers were placed at predetermined borehole locations prior to acquiring radar data. When the UAV was operating, videos were acquired using a camera mounted on the side of the radar. During data processing, the radar waves were compared with the video to obtain instantaneous radar images of the marker. The airborne device emitted radar waves continuously during the flight, meeting the detection requirements of a single site or continuously changing positions. We immediately drilled holes at the sample sites after the aerial measurements to ensure temporal and spatial consistency in the radar ice thickness measurements and the drilled ice thickness measurements (
Figure 4c). After drilling, we measured the thickness of the lake ice using an electronic caliper (
Figure 4d). In total, 35 holes (P1~P35) were drilled after the radar data acquisition and were used for verifying the interpretation of the radar profiles. Using the same methodology, we obtained ice thickness data from a further 5 boreholes (P36~P40) and UAV-borne radar measurements for Gahai Lake in January 2021.
3.4. Radar Data Pre-Processing
The raw radar data were processed using GPR studio 2.0 software. To develop better machine understanding of the targeted characteristics, four processes were completed: zero-offset removal, bandpass filtering, amplitude gain, and layer tracking, respectively. Zero-offset refers to the unavoidable occurrence of low-frequency noise in each radar trace, predominantly arising from inductive coupling effects or electronic saturation along the ground–air interface [
31,
32]. This produces a zero-offset signal, which affects the subsequent process because the low-frequency component blocks the true signal after the amplitude gain. Subtracting the mean signal amplitude is a practical approach to eliminating the zero-offset phenomenon [
33].
where
x′(
t) is the original signal trace,
x*(
t) denotes the signal after zero-offset removal, and
ti (
i = 1, 2,
…,
n) represents the sampling time.
Figure 5 shows a comparison between before and after zero-offset removal.
The bandpass filtering is in the frequency range from 200 to 2000 MHz, and the window function uses the Hamming. During the propagation of radar waves, energy decays with time due to dielectric losses and geometric spreading consumption (strong scattering and refraction) [
34]. Considering the above issues of attenuation, an amplitude gain in radar echoes is required, which is a necessary process when seeking to amplify the return signal to a distinguishable level, and can be achieved by combining the exponential and linear functions. In
Figure 6a, an original image of the radar data is shown.
Figure 6b shows the image after zero-drift removal and bandpass filtering, and
Figure 6c shows the radar echo after the amplitude gain.
Images of thickness measured by the radar typically take the form of pulse waves. The positive and negative peaks of the waveform are represented, respectively, in black and white, or in gray scale and color. In this way, the interfaces of different media and their fluctuations can be clearly inferred from the in-phase axis or equal gray scale and equal color line [
24]. As shown in
Figure 6c, the first interface encountered by the radar wave during propagation is the air–ice interface, where the color change in the in-phase axis is black–white–black and the corresponding change in polarity is + - +, indicating that the radar wave has reached a medium with a higher dielectric constant than that of the original (ε air < ε ice). The second interface is the ice–water interface, where the color change in the in-phase axis is also black–white–black, and the polarity change is also + - +, indicating that the radar wave reaches a medium with a higher dielectric constant than ice (ε ice < ε water). After the above pre-processing of the radar echoes, the air–ice and ice–water interfaces can be clearly identified in the radar echo images, thus allowing us to determine the time taken for the radar waves to propagate through the ice and calculate the ice thickness from the time. We determined the locations of the drilled holes by examining the images taken during drone flight while hovering and flying, corresponding to the range of radar echo image channels in order to calculate the simulated thickness at the drilled hole locations.
Figure 6d shows the results of layer tracking with pre-processed radar echoes.
3.5. Extraction of Amplitudes
In order to obtain the characteristics of radar wave variations at different interfaces, a script program was written to extract the amplitudes from all radar data, and the amplitude variations in the characteristic regions were plotted according to the locations of the sampling sites and the radar wave traces.
3.6. Accuracy Validation
In order to quantify the accuracy of our airborne radar measurements of lake ice thickness, root mean square error (RMSE) and percentage error were used, as follows:
where
H0 is the in situ thickness,
H is the radar thickness,
λ is the percentage error, and
n represents the number of sampling sites.
5. Discussion
When utilizing a radar to calculate ice thickness, it is crucial to accurately determine the velocity of radar wave propagation through the ice. The results of this study for radar wave propagation velocity differ from those of previous studies. In a slight departure from other scholarly studies on sea ice [
16,
17], river ice [
18,
20,
30], glacial ice [
21,
24,
36], lake ice [
26], and reservoir ice [
28], our experiments yielded lower mean velocities and higher dielectric constants than freshwater ice. Experiments by Zhang [
30] undertaken on the Yellow River showed that radar waves propagate at a velocity of 0.163 m/ns during peak ice and at less than 0.150 m/ns during melting. This illustrates that differences in the moisture content on the ice surface or in the ice affect the velocity of radar wave propagation. Therefore, regarding the annual lake ice of Qinghai Lake and Gahai Lake, the moisture content is generally high, causing the radar wave propagation velocity to be low. In addition, considering that Qinghai Lake is a brackish lake located on the Qinghai–Tibet Plateau, the higher salinity of the water causes its freezing point to be lower than that of freshwater ice bodies, and temperature changes in the ice in turn affect the dielectric constant of the ice. In contrast, although Qinghai Lake and Gahai are in the same area, the propagation characteristics of radar waves through the ice on the two lakes are somewhat affected by differences in their salinity. Salt and brine lakes are also present on the Qinghai–Tibet Plateau [
37], and their salinity affects the freezing of lake water, while the dielectric constants of these lakes’ ice differ from those of freshwater ice.
Air bubbles in the ice can increase the velocity of radar wave propagation, and thus affect the accuracy of ice thickness calculations [
35]. When the rate of proportion of air bubbles in the ice is high, the velocity of the radar waves can reach 0.21 m/ns [
38], which would result in a large error if the standard velocity of radar waves in ice (0.17 m/ns) were to be used [
39,
40]. Li [
28] showed that radar waves propagate at a speed of 17.02 m/ns in granular ice and 16.98 in columnar ice. This implies that differences in ice structure and air bubble content affect the propagation characteristics of radar waves in ice. Due to the bubbles in the ice, the GPR waves are scattered and the average velocity increases with the bubble volume ratio; the velocity increases by 0.0014 m/ns for every 1% increase in bubble volume [
35]. Ice thickens via the congelation of ice at the water and ice interface. Congealed ice contains elongated bubbles, which are trapped within the ice, and cause radar waves to be scattered and reflected back [
41]. In general, the low bubble contents in lake ice on the Tibetan Plateau and the unfrozen water in the ice counteract the increase in radar wave velocity in the ice, meaning that the ice thickness can be calculated with high accuracy [
35]. According to field observations, the ice on Qinghai Lake and Gahai Lake during the freeze-up period is columnar. This is because, when the lake water was frozen, salt was entrained in the ice in salt cells, and over time, brine channels formed between the ice, and the residual high concentration of salt was discharged along the channels. The bubble contents of brine channels in the ice are low, meaning that they have a reduced effect on the radar wave propagation velocity. Since the dielectric constant of ice is related to its crystal structure and temperature, it is necessary to determine the in situ thickness of ice using boreholes during each probing, invert the effective dielectric constant of the ice layer based on in situ ice thickness, and then simulate and calculate the final ice thickness based on this value [
30].
The type of lake ice also affects the amplitude of the radar waves. In this experiment, the drilling results show that the majority of the measurement area comprised floating ice. Due to the great difference between the dielectric constants of ice and water, the reflection energy of the ice–water interface is strong, and the amplitude of reflection of radar waves by floating ice is significantly stronger than that of grounded ice—about 3~4 times the average reflection amplitude [
17,
20,
42]. The average reflection amplitude of the floating ice in the experiment was seven times greater than that of touchdown ice, which is higher than values derived in previous studies.
In addition to the propagation velocity within the ice, the accuracy of ice thickness measurements is also related to the frequency of the antenna. The frequency of the transmitting antenna determines the detection depth, while the frequency of the receiving antenna has an effect on the resolution [
43]. High-frequency signals in low-velocity media are characterized by short wavelengths and a high resolution, while low-frequency signals in high-velocity media have a low resolution. The attenuation of radar signals is significantly affected by the antenna frequency, with higher frequencies being more significantly attenuated than lower frequencies. This frequency attenuation behavior therefore explains the greater penetration of radar signals at lower frequencies than at higher frequencies [
44]. In general, longer waves penetrate deeper into the medium, while shorter waves provide a higher resolution. Once the approximate depth of the target is known, it is recommended to select a higher-frequency antenna capable of penetrating to that depth in order to secure the required resolution and penetration depth. In the experiments in this paper, we used 400 and 900 MHz antennae. During data processing, it was found that the images acquired by the 900 MHz frequency radar had a higher resolution and captured data from the interface more accurately compared to 400 MHz, and the results of ice thickness measurements show that the 900 MHz radar also had better accuracy.
Arbitrary underground medium inhomogeneity and undesired measurement noises often disturb radar data interpretation [
45]. A lake with a more homogeneous subsurface, such as Qinghai Lake, has a flatter ice surface, and this ensures better data collection when it enters the stable freeze-up period in winter. However, in the initial freezing period when the ice is thin, or when an ice–water mixture is formed on the lake surface, it is more difficult to distinguish the upper and lower ice interfaces in the radar image; only the air–lake interface can be easily identified. In addition, working with UAV-borne ice-penetrating radar requires control of the UAV’s altitude. We found that when the UAV returned at an altitude of 20 m, ice information was not recognized, even though the radar data had been pre-processed. Therefore, it is advisable to keep UAV flights to a height of less than 10 m.
In conclusion, excluding the above factors that affect the propagation characteristics of radar waves in ice. In order to secure the accuracy of the experimental data, we should obtain the in situ thickness and the propagation time of the radar wave as accurately as possible. Since the velocity of radar wave propagation in ice is affected by the dielectric constant, the critical factor in determining the ice thickness using airborne ice-penetrating radar is to determine the equivalent dielectric constant. As shown in Equation (5), since the propagation velocity of radar waves in air is known, the equivalent dielectric constant can be derived when the in-situ thickness of the ice and the propagation time of the radar waves are obtained. This is the rationale for using two different equivalent dielectric constants for Qinghai Lake and Gahai Lake, two lakes with differences in salinity, temperature, and other properties of the water. When modeling ice thickness using this equipment, it is important to accurately obtain the in situ ice thickness and radar wave propagation time in order to calculate a more accurate equivalent dielectric constant with as many measured points as possible.
6. Conclusions
In this study, UAV-borne ice-penetrating radar was used to detect the ice thickness of Qinghai Lake and Gahai Lake, and the ice thickness at in situ boreholes was used to verify the radar measurements. The echo pattern characteristics of the airborne radar and the propagation characteristics of the radar waves in the lake ice were clarified. The applicability of a UAV-borne ice-penetrating radar for ice thickness detection in winter glacial lakes was explored. The main conclusions are as follows:
After pre-processing the echo images obtained with the UAV ice-penetrating radar, we were able to clearly assess the lake ice developing during the stable freezing period and accurately determine the propagation times of the radar waves through the ice.
Contrasting our in situ measurements of ice thickness, the mean velocity of radar wave propagation in the ice was lower than that in other regions, and the dielectric constant was higher. Compared with in situ ice thickness measurements for Qinghai Lake, measurements taken using a 400 MHz radar showed an RMSE of 0.034 m with a percentage error of 8.6%, while the results for Gahai using a 900 MHz radar showed an RMSE of 0.010 m with a percentage error of 2.8%.
The use of a radar with a higher antenna frequency ensures higher detection accuracy, provided that the detection depth is met. This affirms the need to control the altitude and speed of the UAV during airborne ice-penetrating radar thickness measurements.
In future lake ice surveys, the use of UAVs with an ice-penetrating radar will help quickly obtain accurate information on ice thickness. However, in the alpine region, conditions are vicious, with low temperatures, thin air, and other problems, and the life of a UAV’s single battery is short, meaning that it is easy to lose altitude when flying continuously. These difficulties will need to be overcome in the future in order to complete a large number of lake ice surveys using an airborne radar.