A method and system for ice cover thickness detection based on borehole noise imaging

Abstract

The application provides an ice cover thickness detection method and system based on borehole noise imaging, and relates to the technical field of ice cover geological detection. The method comprises the following steps: arranging a two-dimensional linear survey line composed of multiple geophones around a borehole in an ice cover, collecting continuous noise data composed of mechanical noise generated by borehole operation and auxiliary equipment and environmental background noise; pre-processing the continuous noise data, which comprises the following steps: sequentially performing the following processing on the record of each geophone: removing mean value, removing linear trend, removing strong amplitude events in the time domain, and removing strong amplitude components in the frequency domain, so as to eliminate instrument drift and local strong noise interference, and obtaining pre-processed noise data; selecting one geophone in the two-dimensional linear survey line as a virtual source position, performing mutual coherence calculation on the noise data, and recovering a pseudo-shot gather reflection response corresponding to the virtual source. The application realizes high-resolution imaging of the ice-rock interface.

Description

Technical Field

[0001] This invention relates to the field of ice sheet geological exploration technology, and in particular to a method and system for detecting ice sheet thickness based on drilling noise imaging. Background Technology

[0002] Antarctic ice sheet thickness measurement is a key parameter in studying the structural characteristics and stability of the ice sheet and assessing its impact on global climate. Accurately obtaining ice sheet thickness information is of great significance for building high-precision ice sheet models and their impact on global sea-level change.

[0003] Existing methods for ice sheet thickness detection mainly include geophysical methods such as ice radar and active source seismic exploration. Among them, ice radar is easily affected by changes in shallow snow grains and electromagnetic attenuation at depth, and its detection accuracy is limited under complex conditions. Active source seismic exploration requires special polar seismic sources, which are complex to construct and costly.

[0004] During drilling operations on the Antarctic ice sheet, the interaction between the drill bit and the ice sheet, as well as the drilling equipment, generates high-energy, broadband drilling noise. This noise propagates inside the ice sheet, is reflected at the ice-rock interface, and is received by surface geophones. It can be used for ice-rock interface imaging and ice sheet thickness detection. Based on this, the present invention proposes an ice sheet thickness detection method based on drilling noise, which can effectively detect the thickness of the Antarctic ice sheet without the need for polar-specific seismic source equipment. Summary of the Invention

[0005] The technical problem to be solved by the present invention is to provide a method and system for detecting ice cap thickness based on drilling noise imaging, so as to achieve high-resolution imaging of the ice-rock interface.

[0006] To solve the above-mentioned technical problems, the technical solution of the present invention is as follows: In a first aspect, a method for detecting ice cap thickness based on drilling noise imaging, the method comprising: A two-dimensional linear survey line consisting of multiple geophones is set up around the ice cap drilling site to collect continuous noise data consisting of mechanical noise generated by drilling operations and auxiliary equipment, as well as environmental background noise. The continuous noise data is preprocessed. The preprocessing includes removing the mean, removing the linear trend, removing strong amplitude events in the time domain, and removing strong amplitude components in the frequency domain from the records of each detector in sequence, so as to eliminate instrument drift and local strong noise interference, and obtain the preprocessed noise data. A detector in a two-dimensional linear survey line is selected as the location of a virtual seismic source. Mutual interference calculation is performed on the noise data to recover the pseudo-shot gather reflection response corresponding to the virtual seismic source. The simulated gun gathering reflection response is sequentially subjected to bandpass filtering, top cut-off, amplitude compensation, and frequency wavenumber domain filtering to extract the body wave signal in the preset frequency band, suppress surface wave energy, and enhance deep reflection signal to obtain simulated gun gathering data with improved signal-to-noise ratio. Velocity analysis was performed on the simulated shot set data with improved signal-to-noise ratio to obtain a migration velocity model for migration imaging. The pre-stack time migration method is adopted, and the simulated shot set data with improved signal-to-noise ratio is processed by migration imaging through migration velocity model to obtain time-domain stacked profiles. The time-domain overlay profile is converted to depth to obtain the depth-domain profile, and the ice sheet thickness is determined based on the spatial position of the reflection phase axis of the ice-rock interface in the depth-domain profile.

[0007] Furthermore, the continuous noise data is preprocessed. This preprocessing includes sequentially removing the mean, removing linear trends, removing strong amplitude events in the time domain, and removing strong amplitude components in the frequency domain from the records of each detector. This aims to eliminate instrument drift and localized strong noise interference, resulting in preprocessed noise data, including: The mean-reducing process is performed on each detector record in the continuous noise data to eliminate the DC component of the signal, resulting in the mean-reduced noise data. The noise data is de-linearized to remove the linear drift component in the signal, resulting in de-linearized noise data. The noise data is subjected to time-domain strong amplitude suppression to reduce the impact of sudden strong noise events, resulting in preliminary denoised noise data. The noise data after preliminary denoising is subjected to frequency domain strong amplitude component suppression processing to suppress strong interference energy in specific frequency bands, resulting in the final preprocessed noise data.

[0008] Furthermore, a detector in a two-dimensional linear survey line is selected as the location of a virtual seismic source. Inter-interference calculations are performed on the noise data to reconstruct the pseudo-shot gather reflection response corresponding to the virtual seismic source, including: Based on a two-dimensional linear survey line, a specific detector is selected from all detectors as the virtual seismic source location; The preprocessed noise data is divided into time windows to obtain multiple continuous noise data segments of equal length. For each data segment, the cross-correlation function between the record of the virtual source location geophone and the records of other geophones on the survey line is calculated using the cross-correlation algorithm to extract the phase consistency information between each geophone pair and obtain the corresponding cross-correlation results. The cross-correlation results are superimposed and averaged in the causal and anti-causal parts respectively to enhance signal energy and suppress random noise, thus obtaining a preliminary superimposed cross-correlation sequence. The initially superimposed cross-correlation sequences are symmetrically superimposed to merge causal and anti-causal responses in order to obtain the pseudo-shot gather reflection response corresponding to the virtual source location.

[0009] Furthermore, the simulated gun gathering reflection response is sequentially subjected to bandpass filtering, top cut-off, amplitude compensation, and frequency wavenumber domain filtering to extract the volume wave signal in a preset frequency band, suppress surface wave energy, and enhance deep reflection signals, thereby obtaining simulated gun gathering data with improved signal-to-noise ratio, including: Bandpass filtering is applied to the reflection response of the simulated gun assembly to extract the effective energy within a preset frequency band corresponding to the main frequency of the body wave signal, thus obtaining the bandpass-filtered simulated gun assembly. The top of the bandpass-filtered pseudo-shot collection is cut off to suppress the interference of high-energy noise in the shallow near-surface area, resulting in a pseudo-shot collection with suppressed shallow noise. Amplitude compensation processing is performed on the simulated shot set after shallow noise suppression to enhance the signal energy of the deep reflection waveform, resulting in a simulated shot set with equal amplitude. Frequency wavenumber domain filtering is applied to the simulated gun set. Based on a preset apparent velocity threshold, surface wave energy with low apparent velocity is suppressed, while volume wave reflection events with high apparent velocity are retained, in order to obtain simulated gun set data with improved signal-to-noise ratio.

[0010] Furthermore, velocity analysis is performed on the simulated shot gather data after signal-to-noise ratio enhancement to obtain a migration velocity model for migration imaging, including: On the simulated shot set data with improved signal-to-noise ratio, select common reflection point gathers that contain the main reflection events; A velocity scan is performed on the selected gathers. By calculating the superposition energy at different test velocities, the corresponding velocity that makes the reflection phase axis focus optimal is identified, and a series of discrete velocity control points are obtained. Based on the velocity control points, spatial interpolation is performed along the survey line to obtain a preliminary continuous velocity distribution; The initial continuous velocity distribution is smoothed to eliminate local outliers and maintain the rationality of velocity changes, ultimately yielding a migration velocity model for migration imaging.

[0011] Furthermore, a pre-stack time migration method is employed, and the simulated shot gather data with improved signal-to-noise ratio is processed using a migration velocity model to obtain a time-domain stacked profile, including: Acquire the migration velocity model and the simulated shot set data after signal-to-noise ratio enhancement; based on the geological structural characteristics of the imaging area and the preset migration parameters, set the maximum migration dip angle and migration aperture for pre-stack time migration; The pre-stack time migration method is adopted, and the migration velocity model is used to perform migration calculation on the simulated shot gather data after the signal-to-noise ratio is improved. The data recorded by each receiver point is returned to the underground reflection point position according to the wave field propagation law, and the migrated common reflection point gather is obtained. The common reflection point gathers are stacked, and the energies of all gathers at corresponding positions are summed to obtain the final time-domain stacked profile.

[0012] Furthermore, time-depth conversion is performed on the time-domain overlay profile to obtain the depth-domain profile, and the ice sheet thickness is determined based on the spatial position of the reflection phase axis of the ice-rock interface in the depth-domain profile, including: Obtain the time-domain overlay profile and the migration velocity model; based on the migration velocity model, convert the two-way travel time data in the time-domain overlay profile into the corresponding depth data to obtain a preliminary depth-domain profile. The preliminary depth domain profile is smoothed or interpolated to eliminate depth jumps caused by velocity model errors, thus obtaining a continuous and reliable final depth domain profile. In the final depth domain profile, identify and pick out continuous reflection in-phase axes representing the ice-rock interface; Based on the depth of the reflection phase axis at the ice-rock interface and its variation in the horizontal direction, the ice cover thickness at the corresponding location below the survey line is calculated and output.

