A method for detecting polar subglacial water based on short-time Fourier transform

Abstract

The present invention relates to a method for detecting polar subglacial water bodies based on short-time Fourier transform, comprising the following steps: 1) obtaining radio echo sounding data; 2) obtaining a two-dimensional radar image after heading fuzzy processing; 3) cutting and splicing the two-dimensional radar image after heading fuzzy processing to obtain a cut two-dimensional radar image; 4) waveform reshaping of the interface signal in the cut two-dimensional radar image; 5) short-time Fourier transform processing of the two-dimensional radar image after waveform reshaping to obtain the frequency response and amplitude of each survey line at the interface; 6) combining the normalized frequency response, the frequency response amplitude after short-time Fourier transform and the terrain slope to calculate the ice bottom material detection value, and automatically identify the subglacial water body. Compared with the prior art, the present invention uses the ice bottom material detection value as the basis for judging whether the ice bottom material is a subglacial lake and a subglacial water body, and realizes the high-precision resolution of the ice bottom material of a single-channel radar signal.

Description

Technical Field

[0001] The invention relates to the technical field of earth climate and geodynamics, and in particular to a method for detecting polar subglacial water bodies based on short-time Fourier transform. Background Art

[0002] Polar glaciers are an important part of the Earth's climate and geodynamic systems. For a long time, there has been a lack of sufficient geological and geophysical data from the bottom of the Antarctic ice sheet, making it difficult to support the understanding of Antarctic historical changes and the prediction of global climate change development. To solve this problem, scientific research teams from all over the world have collected geophysical data in the Antarctic region, including gravity, magnetism, altimetry, ice radar, etc. through aerial geophysical mapping, which has played a very important role in understanding the potential changes in the balance of the ice sheet. Among them, subglacial lakes are one of the important components of the Antarctic ice bottom environment. It has an important impact on the material balance of the ice sheet, the stability of the ice sheet, global climate change and sea level change. As of February 2022, a total of 773 subglacial lakes have been discovered worldwide, including 675 in Antarctica and 64 in Greenland.

[0003] There are many technologies for detecting ice sheets and under the ice, but radio echo sounding is the most sensitive to the inside and bottom of the ice sheet, which makes up for the defect that deep ice core drilling cannot be implemented on a large scale in the harsh Antarctic region. At the same time, it also overcomes the shortcomings of geophysical detection technologies such as seismic and gravity exploration that are insensitive to ice sheets, providing an important data source for the detection of subglacial water bodies in Antarctica. In the process of searching for subglacial water bodies, the most commonly used method is visual interpretation, but its limitations are also obvious. For example, it is usually time-consuming, which is extremely unfavorable for the interpretation of a large amount of radio echo sounding data in the future. In addition, the interpretation standards vary with the interpreters and are difficult to unify. Therefore, the development of automatic subglacial water detection methods is essential for polar research.

[0004] Traditional subglacial lake detection methods usually rely on the statistical information of two-dimensional ice radar signal (B-Scan) data, such as echo intensity, reflectivity, roughness and water potential, but ignore the inherent time-frequency characteristics of one-dimensional ice radar signal (A-Scan). The lack of subglacial water detection and identification indicators leads to reduced accuracy of polar subglacial water detection. Summary of the invention

[0005] The purpose of the present invention is to overcome the defects of the above-mentioned prior art and to provide a method for detecting polar subglacial water based on short-time Fourier transform. By utilizing the frequency response characteristics of the substrate signal to process the radio echo sounding data in the polar region, a frequency index that can describe the morphological characteristics of the substrate signal is constructed, and the polar thematic information on the possibility of the existence of water in the substrate is obtained by combining indicators such as slope and amplitude.

[0006] The purpose of the present invention can be achieved by the following technical solutions:

[0007] A method for detecting polar subglacial water based on short-time Fourier transform comprises the following steps:

[0008] 1) Obtain radio echo sounding data;

[0009] 2) Performing data fuzzy processing along the track direction on the original two-dimensional radar image to obtain a two-dimensional radar image after heading fuzzy processing;

[0010] 3) performing cropping and splicing processing on the two-dimensional radar image after the heading blur processing to obtain a cropped two-dimensional radar image;