[0013] Secondly, an ice cap thickness detection system based on drilling noise imaging includes: The acquisition module is used to deploy a two-dimensional linear survey line consisting of multiple geophones around the ice cap drilling site to collect continuous noise data composed of mechanical noise generated by drilling operations and auxiliary equipment, as well as environmental background noise. The calculation module is used to preprocess the continuous noise data. The preprocessing includes removing the mean, removing the linear trend, removing strong amplitude events in the time domain, and removing strong amplitude components in the frequency domain from the records of each detector in sequence, so as to eliminate instrument drift and local strong noise interference and obtain preprocessed noise data. A detector in a two-dimensional linear survey line is selected as the location of a virtual seismic source. Mutual interference calculation is performed on the noise data to recover the pseudo-shot gather reflection response corresponding to the virtual seismic source. The extraction module is used to sequentially perform bandpass filtering, top cut-off, amplitude compensation and frequency wavenumber domain filtering on the simulated gun set reflection response to extract the body wave signal in the preset frequency band, suppress the surface wave energy and enhance the deep reflection signal to obtain simulated gun set data with improved signal-to-noise ratio. The analysis module is used to perform velocity analysis based on the simulated shot set data after the signal-to-noise ratio is improved, and to obtain the migration velocity model for migration imaging. The processing module is used to perform migration imaging processing on the simulated shot gather data after the signal-to-noise ratio has been improved by using the pre-stack time migration method and the migration velocity model to obtain a time-domain stacked profile; the time-domain stacked profile is converted to depth to obtain a depth-domain profile, and the ice cover thickness is determined according to the spatial position of the reflection phase axis of the ice-rock interface in the depth-domain profile.

[0014] Thirdly, a computing device includes: One or more processors; A storage device for storing one or more programs that, when executed by one or more processors, cause the one or more processors to implement the method.

[0015] Fourthly, a computer-readable storage medium storing a program that, when executed by a processor, implements the method.

[0016] The above-described solution of the present invention has at least the following beneficial effects: Because this invention uses drilling noise as a passive source, it collects continuous noise data by deploying two-dimensional linear survey lines around the ice cap drilling site. Interference is eliminated through multi-level preprocessing, including mean removal and delinear trend removal. Then, the pseudo-shot gather reflection response is recovered through time windowing, mutual interference calculation, and averaging and symmetrical superposition of causal and anti-causal responses. Subsequently, the signal-to-noise ratio is optimized through bandpass filtering, top cut-off, amplitude compensation, and frequency-wavenumber domain filtering. A migration velocity model is constructed by combining common reflection point gather selection, velocity scanning, spatial interpolation, and smoothing. A pre-stack time migration method is used to complete signal repositioning and gather superposition to obtain a time-domain superimposed profile. Finally, time-depth conversion, depth-domain profile optimization, ice-rock interface phase axis picking, and... The complete technical system for thickness calculation effectively overcomes the shortcomings of existing ice radar detection, which is susceptible to changes in shallow snow grains and deep electromagnetic attenuation, and active source seismic exploration, which requires a dedicated polar source and is complex and costly to construct. It also solves the technical problems of non-stationary drilling noise, multiple interferences, difficulty in extracting effective body wave signals, inaccurate velocity modeling, blurred imaging, and insufficient accuracy in ice-rock interface positioning and thickness calculation. As a result, it achieves the technical effect of high-resolution imaging of the ice-rock interface and precise and accurate detection of ice sheet thickness below the survey line without the need for a dedicated polar source, simplifying the construction process, and reducing detection costs. This provides reliable data support for ice sheet structural characteristics research, stability analysis, and global climate impact assessment. Attached Figure Description

[0017] Figure 1 This is a schematic flowchart of an ice cover thickness detection method based on drilling noise imaging provided by an embodiment of the present invention.

[0018] Figure 2 This is a schematic diagram of an ice cover thickness detection system based on drilling noise imaging, provided by an embodiment of the present invention.

[0019] Figure 3 This invention provides a passive source volume wave imaging method for detecting ice sheet thickness.

[0020] Figure 4 This invention provides a method for restoring and processing body wave virtual seismic source records based on drilling noise, as provided in an embodiment of the present invention.

[0021] Figure 5 This is a comparison of the migration velocity model and the common reflection gather before and after pre-stack time migration processing provided by the embodiments of the present invention.

[0022] Figure 6 This is an ice sheet internal reflection superposition profile obtained by pre-stack time offset processing provided in an embodiment of the present invention.

[0023] Figure 7 This is the ice sheet depth domain reflection superposition profile obtained after time-depth conversion provided in the embodiments of the present invention. Detailed Implementation

[0024] Exemplary embodiments of the present disclosure will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure may be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete, and will fully convey the scope of the disclosure to those skilled in the art.

[0025] like Figure 1 As shown, an embodiment of the present invention proposes a method for detecting ice cap thickness based on drilling noise imaging, the method comprising the following steps: Step 1: Deploy a two-dimensional linear survey line consisting of multiple geophones around the ice cap drilling site to collect continuous noise data consisting of mechanical noise generated by drilling operations and auxiliary equipment, as well as environmental background noise. Step 2: Preprocess the continuous noise data. The preprocessing includes removing the mean, removing the linear trend, removing strong amplitude events in the time domain, and removing strong amplitude components in the frequency domain from the records of each detector in sequence, so as to eliminate instrument drift and local strong noise interference, and obtain the preprocessed noise data. Step 3: Select a detector in the two-dimensional linear survey line as the virtual source location, perform mutual interference calculation on the noise data, and recover the pseudo-shot gather reflection response corresponding to the virtual source. Step 4: The simulated gun gathering reflection response is sequentially subjected to bandpass filtering, top cut-off, amplitude compensation, and frequency wavenumber domain filtering to extract the body wave signal in the preset frequency band, suppress the surface wave energy, and enhance the deep reflection signal to obtain simulated gun gathering data with improved signal-to-noise ratio. Step 5: Perform velocity analysis based on the simulated shot set data with improved signal-to-noise ratio to obtain the migration velocity model for migration imaging; Step 6: Using the pre-stack time migration method, the simulated shot set data with improved signal-to-noise ratio is processed by migration imaging through the migration velocity model to obtain the time-domain stacked profile. Step 7: Perform time-depth conversion on the time-domain overlay profile to obtain the depth-domain profile, and determine the ice sheet thickness based on the spatial position of the reflection phase axis of the ice-rock interface in the depth-domain profile.

[0026] In this embodiment of the invention, because the invention uses two-dimensional linear survey lines to collect continuous noise data consisting of mechanical noise generated by drilling operations and auxiliary equipment and environmental background noise around the ice cap drilling site, and eliminates instrument drift and local strong noise interference through multi-level preprocessing such as mean removal and linear trend removal, and combines mutual interference calculation to recover the pseudo-shot gather reflection response, and then optimizes the signal-to-noise ratio through bandpass filtering, top cut-off, amplitude compensation, and frequency wavenumber domain filtering, combined with velocity analysis to construct a migration velocity model, pre-stack time migration imaging, and time-depth conversion, the invention effectively overcomes the technical defects of existing ice radar detection which is easily affected by shallow snow grain changes and deep electromagnetic attenuation, and active source seismic exploration which requires a special polar source and is complex and costly to construct. It solves the problems of non-stationary drilling noise, uneven energy distribution, and strong amplitude events that can easily mask effective reflection information. Thus, it achieves the technical effect of high-resolution imaging of ice-rock interface without the need for additional special polar sources, reducing detection costs and construction complexity, ensuring the stability and high accuracy of ice cap thickness detection, and completing fine detection in complex ice cap environments.

[0027] In a preferred embodiment of the present invention, step 1 above may include: Step 1: Deploy a two-dimensional linear survey line consisting of multiple geophones around the drilling site on the ice sheet to collect continuous noise data composed of mechanical noise generated by drilling operations and auxiliary equipment, as well as environmental background noise. Specifically, this includes planning and deploying a two-dimensional linear survey line in the areas on both sides of the drilling site where drilling operations have commenced on the Antarctic ice sheet. The deployed survey line consists of multiple single-component nodal seismographs, which are used as geophones. These geophones have a natural frequency of 5 Hz and are placed sequentially at a standard interval of 20 meters, ultimately resulting in a complete two-dimensional linear survey line with a total length of approximately 2,100 meters, ensuring the accuracy of the measurements. The line can fully cover the area around the well that needs to be detected. All geophones are activated and put into continuous data acquisition mode. The noise data collected mainly includes three parts: first, the strong vibration mechanical noise generated by the interaction between the drill bit and the ice sheet during the drilling operation; second, the stable vibration noise generated during the continuous operation of the drilling auxiliary equipment; and third, the environmental background noise of the Antarctic ice sheet. During the acquisition process, the sampling interval is set to four milliseconds, and the continuous acquisition time is guaranteed to be no less than twenty days, so as to obtain complete and effective continuous noise data, providing sufficient data support for ice-rock interface imaging and ice sheet thickness detection.

[0028] In this embodiment of the invention, because the invention employs a two-dimensional linear survey line composed of multiple geophones deployed around the ice cap drilling site to collect continuous noise data consisting of mechanical noise generated by drilling operations and auxiliary equipment, as well as environmental background noise, it effectively overcomes the technical difficulties of existing active source seismic exploration, which requires the use of polar-specific seismic sources, is complex to construct, and is costly. At the same time, it solves the problem of the ineffective utilization of naturally generated high-energy broadband noise during drilling. Thus, it achieves the technical effect of not needing to set up additional artificial seismic sources, simplifying on-site construction procedures, reducing detection costs, and obtaining continuous, comprehensive, and high-energy passive source noise data, providing a reliable data foundation for reflection response recovery and ice cap thickness detection.

[0029] In a preferred embodiment of the present invention, step 2 above may include: Step 2.1 involves performing a mean-reduction process on each detector record in the continuous noise data to eliminate the DC component of the signal, resulting in mean-reduced noise data. Specifically, this includes: performing a comprehensive statistical analysis on each complete signal sequence recorded by each detector in the continuous noise data, traversing every data point in the signal sequence, collecting the numerical information of all valid records, and then calculating the arithmetic mean of all data points in the sequence using an arithmetic mean algorithm, ensuring that the average calculation covers the entire sequence without omitting any data point. Each data point in the signal sequence is then subtracted from the previously calculated arithmetic mean. This point-by-point cancellation method completely eliminates the prevalent DC component in the signal. This DC component causes the signal baseline to deviate from the normal level, interfering with the subsequent identification and analysis of the valid noise signal. After this processing, the signal fluctuates smoothly around the zero baseline, ultimately yielding the mean-reduced noise data, thus eliminating the interference caused by the DC component for further signal processing.

[0030] Step 2.2 involves de-linearizing the noise data to remove linear drift components from the signal, resulting in de-trended noise data. Specifically, this includes: performing professional linear fitting analysis on the signal sequence corresponding to each detector, considering the detector's operating characteristics in polar environments and the potential for instrument zero drift to cause linear drift. Using an appropriate fitting algorithm, a linear trend curve is constructed based on the time distribution and numerical variation characteristics of the data points in the signal sequence, accurately reflecting the overall trend of the sequence. This ensures the curve closely matches the long-term drift pattern of the signal. Then, for each data point in the signal sequence, the value at the same time point on the linear trend curve is subtracted according to its corresponding time position, correcting the linear drift components in the signal point by point. This correction process removes interference from factors such as instrument zero drift, allowing the signal to more accurately reflect the propagation characteristics of drilling noise within the ice cap, ultimately yielding de-trended noise data and further improving signal quality.

[0031] Step 2.3 involves performing time-domain strong amplitude suppression processing on the noise data to reduce the impact of sudden strong noise events, resulting in preliminary denoised noise data. Specifically, this includes: focusing on the energy distribution of the signal recorded by each detector in the time dimension, monitoring the dynamic changes in signal energy in real time, and considering the characteristics of scenarios that may generate strong amplitude events during drilling operations, such as drill bit impacting ice or sudden equipment vibrations, to comprehensively judge and set a reasonable energy threshold. This threshold setting requires multiple adjustments to accurately identify data points exceeding the threshold corresponding to sudden strong noise events while avoiding misjudging normal-intensity drilling noise as interference signals. For strong amplitude data points in the signal sequence that exceed the energy threshold, a progressive smoothing technique is used to gradually reduce the energy amplitude of these abnormal data points, allowing their energy levels to slowly transition to the energy range of surrounding normal noise data. This prevents strong amplitude data points from dominating the signal processing process due to excessive energy, effectively weakening the masking effect of sudden strong noise events on the effective signal, thus obtaining preliminary denoised noise data.

[0032] Step 2.4 involves suppressing the strong amplitude components in the frequency domain of the noise data after initial denoising to suppress strong interference energy in specific frequency bands, resulting in the final pre-processed noise data. Specifically, this includes conducting a comprehensive and detailed spectrum analysis on the signals recorded by each detector for the noise data after initial denoising. By utilizing spectrum analysis technology, the time-domain signal is converted into a frequency-domain signal, clearly revealing the energy distribution of the signal in different frequency bands. Specific interference frequency bands with abnormally concentrated energy are accurately located and identified. These interference bands often originate from the continuous operation of machinery such as power generation equipment, and their strong energy severely masks the effective signals related to the propagation of body waves within the ice sheet. Based on the frequency range, energy intensity, and other characteristics of the identified interference bands, a targeted filtering scheme is designed. This scheme employs frequency-selective attenuation technology to precisely locate specific interference frequency bands and effectively attenuate the strong interference energy within them, while maximizing the preservation of the effective frequency band energy related to the propagation of body waves within the ice sheet. This ensures that the effective signal is not over-filtered. Through this targeted processing, the adverse effects of strong interference in specific frequency bands on the recovery of the reflection response are completely eliminated. Finally, purified pre-processed noise data is obtained, providing high-quality and highly reliable data support for conducting mutual interference calculations and recovering the reflection response of virtual seismic sources.

[0033] In this embodiment of the invention, because the invention employs a progressive technical approach of sequentially performing mean-reduction processing, linear trend removal processing, time-domain strong amplitude suppression processing, and frequency-domain strong amplitude component suppression processing on each detector record in continuous noise data, it effectively overcomes problems such as DC signal components, linear drift components caused by instrument zero drift, sudden strong noise events, and strong interference energy in specific frequency bands present in drilling noise. Simultaneously, it solves the technical difficulties of unstable subsequent reflection response recovery and easy masking of effective reflection information caused by these interferences. This achieves the effect of completely eliminating instrument drift and local strong noise interference, and purifying noise data, laying a reliable data foundation for mutual interference calculation to recover the pseudo-shot gather reflection response and for accurate imaging, thus ensuring the stability and accuracy of ice cover thickness detection.

[0034] In a preferred embodiment of the present invention, step 3 above may include: Step 3.1: Based on the two-dimensional linear survey line, select a specific geophone from all geophones as the virtual seismic source location. This includes: comprehensively checking and evaluating the working status and data quality of all geophones on the completed two-dimensional linear survey line. The check includes the geophone's signal reception stability, data integrity, and the presence of obvious faults or abnormal interference. Select geophones with continuous data recording, no significant signal distortion, and no severe environmental interference. Considering the propagation characteristics of drilling noise in the Antarctic ice sheet and the need for subsequent reflection response recovery, prioritize geophones located in the central region of the survey line with the best data quality as candidates. Geophones in the central region can receive noise signals from different directions around the well more evenly, which is beneficial for comprehensively capturing reflection information from the ice-rock interface. After multiple rounds of data quality verification, determine one geophone as the virtual seismic source location. For example, select the 85th geophone on the survey line. This geophone records the noise signal with the highest signal-to-noise ratio and the richest effective information, providing a reliable reference signal for accurate recovery of the simulated shot gather reflection response.

[0035] Step 3.2 involves dividing the preprocessed noise data into time windows to obtain multiple continuous noise data segments of equal length. Specifically, this includes determining a reasonable time window length based on the strong energy and wideband characteristics of drilling noise and the continuous pattern of effective reflected signals. Considering the duration of effective signals and the completeness of noise events in drilling noise, the time window length is set to thirty seconds. This length can fully encompass the noise events generated by a single drilling operation while avoiding signal mixing due to an excessively long time window or loss of effective information due to an excessively short time window. Starting from the beginning of the preprocessed noise data, the entire noise data sequence is divided into continuous and equal-length segments according to the set thirty-second duration. During the division process, seamless connection, no overlap, and no gaps are ensured between each data segment to avoid data omission or duplicate processing. Through this division method, the complete continuous noise data is decomposed into multiple continuous noise data segments with uniform structure and consistent length. Each data segment contains relatively independent and complete noise signal information, laying the foundation for subsequent segment-by-segment mutual interference calculations to improve calculation accuracy and efficiency.

[0036] Step 3.3: For each data segment, the cross-correlation algorithm is used to calculate the cross-correlation function between the record of the virtual source location geophone and the records of other geophones on the survey line, in order to extract the phase consistency information between each geophone pair and obtain the corresponding cross-correlation results. Specifically, this includes: for each segment of continuous noise data after division, clarifying the correspondence between the virtual source location geophone and all other geophones on the survey line, establishing a one-to-one geophone pairing combination, and for each pairing combination, using the cross-correlation algorithm to perform calculations. Before the calculation, the amplitude of the recorded signal of the virtual source geophone and the recorded signals of the paired other geophones is normalized to weaken strong noise. The dominant role of amplitude anomalies in the calculation results focuses the calculation on the phase relationship between signals. In the application of the cross-correlation algorithm, the phase changes of the two sets of signals are compared point by point in the time dimension to accurately capture the phase consistency feature information between them. The phase deviation caused by random noise or interference signals is eliminated. Through this calculation process, stable phase consistency information between each pair of detectors is extracted. This information can reflect the correlation characteristics of noise signals propagating inside the ice sheet and reflected by the ice-rock interface. Finally, the cross-correlation results of each pair of detectors under each data segment are obtained, providing raw data support for signal energy enhancement and noise suppression.

[0037] Step 3.4 involves superimposing and averaging the cross-correlation results in both the causal and anti-causal components to enhance signal energy and suppress random noise, resulting in a preliminary superimposed cross-correlation sequence. Specifically, for each data segment's cross-correlation result, it is first divided into a causal component and an anti-causal component based on the time logic of signal propagation. The causal component corresponds to the normal time-sequence propagation signal of the noise signal originating from the virtual seismic source location, propagating through the ice sheet, reflecting, and then being received by other detectors. The anti-causal component corresponds to the characteristics of the reverse propagation signal after the signal propagation time axis is reversed. For the cross-correlation results corresponding to all data segments, the following steps are performed: Instead of performing centralized superposition of the causal components, the causal signals of multiple data segments are aligned along the time axis and their energy is accumulated, and then the average value is calculated. Using the same method, the inverse causal components of the cross-correlation results of all data segments are also superimposed and averaged. Through this multi-data-segment superposition and averaging operation, the characteristic of the effective reflected signal being in phase is utilized, so that the energy of the effective signal is mutually enhanced during the superposition process, while random noise, due to its disordered phase, cancels each other out during the superposition process, thereby significantly enhancing the energy intensity of the effective signal and suppressing the interference of random noise, resulting in a preliminary superposition cross-correlation sequence.

[0038] Step 3.5 involves symmetrically superimposing the initially superimposed cross-correlation sequences to merge causal and anti-causal responses, thereby obtaining the pseudo-shot gather reflection response corresponding to the virtual source location. Specifically, this includes: symmetrically aligning the causal and anti-causal parts of the initially superimposed cross-correlation sequences with time zero as the center of symmetry; reversing the anti-causal part along time zero to ensure a symmetrical distribution of the reversed anti-causal part and the causal part on the time axis, guaranteeing precise phase matching of the corresponding effective reflection signals; superimposing and merging the symmetrically aligned causal part and the reversed anti-causal part to further enhance the energy of the corresponding effective reflection signals in both parts, while further canceling residual random noise and local anomalous interference signals. Through this symmetrical superimposition process, the effective reflection information contained in the causal and anti-causal responses is fully integrated, compensating for potential energy deficiencies or noise residues in a single signal, ultimately obtaining the pseudo-shot gather reflection response corresponding to the virtual source location. This response signal features concentrated energy, low noise interference, and clear reflection characteristics, accurately reflecting the reflection situation of the ice-rock interface inside the ice sheet.