[0011] 4) reshaping the waveform of the interface signal in the cropped two-dimensional radar image;

[0012] 5) Perform short-time Fourier transform processing on the two-dimensional radar image after waveform reorganization to obtain the frequency response and amplitude of each measuring line at the interface;

[0013] 6) The detection value of ice bottom material is calculated by combining the normalized frequency response, the frequency response amplitude after short-time Fourier transform and the terrain slope to automatically identify the subglacial water body.

[0014] In the step 1), the radio echo sounding data are specifically different coherent radar sounder data provided by the Ice Sheet Remote Sensing Center.

[0015] In the step 2), the expression of the two-dimensional radar image after heading fuzzy processing is:

[0016]

[0017] Among them, BScan represents the original two-dimensional radar image, w 1 is the smoothing window width, BScan′ is the 2D radar image after heading blur processing, m represents the total number of rows in the 2D radar image, n represents the total number of columns in the 2D radar image, and subscript j 1 Indicates the row number of the two-dimensional radar image, subscript k 1 Indicates the column number of the two-dimensional radar image, subscript i 1 Indicates the column number within the smoothing window.

[0018] In the step 3), the expression for data clipping is:

[0019]

[0020] Where BScan″ represents the cropped two-dimensional radar image, is the kth in the two-dimensional radar image after heading fuzzy processing2 List the sampling points corresponding to the ice-bottom interface, w 2 is the half window width of the data clipping, and each column signal starts from the The sampling points extend upward and downward w 2 The number of sampling points for cutting is 2w, and the number of sampling points for each A-Scan after cutting is 2w 2 +1.

[0021] The step 4) specifically comprises the following steps:

[0022] 41) For each track of the cropped 2D radar image, the size is 2w 2 The +1×1 A-Scan is used to remove DC drift, and then:

[0023]

[0024] Among them, AScan represents any A-Scan in the cropped two-dimensional radar image, 2w 2 +1 is the number of sampling points for each A-Scan, is the jth 2 sampling points, AScan′ represents the A-Scan after removing DC drift;

[0025] 42) Reshape the interface waveform in the two-dimensional radar image after the DC drift removal process.

[0026] In the step 42), the waveform reshaping specifically includes:

[0027] 421) setting the left and right endpoints of the main peak, and flipping the main peak and splicing them to the left and right endpoints of the main peak respectively;

[0028] 422) Set the values ​​of other sampling points to 0, then:

[0029]

[0030] in, is the jth A-Scan after waveform reorganization 3 sampling point values, l represents the sampling point number corresponding to the position of the left endpoint of the main peak before the waveform is reshaped, r represents the sampling point number corresponding to the position of the right endpoint of the main peak before the waveform is reshaped, and Peak represents the main peak signal.

[0031] The left and right endpoints of the main peak are defined as:

[0032] One sixth of the peak intensity of the interface is set as the threshold, and the two intersection points with the threshold and the closest peak positions on the left and right sides of the A-Scan peak are respectively used as the left and right endpoints of the main peak of the interface.

[0033] In the step 5), the short-time Fourier transform formula is:

[0034]

[0035] Among them, x(τ) is the signal to be processed, w(*) is the window function, the window function for short-time Fourier transform processing at time t is w(τ-t), and ω is the angular velocity.

[0036] In the step 6), the calculation formula of the ice bottom material detection value C is:

[0037]

[0038]

[0039]

[0040] in, For the kth 3 The initial value of A-Scan's ice bottom material detection, For the kth 3 The ice bottom material detection value after A-Scan smoothing, w 3 is the smoothing window width, i 2 is the column number in the smoothing window, F is the normalized frequency response, A is the amplitude of the frequency response after short-time Fourier transform, slope is the terrain slope, α is the influence weight of the terrain slope, e is the natural logarithm, Δz is the elevation difference between the main peaks of two adjacent A-Scans, and Δx is the horizontal distance between two adjacent A-Scans.

[0041] In step 6), the automatic identification of subglacial water bodies is specifically performed as follows:

[0042] The ice bottom material detection value C is used as the basis for the final judgment of whether the ice bottom material is a subglacial lake or subglacial water body. The value of is greater than the set threshold, then the corresponding kth 3 If column A-Scan is a subglacial lake or subglacial water body, then:

[0043] A-Scan for subglacial water bodies, A-Scan for subglacial sediments, A-Scan for bedrock.

[0044] Compared with the prior art, the present invention has the following advantages:

[0045] In the existing subglacial lake detection technology, the absolute and relative sizes of continuous multi-channel echo intensities are often relied upon, supplemented by parameters such as terrain, water potential, and roughness for detection. However, this detection method ignores the inherent characteristics of the A-Scan waveform and lacks the ability to detect and interpret single-channel radar signals. The waveform reshaping method based on short-time Fourier transform proposed in the present invention extracts additional features, and uses the ice bottom material detection value as the basis for the final judgment of whether the ice bottom material is a subglacial lake or subglacial water body, thereby achieving high-precision resolution of ice bottom materials for single-channel radar signals. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] Figure 1 This is a flowchart of polar subglacial lake detection based on short-time Fourier transform.

[0047] Figure 2a It is a two-dimensional radar profile image.

[0048] Figure 2b The waveforms are of different interfaces.

[0049] Figure 3a It is a two-dimensional radar profile after heading blur. The curve represents the position of the sampling points 150 degrees on both sides of the ice bottom interface.

[0050] Figure 3b This is the two-dimensional radar profile after data cropping.

[0051] Figure 4 The A-Scan images after data cropping, where Figure (4a) is the ice-water interface, Figure (4b) is the ice-rock interface, the wavy curve on the left represents the waveform after removing the DC drift, and the straight curve on the left represents the waveform after waveform reshaping.

[0052] Figure (5a) is the short-time Fourier transform result of the ice-water interface. In the figure, (a) represents the A-Scan position of the ice-water interface, (b) represents the reorganized A-Scan waveform of the ice-water interface, and (c) represents the result of short-time Fourier transform of Figure (b). The circular position is the frequency response.

[0053] Figure (5b) is the short-time Fourier transform result of the ice-rock interface. In the figure, (d) represents the A-Scan position of the ice-rock interface, (e) represents the reorganized A-Scan waveform of the ice-rock interface, and (f) represents the result of short-time Fourier transform of Figure (e). The circular position is the frequency response.

[0054] Figure 6The ice radar image and detection results corresponding to the survey line 20081225_02_024, where Figure (6a) is the two-dimensional radar image and detection results (color band), Figure (6b) is the frequency response and frequency corresponding amplitude after short-time Fourier transform, Figure (6c) is the slope of the ice bottom terrain and the echo intensity at the ice bottom interface, and Figure (6d) is the ice thickness and water potential. In Figure (6a), the dark semi-transparent rectangle is the boundary of the subglacial lake of Livingstone et al., and the light semi-transparent rectangle is the boundary of the subglacial lake of Wolovick et al.

[0055] Figure 7 The results of subglacial lake and subglacial water body detection in different radar profiles in the AGAP region are shown in Figure (7a) for survey line 20081228_03_006, Figure (7b) for survey line 20081225_02_024, Figure (7c) for survey line 20090109_01_023, Figure (7d) for survey lines 20090107_03_015 and 20090107_03_016, and Figure (7e) is the survey line 20081223_01_015, (7f) is the survey line 20090102_03_001, (7g) is the survey line 20090108_01_016 and the survey line 20090108_01_017, (7h) is the survey line 20090106_05_013, (7i) is the survey line 20090108_03_021, and (7j) is the survey line 2 0081228_01_020, Figure (7k) is the survey line 20081226_02_017, Figure (7l) is the survey line 20081225_02_018, Figure (7m) is the survey line 20090106_05_020 and the survey line 20090106_05_021, Figure (7n) is the survey line 20090107_02_014 and the survey line 20090107_02_015 , line 20090107_02_016, line 20090107_02_017 and line 20090107_02_018. The dark translucent rectangle is the boundary of the subglacial lake on the CReSIS line published by Livingstone et al. in 2022, and the light translucent rectangle is the boundary of the subglacial lake on the CReSIS line published by Wolovick et al. in 2013.