[0039] In this embodiment of the invention, because it employs the following techniques—selecting specific detectors from all detectors along a two-dimensional linear survey line as virtual seismic source locations, dividing the pre-processed noise data into equal-length time windows to obtain multiple continuous noise data segments, calculating the cross-correlation function between the virtual seismic source location detector and other detectors along the survey line using a cross-coherence algorithm for each data segment to extract phase consistency information, averaging the cross-correlation results in the causal and anti-causal parts respectively, and performing symmetrical superposition processing on the initially superimposed cross-correlation sequences to merge causal and anti-causal responses—it effectively overcomes the limitations of drilling. The non-stationarity, uneven energy distribution, and dominance of high-amplitude events inherent in well noise make traditional cross-correlation methods susceptible to the dominance of high-energy events, resulting in unstable reflection response recovery results. Furthermore, random noise can easily mask effective phase information. This approach achieves the effects of accurately extracting phase consistency information between each detector pair, significantly enhancing effective signal energy, and efficiently suppressing random noise interference. It successfully obtains stable and high-quality reflection responses of the virtual source corresponding to the pseudo-shot gather, providing reliable reflection response data support for subsequent pseudo-shot gather signal optimization, velocity modeling, and migration imaging, thus ensuring the accuracy of ice cover thickness detection.

[0040] In a preferred embodiment of the present invention, step 4 above may include: Step 4.1 involves bandpass filtering the simulated shot gathering reflection response to extract the effective energy within a preset frequency band corresponding to the dominant frequency of the body wave signal, resulting in the bandpass-filtered simulated shot gathering. Specifically, this includes: performing a comprehensive spectral analysis of the simulated shot gathering reflection response; considering the strong energy and wideband characteristics of Antarctic ice sheet drilling noise, the inherent operating parameters of the detector, and the propagation characteristics of the body wave signal; clarifying the dominant frequency range of the effective body wave signal; based on the background information that the detector's inherent frequency is 5 Hz and the dominant frequency of the body wave event is below 10 Hz; and considering the potential high-frequency interference and low-frequency redundant signals in the noise, spectral analysis reveals that the simulated shot gathering reflection response is related to the reflection from the ice-rock interface. The effective body wave signal energy is mainly concentrated in the frequency band from 4 Hz to 40 Hz. Signals outside this band are mostly irrelevant interference energy. Based on this analysis, a targeted bandpass filter is designed. This filter allows signals in the 4 Hz to 40 Hz band to pass smoothly, while effectively attenuating low-frequency interference signals below 4 Hz and high-frequency noise signals above 40 Hz. The designed bandpass filter is applied to the simulated gun gathering reflection response, and the signal is filtered point by point to accurately extract the effective energy in the preset frequency band corresponding to the main frequency of the body wave signal. Interference in irrelevant frequency bands is filtered out, and finally, the simulated gun gathering after bandpass filtering is obtained, making the preliminary outline of the body wave signal clearer.

[0041] Step 4.2 involves top-cutting the bandpass-filtered simulated shot gather to suppress interference from near-surface shallow high-energy noise, resulting in a simulated shot gather with suppressed shallow noise. Specifically, this includes: firstly, a detailed analysis of the time-domain signal characteristics of the bandpass-filtered simulated shot gather. Since the noise energy in the near-surface region of the Antarctic ice sheet is relatively strong, and this type of shallow noise is mainly concentrated in the early period of signal propagation, the high-energy noise in this period will mask the subsequent deep body wave reflection signals, affecting the identification of the effective signal. By analyzing the time axis distribution of the simulated shot gather, the time range corresponding to shallow high-energy noise is determined, which is usually a short period of time after the signal initiation. Based on this time range, reasonable top cut-off parameters are set to clarify the time interval to be cut off, ensuring that this interval only covers the area where shallow noise is concentrated and does not involve the effective deep body wave signal. Top cut-off processing is applied to the bandpass-filtered simulated shot gather, and the signal amplitude within the set time interval is significantly attenuated or reduced to zero, thereby suppressing the interference of near-surface shallow high-energy noise and preventing it from masking the effective deep reflection signal. After processing, a simulated shot gather with suppressed shallow noise is obtained, and the relative intensity of the deep body wave signal is improved, laying the foundation for amplitude compensation and filtering.

[0042] Step 4.3 involves amplitude compensation processing of the simulated shot gather after shallow noise suppression to enhance the signal energy of the deep reflection waveform, resulting in an amplitude-equalized simulated shot gather. Specifically, this includes amplitude characteristic analysis of the simulated shot gather after shallow noise suppression. It was found that as the body wave signal propagates within the ice sheet, its energy gradually attenuates with increasing propagation distance and depth, leading to significantly weaker signal energy in the deep reflection waveform compared to the shallow signal. This makes it difficult to identify potential deep body wave reflection events. To address this issue, amplitude compensation processing is necessary. Based on the propagation characteristics of the ice sheet medium and the amplitude attenuation law of the simulated shot gather signal, an amplitude attenuation model is established. The model can accurately reflect the attenuation trend of signal amplitude with propagation time. Based on this model, the amplitude compensation coefficient corresponding to different propagation time periods is calculated. The longer the propagation time, the larger the compensation coefficient is to offset the influence of energy attenuation. The calculated compensation coefficient is applied point by point to each signal data point of the simulated gun set to enhance the signal energy of the deep reflection waveform in a targeted manner, while maintaining the relative stability of the shallow signal energy to avoid over-compensation leading to new interference. Through this processing, the overall amplitude distribution of the simulated gun set is more balanced, the signal energy of the deep reflection waveform is significantly improved, and potential body wave reflection events are more obvious, finally obtaining the simulated gun set with balanced amplitude.

[0043] Step 4.4 involves frequency-wavenumber domain filtering of the simulated shot gather. Based on a preset apparent velocity threshold, surface wave energy with low apparent velocities is suppressed, while volume wave reflection events with high apparent velocities are retained to obtain simulated shot gather data with improved signal-to-noise ratio. Specifically, this includes performing frequency-wavenumber domain spectral analysis on the amplitude-equalized simulated shot gather, converting the signals in the time and spatial domains into frequency-wavenumber domain signals, clearly presenting the signal energy distribution corresponding to different frequencies and wavenumbers. Spectral analysis reveals that the energy of surface wave signals is mainly concentrated in the low apparent velocity region, while the apparent velocities corresponding to volume wave reflection events are significantly higher. This is consistent with the characteristics of drilling noise propagation in the background. Combining the analysis results of actual detection data, 2500 m / s is set as the apparent velocity threshold. This threshold can effectively distinguish between surface waves with low apparent velocities and volume waves with high apparent velocities. Subsequently, a targeted filtering operator is constructed in the frequency-wavenumber domain. This operator significantly suppresses signal energy with apparent velocities below 2500 m / s while retaining signal energy with apparent velocities above this threshold to the maximum extent possible. The filtering operator is applied to the frequency-wavenumber domain signal to accurately filter out surface wave energy with low apparent velocity, thus avoiding its interference with body wave reflection events. After filtering, the signal is converted back from the frequency-wavenumber domain to the time and spatial domains, resulting in pseudo-shot gather data with a significantly improved signal-to-noise ratio. The characteristics of body wave reflection events are clearer and more continuous, providing high-quality and effective data support for velocity analysis and migration imaging.

[0044] In this embodiment of the invention, because a progressive technical approach is adopted—namely, performing bandpass filtering on the simulated gun assembly reflection response to extract the effective energy of the preset frequency band corresponding to the main frequency of the body wave signal, performing top cut-off on the bandpass-filtered simulated gun assembly to suppress high-energy noise in the shallow near-surface region, performing amplitude compensation on the simulated gun assembly after suppressing shallow noise to enhance the energy of the deep reflection waveform signal, and performing frequency-wavenumber domain filtering on the simulated gun assembly and suppressing low-apparent-velocity surface wave energy while retaining high-apparent-velocity body wave reflection events based on a preset apparent velocity threshold—it effectively overcomes the interference of irrelevant frequency bands and the masking effect of high-energy noise in the shallow near-surface region in the simulated gun assembly reflection response. The study addresses the technical challenges of identifying volume wave reflection events caused by weak effective signals, weak energy in deep reflection waveforms, and interference from low apparent velocity surface wave energy. It also solves the problems of low signal-to-noise ratio in simulated shot assemblies and difficulty in distinguishing effective volume wave signals due to drilling noise non-stationarity and uneven energy distribution. This results in the accurate extraction of effective volume wave signals related to ice cover thickness detection, a significant improvement in the signal-to-noise ratio of simulated shot assemblies data, and clearer characteristics of volume wave reflection events. This provides high-quality data support for subsequent migration velocity model construction and pre-stack time migration imaging, further ensuring the accuracy and reliability of ice-rock interface imaging and ice cover thickness detection.

[0045] In a preferred embodiment of the present invention, step 5 above may include: Step 5.1: On the simulated shot gather data with improved signal-to-noise ratio (SNR), select common reflection point gathers containing major reflection events. Specifically, this includes: performing a comprehensive temporal and spatial domain feature analysis on the simulated shot gather data with improved SNR, focusing on observing the reflection phase axis morphology of signals in the simulated shot gather, and combining the propagation law of body waves inside the Antarctic ice sheet and the characteristics of reflection signals from the ice-rock interface to identify major reflection events with concentrated energy, strong continuity, and the ability to reflect the internal structure of the ice sheet. Reflection events usually manifest as hyperbolic or approximately horizontal signal phase axes with fixed propagation laws, and their energy is significantly higher than the surrounding random noise. According to the common reflection point principle, signals from the same underground reflection point and corresponding to different shot-receiver distances in the simulated shot gather data are grouped and classified to form multiple common reflection point gathers. During the classification process, the location of the underground reflection point corresponding to each signal is accurately matched to ensure that all signals in the same gather correspond to the same underground reflection interface. At the same time, invalid gathers containing messy interference signals and unclear reflection phase axes are removed. Finally, common reflection point gathers containing only major reflection events and with high signal quality are selected to provide accurate target data for subsequent velocity analysis.