[0056] Figure 8These are the detection results of the 2009_Antarctica_TO_Gambit directory, where the AGAP area is located in the middle rectangular survey line area, the gray dots are subglacial lakes on the CReSIS survey line published by Livingstone et al. in 2022, the white dots are subglacial lakes on the CReSIS survey line published by Wolovick et al. in 2013, and the dark dots are subglacial lakes calculated by the method of the present invention. The size of the dots represents the possibility that this is a subglacial lake or subglacial water body.

[0057] Fig. 9 Figure 1 is the detection result of the CReSIS Antarctic regional dataset, where the blue solid dots are subglacial lakes and subglacial water bodies detected by the method of the present invention, the red solid dots are the 675 subglacial lakes in the Antarctic region summarized by Livingstone et al., and the white solid line is the ice radar survey line. In the figure, (a) represents the detection result of the entire Antarctic, (b) represents the detection result of the AGAP region, (c) represents the detection result of the Vostok Lake region, (d) represents the detection result of the Byrd Glacier region, (e) represents the detection result of the West Antarctic, (f) represents the detection result of the west side of the Ross Ice Shelf, and (g) represents the detection result of the Coats Land region.

[0058] Fig.10 : This is a comparison between the detection results of the present invention and the water potential prediction results, wherein (a) represents the comparison of detection results in the entire Antarctic, (b) represents the comparison of detection results in the AGAP region, (c) represents the comparison of detection results in the Vostok Lake region, (d) represents the comparison of detection results in the Byrd Glacier region, (e) represents the comparison of detection results in the West Antarctic, (f) represents the comparison of detection results on the west side of the Ross Ice Shelf, and (g) represents the comparison of detection results in the Coats Land region. DETAILED DESCRIPTION

[0059] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0060] like Figure 1 As shown, the present invention aims to address the lack of existing subglacial water detection and identification indicators in subglacial water detection methods, and provides a polar subglacial water detection method based on Short-Time Fourier Transform (STFT). Based on the public ice radar data collated by the Center for Remote Sensing of Ice Shelves (CReSIS) project, a MATLAB algorithm program was written to process airborne ice radar signals in the entire Antarctic and Greenland, and the ice bottom material was classified into bedrock-water, and the geographical locations of all possible subglacial water bodies and subglacial channels in the survey area were found and mapped.

[0061] When electromagnetic waves reach the ice-water interface, most of the electromagnetic waves are reflected due to the extremely high dielectric constant of water (∈w≈81) and its smooth surface, and only very little energy passes through the interface into the water, resulting in a sharply changing feature on the A-Scan. On the contrary, the dielectric constant of rock or other non-aqueous materials is not high (such as rock, ∈r≈4~9), and more electromagnetic waves can enter the ice-rock interface and produce multiple reflections in the rock layer due to the material inhomogeneity. At the same time, the multiple reflections produced by the rough ice-rock interface will coherently cancel each other, which makes the echo intensity at the interface in the A-Scan gradually decrease.

[0062] The change in the relative power of electromagnetic waves under the ice sheet depends only on the change in reflectivity. According to the reflectivity calculation formula (1), the greater the difference in the dielectric constants of different materials on both sides of the interface, the higher the reflectivity, the more drastic the change in echo power, and the narrower and sharper the features at the interface on the A-Scan. Figure 2b The A-Scan waveform at the interface of different materials is shown (the upper black line is the ice-rock interface, and the lower black line is the ice-water interface).

[0063]

[0064] Among them, R c is the reflectivity, ∈ b ,∈ i are the dielectric constants of the materials on both sides of the interface.

[0065] As can be seen from the above, the reflection coefficients of radar signals at the ice-rock interface and the ice-water interface are different, and the main peaks of the ice-water and ice-rock interfaces will also show different characteristics - the main peak of the ice-water interface is narrow and sharp, while the main peak of the ice-rock interface is wider, which makes the subglacial lake show obvious high-frequency characteristics on the radar map, while the subglacial bedrock shows low-frequency characteristics. In the existing subglacial lake detection, it often relies on the absolute and relative size of the continuous multi-channel echo intensity, supplemented by parameters such as terrain, water potential, and roughness for detection. This detection method ignores the inherent characteristics of the A-Scan waveform and lacks the detection and interpretation capabilities of single-channel radar signals. Based on this, the present invention proposes a subglacial water body detection method based on short-time Fourier transform to achieve ice bottom material resolution through single-channel radar signals. The specific implementation steps are as follows:

[0066] Step 1: Preparation of radio echo sounding data;

[0067] The radio echo sounding data used are data products of different coherent radar sounders provided by the Center for Remote Sensing of Ice Sheets (CReSIS), which can be downloaded free of charge. The downloaded product data uses the coverage of the survey line as the data location, and provides information including the longitude and latitude on the survey line, aircraft flight altitude, distance from the aircraft to the ice bed, ice thickness, ice surface elevation, ice bed elevation, radar echo map, and the location of the ice surface and ice bed on the echo map.

[0068] Step 2: Use a data blurring method along the track direction to reduce noise and enhance features of the original data;

[0069] Each A-Scan measured by airborne radar is independent of each other, but adjacent A-Scans often have some common features:

[0070] 1) Interfaces with similar geographical locations often have similar altitudes, that is, the peaks of the interfaces in the A-Scan waveform are located at similar sampling points;

[0071] 2) The noise in each A-Scan is random noise.

[0072] Therefore, the B-Scan is processed for heading ambiguity, that is, the average of several consecutive A-Scans is taken to highlight the similar parts of the A-Scan, and the random noise is filtered out without changing the natural frequency of the strongest peak of the signal. At the same time, since most of the subglacial lake is horizontal and the bedrock often has a certain slope, this will make the peak of the ice-water interface of the reconstructed A-Scan more prominent, while the peak at the ice-rock interface is more gentle, so that the two can be further distinguished by frequency response. The formula for heading ambiguity is:

[0073]

[0074] Among them, BScan represents the original two-dimensional radar image, w 1 is the smoothing window width, BScan′ is the 2D radar image after heading blur processing, m represents the total number of rows in the 2D radar image, n represents the total number of columns in the 2D radar image, and subscript j 1 Indicates the row number of the two-dimensional radar image, subscript k 1 Indicates the column number of the two-dimensional radar image, subscript i 1 Indicates the column number in the smoothing window. The smoothing window width in the present invention is fixed to the number of A-Scans corresponding to the horizontal distance of 300m.

[0075] Step 3: crop and splice the two-dimensional radar image after heading blur processing;

[0076] A complete A-Scan usually includes direct waves, surface reflection waves, isochronous layers, ice bottom reflection waves (bedrock or water body), and some also include isochronous layers caused by ice flow (internal layer, flow disrupted layer), ice (frozen-on ice) that refreezes after melting at the bottom of the ice, and interference stripes similar to barcodes generated by radar. Usually, there is no need to process the entire A-Scan from the surface of the ice sheet to the bottom of the ice sheet. Only some sampling points near the ice bottom interface need to be processed. The present invention uses the interface extracted by CReSIS, but selects the maximum value position within the range of 50 sampling points above and below the interface as the real interface position, and then takes a sampling point above and below the new interface to form a new strip-shaped B-Scan. At this time, the B-Scan usually only contains ice-rock, ice-water interface and a small amount of isochronous layers. The formula for data clipping is:

[0077]

[0078] Among them, BScan″ is the cropped two-dimensional radar image, is the kth in the two-dimensional radar image after heading fuzzy processing 2 List the sampling points corresponding to the ice-bottom interface, w 2 is the half window width of the data clipping, and each column signal starts from the The sampling points extend upward and downward w 2 The number of sampling points for cutting is 2w after cutting. 2 +1, in actual calculation, m is usually 150, and the specific image after cropping is as follows Figure 3b shown.

[0079] Step 4: reshape the waveform of the interface signal;

[0080] The absolute value corresponding to each sampling point has no physical meaning. The difference in echo intensity between bedrock and water is reflected in the change of relative values ​​within the range. In addition, when Fourier transform is performed on A-Scan, it can be found that a large absolute value of the average value of a signal will cause a very strong zero-frequency component in its spectrum. Therefore, each channel in the cropped B-Scan is 2w in size. 2 +1×1 A-Scan to remove DC drift:

[0081]

[0082] Among them, AScan represents any A-Scan in the cropped two-dimensional radar image, 2w 2 +1 is the number of sampling points for each A-Scan, is the jth 2sampling points, AScan′ is the A-Scan after removing DC drift.

[0083] One sixth of the peak intensity of the interface after DC drift processing is set as the threshold, and the two intersection points with the threshold and the closest peak positions on the left and right sides of the A-Scan peak are taken as the left and right endpoints of the main peak of the interface, such as points A and B in Figures (4a) and (4b). The part where the main peak intensity exceeds the threshold is the peak height, and the distance between the left and right endpoints is the peak width.

[0084]

[0085] Among them, l is the sampling point number corresponding to the left intersection of the main peak and the threshold, r is the sampling point number corresponding to the right intersection of the main peak and the threshold, Peak is the part of the main peak above the threshold, The jth peak 4 The value at the sampling point, i 3 is the column number of AScan′.

[0086] Figure 5 The frequency value of the signal between points A and B can be used as the natural frequency value of the interface signal waveform. Since the mean value of the main peak part is still large, the main peak is directly subjected to Fourier transform. The proportion of its zero-frequency component is very high, resulting in its frequency response of 0MHz. Based on this, the present invention performs waveform reshaping on the interface waveform. The specific operation is as follows:

[0087] (1) Flip the main peak and splice it to the left and right ends of the main peak;

[0088] (2) Set the values ​​of other sampling points to 0.

[0089] The specific formula is:

[0090]

[0091] in, is the jth A-Scan after waveform reorganization 3 sampling point values, l represents the sampling point number corresponding to the position of the left endpoint of the main peak before the waveform is reshaped, r represents the sampling point number corresponding to the position of the right endpoint of the main peak before the waveform is reshaped, and Peak represents the main peak signal.

[0092] The comparison chart before and after waveform reshaping is as follows Figure 4 shown.

[0093] Step 5: Calculate the frequency response and amplitude of each measuring line at the interface through short-time Fourier transform, and perform short-time Fourier transform on the signal after waveform reshaping. The formula is:

[0094]

[0095] Among them, x(τ) is the signal to be processed, w(*) is the window function, the window function for short-time Fourier transform processing at time t is w(τ-t), and the short-time signal x intercepted in the window i (τ) = x(τ)w(τ-t), STFT(t,ω) is the time-frequency joint function obtained by short Fourier transform.

[0096] Through the analysis of the short-time Fourier transform method, the two-dimensional time domain-frequency domain matrix of each target A-Scan can be obtained, and the different characteristics of the ice-water interface and the ice-rock interface can be summarized. The results of the short-time Fourier transform analysis of different material interfaces are shown in Figure 2. Figure 5 As shown, by performing a short-time Fourier transform on the reconstructed waveform, it can be seen that the main peak frequency response of the reconstructed ice-water interface is around 1 MHz, while the main peak frequency response of the ice-rock interface is still around 0 MHz. Therefore, a preliminary distinction can be made between bedrock and subglacial lakes based on the frequency response.

[0097] Step 6: Combine normalized frequency value, frequency amplitude, ice bottom slope and other parameters to calculate the detection value of ice bottom material and automatically identify subglacial water bodies.

[0098] The frequency response of some subglacial water-containing sediments is not zero. In order to further distinguish the interfaces of different materials, the results of short-time Fourier transform are input as parameter values ​​into the final judgment function:

[0099]

[0100]

[0101] in, is the initial value of the ice bottom material detection value, For the kth 3 The detection value after A-Scan smoothing, w 3 is the smoothing window width, i 2 is the column number in the smoothing window, F is the normalized frequency response (frequency response divided by the sampling frequency of the radar image), A is the amplitude of the frequency response after short-time Fourier transform; slope is the terrain slope, α is the influence weight of the terrain slope, the larger α is, the greater the influence of the terrain slope on the result, in the present invention, the value of α is 5, and the terrain slope calculation formula is:

[0102]