[0046] Step 5.2 involves velocity scanning of the selected gathers. By calculating the superposition energy at different test velocities, the corresponding velocities that optimally focus the reflection phase axis are identified, resulting in a series of discrete velocity control points. Specifically, for each common reflection point gather, a reasonable test velocity range is determined by combining the physical properties of the Antarctic ice sheet medium, the apparent velocity threshold of the body wave determined in the previous FK filtering, and relevant geological exploration experience. This range needs to cover the possible propagation velocity range of the ice sheet body wave, avoiding both omission of effective velocity values ​​due to an overly narrow range and increased computational load and introduction of irrelevant interference due to an overly wide range. Within the set test velocity range, a series of continuous test velocity values ​​are divided according to fixed velocity intervals. For each test velocity, all signals in the corresponding common reflection point gather are processed according to that velocity. Dynamic correction is performed to eliminate the propagation time differences caused by different shot-receiver distances, making the originally dispersed reflection phase axes tend to align. The superposition energy of all signals in the gather after dynamic correction is calculated. The magnitude of the superposition energy directly reflects the focusing degree of the reflection phase axes. The higher the energy, the closer the test speed is to the actual propagation speed of the signal, and the better the focusing effect of the reflection phase axes. The above dynamic correction and superposition energy calculation operations are repeated for all test speeds. By comparing the peak superposition energy under different test speeds, the corresponding speed value that makes the reflection phase axes focus best is accurately identified. This speed value is used as the effective speed parameter corresponding to the common reflection point gather. This operation is performed on all selected common reflection point gathers, and finally a series of discrete speed control points corresponding to different underground reflection point locations are obtained.

[0047] Step 5.3: Based on velocity control points, perform spatial interpolation along the survey line to obtain a preliminary continuous velocity distribution. Specifically, this includes: based on discrete velocity control points, determining the horizontal coordinates of the corresponding subsurface co-reflection points on the two-dimensional linear survey line, and establishing a one-to-one correspondence between velocity values ​​and the horizontal positions of the survey line. Considering the continuity of the Antarctic ice sheet medium, the propagation velocity of its internal body waves typically exhibits a gradual trend and does not show irregular abrupt changes. Therefore, a spatial interpolation method suitable for continuous media is selected. This type of method reasonably infers the velocity value at any position between two points based on the velocity information of adjacent control points. To ensure the continuity and rationality of velocity changes, interpolation calculations are performed point-by-point in the blank areas between discrete velocity control points along the horizontal direction of the two-dimensional linear survey line, using the CMP point spacing of the survey line as the unit. During the interpolation process, the magnitude and trend of adjacent velocity control points are fully considered to ensure that the interpolated velocity values ​​can transition smoothly, neither deviating from the overall velocity distribution pattern nor failing to accurately reflect the local velocity differences in the ice cover medium. Through this interpolation operation, the originally discrete velocity control points are transformed into continuously distributed velocity data covering the entire survey line range, obtaining a preliminary continuous velocity distribution and providing a foundation for smoothing processing.

[0048] Step 5.4 involves smoothing the initial continuous velocity distribution to eliminate local outliers and maintain the reasonableness of velocity variations, ultimately obtaining the migration velocity model for migration imaging. This process includes: outlier detection and identification of the initial continuous velocity distribution; statistical analysis of the overall range, mean, and standard deviation of the velocity distribution to identify local anomalous velocity values ​​exceeding the reasonable velocity range. These outliers are often caused by calculation errors during velocity scanning, data deviations during interpolation, or signal interference from individual gathers. Failure to address these outliers will affect the accuracy of subsequent migration imaging. For the identified local outliers, local correction or replacement methods are used, i.e., replacing them with the mean or interpolation result of adjacent reasonable velocity values. To avoid outliers interfering with the overall velocity distribution, a smoothing algorithm suitable for the velocity characteristics of the ice sheet medium is selected, and a reasonable smoothing window size is set. Along the two-dimensional linear survey line, the initial continuous velocity distribution is smoothed point-by-point. By averaging or weighting the velocity values ​​within the window, local random fluctuations are further weakened, making the velocity distribution smoother and more continuous. After smoothing, the rationality of the velocity distribution is verified again: whether the velocity range conforms to the physical laws of ice sheet bulk wave propagation, and whether the velocity changes are consistent with the gradual characteristics of the ice sheet's internal structure, ensuring no obvious velocity jumps or unreasonable fluctuations. Finally, a stable, reliable migration velocity model for migration imaging that conforms to the propagation characteristics of the Antarctic ice sheet medium is obtained.

[0049] In this embodiment of the invention, a progressive technical approach is adopted. This approach involves selecting common reflection point gathers containing major reflection events from the simulated shot gather data after signal-to-noise ratio improvement, performing velocity scanning on the selected gathers, identifying the optimal corresponding velocity for focusing the reflection phase axis by calculating the superposition energy at different test velocities to obtain discrete velocity control points, performing spatial interpolation along the survey line based on the velocity control points to obtain a preliminary continuous velocity distribution, and smoothing the preliminary continuous velocity distribution to eliminate local outliers and maintain the rationality of velocity changes. This approach effectively overcomes the technical problems of complex velocity distribution in the ice sheet medium, the difficulty of accurately capturing its variation patterns in traditional velocity modeling, the dispersion of velocity information corresponding to effective reflection events in the simulated shot gather, and the tendency of local outliers to distort the velocity model and thus affect the accuracy of migration imaging. It also solves the problem of not being able to accurately locate reflection signals due to inaccurate velocity models. This achieves the effect of constructing an accurate, continuous, and reasonable migration velocity model, providing a reliable velocity foundation for pre-stack time migration imaging, ensuring that the reflection phase axis can be accurately focused, and further improving the clarity of ice-rock interface imaging and the accuracy of ice sheet thickness detection.

[0050] In a preferred embodiment of the present invention, step 6 above may include: Step 6.1: Obtain the migration velocity model and the signal-to-noise ratio (SNR) enhanced simulated shot gather data. Based on the geological structural characteristics of the imaging area and the preset migration parameters, set the maximum migration dip angle and migration aperture for pre-stack time migration. Specifically, this includes: accurately retrieving two types of core data from the preprocessing workflow: the completed migration velocity model, which has undergone velocity scanning, spatial interpolation, and smoothing to accurately reflect the distribution of body wave propagation velocity in the Antarctic ice sheet medium; and the optimized SNR enhanced simulated shot gather data, which has eliminated irrelevant interference, enhanced body wave signals, and possesses a high-quality reflection information foundation. Subsequently, a comprehensive analysis of the geological structure of the imaging area was conducted. The internal structure of the Antarctic ice sheet is relatively homogeneous, and the ice-rock interface generally exhibits a gentle undulation without severe folds or faults, with only minor local topographical changes. Based on this geological characteristic, and referring to engineering experience in passive source seismic detection and previous experimental data, key parameters for pre-stack time migration were set: Considering the gentle characteristics of the ice-rock interface, the maximum migration dip angle was set to 40 degrees. This value can cover the maximum possible undulation angle of the interface while avoiding the introduction of irrelevant interference energy due to an excessively large dip angle. For the 300-millisecond time window where the main reflection events are concentrated, the propagation range of the effective reflection signal was calculated by combining the propagation speed of body waves in the ice sheet. The migration aperture was set to 500 meters to ensure that the aperture can completely cover the propagation path of the main reflection events without missing any effective reflection energy, while avoiding the inclusion of additional noise interference due to an excessively large aperture.

[0051] Step 6.2 employs the pre-stack time migration method, using a migration velocity model to perform migration calculations on the signal-to-noise ratio (SNR) improved simulated shot gather data. Data recorded at each receiver point is repositioned to the underground reflection point location based on wavefield propagation laws, resulting in the migrated common reflection point gather. Specifically, this includes: defining the effective data range for this migration calculation: based on previous simulated shot gather analysis results, a data segment with high body wave SNR between zero meters and 700 meters is selected. This data segment has concentrated body wave reflection signal energy and less interference, providing reliable input for accurate imaging. The Kirchhoff pre-stack time migration method adapted to passive source seismic data is adopted, using the acquired migration velocity model as the core calculation basis, incorporating the body wave propagation characteristics of the Antarctic ice sheet medium, and targeting the SNR improved simulated shot gather data. For each geophone record, based on the geometric laws and dynamic characteristics of wave field propagation, and combined with the body wave propagation velocity at the corresponding position in the migration velocity model, the coordinates of the underground true reflection point corresponding to the recorded signal are calculated in reverse. That is, by calculating the propagation time and path of the signal from the underground reflection point to the geophone, the reflected signals that were originally scattered at different geophones due to differences in propagation paths are accurately returned to their corresponding underground reflection point positions. During the migration calculation process, the migration velocity parameters are iteratively adjusted multiple times. After each adjustment, the shape of the reflection phase axis in the common reflection point gather is observed until the reflection phase axis tends to be horizontal in the time-distance domain, ensuring the accuracy of the reflection signal return. Finally, the migrated common reflection point gather is obtained, and the effective reflected signal in the gather has been accurately focused, while the messy interference signals are suppressed.