[0103] Among them, Δz is the elevation difference between the main peaks of two adjacent A-Scans, and Δx is the horizontal distance between two adjacent A-Scans, both in meters.

[0104] In this example, C is the ice bottom material detection value before smoothing, which is a vector, and C′ is the ice bottom material detection value after smoothing, which is a vector; is the kth smoothed 3 The ice bottom material detection value in the column is a numerical value.

[0105] The present invention uses the ice bottom material detection value C as the basis for the final judgment of whether the ice bottom material is a subglacial lake or a subglacial water body. The value of is greater than the set threshold, then the corresponding kth 3 Column A-Scan is a subglacial lake or subglacial water body. The threshold is selected based on the known lakes on the 20081225_02_024 survey line. The minimum value of the C value of the three known lakes is set as the threshold, that is, A-Scan for subglacial water bodies, A-Scan for subglacial sediments, The A-Scan of the detection result is bedrock. The detection result range is represented by the color bar. The closer the color is to blue, the more likely it is subglacial water, and the closer it is to red, the more likely it is bedrock.

[0106] Example

[0107] In the CReSIS dataset, the survey lines in the AGAP region are located in the 2009_Antarctica_TO_Gambit directory, which also includes some survey lines extending to Dome Fuji and Vostok, and a survey line extending to Coats Land. There are 889 valid B-Scans and 1,471,608 valid A-Scans in the directory. B-Scans connected to the same survey line are distinguished by thick black and white lines.

[0108] The detection results of the AGAP area are as follows Figure 8 As shown, the grey dots are subglacial lakes on the CReSIS survey line published by Livingstone et al. in 2022; the white dots are subglacial lakes on the CReSIS survey line published by Wolovick et al. in 2013; the dark dots are subglacial lakes calculated by the method of the present invention, and the size of the dots represents the possibility that this is a subglacial lake or subglacial water body.

[0109] Figure 8 The white line in the middle is Figure 7 At the location of the typical B-Scan example, some terrains span across multiple B-Scans. The present invention splices them and displays them, as shown in Figure (7m). The elements in each sub-image are B-Scan radar image, detection value, frequency response, frequency amplitude, slope, relative echo intensity after geometric correction, ice thickness and water potential from top to bottom.

[0110] The detection method proposed in this invention judges the measurement line of the CReSIS Antarctic dataset as follows: Fig. 9 As shown. Among them, the dark solid dots are subglacial water bodies detected by the method of the present invention. The size of the dark solid dots represents the possibility of being subglacial water bodies. The display threshold of the dots is the threshold calculated in the previous text. The present invention sets the threshold for the Antarctic region to 9. Detection values ​​exceeding the threshold will be displayed on the map. The light solid dots are 675 subglacial lakes in the entire Antarctic region summarized and sorted by Livingstone et al. The latest version of the subglacial lake list summarizes all subglacial lakes on the survey lines of various countries and projects in previous years. It has a wider range than the survey lines in the CReSIS database. Therefore, many lakes that are not on the survey lines in the CReSIS database cannot be cross-compared and verified. The comparison results between the method of the present invention and the water potential prediction are shown in the figure. Fig.10 shown.

[0111] The above is only a preferred embodiment of the present invention, and the protection scope of the present invention is not limited to the above embodiment. All technical solutions belonging to the principle of the present invention belong to the protection scope of the present invention. For those skilled in the art, some improvements and modifications made without departing from the principle of the present invention should also be regarded as the protection scope of the present invention.

Claims

1. A method for detecting polar subglacial water based on short-time Fourier transform, characterized in that: The following steps are involved: 1) Obtain radio echo sounding data; 2) Performing data fuzzy processing along the track direction on the original two-dimensional radar image to obtain a two-dimensional radar image after heading fuzzy processing; 3) performing cropping and splicing processing on the two-dimensional radar image after the heading blur processing to obtain a cropped two-dimensional radar image; 4) reshaping the waveform of the interface signal in the cropped two-dimensional radar image; 5) Perform short-time Fourier transform processing on the two-dimensional radar image after waveform reorganization to obtain the frequency response and amplitude of each measuring line at the interface; 6) Calculate the ice bottom material detection value by combining the normalized frequency response, the frequency response amplitude after short-time Fourier transform, and the terrain slope to automatically identify the subglacial water body; In the step 2), the expression of the two-dimensional radar image after heading fuzzy processing is: Wherein, BScan represents the original 2D radar image, w1 is the smoothing window width, BScan′ is the 2D radar image after heading blur processing, m represents the total number of rows of the 2D radar image, n represents the total number of columns of the 2D radar image, subscript j1 represents the row number of the 2D radar image, subscript k1 represents the column number of the 2D radar image, and subscript i1 represents the column number within the smoothing window; In the step 3), the expression for data clipping is: Where BScan″ represents the cropped two-dimensional radar image, is the sampling point corresponding to the ice bottom interface in the k2th column of the two-dimensional radar image after heading fuzzification, w2 is the half window width of data cutting, and each column signal is obtained from the k2th column. The A-Scan is cut from w2 sampling points upward and downward, and the number of sampling points of each A-Scan after cutting is 2w2+1; The step 4) specifically comprises the following steps: 41) After the cropped two-dimensional radar image, each A-Scan with a size of 2w2+1×1 is subjected to DC drift removal processing, then: Among them, AScan represents any A-Scan in the cropped two-dimensional radar image, 2w2+1 is the number of sampling points of each A-Scan, is the j2th sampling point in A-Scan, and AScan′ represents the A-Scan after removing DC drift; 42) reshaping the interface waveform in the two-dimensional radar image after the DC drift removal process; In the step 42), the waveform reshaping specifically includes: 421) setting the left and right endpoints of the main peak, and flipping the main peak and splicing them to the left and right endpoints of the main peak respectively; 422) Set the values ​​of other sampling points to 0, then: in, is the j3rd sampling point value of A-Scan after waveform reorganization, l represents the sampling point number corresponding to the position of the left endpoint of the main peak before waveform reorganization, r represents the sampling point number corresponding to the position of the right endpoint of the main peak before waveform reorganization, and Peak represents the main peak signal.

2. The method for detecting polar subglacial water based on short-time Fourier transform according to claim 1, characterized in that: In the step 1), the radio echo sounding data are specifically different coherent radar sounder data provided by the Ice Sheet Remote Sensing Center.

3. The method for detecting polar subglacial water based on short-time Fourier transform according to claim 1, characterized in that: The left and right endpoints of the main peak are defined as: One sixth of the peak intensity of the interface is set as the threshold, and the two intersection points with the threshold and the closest peak positions on the left and right sides of the A-Scan peak are respectively used as the left and right endpoints of the main peak of the interface.

4. The method for detecting polar subglacial water based on short-time Fourier transform according to claim 1, characterized in that: In the step 5), the short-time Fourier transform formula is: Among them, x(τ) is the signal to be processed, w(*) is the window function, the window function for short-time Fourier transform processing at time t is w(τ-t), and ω is the angular velocity.

5. The method for detecting polar subglacial water based on short-time Fourier transform according to claim 1, characterized in that: In the step 6), the calculation formula of the ice bottom material detection value C is: in, is the initial value of ice bottom material detection of the k3th A-Scan, is the ice bottom material detection value after smoothing of the k3th A-Scan, w3 is the smoothing window width, i2 is the column number in the smoothing window, F is the normalized frequency response, A is the amplitude of the frequency response after short-time Fourier transform, slope is the terrain slope, α is the influence weight of the terrain slope, e is the natural logarithm, Δz is the elevation difference between the main peaks of two adjacent A-Scans, and Δx is the horizontal distance between two adjacent A-Scans.

6. The method for detecting polar subglacial water based on short-time Fourier transform according to claim 5, characterized in that: In step 6), the automatic identification of subglacial water bodies is specifically performed as follows: The ice bottom material detection value C is used as the basis for the final judgment of whether the ice bottom material is a subglacial lake or subglacial water body. If the value of is greater than the set threshold, the corresponding A-Scan in the k3th column is a subglacial lake or subglacial water body, then: A-Scan for subglacial water bodies, A-Scan for subglacial sediments, A-Scan for bedrock.

Images