[0052] Step 6.3 involves stacking the common reflection point gathers, summing the energy at corresponding locations in all gathers to obtain the final time-domain stacked profile. Specifically, this includes: preprocessing the offset common reflection point gathers by applying a stretching cut operation to remove signal stretching distortion caused by differences in shot-receiver distances or minor deviations in the velocity model during the offset calculation, thus preventing such distorted signals from affecting the stacking effect; and performing gather stacking processing: according to the spatial location of the underground reflection points, summing the energy of signals corresponding to the same reflection location in all common reflection point gathers point-by-point; for each time point in the time domain and each corresponding location in the spatial domain, summing the signal amplitude at that location in all gathers; and enhancing the intensity of the effective reflection signal through energy accumulation. Residual random noise, due to its chaotic phase, cancels each other out during accumulation, further reducing noise interference. During the stacking process, the spatial location of the gathers is strictly aligned with the time scale to ensure accurate stacking and energy convergence of signals from the same reflection interface. Finally, through this superposition process, a time-domain superimposed profile with clear reflection interface features and strong continuity is obtained. The reflection phase axis of the ice-rock interface and other possible structural interfaces inside the ice sheet can be clearly identified in the profile, providing a high-quality imaging basis for time-depth conversion and ice sheet thickness calculation.

[0053] In this embodiment of the invention, by first acquiring the migration velocity model and the simulated shot gather data with improved signal-to-noise ratio, combining the geological structural features of the imaging area and preset migration parameters to set the maximum migration dip angle and migration aperture for pre-stack time migration, and then using the pre-stack time migration method to relocate the recorded data of each receiver point to the corresponding underground reflection point position according to the wave field propagation law, and finally performing energy summation and superposition processing on the common reflection point gathers, the technical problems of the complex geological structure of the Antarctic ice sheet leading to the dispersion of the reflection signal propagation path, the difficulty in accurately matching the recorded data of the receiver point with the underground reflection point, and the easy occurrence of imaging blur and poor continuity of the reflection interface due to improper migration parameter settings are effectively overcome. At the same time, it solves the problem that traditional migration methods cannot effectively integrate dispersed reflection energy and are difficult to highlight the characteristics of the ice-rock interface, thereby achieving the effect of accurately relocating the reflection signal and effectively converging the reflection energy, obtaining a clear and continuous time-domain superimposed profile of the reflection interface, providing a high-quality imaging foundation for subsequent time-depth conversion, ice-rock interface identification, and ice sheet thickness calculation, and further ensuring the accuracy and reliability of ice sheet thickness detection.

[0054] In a preferred embodiment of the present invention, step 7 above may include: Step 7.1: Obtain the time-domain stacked profile and migration velocity model. Based on the migration velocity model, convert the two-way travel time data in the time-domain stacked profile into corresponding depth data to obtain a preliminary depth-domain profile. Specifically, this includes: obtaining two types of core data from the pre-processing workflow: first, the time-domain stacked profile, which has undergone pre-stack time migration and stacking processing, clearly presenting the temporal distribution characteristics of the internal reflection interface of the ice sheet, with continuous and concentrated reflection phase axes; second, the migration velocity model, which has undergone velocity scanning, spatial interpolation, and smoothing processing to realistically reflect the volume wave propagation velocity law of the Antarctic ice sheet medium, with continuous and reasonable velocity changes. Based on the basic principle of time-depth conversion, combined with the Antarctic ice... Based on the characteristics of the ice sheet as a continuous and homogeneous medium, and using the migration velocity model as the basis for calculation, a correspondence between two-way travel time and depth is established. For each common reflection point in the time-domain overlay profile, the corresponding two-way travel time data is extracted at each time scale. According to the body wave propagation velocity corresponding to the CMP point in the migration velocity model, the two-way travel time data at each time point is converted into the corresponding subsurface depth data point by point through the calculation logic of multiplying the two-way travel time and velocity and then dividing by two. During the conversion process, the one-to-one correspondence between the spatial position of the time-domain overlay profile and the depth data is strictly maintained to ensure that the spatial position of each reflection phase axis is accurately converted. Finally, a preliminary depth-domain profile that can initially reflect the depth distribution of the internal structure of the ice sheet is obtained.

[0055] Step 7.2 involves smoothing or interpolating the preliminary depth domain profile to eliminate depth jumps caused by velocity model errors, resulting in a continuous and reliable final depth domain profile. Specifically, this includes: firstly, performing a depth continuity analysis on the preliminary depth domain profile, comparing the depth data between adjacent common reflection points one by one to identify depth jumps caused by minor errors in the migration velocity model, time-depth conversion calculation deviations, etc. These jumps manifest as depth values ​​at adjacent CMP points at the same reflection interface exceeding a reasonable range, inconsistent with the gently undulating geological characteristics of the Antarctic ice sheet's ice-rock interface. To eliminate these anomalies, an appropriate processing method is selected based on the distribution characteristics of the depth data: if the depth... The depth jumps are relatively scattered and small in amplitude. A moving average method is used for smoothing, and a reasonable smoothing window is set. The average depth data within the window is taken to weaken local jumps and make the depth changes more gradual. If the depth jumps are concentrated in local areas and the data is missing, a linear interpolation method is used. Based on the reasonable depth data on both sides of the jump area, the depth values ​​of each CMP point in the jump area are calculated to fill the data gaps. Throughout the process, the geological structure characteristics of the Antarctic ice sheet are always used as constraints to ensure that the depth changes conform to the natural gradual change of the ice-rock interface and do not destroy the true reflection interface morphology. Finally, a depth domain profile that is continuous in depth, without abnormal jumps, and reflects the depth structure inside the ice sheet is obtained.

[0056] Step 7.3: In the final depth domain profile, identify and pick out continuous reflection phase axes representing the ice-rock interface. Specifically, this includes: conducting a comprehensive characteristic analysis of the reflection phase axes in the final depth domain profile; combining the structural composition of the Antarctic ice sheet and the propagation law of reflected signals, clarifying the typical characteristics of the reflection phase axes at the ice-rock interface: compared to the reflection phase axes at the snow-ice interface, the ice-rock interface, as a rigid interface between the ice sheet and the bedrock, has stronger reflected signal energy, higher amplitude, and better continuity of the phase axis, which can be stably presented at most CMP points along the survey line. Based on these characteristics, identify and trace the reflection phase axes point by point in the final depth domain profile. First, locate the strong reflection phase axis suspected to be the ice-rock interface in the area where the reflected energy is concentrated and the phase axis is clear. Then, continuously track the phase axis along the horizontal direction of the survey line to ensure its continuity throughout the entire survey line. For areas with weaker local energy but still identifiable, make reasonable supplements by combining the phase axis positions of adjacent CMP points to avoid phase axis breakage due to local energy fluctuations. Eliminate interference from other secondary reflection phase axes such as the snow-ice interface. By comparing the energy intensity, continuity and depth distribution range of the phase axis, finally determine the unique continuous reflection phase axis that represents the ice-rock interface and complete the phase axis picking work throughout the entire survey line.

[0057] Step 7.4: Based on the depth of the reflection phase axis at the ice-rock interface and its variation in the horizontal direction, calculate and output the ice cap thickness at the corresponding location below the survey line. Specifically, this includes: defining ice cap thickness as the vertical distance from the ice cap surface to the ice-rock interface; since the survey line is laid on the ice cap surface, the surface depth is considered zero; therefore, the ice cap thickness corresponding to each co-reflection point is the depth value of the reflection phase axis at the ice-rock interface picked up at that point; compiling the ice-rock interface depth data for all CMP points along the entire survey line; and establishing the correspondence between depth values ​​and the horizontal position of the CMP points: for example, a drilling location at a horizontal distance of 1070 meters... At point CMP1, the ice-rock interface depth is 540 meters, which represents the ice sheet thickness at that location. At point CMP1, the ice-rock interface depth is 750 meters, corresponding to an ice sheet thickness of 750 meters. At point CMP199, the ice-rock interface depth is 650 meters, corresponding to an ice sheet thickness of 650 meters. During the data processing, the data were validated for reasonableness to ensure that the ice sheet thickness at each location conforms to the overall topographic features of the Antarctic ice sheet. Abnormal data were removed, and the ice sheet thickness data at each CMP point were output in a visual format, clearly showing the horizontal variation pattern of the ice sheet thickness below the survey line, thus achieving precise detection and data presentation of the Antarctic ice sheet thickness.

[0058] In this embodiment of the invention, a progressive technique is employed. First, a time-domain overlay profile and a migration velocity model are acquired. Based on this model, two-way travel time data are converted into depth data to obtain a preliminary depth-domain profile. Then, the preliminary depth-domain profile is smoothed or interpolated to eliminate depth jumps. Subsequently, continuous reflection phase axes representing the ice-rock interface are identified and picked up in the final depth-domain profile. Finally, the ice cover thickness at the corresponding location is calculated and output based on the depth and horizontal changes of these phase axes. This technique effectively overcomes the technical problems of time-domain data not directly reflecting the true underground depth, velocity model errors easily causing depth jumps, difficulty in accurately identifying the ice-rock interface in complex noise environments, and the inability of traditional methods to achieve refined calculation of ice cover thickness. It also solves the problem of insufficient ice cover thickness detection accuracy caused by discontinuous depth data and inaccurate interface positioning. Thus, it achieves the technical effects of obtaining continuous and reliable depth-domain profiles, accurately locking the ice-rock interface location, and realizing refined calculation and output of ice cover thickness at various locations below the survey line. This ensures that the detection results highly match the actual ice cover thickness, providing accurate and reliable thickness data support for ice cover structure research and global climate impact assessment.

[0059] like Figure 2 As shown, embodiments of the present invention also provide an ice cap thickness detection system based on drilling noise imaging, comprising: The acquisition module is used to deploy a two-dimensional linear survey line consisting of multiple geophones around the ice cap drilling site to collect continuous noise data composed of mechanical noise generated by drilling operations and auxiliary equipment, as well as environmental background noise. The calculation module is used to preprocess the continuous noise data. The preprocessing includes removing the mean, removing the linear trend, removing strong amplitude events in the time domain, and removing strong amplitude components in the frequency domain from the records of each detector in sequence, so as to eliminate instrument drift and local strong noise interference, and obtain the preprocessed noise data. Then, a detector in the two-dimensional linear survey line is selected as the virtual source location, and the noise data is subjected to mutual interference calculation to recover the pseudo-shot gather reflection response corresponding to the virtual source. The extraction module is used to sequentially perform bandpass filtering, top cut-off, amplitude compensation and frequency wavenumber domain filtering on the simulated gun set reflection response to extract the body wave signal in the preset frequency band, suppress the surface wave energy and enhance the deep reflection signal to obtain simulated gun set data with improved signal-to-noise ratio. The analysis module is used to perform velocity analysis based on the simulated shot set data after the signal-to-noise ratio is improved, and to obtain the migration velocity model for migration imaging. The processing module is used to perform migration imaging processing on the simulated shot gather data after the signal-to-noise ratio has been improved by using the pre-stack time migration method and the migration velocity model to obtain a time-domain stacked profile; the time-domain stacked profile is converted to depth to obtain a depth-domain profile, and the ice cover thickness is determined according to the spatial position of the reflection phase axis of the ice-rock interface in the depth-domain profile.

[0060] In a preferred embodiment of the present invention, drilling noise imaging is based on the passive source reflected wave seismic interferometry theory. This theory eliminates propagation effects on the common propagation path and recovers the reflection response between receiving points by performing autocorrelation or cross-correlation operations on the wave field generated by the noise source recorded at the receiving point. In an ideal one-dimensional medium model, the medium structure can be simplified as a free surface, a single layer of isotropic and lossless medium overlying a uniform half-space. When a noise source exists in the lower half-space and generates an upward-propagating plane wave, the wave field propagates in the medium and is reflected at the free surface and layer interface. Under this condition, by performing autocorrelation operations on the transmitted wave field received at the surface, the reflection response received at the surface of the surface-excited source in this medium can be recovered. The one-dimensional reflected wave seismic interferometry theory can be written as follows: ,in The reflection response received at the surface of the earthquake source is generated by the surface. The autocorrelation of the source function. The above shows that, in the one-dimensional case, autocorrelation of the transmitted response received at the surface can recover the reflected response received at the surface from the surface-excited seismic source. In the two-dimensional and three-dimensional cases, the wave propagation path is more complex, but the reflected response can still be approximately recovered. When the noise sources are spatially uniformly distributed and the source functions are uncorrelated, cross-correlation of the wave fields received at different locations on the surface can eliminate common propagation paths, thereby approximately recovering the reflection response between receiving points. The seismic interferometry theory of reflected waves under two-dimensional and three-dimensional conditions can be written as follows: ; in for Excitation by pulse source The reflected response received by the detector, Represents the Dirac function, The autocorrelation of the source function, and For respectively and The horizontal component of the detector, and for and The above formula shows that, under the condition that the two-dimensional and three-dimensional noise sources satisfy spatial random distribution, cross-correlation of the records received by the detectors can approximately recover the reflection response between the detectors. In actual drilling environments, noise signals have complex sources, including both strong vibrations generated during drilling operations and stable vibration signals generated by the continuous operation of machinery such as power generation equipment. Noise has obvious characteristics such as non-stationarity, uneven energy distribution, and dominance of strong amplitude events. Directly using the cross-correlation method is easily dominated by strong energy events, leading to unstable reflection response recovery results and even masking effective reflection information. To recover the reflection signal more stably, mutual coherent seismic interferometry is used: ; in, and For located and The record received by the detector For the mutual interference results of the two records, As a regularization parameter, the mutual coherence method normalizes the signal amplitude information before correlation operations, ensuring that the correlation results primarily reflect the phase consistency between signals, while weakening the influence of amplitude magnitude on the results. This approach effectively suppresses interference from strong impulse events and local anomalous energy on the interferometric results, improving the adaptability of noise interferometric imaging to different types of noise sources.

[0061] like Figure 3As shown, embodiments of the present invention also provide the application of dynamic source reflection wave seismic interferometry to the drilling scenario of the Antarctic ice sheet. The drilling noise generated during the interaction between the drill bit and the ice body, as well as the mechanical vibration generated by the drilling auxiliary equipment, can be regarded as a passive noise source. This noise propagates inside the ice sheet and is reflected at the ice bottom interface and the internal structural interface. By deploying geophones around the drilling site to obtain continuous noise records, the thickness of the ice sheet below can be detected by passive source reflection imaging technology. Since the drilling noise signal has the characteristics of non-stationarity, uneven frequency band distribution and multi-path propagation, the original signal needs to be preprocessed before the seismic interferometry recovers the volume wave response. This includes operations such as removing the mean, removing the trend, and removing strong amplitudes to enhance the reflected wave energy, suppress strong energy surface waves and random noise interference, thereby improving the stability and imaging quality of the subsequent imaging results.

[0062] After preprocessing, drilling noise signals recorded at different receiving locations or within different time windows are subjected to mutual interference operations to recover the reflection response of the virtual seismic source. After denoising, deconvolution, velocity analysis, migration, and time-depth conversion, a depth domain imaging profile of the ice-rock interface is obtained, thereby realizing the imaging and detection of ice sheet thickness and ice-bottom interface structure. To verify the feasibility and operability of the above-mentioned drilling noise-based ice thickness detection method and its data processing flow in the polar ice sheet environment, an actual drilling operation area on the northwestern edge of Princess Elizabeth Land in East Antarctica was selected as the application scenario. Passive source seismic survey lines were deployed in this area to collect continuous noise data generated by drilling operations and environmental background. Following the aforementioned data preprocessing, noise signal extraction, correlation calculation, and imaging processing flow, the snow-ice interface and ice-rock interface structure below the drilling site were detected and imaged for analysis, in order to verify the application effect of the method under complex ice sheet conditions.

[0063] In this application example, during the 2023–2024 Antarctic expedition season, a drilling operation was conducted in the northwestern margin of Princess Elizabeth Land in East Antarctica, penetrating the ice sheet to the bedrock, while simultaneously carrying out passive source seismic surveys based on drilling noise. The drilling site was located approximately 28 km south of the coast of Prydz Bay, with geographical coordinates of 69.585591°S, 76.385165°E, an elevation of approximately 680 m, and approximately 23.5 km from Zhongshan Station. Drilling was completed between January and February 2024, penetrating the ice layer and reaching the bedrock at a depth of 541.12 m. A two-dimensional linear passive source seismic survey line was deployed around both sides of the drilling site, consisting of 100 single-component nodal seismographs with a natural frequency of 5 Hz and a conventional spacing of 20 m, for a total survey line length of approximately 2100 m. Passive source seismic data were continuously acquired for 20 days with a sampling interval of 4 ms. The recorded signals mainly included mechanical noise generated by drilling operations and auxiliary equipment operation, as well as environmental background noise. Based on the acquired data, the ice-rock interface below the well was imaged. Before performing passive source seismic interferometry, the acquired passive source noise data was preprocessed. The data from each geophone were processed to remove the mean and linear trend to eliminate the influence of instrument zero drift and long-term trends. Subsequently, the data were processed to remove strong amplitudes in the time domain and frequency domain to reduce the impact of local strong noise events on seismic interferometry calculations.

[0064] like Figure 4 As shown, after preprocessing, the 85th detector was selected as the virtual source location. Seismic interferometry was used to recover the body wave signal from the noise data of the first day. The pseudo-shot set was recovered from the noise record using the mutual coherence method. Causal and anti-causal relationships were superimposed with a mutual coherence time window of 30 s to enhance the signal energy. Significant low-frequency energy was observed in the recovered pseudo-shot set in the time period of 0 to 250 ms.

[0065] Considering the instrument's natural frequency of approximately 5 Hz and the dominant frequency of body wave events below 10 Hz, a 4–40 Hz bandpass filter was used to process the simulated shot gather. A clear hyperbolic body wave event was observed at approximately 370 ms. To further reduce shallow noise interference, a top cut was applied to the filtered simulated shot gather to reduce the impact of high-energy near-surface noise. After processing, amplitude compensation was performed to enhance the signal energy of the deep reflection waveform, making the potential body wave energy more apparent. Based on this, f–k analysis and filtering were performed to further suppress strong surface wave components. The f–k spectrum analysis of the simulated shot gather showed that its energy was mainly concentrated at low apparent velocities of surface waves. During the f–k filtering stage, an apparent velocity limit of 2500 m / s was set, retaining only energy above this apparent velocity, significantly weakening the surface wave component while preserving continuous body wave reflection events. After processing, the f–k filtered simulated shot gather was obtained, exhibiting a high body wave signal-to-noise ratio, suitable for velocity analysis and migration imaging.

[0066] like Figure 5 As shown, velocity analysis is performed on the virtual source shot set data after denoising and amplitude compensation. By scanning and analyzing the velocity of reflection events, the corresponding migration velocity is picked up, and the picked velocity is spatially interpolated and smoothed to obtain a time migration velocity model for imaging.

[0067] Based on the obtained time migration velocity model, data with high volume wave signal-to-noise ratio within the migration range of 0–700 m were selected. The Kirchhoff pre-stack time migration method was used to process the simulated shot gather data. During the migration, the maximum migration tilt angle was set to 40°. For the 300 ms time window where the main reflection events are located, a migration aperture range of 500 m was selected. By iteratively adjusting the migration velocity parameters multiple times, the common reflection channel concentrated reflection phase axis obtained by pre-stack time migration tends to be horizontal in the time-distance domain to obtain stable imaging results.

[0068] like Figure 6 As shown, after stretching and stacking the migrated gathers, a migration stacked profile is obtained. The reflection continuity on the stacked profile is good. Two clear and relatively continuous reflection interfaces are visible at CMP1–81 at 300ms and 400ms. A very clear and continuous reflection interface is visible near 350ms at CMP81–199. The clear and continuous reflection interfaces indicate that the migration parameters and velocity are appropriate, providing a reliable root mean square velocity model for time-depth conversion.

[0069] Will Figure 6 The time-domain profile is transformed into a depth-domain profile as follows: Figure 7 As shown, two clear and relatively continuous reflection interfaces can be identified within a depth range of 450–850m between CMP50–90 and CMP50. The stronger in-phase axis in this range is the interface between the ice sheet and the bedrock. The drilling location at a horizontal distance of 1070m is located at the highest point of the bedrock, where the ice sheet thickness is approximately 540m. The ice thickness at CMP=1 is approximately 750m, and at CMP=199 it is approximately 650m, enabling precise detection of the thickness of the ice sheet below.

[0070] In this application example, the ice thickness below the drilling point is detected based on passive source volume wave imaging results. The results are compared and analyzed with ice thickness data obtained from actual drilling. The results show that by utilizing environmental noise and passive source signals generated during drilling operations, the reflection response of the ice-bedrock interface below the drilling point can be effectively obtained without artificial seismic sources, achieving reliable ice thickness detection. This method provides stable imaging, has low requirements for on-site construction conditions, and can perform detailed ice thickness detection near the drilling site without significantly increasing equipment investment and personnel workload.

[0071] The above description represents the preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A method for detecting ice cap thickness based on drilling noise imaging, characterized in that, The method includes: A two-dimensional linear survey line consisting of multiple geophones is set up around the ice cap drilling site to collect continuous noise data consisting of mechanical noise generated by drilling operations and auxiliary equipment, as well as environmental background noise. The continuous noise data is preprocessed. The preprocessing includes removing the mean, removing the linear trend, removing strong amplitude events in the time domain, and removing strong amplitude components in the frequency domain from the records of each detector in sequence, so as to eliminate instrument drift and local strong noise interference, and obtain the preprocessed noise data. A detector in a two-dimensional linear survey line is selected as the location of a virtual seismic source. Mutual interference calculation is performed on the noise data to recover the pseudo-shot gather reflection response corresponding to the virtual seismic source. The simulated gun gathering reflection response is sequentially subjected to bandpass filtering, top cut-off, amplitude compensation, and frequency wavenumber domain filtering to extract the body wave signal in the preset frequency band, suppress surface wave energy, and enhance deep reflection signal to obtain simulated gun gathering data with improved signal-to-noise ratio. Velocity analysis was performed on the simulated shot set data with improved signal-to-noise ratio to obtain a migration velocity model for migration imaging. The pre-stack time migration method is adopted, and the simulated shot set data with improved signal-to-noise ratio is processed by migration imaging through migration velocity model to obtain time-domain stacked profiles. The time-domain overlay profile is converted to depth to obtain the depth-domain profile, and the ice sheet thickness is determined based on the spatial position of the reflection phase axis of the ice-rock interface in the depth-domain profile.

2. The ice cap thickness detection method based on drilling noise imaging according to claim 1, characterized in that, The continuous noise data is preprocessed, including removing the mean, linear trend, strong time-domain amplitude events, and strong frequency-domain amplitude components from the records of each detector, in order to eliminate instrument drift and local strong noise interference, resulting in preprocessed noise data, including: The mean-reducing process is performed on each detector record in the continuous noise data to eliminate the DC component of the signal, resulting in the mean-reduced noise data. The noise data is de-linearized to remove the linear drift component in the signal, resulting in de-linearized noise data. The noise data is subjected to time-domain strong amplitude suppression to reduce the impact of sudden strong noise events, resulting in preliminary denoised noise data. The noise data after preliminary denoising is subjected to frequency domain strong amplitude component suppression processing to suppress strong interference energy in specific frequency bands, resulting in the final preprocessed noise data.

3. The ice cap thickness detection method based on drilling noise imaging according to claim 2, characterized in that, A geophone in a two-dimensional linear survey line is selected as the location of a virtual seismic source. Inter-interference calculations are performed on the noise data to reconstruct the pseudo-shot gather reflection response corresponding to the virtual seismic source, including: Based on a two-dimensional linear survey line, a specific detector is selected from all detectors as the virtual seismic source location; The preprocessed noise data is divided into time windows to obtain multiple continuous noise data segments of equal length. For each data segment, the cross-correlation algorithm is used to calculate the Green's function between the record of the virtual source location geophone and the records of other geophones on the survey line, so as to extract the phase consistency information between each geophone pair and obtain the corresponding cross-correlation results. The cross-correlation results are superimposed and averaged in the causal and anti-causal parts respectively to enhance signal energy and suppress random noise, thus obtaining a preliminary superimposed cross-correlation sequence. The initially superimposed cross-correlation sequences are symmetrically superimposed to merge causal and anti-causal responses in order to obtain the pseudo-shot gather reflection response corresponding to the virtual source location.

4. The ice cap thickness detection method based on drilling noise imaging according to claim 3, characterized in that, The simulated gun gathering reflection response is sequentially subjected to bandpass filtering, top cut-off, amplitude compensation, and frequency-wavenumber domain filtering to extract the body wave signal in a preset frequency band, suppress surface wave energy, and enhance deep reflection signals, thereby obtaining simulated gun gathering data with improved signal-to-noise ratio, including: Bandpass filtering is applied to the reflection response of the simulated gun assembly to extract the effective energy within a preset frequency band corresponding to the main frequency of the body wave signal, thus obtaining the bandpass-filtered simulated gun assembly. The top of the bandpass-filtered pseudo-shot collection is cut off to suppress the interference of high-energy noise in the shallow near-surface area, resulting in a pseudo-shot collection with suppressed shallow noise. Amplitude compensation processing is performed on the simulated shot set after shallow noise suppression to enhance the signal energy of the deep reflection waveform, resulting in a simulated shot set with equal amplitude. Frequency wavenumber domain filtering is applied to the simulated gun set. Based on a preset apparent velocity threshold, surface wave energy with low apparent velocity is suppressed, while volume wave reflection events with high apparent velocity are retained, in order to obtain simulated gun set data with improved signal-to-noise ratio.

5. The ice cap thickness detection method based on drilling noise imaging according to claim 4, characterized in that, Velocity analysis is performed on the simulated shot gather data with improved signal-to-noise ratio to obtain a migration velocity model for migration imaging, including: On the simulated shot set data with improved signal-to-noise ratio, select common reflection point gathers that contain the main reflection events; A velocity scan is performed on the selected gathers. By calculating the superposition energy at different test velocities, the corresponding velocity that makes the reflection phase axis focus optimal is identified, and a series of discrete velocity control points are obtained. Based on the velocity control points, spatial interpolation is performed along the survey line to obtain a preliminary continuous velocity distribution; The initial continuous velocity distribution is smoothed to eliminate local outliers and maintain the rationality of velocity changes, ultimately yielding a migration velocity model for migration imaging.

6. The ice cap thickness detection method based on drilling noise imaging according to claim 5, characterized in that, Using a pre-stack time migration method, the simulated shot gather data with improved signal-to-noise ratio is processed for migration imaging through a migration velocity model to obtain a time-domain stacked profile, including: Acquire the migration velocity model and the simulated shot set data after signal-to-noise ratio enhancement; based on the geological structural characteristics of the imaging area and the preset migration parameters, set the maximum migration dip angle and migration aperture for pre-stack time migration; The pre-stack time migration method is adopted, and the migration velocity model is used to perform migration calculation on the simulated shot gather data after the signal-to-noise ratio is improved. The data recorded by each receiver point is returned to the underground reflection point position according to the wave field propagation law, and the migrated common reflection point gather is obtained. The common reflection point gathers are stacked, and the energies of all gathers at corresponding positions are summed to obtain the final time-domain stacked profile.

7. The ice cap thickness detection method based on drilling noise imaging according to claim 6, characterized in that, Time-depth transformation is performed on the time-domain overlay profile to obtain the depth-domain profile. The ice sheet thickness is then determined based on the spatial location of the reflection phase axis of the ice-rock interface within the depth-domain profile, including: Obtain the time-domain overlay profile and the migration velocity model; based on the migration velocity model, convert the two-way travel time data in the time-domain overlay profile into the corresponding depth data to obtain a preliminary depth-domain profile. The preliminary depth domain profile is smoothed or interpolated to eliminate depth jumps caused by velocity model errors, thus obtaining a continuous and reliable final depth domain profile. In the final depth domain profile, identify and pick out continuous reflection in-phase axes representing the ice-rock interface; Based on the depth of the reflection phase axis at the ice-rock interface and its variation in the horizontal direction, the ice cover thickness at the corresponding location below the survey line is calculated and output.

8. An ice cap thickness detection system based on drilling noise imaging, wherein the system implements the method as described in any one of claims 1 to 7, characterized in that, include: The acquisition module is used to deploy a two-dimensional linear survey line consisting of multiple geophones around the ice cap drilling site to collect continuous noise data composed of mechanical noise generated by drilling operations and auxiliary equipment, as well as environmental background noise. The calculation module is used to preprocess the continuous noise data. The preprocessing includes removing the mean, removing the linear trend, removing strong amplitude events in the time domain, and removing strong amplitude components in the frequency domain from the records of each detector in sequence, so as to eliminate instrument drift and local strong noise interference and obtain preprocessed noise data. A detector in a two-dimensional linear survey line is selected as the location of a virtual seismic source. Mutual interference calculation is performed on the noise data to recover the pseudo-shot gather reflection response corresponding to the virtual seismic source. The extraction module is used to sequentially perform bandpass filtering, top cut-off, amplitude compensation and frequency wavenumber domain filtering on the simulated gun set reflection response to extract the body wave signal in the preset frequency band, suppress the surface wave energy and enhance the deep reflection signal to obtain simulated gun set data with improved signal-to-noise ratio. The analysis module is used to perform velocity analysis based on the simulated shot set data with improved signal-to-noise ratio to obtain a migration velocity model for migration imaging. The processing module is used to perform migration imaging processing on the simulated shot gather data after the signal-to-noise ratio has been improved by using the pre-stack time migration method and the migration velocity model to obtain a time-domain stacked profile; the time-domain stacked profile is converted to depth to obtain a depth-domain profile, and the ice cover thickness is determined according to the spatial position of the reflection phase axis of the ice-rock interface in the depth-domain profile.

9. A computing device, characterized in that, include: One or more processors; A storage device for storing one or more programs, which, when executed by one or more processors, cause the one or more processors to implement the method as described in any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores a program that, when executed by a processor, implements the method as described in any one of claims 1 to 7.

Images