Stratum water content layering inversion method based on shallow earthquake full waveform

Abstract

The invention discloses a stratigraphic water content layering inversion method based on a shallow seismic full waveform, and relates to the technical field of geophysical exploration data processing, and the method comprises the following steps: 1, under the same survey line condition, synchronously obtaining a seismic wave full waveform, complex impedance and an induced polarization frequency spectrum; 2, constructing a layered initial model containing a layer boundary line and an in-layer water content initial value; based on the layered initial model, executing dual-domain quantile alignment interface keeping iteration until convergence, and obtaining a final layer boundary line and an in-layer water content field; the double-domain quantile alignment interface maintaining iteration comprises the following steps: updating an in-layer moisture content field and a layer boundary line according to a seismic wave full waveform; and performing total variation regularization processing on the updated in-layer water content field to maintain the transition characteristic of the layer boundary line. According to the method, an efficient and reliable solution with actual guidance value is provided for mineral resource exploration, geological structure detection, groundwater survey and engineering geology evaluation.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration data processing technology, and particularly to a method for layered inversion of formation water content based on the full waveform of seismic waves, applicable to mineral resource exploration, engineering geological investigation, and groundwater detection. It is an intelligent sensing system for mineral geological exploration services utilizing advanced technology. Background Technology

[0002] In mineral resource exploration, oil and gas exploration, engineering geological surveys, and groundwater resource assessment, formation water cut is a key parameter affecting reservoir evaluation, mineralization prediction, engineering stability assessment, and groundwater resource evaluation. Changes in water cut directly lead to changes in the formation's dielectric constant, electrical conductivity, elastic wave velocity, and mechanical properties, thereby affecting the accuracy of exploration target identification, reservoir quality evaluation, and engineering safety assessment. Therefore, accurately obtaining the water cut distribution of different formation layers is of great significance for improving mineral resource exploration efficiency, optimizing drilling layout, and ensuring engineering safety.

[0003] Existing methods for detecting formation water cut mainly include borehole sampling, well logging interpretation, seismic wave inversion, time-domain reflectometry, frequency-domain reflectometry, and resistivity or induced polarization measurements. While borehole sampling offers high accuracy, it is costly, time-consuming, and only provides discrete point information, failing to reflect the continuous distribution characteristics of the formation. Well logging methods, although providing high-resolution data near the wellbore, have limited lateral detection range, unable to meet the needs of large-area exploration. Seismic wave exploration, as the most important geophysical exploration method, inverts formation physical parameters by analyzing the amplitude, arrival time, and phase information of reflected waves. It has advantages such as large detection depth, high lateral resolution, and relatively low cost, and is widely used in mineral exploration and structural exploration. Traditional seismic wave inversion methods are mainly based on wave impedance inversion or velocity inversion, converting wave velocity or wave impedance into water cut through empirical formulas. However, in multi-layered structures, reflected waves from different layer interfaces superimpose on each other, resulting in complex waveforms. Traditional processing methods based on envelope or time delay are difficult to distinguish between adjacent thin layers, especially under complex geological conditions where aquifers and dry layers alternate, resulting in large errors in layer interface identification. Summary of the Invention

[0004] The purpose of this invention is to provide a layered inversion method for formation water content based on shallow seismic full waveform. By organically combining dual-domain joint constraints of seismic full waveform and electrochemical spectrum, quantile alignment strategy, total variation regularization processing, and interface-preserving iterative mechanism, significant technological advancements have been achieved in areas such as clear identification of layer interfaces, stable inversion of water content within layers, non-uniqueness suppression, automated processing, and adaptability to complex geological conditions. This provides an efficient, reliable, and practically valuable solution for mineral resource exploration, geological structure detection, groundwater investigation, and engineering geological assessment, and has significant theoretical innovation significance and broad application and promotion value.

[0005] To address the aforementioned technical problems, this invention provides a method for layered inversion of formation water content based on shallow seismic full waveforms, comprising the following steps: Step 1: Under the same survey line conditions, simultaneously acquire the full waveform, complex impedance, and induced polarization spectrum of the seismic wave; Step 2: Construct a layered initial model containing the initial values ​​of the layer boundaries and the intralayer water content; and based on the layered initial model, perform a two-domain quantile alignment interface maintenance iteration until convergence to obtain the final layer boundaries and the intralayer water content field; the two-domain quantile alignment interface maintenance iteration includes: updating the intralayer water content field and layer boundaries according to the full waveform of the seismic wave; updating the intralayer water content field according to the complex impedance and induced polarization spectrum through quantile alignment of spectral fingerprint and envelope fingerprint; and performing total variation regularization on the updated intralayer water content field to maintain the transition characteristics of the layer boundaries.

[0006] Furthermore, the steps for constructing a layered initial model that includes layer boundaries and initial values ​​of water content within the layers include: obtaining candidate reflection zones by peak-valley tracking based on the envelope peak, phase plateau, and instantaneous frequency steps of the full seismic wave waveform; obtaining strong response zones by frequency band clustering based on the arc shape and amplitude-phase peak shape of the complex impedance and induced polarization spectrum; performing voxel-level registration between the candidate reflection zones and the strong response zones; extracting connected ridges through connected domain analysis; and generating a set of continuous layer boundaries as layer boundaries using shortest path search.

[0007] Furthermore, before constructing the hierarchical initial model, a data standardization step is also included, including: performing DC removal, bandpass filtering, first arrival template subtraction, inter-channel amplitude equalization, and time zero-point correction sequentially on the full waveform of the seismic wave to obtain the standardized full waveform; and performing electrode geometry correction, amplitude and phase consistency, and frequency resampling on the complex impedance and induced polarization spectrum to obtain the standardized spectrum.

[0008] Furthermore, the process of updating the intralayer water content field and layer boundary based on the full waveform of seismic waves specifically includes: generating a one-dimensional reflection sequence under the current layer boundary and intralayer water content field, and using zero-phase source wavelets to convolve and synthesize the full waveform; cross-correlating and aligning the synthesized full waveform with the standardized full waveform to obtain the envelope difference map, phase difference map, and instantaneous frequency difference map; constructing a pyramid sequence from the low-frequency sub-band to the high-frequency sub-band, and matching and tracking according to the sub-band order to locate the deviation area in the difference map, and locally updating the intralayer water content field in the corresponding intralayer region.

[0009] Furthermore, after the process of locally updating the water content field within the corresponding layer region, the process also includes: calculating the quantile distribution of the envelope time shift within the deviation region, pairing it with the quantile distribution of the phase difference to generate a time shift-phase pairing table, adjusting the water content field within the layer based on the time shift-phase pairing table, and straightening the time-to-depth conversion relationship based on the phase anchoring results to update the layer boundary.

[0010] Furthermore, based on the complex impedance and induced polarization spectrum, the process of updating the intralayer water content field through quantile alignment of spectral fingerprint and envelope fingerprint specifically includes: dividing the normalized spectrum into several frequency bands according to the logarithmic frequency; performing circular arc fitting on the Nyquist trajectory within each frequency band and removing the relaxation component; within each intralayer region, statistically analyzing the arc radius obtained from the circular arc fitting and the area enclosed by the trajectory as the spectral fingerprint, and statistically analyzing the empirical quantile curve of the spectral fingerprint; quantile alignment of the empirical quantile curve of the spectral fingerprint with the quantile curve of the envelope intensity to form a monotonic mapping table of spectral fingerprint-envelope fingerprint; and based on the monotonic mapping table of spectral fingerprint-envelope fingerprint, converting the spectral fingerprint of each voxel into a water content increment and updating the intralayer water content field.

[0011] Furthermore, the process of performing total variation regularization on the updated intralayer water cut field specifically includes: calculating the first-order difference map in the mileage direction and the depth direction within each intralayer region; using the upper quartile of each first-order difference map as the shrinkage threshold, setting the difference below the threshold to zero, and reducing the amplitude of the difference above the threshold to obtain the shrunken difference map; performing integral reconstruction on the shrunken difference map, and using the water cut value at the boundary of the intralayer region as the reconstruction anchor point to restore the intralayer water cut field after total variation regularization.

[0012] Furthermore, after restoring the intralayer water content field after total variation regularization, the process also includes: extracting the intersection points of the envelope extreme value zone and the water content gradient extreme value zone on the reconstructed water content field to form a set of interface candidate points; performing dynamic programming on the set of interface candidate points and constraining the mileage difference and curvature between adjacent points to obtain a continuous and smooth optimal layer boundary, and replacing the original layer boundary.

[0013] Furthermore, the convergence criterion for iteration until convergence is: the maximum displacement of the layer boundary between two cycles is less than the displacement threshold determined by the grid resolution, and the maximum change amplitude in the envelope difference map and the phase difference map is less than the amplitude threshold determined by the difference map distribution.

[0014] Furthermore, after iterative convergence, the process also includes a result stitching step, which includes: stitching the layer boundaries of adjacent data blocks with the intralayer water content field in mileage order; aligning the overlapping areas by endpoints and performing a linear transition; and performing median smoothing of a sliding window on the intralayer water content field to form a continuous profile.

[0015] The formation water content inversion method based on shallow seismic full waveform of the present invention has the following beneficial effects: This invention establishes a statistical mapping relationship between the full waveform envelope characteristics of seismic waves and the complex impedance-induced polarization spectral fingerprint through a dual-domain quantile alignment strategy, achieving deep fusion of these two types of physical field data in water content inversion. The full waveform of seismic waves provides high spatial resolution layer interface information and vertical structural features, while the complex impedance spectrum provides electrochemical response information sensitive to water content changes. The alignment of the two in the quantile space effectively avoids the empirical and subjective nature of weight parameter selection in traditional weighted combination methods. This results in inversion results that inherit the structural clarity of seismic wave data while incorporating the sensitivity of electrical resistivity data to fluid saturation, significantly improving the accuracy of water content inversion in complex geological environments such as mineralized water layers, clay interlayers, and metallic mineralization zones.

[0016] The interface-preserving iterative mechanism proposed in this invention maintains the transition characteristics of layer boundaries through total variation regularization, while allowing for the smooth evolution of the water content field within the layers. In each iteration, multi-dimensional matching tracking is first performed based on the envelope difference, phase difference, and instantaneous frequency difference of the entire seismic wave waveform to accurately locate the layer interface; then, the circular fitting characteristics of the complex impedance spectrum are used to fine-tune the water content within the layers; finally, through the contraction-reconstruction process of total variation regularization, intra-layer noise is suppressed while preserving the inter-layer transitions. This mechanism ensures that the accuracy of layer boundary positioning and the smoothness of intra-layer water content no longer compromise, and the inversion results show clear transitions at the layer boundaries while maintaining physically reasonable continuity within the layers. It is particularly suitable for fine-grained exploration of complex geological structures such as thin interbedded layers and lenticular bodies.

[0017] The high-precision layered water content profiles output by the method of this invention can be directly used for 3D geological modeling, reservoir parameter field establishment, and groundwater flow numerical simulation, providing quantitative basis for mineral resource target area selection, borehole location design, and engineering risk assessment. The inversion results are output in a standard format, seamlessly integrating with GIS systems, supporting automatic geological profile mapping, water content contour plotting, and 3D visualization, significantly improving the informatization level and decision support capabilities of geological exploration results. In projects such as regional geological surveys, mineral prospect evaluation, and groundwater resource exploration, the method of this invention has become an important technical means, promoting the development of geophysical exploration technology towards refinement and intelligence. Attached Figure Description

[0018] Figure 1 A schematic diagram of the mileage-depth profile of the layer boundary and water content field provided for an embodiment of the present invention; Figure 2 This is a schematic diagram of the complex impedance Nyquist trajectory and circular arc fitting features provided in an embodiment of the present invention; Figure 3 This is a schematic diagram of the quantile alignment mapping between spectral fingerprints and envelope fingerprints provided in an embodiment of the present invention; Figure 4 This is a schematic diagram illustrating how total variation regularization preserves the characteristics of layer boundary transitions, as provided in an embodiment of the present invention. Figure 5 This is a schematic diagram of the quantile pairing update mechanism for envelope peak time shift and phase difference provided in an embodiment of the present invention. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0020] Example 1: A method for layered inversion of formation water content based on shallow seismic full waveform, comprising the following steps: Step 1: Under the same survey line conditions, simultaneously acquire the full waveform, complex impedance, and induced polarization spectrum of the seismic wave; Step 2: Construct a layered initial model containing the initial values ​​of the layer boundaries and the intralayer water content; and based on the layered initial model, perform a two-domain quantile alignment interface maintenance iteration until convergence to obtain the final layer boundaries and the intralayer water content field; the two-domain quantile alignment interface maintenance iteration includes: updating the intralayer water content field and layer boundaries according to the full waveform of the seismic wave; updating the intralayer water content field according to the complex impedance and induced polarization spectrum through quantile alignment of spectral fingerprint and envelope fingerprint; and performing total variation regularization on the updated intralayer water content field to maintain the transition characteristics of the layer boundaries.

[0021] refer to Figure 1 This figure illustrates the distribution of layer boundaries and water content fields on the mileage-depth profile after iterative convergence of the bi-domain quantile alignment interface. The figure uses a mileage-depth grid as the unified display medium. The horizontal axis represents the mileage direction in meters (m), ranging from 0 to 50 meters with a step size of 0.02 meters. The vertical axis represents the depth direction in meters (m), ranging from 0 to 1.0 meters with a step size of 0.01 meters. The origin is located in the upper left corner of the figure, with mileage to the right being the positive direction and depth downwards being the positive direction. Within this mileage-depth grid, the water content field is represented by grayscale filling. The grayscale value is negatively correlated with the volumetric water content; that is, a lighter color indicates a lower water content, and a darker color indicates a higher water content. Figure 4 A moisture content legend is displayed on the right, using vertical color bars with gray levels ranging from light to dark, corresponding to a continuous change in moisture content from 5% to 60%. Several scale values ​​are labeled next to the legend for easy quantitative reading of moisture content at any location. Three continuous black solid lines are superimposed on the moisture content field background, representing the first, second, and third interfaces. The first interface is located at a depth of approximately 0.3 meters, exhibiting slight undulations along the mileage direction with an amplitude of approximately ±0.02 meters and a period of 5 to 10 meters. The second interface is located at a depth of approximately 0.55 meters, also showing periodic undulations along the mileage direction with an amplitude of approximately ±0.025 meters.

[0022] The third interface is located at a depth of approximately 0.8 meters, exhibiting undulations similar to the first two interfaces but with a longer period. All three layer boundaries remain continuous and smooth, with curvature variations between adjacent points controlled to no more than 5 degrees per 0.02 meters, meeting the requirements for maintaining transition characteristics after total variation regularization. Four intralayer regions are divided by the three layer boundaries. The first region lies between the surface and the first interface, with a depth ranging from 0 to approximately 0.3 meters. This region has a generally low water content and a shallow grayscale value, with volumetric water content mainly distributed between 5% and 15%, exhibiting slight random fluctuations along the mileage direction. The second region lies between the first and second interfaces, with a depth ranging from approximately 0.3 to 0.55 meters. This region has a significantly increased water content and a markedly deeper grayscale value, with volumetric water content mainly distributed between 20% and 30%. The third layer, located between the second and third interfaces, ranges in depth from approximately 0.55 to 0.8 meters, exhibiting a further increase in water content, with volumetric water content primarily between 30% and 40%. The fourth layer, situated below the third interface, ranges in depth from 0.8 to 1.0 meter. This region has the highest water content and the deepest grayscale, with a volumetric water content reaching 40% to 50%. At the layer interfaces, the water content field exhibits a clear transition characteristic. When crossing layer boundaries, the grayscale value abruptly changes in the depth direction, indicating a significant difference in water content between adjacent layers. This transition characteristic is preserved through total variation regularization, ensuring a gradient transition of one depth unit on each side of the layer boundary, while maintaining a relatively uniform water content distribution within the layer, effectively suppressing minor fluctuations within the layer. Within the same layer, the water content field shows a continuous gradual change in the kilometer direction, without any isolated gradient peaks exceeding 0.10 meters, reflecting the kilometer continuity after median smoothing. The intersection of the layer boundary and the envelope extreme value zone forms a continuous curve in the mileage direction, indicating that the location of the layer boundary is based on the overlapping area of ​​the local maximum of the envelope peak of the full seismic wave waveform and the water content gradient extreme value zone, and the optimal path for the entire segment is obtained through dynamic programming.

[0023] In the specific implementation process, regarding the establishment of the same survey line conditions, a baseline is first laid out on the road surface using wear-resistant marking paint along the planned route. The baseline is parallel to the road centerline, and the lateral distance between the baseline and the road boundary is fixed at 1.50 meters. Start and end markers are placed at both ends of the baseline, and a positioning marker is set every 50 meters. A mileage scale sequence is generated along the baseline direction in 0.02-meter increments, using a ground-based laser displacement sensor as the mileage trigger source, with a trigger error not exceeding 0.002 meters. The allowable error for lateral offset of the survey line is controlled within 0.05 meters. Longitudinal and cross slope information is recorded in real time by inertial measurement equipment with a resolution better than 0.01 degrees. This ensures repeatability under the same survey line conditions, making the surface coupling state and geometric conditions consistent at the same mileage point, thereby reducing the difference in reflection path and electric field distribution caused by path offset.

[0024] The system employs a ground-level pushing or vehicle-mounted ground-level method to maintain the transmitting and receiving antennas at a height between 0.12 meters and 0.18 meters above the ground, with a variation not exceeding 0.01 meters throughout the entire measurement range. The center frequency is selected based on the actual layer thickness resolution: a 900 MHz channel is used when the roadbed thickness is between 0.15 meters and 0.40 meters; a 400 MHz channel is added when the roadbed thickness is between 0.30 meters and 0.80 meters, and dual-channel simultaneous acquisition is performed. The sampling rate is set to 2.56 gigabits per second, the recording window is 1000 nanoseconds, and the quantization precision is 16 bits. Distance triggering is used along the mileage direction, triggering a sampling channel every 0.02 meters, with trigger jitter less than 1 microsecond. Superimposed sampling is performed at each mileage point, with a default superposition count of 8 times and a superposition interval of 2 milliseconds, using energy superposition to increase the square root of the signal-to-noise ratio. To ensure phase stability across the entire waveform, an aluminum alloy reflector, 10 mm thick and 1 m x 1 m in size, was placed at the starting point and every 200 meters. These were used to record a stable reference for the first arrival wave and a strong near-surface reflection reference. During acquisition, the earliest arriving and largest amplitude paired positive and negative peaks were used as the first arrival wave anchor points. This is because these peaks are less sensitive to changes in the near-surface medium and have stable arrival times under the same antenna geometry. Anchoring these peaks aligns all channels to the same time reference, reducing subsequent time-to-depth conversion errors. To reduce inter-channel geometric stretching caused by variations in travel speed, the travel speed was controlled between 40 km / h and 50 km / h, with speed fluctuations not exceeding ±2 km / h. At strong scattering locations such as bridgeheads, culverts, and metal covers, the speed was reduced to 20 km / h, and the number of superpositions was automatically increased to 16 to maintain consistent inter-channel signal-to-noise ratios.

[0025] A multi-electrode wheel-type contact array is used, rolling along the same measurement line. The array contains 32 contact electrodes with a spacing of 0.50 meters, and a total unfolded length of 15.5 meters. Current injection and potential measurement are performed in pairs. The current injection range is 1 mA to 50 mA, automatically limited according to the contact resistance, with a maximum contact resistance threshold of 1000 ohms. The frequency scanning range is 0.1 Hz to 10000 Hz, with 10 frequency points set every ten octaves of the logarithmic frequency, for a total of 41 frequency points. At each frequency point, stable injection is maintained for at least 5 cycles before amplitude and phase sampling is performed. The sampling duration is not less than 3 times the measurement period to allow the electric field distribution to reach a steady state at that frequency point, thereby forming the amplitude and phase curves and Nyquist trajectory of the complex impedance and induced polarization spectrum. Considering the need for longer settling time in the low-frequency band, a 20-second dwell time is set for each point from 0.1 Hz to 1 Hz; a 5-second dwell time is set for each point from 1 Hz to 10 Hz; and a 1-second dwell time is set for each point from 10 Hz to 10000 Hz. To avoid current distribution disturbances caused by travel, a brief pause is made at every 5-meter distance, with a pause duration of 3 to 6 seconds. Low-frequency measurements are completed within the pause window; mid-to-high frequency measurements are completed within the travel window between two pause points. This segmented measurement arrangement ensures sufficient steady-state conditions at low-frequency points while maintaining high mileage resolution at high-frequency points, both of which together guarantee the geometric integrity of the Nyquist trajectory. Each electrode wheel is equipped with a pressure sensor, and the single-electrode contact pressure is maintained between 100 N and 150 N. Pressure closed-loop control reduces and stabilizes the contact resistance, thereby reducing jumps in the amplitude-phase curve and dispersion of the Nyquist trajectory.

[0026] A single master time base is used for unified time stamping. The full-wavelength seismic wave acquisition equipment and the complex impedance and induced polarization spectrum acquisition equipment are connected to the same time source, a combination of satellite timing pulses and a 10 MHz steady-state reference, with a timestamp resolution of 1 microsecond. All mileage trigger events are generated by ground-based laser displacement sensors and simultaneously delivered to both acquisition devices via a hardware distributor. The event arrival time difference is maintained within 5 microseconds through laboratory testing. Each mileage trigger event generates a data block identifier, containing a full-wavelength seismic wave channel and a set of complex impedance and induced polarization spectrum data points completed within that mileage window, along with recorded values ​​for attitude angle, temperature, and contact pressure. This alignment aims to ensure that the envelope peaks, phase, and instantaneous frequency characteristics of the full-wavelength seismic wave, along with the amplitude-phase curves and Nyquist trajectories at the same location, are entered into the same mileage-depth grid cell, reducing interpolation errors during cross-domain feature registration.

[0027] The longitudinal spacing between the GPR antenna and the electrode array is set to 2.0 meters, and the lateral offset is maintained at 0.20 meters. A continuous shielding curtain, made of composite fabric containing metal wires, is installed between them, with a height of 0.30 meters, dragging close to the ground to reduce mutual electromagnetic coupling. Before the start of the test, current injection is turned off, and 100 channels of the full seismic waveform are independently acquired to calculate the root mean square value of the background noise. Then, current injection is turned on, and 100 channels are acquired again. If the root mean square value increases by more than 20%, the spacing between the antenna and the electrode is increased to 2.5 meters, and the test is repeated until the increase does not exceed 20%. To address external power frequency interference, the 50 Hz and 100 Hz components are suppressed using a narrowband notch filter, and the harmonics of these frequencies are avoided in the frequency planning of the complex impedance and induced polarization spectrum. For areas with metal manhole covers, expansion joints, and dense steel reinforcement, an automatic avoidance strategy is set up: when the early strong amplitude pulses of the full waveform of the seismic wave exceed 10 consecutive channels, the low-frequency measurement of complex impedance and induced polarization spectrum is suspended in the same mileage interval, and only the high-frequency measurement is retained. The low-frequency measurement is then resumed after crossing the interval, in order to avoid abnormal stretching of the Nyquist trajectory caused by local short circuits in the metal.

[0028] Before the initial measurement, two dielectric plates with thicknesses of 0.05 meters and 0.10 meters, respectively, were placed sequentially on the road surface. The dielectric properties of these materials were stable. The arrival time difference from the first arrival wave anchor point to the first main reflection peak was measured by recording the full waveforms of the seismic waves on and below the plates. Since the thicknesses of the two dielectric plates were known, two-point constraints were established using these two sets of arrival time differences to create a time-to-depth conversion table. Combined with the stability of the first arrival wave anchor points along the route, a time-to-depth conversion table was generated for each data block in subsequent processing. The advantage of this approach is that the time reference and near-surface propagation characteristics are anchored by actual measurements, thereby reducing depth-scale drift caused by changes in near-surface moisture, temperature, and roughness.

[0029] The signal-to-noise ratio (SNR) threshold for the full seismic wave waveform is set to 12 dB. If more than 20 consecutive channels fall below the threshold, the stacking order is automatically increased to 16 and the travel speed is reduced to 20 km / h. The Nyquist trajectory of the complex impedance and induced polarization spectrum should be approximately a piecewise circular arc. If the coefficient of variation of the Euclidean distance between points in the trajectory exceeds 0.30 within the same frequency band, the frequency band is considered not to have reached steady state, and the dwell time at that frequency point is automatically extended until the coefficient of variation drops below 0.20. All real-time correction actions are recorded in the metadata of the data block for subsequent identification of anomalous intervals and arrangement of supplementary measurements. Each 0.02-meter mileage step constitutes a data block, which stores one original full seismic wave waveform channel, the corresponding amplitude and phase curves and Nyquist trajectory of the complex impedance and induced polarization spectrum, attitude angle and temperature records, contact pressure records, and time and mileage markers. A segment file is generated every 50 meters. The segment file contains 2500 consecutive data blocks, which facilitates the construction of the initial model by segment and the iterative maintenance of the dual-domain stratification alignment interface in subsequent steps. To adapt to the subsequent mileage-depth grid, the paragraph file includes longitudinal and transverse slope curves along the route, ensuring that registration and interpolation are performed on the same grid in the future.

[0030] In one alternative implementation, a pure push-type method can be used instead of a vehicle-mounted method, with a travel speed between 1 km / h and 3 km / h. In this case, the low-frequency band of the complex impedance and induced polarization spectrum can be completed while in motion without a short stop, and the dwell time can be set according to the same rules. A fixed station method can be used, stopping at points every 5 meters. The full waveform of the seismic wave is superimposed 32 times during the stop, and the complex impedance and induced polarization spectrum are scanned across the entire frequency band. This method is suitable for high-precision measurement sections under closed traffic conditions. The number of electrodes can be expanded to 48, and the electrode spacing can be reduced to 0.30 meters for areas sensitive to changes in thin-layer moisture. In cold seasons, low-salinity water can be sprayed on the electrode contact area to reduce contact resistance. The spraying amount is controlled to within 10 ml per electrode per spray to avoid significant disturbance to the moisture field of the mileage-depth grid.

[0031] In the initial model construction phase, a mileage-depth grid was used as the unified carrier. The mileage step size was set to 0.02 meters, and the depth step size was set to 0.01 meters. Envelope maps, phase maps, and instantaneous frequency maps were calculated for the entire seismic wave waveform. The envelope map was obtained using analytical signal methods. Subsequently, a local median and local absolute deviation were determined along the mileage using a sliding window of 0.20 meters in length for each track. The peak higher than the local median plus three times the local absolute deviation was selected as the envelope peak. The advantage of this selection is that it can highlight the concentrated area of ​​reflected energy against the background of slow changes in road material, and the position of the envelope peak is relatively stable under the same antenna geometry. The phase map was constructed using phase expansion and second-order difference smoothing along the mileage. Continuous bands with a phase change rate of less than 5 degrees per 0.02 meters were marked as phase-stable bands. Instantaneous frequency maps calculate local frequencies within each channel with a depth step of 0.01 meters. Steps are then detected along the mileage direction with a 0.20-meter window. A step is identified when the frequency difference between two adjacent depth layers continuously exceeds 10 Hz and its duration exceeds 0.10 meters within the same window. The intersection and union of the envelope peak, phase plateau, and instantaneous frequency steps on the mileage-depth grid are used to obtain candidate reflection bands. Voxel elements within a candidate reflection band with a continuous length of at least 0.50 meters are considered valid; isolated segments with insufficient length are discarded.

[0032] On the complex impedance and induced polarization spectrum side, strong response bands are extracted using Nyquist trajectories and amplitude-phase curves. Frequency points are logarithmically divided into five bands: 0.1 Hz to 1 Hz, 1 Hz to 10 Hz, 10 Hz to 100 Hz, 100 Hz to 1000 Hz, and 1000 Hz to 10000 Hz. For each band, a circular arc is fitted in the Nyquist plane, with at least 10 fitting points and a residual standard deviation not exceeding 3% of the fitted arc length. The arc radius and the area enclosed by the trajectory are recorded for each arc, and the peak position and peak width of the same frequency band are recorded in the amplitude-phase curve. When a large arc radius and a large trajectory enclosed area appear at the same mileage-depth location in the low-frequency band, and a distinct peak position and a wide peak width appear in the amplitude-phase curve, this location is marked as a strong response band. This approach combines the low-frequency relaxation characteristics, which are more sensitive to pore liquid phase behavior, with the amplitude-phase peak shape, prioritizing the marking of regions where water content significantly affects the electrochemical response. To suppress occasional spikes, isolated outliers are removed by three-point median filtering, and the strong response band is required to have a continuous length of not less than 0.50 meters in the mileage direction.

[0033] During the voxel-level registration stage, candidate reflection bands and strong response bands are mapped to the same mileage-depth grid cells. If a cell is covered by both, it is marked as a high-confidence cell; if it is covered by only one, it is marked as a general cell. Connectivity analysis is performed on high-confidence cells, and regions with a connected area of ​​less than 50 cells are deleted. The remaining high-confidence cells are skeletonized to obtain connected ridges. Shortest path search is performed on each connected ridge, extending progressively from the left end mileage to the right end mileage, with a step size of 1 cell. The allowed vertical step size is 1 cell per step, and the allowed curvature is limited to no more than 5 degrees per 0.02 meters by the angle change between adjacent three points. This setting can suppress unreasonable high-frequency jitter while preserving the true curvature characteristics of the layer boundary. Multiple shortest paths constitute a continuous set of layer boundaries, serving as the initial values ​​for the layer boundaries. Intra-layer regions are formed by being surrounded by adjacent layer boundaries. To provide an initial value for the water content within the layer, the envelope intensity is binned according to quantiles within each layer region. The quantile boundaries are set at five positions: 0.10, 0.30, 0.50, 0.70, and 0.90. These five quantile intervals are then mapped to volume fractions of 0.05, 0.15, 0.25, 0.35, and 0.45, respectively. If a position also belongs to a strong response band, the corresponding volume fraction is increased by 0.05, with an upper limit of 0.60.

[0034] The process employs a two-domain quantile alignment interface for iterative updates, using a fixed sequence: first, updating the intralayer water brine field and layer boundaries based on the full seismic waveform; then, updating the intralayer water brine field based on complex impedance and induced polarization spectra using quantile alignment of spectral and envelope fingerprints; finally, applying total variation regularization to the updated intralayer water brine field to preserve the transition characteristics of the layer boundaries. A convergence check is performed after each iteration. If the check condition is not met, the process proceeds to the next iteration.

[0035] Based on the full waveform of the seismic wave, the water content field and layer boundaries within the layer are updated. First, a one-dimensional reflection sequence is generated under the current layer boundary and water content field, and aligned with the standardized full waveform to form an envelope difference map, phase difference map, and instantaneous frequency difference map. Then, a pyramid sequence is constructed from low-frequency sub-bands to high-frequency sub-bands, with four sub-bands. Each sub-band covers a dominant frequency range of 100 MHz to 200 MHz, 200 MHz to 400 MHz, 400 MHz to 800 MHz, and 800 MHz to 1200 MHz, respectively. Deviation zones are identified within each sub-band. A deviation zone is defined as an area within the same mileage window where the envelope difference map continuously exceeds the local median plus twice the local absolute deviation, with a duration of not less than 0.20 meters. Simultaneously, the median absolute value of the phase difference map within the same area must be greater than 5 degrees. For each deviation zone, an update block is delineated with a length of 0.20 meters in the mileage direction and a thickness of 0.05 meters in the depth direction. Calculate quantile curves for the envelope peak time shift and phase difference within the update block, with quantile positions set at 0.20, 0.50, and 0.80. Pair the two quantile curves to form a paired table of time shift and phase. If, at the same quantile position, the envelope peak time shift is positive and the phase difference is in the same direction, it indicates that the propagation delay and phase lag within the update block point in the same direction. In this case, uniformly increase the intralayer water content field of the update block by 0.01 across all cells. If the envelope peak time shift is negative and the phase difference is in the same direction, uniformly decrease by 0.01. If the two are opposite, adjust by 0.01 only on the top 30 percentile cells with the highest envelope peak intensity, leaving the remaining cells unchanged. After adjustment, recalculate the envelope difference map and phase difference map. If the absolute value of the median of both decreases within the update block, accept the adjustment; otherwise, roll back. A maximum of 12 consecutive adjustments can be performed on the same update block. After completing the subband processing, accumulate the subband results into the full-frequency results. Under full-frequency results, the same scan line is straightened along its entire length according to the time-to-depth conversion table. The waveform morphology along the mileage direction is aligned to the same first-arrival anchor point. The layer boundary is then repositioned at the location where it crosses the envelope extremum zone, prioritizing the passage of the layer boundary through the local maximum of the envelope peak. The reason for locating the layer boundary at the local maximum of the envelope peak is that the reflection intensity is most concentrated at the interface, resulting in a more stable energy distribution of the cross-boundary unit and improving the repeatability accuracy of the layer boundary.

[0036] In updating the intralayer water content field through quantile alignment of spectral fingerprints and envelope fingerprints, the Nyquist trajectory is first fitted with arcs segment by segment according to five fixed frequency bands. The arc fitting for each segment is solved using least-squares geometric distance, and the standard deviation of the fitting residual does not exceed three percent of the diagonal of the circumscribed rectangle of that segment's point cloud. The obtained arc radius and trajectory enclosed area are used as the basic quantities of the spectral fingerprint. Simultaneously, the peak position and peak width are calculated on the amplitude-phase curve as auxiliary entries. To establish a correspondence with the full waveform of the seismic wave, quantile curves are calculated for the envelope intensity within each layer region, with quantile positions of 0.10, 0.30, 0.50, 0.70, and 0.90, serving as the envelope fingerprint. The quantile curves of the spectral fingerprint and the envelope fingerprint are then aligned one-to-one to form a monotonic mapping table between the spectral and envelope fingerprints. The quantile alignment method is as follows: at the same quantile position, if the arc radius is larger and the trajectory enclosed area is larger, the corresponding envelope intensity is higher; if the arc radius is smaller and the trajectory enclosed area is smaller, the corresponding envelope intensity is lower. This correspondence reflects the co-enhancing trend of electrochemical relaxation response and electromagnetic reflection intensity as water content increases. Using a monotonic mapping table, the spectral fingerprint entry for each voxel is converted into a water content increment. Specifically, the increment is 0.02 when the voxel's spectral fingerprint is in the high quantile, 0.02 when it is in the low quantile, and no adjustment when it is in the middle quantile. If a voxel has a frequency band with abnormal fitting residuals in its neighborhood along the mileage direction, this frequency band information is temporarily not used; only the spectral fingerprint of a reliable frequency band is used for conversion to avoid interference from abnormal frequencies on the update direction. After conversion, the gradient magnitudes of water content along the mileage and depth directions are calculated within the layer region. If the local gradient magnitude continuously exceeds a preset upper limit of 0.10 meters in thickness, the increment is limited to 0.01 at that point to prevent excessive spatial diffusion of strong local responses on the spectral side.

[0037] To preserve the transition characteristics at layer boundaries, total variation regularization was applied to the updated intralayer water content field, performed separately for each intralayer region. First, first-order difference maps in the mileage and depth directions were calculated. The upper quartiles of each difference map were calculated, setting differences below the upper quartile to zero and proportionally reducing differences above the upper quartile to the midpoint between the upper quartile and the maximum value. This process suppresses minor fluctuations within the layer while preserving significant transitions at layer boundaries. Subsequently, integral reconstructions were performed in both the mileage and depth directions. The integration starting point was selected as the water content value at the boundary of the intralayer region as the anchor point to ensure continuity between the reconstructed values ​​and the surrounding area. After reconstruction, candidate point sets were extracted at the intersection of the reconstructed results and the envelope extremum zone of the full seismic waveform. The mileage interval between adjacent candidate points was required to be one unit, and the depth interval was required to be no more than one unit. Dynamic programming is performed on the candidate point set, progressively selecting points with the maximum envelope strength and the most significant water content gradient transition from the left end to the right end, while constraining the curvature of adjacent points to no more than 5 degrees per 0.02 meters, to obtain a continuous and smooth optimal layer boundary. This optimal layer boundary replaces the original layer boundary, and the next iteration begins.

[0038] Convergence is determined by simultaneously meeting two conditions as the stopping criteria. The first condition is that the maximum displacement of the layer boundary between two cycles is less than a displacement threshold determined by the depth step size, which is 0.02 meters. The second condition is that the maximum variation amplitude in both the envelope difference map and the phase difference map is less than an amplitude threshold determined by the difference map distribution, which is the tenths of the corresponding values ​​in the envelope difference map and the phase difference map, typically falling between 0.05 and 0.10 of the normalized amplitude. If either condition is not met, the cycle continues to the next iteration. Convergence is usually achieved between 3 and 7 iterations. If convergence is not achieved after 10 consecutive iterations, the mileage interval most prone to misregistration is locally reinitialized, and the above process is repeated only within that interval.

[0039] In one alternative implementation, the mileage step size can be set to 0.01 meters to improve the resolution along the line, while the depth step size can be set to 0.005 meters. When generating the initial layer boundary, a candidate path pre-screening based on the instantaneous frequency step connectivity can be added, with the pre-screening threshold being a continuous length of not less than 0.30 meters. When converting the spectral fingerprint into a water content increment, a three-level increment strategy can be used, namely 0.01, 0.02, and 0.03, corresponding to the quantile intervals of 0.10 to 0.30, 0.30 to 0.70, and 0.70 to 0.90, to refine the update intensity under different response intensities. After the total variation regularization processing, an additional median smoothing with a mileage direction length of 0.20 meters can be performed to further improve intra-layer consistency without weakening the layer boundary transition.

[0040] refer to Figure 4This figure compares the effect of total variation regularization on preserving the boundary transition characteristics. The figure is divided into two parts: the upper part shows the water cut profile before regularization, and the lower part shows the water cut profile after regularization. In the upper figure, the horizontal axis represents mileage in meters (m), ranging from 0.0 to 2.0; the vertical axis represents water cut, ranging from 0.00 to 0.40. The water cut curve exhibits a clear stratified structure, transitioning from 0.20 to 0.15 near mileage 0.5m, from 0.15 to 0.08 near mileage 1.0m, from 0.08 to 0.10 near mileage 1.5m, and from 0.10 back to 0.20 near mileage 2.0m. However, the curve contains many minor fluctuations and noise within each layer; two typical intra-layer noise regions are marked with black dashed ellipses. Near mileage 0.3m, the water cut oscillates repeatedly between 0.19 and 0.20; near mileage 1.2m, the water cut exhibits minor fluctuations between 0.07 and 0.09. This intra-layer noise interferes with the accurate positioning of layer boundaries, reducing the reliability of the inversion results. In the lower figure, the coordinate axis settings are the same as in the upper figure. After total variation regularization, the water cut curves become smoother within each layer, and intra-layer noise is effectively suppressed. In the mileage 0.0 to 0.5m range, the water cut remains stable at 0.10; in the mileage 0.5 to 1.0m range, the water cut remains stable at 0.15; in the mileage 1.0 to 1.5m range, the water cut remains stable at 0.08; and in the mileage 1.5 to 2.0m range, the water cut initially remains at 0.10 before rising back to 0.10. More importantly, the water cut transition characteristics at each layer boundary are well preserved. Three key locations are marked with black dashed lines: at mileage 0.5m, the water cut jumps significantly and clearly from 0.10 to 0.15; at mileage 1.0m, the water cut jumps clearly from 0.15 to 0.08; and within the mileage 1.5m range, the intralayer smoothness is good. The right side of the figure highlights the three key functions of total variation regularization: suppressing minor intralayer fluctuations, preserving significant interlayer boundary transitions, and reducing gradients to the upper quartile midpoint. This method calculates the first-order differences along the mileage and depth directions, setting differences below the upper quartile to zero and proportionally reducing differences above the upper quartile, thereby suppressing minor intralayer variations while preserving significant gradient transitions at interlayer boundaries.

[0041] Example 2: When constructing a layered initial model including layer boundaries and initial values ​​of water content within layers, envelope maps, phase maps, and instantaneous frequency maps of the seismic wave full waveform are first generated on a kilometer-depth grid. The envelope map is calculated using analytical signal methods, and then the local median and local absolute deviation are calculated within a 0.20-meter window along the kilometer direction. Peaks exceeding the local median plus three times the local absolute deviation are designated as envelope peaks. The phase map, after phase expansion, calculates the phase change rate along the kilometer direction within a 0.20-meter window. Continuous zones with a rate lower than 5 degrees per 0.02 meters are designated as phase plateaus. The instantaneous frequency map is calculated with a depth step of 0.01 meters. Within the same 0.20-meter window along the kilometer direction, when the local frequency difference between adjacent depth layers continuously exceeds 10 Hz and its length exceeds 0.10 meters, it is marked as an instantaneous frequency step. The intersection and union of these three types of markers on the grid are used to form candidate reflection zones. Subsequently, on the complex impedance and induced polarization spectrum side, the frequency was logarithmically divided into five bands: 0.1 Hz to 1 Hz, 1 Hz to 10 Hz, 10 Hz to 100 Hz, 100 Hz to 1000 Hz, and 1000 Hz to 10000 Hz. Circular fitting was performed on the Nyquist plane, with the residual standard deviation controlled within 3% of the fitted arc length. Simultaneously, peak positions and peak widths for the same frequency band were extracted from the amplitude-phase curves. When a large arc radius and a large trajectory enclosed area appeared in the low-frequency band, and the corresponding amplitude-phase peak position was obvious with a wide peak width, this location was marked as a strong response band. Candidate reflection bands and strong response bands were mapped to the same voxel element; elements covered by both were designated as high-confidence elements, while elements covered by only one were designated as general elements. Connectivity analysis was performed on the high-confidence elements, regions with an area less than 50 elements were removed, and the remaining regions were skeletonized to extract connected ridges. On each connected ridge, a shortest path search is performed from the left end mileage to the right end mileage in one-unit steps, with a permissible vertical step of one unit and a curvature limit of no more than 5 degrees per 0.02 meters, thus obtaining a set of continuous layer boundaries. After adjacent layer boundaries enclose the intralayer region, the envelope intensity within each intralayer region is binned according to five quantiles: 0.10, 0.30, 0.50, 0.70, and 0.90. These bins are then mapped to volume fractions of 0.05, 0.15, 0.25, 0.35, and 0.45, respectively. If the location also belongs to a strong response band, the corresponding volume fraction is increased by 0.05, not exceeding 0.60, forming the initial water content field within the layer. This process prioritizes locations where energy concentration, phase stability, and strong spectral relaxation occur simultaneously as the basis for interfaces and high water content initial values. This ensures continuity in the mileage direction and avoids excessive fluctuations in the depth direction, providing a stable starting point for subsequent iterations.

[0042] Example 3: Before constructing the layered initial model, data standardization is performed. For the full waveform of the seismic wave, DC is first linearly removed from each channel according to the recording time window, and then out-of-band noise is suppressed by bandpass filtering. The typical passband is 100 MHz to 1200 MHz. Then, the first arrival template is generated by the reference channel of the metal plate collected in the field. The earliest and largest paired peaks are located by sliding correlation. After alignment by channel, the first arrival template is subtracted to improve the contrast of shallow reflection. Next, the amplitude equalization between channels is performed. The scaling factor is calculated according to the 95th percentile amplitude of each channel so that the difference of the percentile amplitude of adjacent channels is less than 10%. Finally, the arrival time of the fixed trigger mark is used as the zero point of time, with an error of no more than 1 microsecond, to obtain the standardized full waveform. For the complex impedance and induced polarization spectrum, the geometry was first corrected based on the actual electrode spacing and contact posture, with the spacing error corrected to within 0.01 meters. Then, amplitude and phase uniformity was performed, using the amplitude and phase at the 0.1 Hz and 10000 Hz endpoints as a reference. Linear interpolation was used to align the amplitude and phase curves within the same mileage window, ensuring that the endpoint difference did not exceed 1 degree and 1 / 100th of a degree. Next, the unequally spaced original frequency points were resampled to a logarithmically uniform grid, resulting in 41 frequency points. The frequency point spacing was distributed at 10 points per decibel, yielding a standardized spectrum. To reduce power frequency interference, a notch width of 2 Hz was set in the 50 Hz and 100 Hz neighborhoods, and three-point median filtering was performed on three adjacent frequency points to preserve the geometry of the Nyquist trajectory.

[0043] The entire standardization process is triggered at a mileage step of 0.02 meters. The same trigger number corresponds to a standardized full waveform and a set of standardized spectral data simultaneously, with a timestamp resolution of 1 microsecond. Through the above processing, the seismic wave full waveform is comparable in phase and amplitude, and the complex impedance and induced polarization spectrum are superimposed in amplitude, phase, and frequency. The two can be directly registered on the mileage-depth grid and proceed to subsequent steps.

[0044] Example 4: When updating the intralayer water content field and layer boundaries based on the full waveform of seismic waves, a one-dimensional reflection sequence is first generated under the current layer boundary and intralayer water content field. The generation method is as follows: a strong reflection strip is established at each layer boundary location, and a weak reflection background strip is established for the intralayer region. The amplitude ratio of the strips is determined by the relative strength of the difference between the current intralayer water content field and the adjacent layer, so that the energy is concentrated at the interface and uniform within the layer. Subsequently, the zero-phase source wavelet is estimated from the standardized full waveform through homomorphic deconvolution. The zero-phase source wavelet is convolved with the one-dimensional reflection sequence to obtain the synthetic full waveform. The synthetic full waveform and the standardized full waveform are cross-correlated and aligned channel by channel, outputting the envelope difference map, phase difference map, and instantaneous frequency difference map. A pyramid sequence from low-frequency sub-bands to high-frequency sub-bands is constructed, with four sub-bands, and the center frequencies successively covering 100 MHz to 1200 MHz. Within each sub-band, deviation zones are detected using windows of 0.20 meters in the mileage direction and 0.05 meters in the depth direction: when the envelope difference map is continuously higher than the local median plus twice the local absolute deviation within this window and the median absolute deviation of the phase difference map is greater than 5 degrees, this window is marked as a deviation zone.

[0045] For each deviation zone, update blocks are defined at intervals of 0.20 meters along the mileage direction and 0.05 meters along the depth direction. Within each update block, the top 30 percentile cells with the highest envelope peak intensity are processed first. The water content field of these cells is uniformly adjusted by a step size of 0.01, and the difference map is recalculated. If the median absolute values ​​of the envelope difference and phase difference both decrease within the update block, they are retained; otherwise, fine-tuning is performed in the top 50 percentile cells with a step size of 0.005. Each update block undergoes a maximum of 12 consecutive fine-tunings, with sub-bands completed sequentially from low to high, ensuring that large-scale time shifts are corrected first, followed by compensation for small-scale details. After the sub-bands are completed, the full-frequency results are used for straightening the entire segment: using the earliest and largest paired peaks as the first arrival anchor points, all channels within the same scan line are aligned to the same time reference. Then, the layer boundary is repositioned at the location where the layer boundary crosses the envelope extreme value zone, prioritizing the passage of the layer boundary through the local maximum of the envelope peak. The reason for choosing the local maximum is that reflected energy accumulates at the interface, and the energy distribution of the cross-boundary unit is more stable, which can improve the repeatability accuracy and mileage continuity of the layer boundary.

[0046] Example 5: After locally updating the intralayer water content field within the corresponding intralayer region, the quantile distributions of the envelope time shift and phase difference within the deviation region are further calculated, with quantile positions set at 0.20, 0.50, and 0.80. The two quantile distributions are paired to form a time shift-phase pairing table. The meaning of pairing is as follows: at the same quantile position, if the envelope time shift and phase difference are in the same direction and their amplitudes are both within the corresponding quantile interval, it indicates that the propagation delay and phase lag within the update block have a consistent trend. In this case, the intralayer water content field of all units within the update block is adjusted in the same direction with a step size of 0.01. If the envelope time shift and phase difference are in the same direction but their amplitudes are inconsistent, then only the top 30 percentile units with the highest envelope peak intensity are adjusted with a step size of 0.01. If the envelope time shift and phase difference are in opposite directions, then the top 50 percentile units within the update block are finely adjusted in the opposite direction with a step size of 0.005. After each adjustment, immediately recalculate the envelope difference map and phase difference map. If the absolute value of the median of both decreases, retain them; otherwise, roll back to the previous state and reduce the step size to 0.005 and try again, up to a maximum of 3 times.

[0047] After completing the pairwise adjustment, phase anchoring-guided straightening is performed: within the same scan line, the earliest and largest paired peaks are used as the reference, aligned to a unified time reference, and the time-to-depth conversion relationship is straightened based on the mileage distribution of the phase difference, making the phase horizontal lines of adjacent channels tend to be parallel in the mileage direction. After straightening, the layer boundary is repositioned so that the intersection of the layer boundary and the envelope extreme value zone forms a continuous curve in the mileage direction. If the depth difference between two adjacent intersection points exceeds 0.02 meters, linear interpolation is used to supplement the points between the two points to maintain the curvature not exceeding 5 degrees per 0.02 meters. This process transforms the ordering relationship of propagation time shift and phase lag into an executable step rule through quantization pairing, avoiding the uncertainty caused by subjective thresholds. At the same time, by unifying the time-to-depth scale through phase anchoring, the geometric stretching in the mileage direction is reduced, ultimately obtaining layer boundary and intralayer water cut field update results that are consistent with the full waveform morphology and continuous in mileage.

[0048] refer to Figure 5This diagram illustrates the quantile pairing update mechanism for envelope peak time shift and phase difference. The left side of the diagram shows the envelope peak time shift quantile curve, and the right side shows the phase difference quantile curve. The two are connected by a dashed arrow to indicate the quantile pairing relationship. In the coordinate system on the left, the horizontal axis represents the quantile position, with values ​​at three key quantile points: 0.2, 0.5, and 0.8. The vertical axis represents the time shift, in nanoseconds (ns), ranging from -2.0 to +2.0. The time shift quantile curves are obtained by statistically analyzing the time shift of all envelope peaks within the deviation region. Quantile values ​​are calculated at the three quantile positions and marked with solid circles. The curve starts at the lower left (100, 362), corresponding to the 0.2 quantile and a time shift Q0.2 = -1.0 ns, indicating that the envelope peak at this quantile arrives earlier than the synthesized waveform. It passes through the midpoint (225, 200), corresponding to the 0.5 quantile and a time shift Q0.5 = +0.6 ns; extending to the upper right (350, 140), corresponding to the 0.8 quantile and a time shift Q0.8 = +1.2 ns, indicating that the envelope peak at higher quantiles arrives relatively later. In the coordinate system on the right, the horizontal axis also represents the quantile positions of 0.2, 0.5, and 0.8, and the vertical axis represents the phase difference in degrees, ranging from -10° to +10°. The phase difference quantile curve is obtained by statistically analyzing all phase difference values ​​within the same deviation region, calculating the quantile values ​​at each of the three quantile positions, and marking them with solid circles. The curve starts at the lower left point (450, 362), corresponding to the 0.2 quantile and a phase difference P0.2 = -5.2°, indicating that the phase at this quantile is relatively ahead. Passing through the midpoint (575, 195), it corresponds to the 0.5 quantile and a phase difference P0.5 = +3.8°. Extending to the upper right point (700, 135), it corresponds to the 0.8 quantile and a phase difference P0.8 = +7.5°, indicating that the phase at higher quantiles is relatively lagging. Dashed arrows connect two feature points at the same quantile, forming quantile pairings. At the 0.2 quantile, the time shift and phase difference are negative, indicating that propagation is advanced and the phase is ahead within this interval, belonging to the same direction and amplitude case. At the 0.5 and 0.8 quantiles, the time shift and phase difference are positive, indicating that propagation is delayed and the phase is lagging, also belonging to the same direction and amplitude case. The figure below illustrates three rules for the quantile pairing update strategy: When the same direction and amplitude appear at the 0.2 quantile, it indicates that the water content in the lower quantile range is overestimated, and the water content of all cells within the update block is reduced by 0.01; when the same direction and amplitude appear at the 0.5 and 0.8 quantiles, it indicates that the water content in the higher quantile range is underestimated, and the water content of all cells is increased by 0.01; when the same direction and different amplitudes appear, adjustments are only made on the top 30% of high-intensity cells. This mechanism transforms the time shift and phase ordering relationship into an executable water content update step rule.

[0049] Example 6: When updating the water content field within a layer based on the quantile alignment of spectral fingerprints and envelope fingerprints using complex impedance and induced polarization spectra, the normalized spectrum is first divided into five fixed frequency bands according to logarithmic frequencies, ranging from 0.1 Hz to 1 Hz, 1 Hz to 10 Hz, 10 Hz to 100 Hz, 100 Hz to 1000 Hz, and 1000 Hz to 10000 Hz. Circular fitting is performed on the Nyquist trajectory of each frequency band, with at least 10 fitting points. After fitting, relaxation components are stripped using a distance threshold set to 5% of the radius of the fitted circular arc for that segment. The fitting is repeated on the remaining point set until the number of residual points is less than 20% of the total number of points in that segment or the standard deviation of the residuals is less than 3% of the fitted arc length. The radius of the circular arc and the area enclosed by the trajectory are recorded for each stripped frequency band as two fundamental quantities of the spectral fingerprint. Simultaneously, the peak position and peak width of the same frequency band are recorded on the amplitude-phase curve as auxiliary quantities of the spectral fingerprint.

[0050] Within each layer region, empirical quantile curves for the arc radius and the area enclosed by the trajectory are calculated, with quantile positions uniformly set to 0.10, 0.30, 0.50, 0.70, and 0.90. Within the same region, quantile curves for the envelope intensity are also calculated, with consistent quantile positions to ensure direct correspondence between the two domains. The two sets of quantile curves are aligned to form a monotonic mapping table of spectral fingerprint-envelope fingerprint: when both the arc radius and the area enclosed by the trajectory at a certain quantile position are at high quantiles, the corresponding envelope intensity quantile is also set to a high quantile; when both are at low quantiles, the corresponding envelope intensity quantile is set to a low quantile; when one is high and the other low, the nearest quantile principle is used to ensure monotonicity. Based on this table, the spectral fingerprint of each voxel is converted into a water content increment: increasing by 0.02 when located in the high quantile range, decreasing by 0.02 when located in the low quantile range, and remaining unchanged when located in the middle quantile range.

[0051] To avoid the spread of isolated anomalies, increments are applied only within connected segments with a continuous length of at least 0.20 meters in the mileage direction and a thickness of at least 0.03 meters in the depth direction. If the fitting residual of a certain frequency band exceeds a threshold, the entry for that frequency band is removed, and only the entries for the remaining reliable frequency bands participate in the conversion. After each conversion, the ranking of envelope intensity within the layer is checked to ensure it is aligned with the updated water content ranking. If a reverse segment appears and its length exceeds 0.10 meters, the increment of that segment is halved to 0.01 and checked again until the ranking is aligned or the segment length is reduced to no more than 0.10 meters. Through this process, the change in the water content field within the layer progresses along the common enhancement trend of the two domains and does not depend on external weights or historical data.

[0052] refer to Figure 2 ,like Figure 2The figure shows a schematic diagram of the Nyquist trajectory and circular arc fitting characteristics of complex impedance. The horizontal axis represents the real part Z' (Ω·m), and the vertical axis represents the imaginary part -Z'' (Ω·m). The figure shows two typical Nyquist trajectories, corresponding to roadbed materials with different moisture contents. The solid line trajectory represents a high moisture content sample with a moisture content θ=0.35. This trajectory exhibits a large circular arc shape, with the starting point near (540, 495) at a frequency of 0.1Hz, marked with a solid circle; the ending point is near (150, 495) at a frequency of 10000Hz, marked with a solid triangle. During the frequency scan from 0.1Hz to 10000Hz, the trajectory extends from right to left, with the direction of frequency increase indicated by arrows. A circular arc fitting is performed on this solid line trajectory, with the fitting center located at (345, 500), a fitting radius R1 of 160Ω·m, and a trajectory envelope area A1 of 8050Ω²·m². The dashed trajectory represents a low-moisture sample with a water content of θ=0.15. This trajectory exhibits a small arc shape, with the starting point near (370, 492), corresponding to a frequency of 0.1 Hz, marked by a solid gray circle; the ending point is near (180, 497), corresponding to a frequency of 10000 Hz, marked by a solid gray triangle. A circular arc fitting is performed on this dashed trajectory, with the fitting center at (275, 500), a fitting radius R² of 100 Ω·m, and a trajectory envelope area A² of 3140 Ω²·m². Comparing the two trajectories reveals that the sample with higher water content has a larger arc radius and a larger trajectory envelope area. This is because increased water content leads to an increase in the porous liquid phase, resulting in a more significant electrochemical relaxation process, manifested as a larger arc response on the Nyquist plane. The arc radius and trajectory envelope area, as fundamental quantities in spectral fingerprinting, can quantitatively characterize the water content level of a material. In the dual-domain quantization alignment method of the present invention, circular arc fitting is performed on five logarithmic frequency bands (0.1-1Hz, 1-10Hz, 10-100Hz, 100-1000Hz, 1000-10000Hz) respectively, and the radius of the circular arc and the area of ​​the trajectory envelope of each frequency band are extracted to form a complete spectral fingerprint feature, which is used for subsequent quantization alignment with the GPR envelope fingerprint.

[0053] refer to Figure 3 ,like Figure 3The figure shows a schematic diagram of the quantile alignment mapping relationship between spectral fingerprint and envelope fingerprint. The left side of the figure shows the envelope intensity quantile curve in the GPR domain, and the right side shows the arc radius quantile curve in the spectral domain. The two are connected by a dashed line to represent the quantile alignment mapping relationship. In the coordinate system on the left, the horizontal axis represents the quantile, ranging from 0.1 to 0.9, and the vertical axis represents the envelope intensity, ranging from 0.0 to 0.8. The envelope intensity quantile curve is obtained by statistically analyzing the full waveform envelope intensity of seismic waves within the layer region. Quantile values ​​are calculated at quantile positions of 0.1, 0.3, 0.5, 0.7, and 0.9, forming five feature points marked with solid circles. The curve starts from the lower left point (100, 490), corresponding to low quantile and low envelope intensity, and extends to the upper right to the point (350, 130), corresponding to high quantile and high envelope intensity, showing a monotonically increasing trend. In the coordinate system on the right, the horizontal axis also represents the quantiles, ranging from 0.1 to 0.9, and the vertical axis represents the radius of the arc, in Ω·m, ranging from 50 to 250. The arc radius quantile curve is obtained by statistically analyzing the arc radii obtained by fitting the complex impedance Nyquist trajectory within the same layer. Quantile values ​​are calculated at the same quantile positions of 0.1, 0.3, 0.5, 0.7, and 0.9, forming five feature points marked with solid circles. The curve starts from the lower left point (450, 485), corresponding to the low quantile and small arc radius, and extends to the upper right to the point (700, 125), corresponding to the high quantile and large arc radius, also showing a monotonically increasing trend. The two sets of quantile curves are connected by dashed arrows, indicating the mapping relationship of quantile alignment. At the same quantile position, a one-to-one correspondence is established between the envelope strength quantile on the left and the arc radius quantile on the right. For example, at the 0.1 quantile, a low envelope strength corresponds to a small radius of curvature; at the 0.9 quantile, a high envelope strength corresponds to a large radius of curvature. This correspondence reflects the common enhancing trend of electromagnetic reflection intensity and electrochemical relaxation response as water content increases.

[0054] Example 7: When performing total variation regularization on the updated intralayer water content field to maintain the transition characteristics of the layer boundary, the process is performed separately for each intralayer region. First, a first-order difference map is calculated in the mileage direction, and then a first-order difference map is calculated in the depth direction. For each difference map, the upper quartile is counted, and the differences below the upper quartile are set to zero. The differences above the upper quartile are linearly reduced to the midpoint between the upper quartile and the maximum value of the map, resulting in a shrunken difference map.

[0055] This contraction strategy can suppress fine undulations within the layer while preserving obvious transitions near the layer boundaries. Subsequently, the contraction difference map along the mileage direction is reconstructed by integration, with the water content value at the left boundary of the region within the layer selected as the anchor point for integration, and the data is accumulated layer by layer along the mileage direction. Then, the contraction difference map along the depth direction is reconstructed by integration, with the water content value at the upper boundary of the region within the layer selected as the anchor point for integration, and the data is accumulated layer by layer along the depth direction. To reduce drift caused by unidirectional integration, the two reconstruction results are arithmetically averaged point-to-point to obtain a preliminary draft of the intralayer water content field after total variation regularization. A small-scale correction is performed on the preliminary draft within each sub-block, with the sub-block size being 0.10 meters along the mileage direction and 0.02 meters along the depth direction. The correction rule is: if the difference between the maximum and minimum values ​​within a sub-block is less than 5% of the upper quartile of that region, it remains unchanged; if it exceeds this threshold, the extreme values ​​within the sub-block are contracted towards the sub-block mean by 20%.

[0056] After correction, the transitions on both sides of the layer boundary are checked again to see if they have been weakened. The method is to take an adjacent cell in the depth direction on each side of the layer boundary and calculate the difference between the two cells. If the average difference decreases by more than 20% compared to before correction, the correction of that sub-block is cancelled. After completing the above process, the final intralayer water content field after total variation regularization is obtained, and this result is sent to the next step of interface candidate point extraction.

[0057] Example 8: After restoring the intralayer water content field after total variation normalization, the candidate interface points and continuous layer boundaries are determined. First, the gradient magnitude map in the depth direction is calculated on the reconstructed intralayer water content field, and the local upper quartiles are statistically analyzed using a 0.10-meter window in the mileage direction. Connected strips not lower than the local upper quartiles are taken as the extreme value zones of the water content gradient.

[0058] Simultaneously, on the envelope diagram of the full seismic wave waveform, envelope extreme value zones are formed at the top three local maxima of each trace, and the highest energy trace is retained within a 0.10-meter window along the mileage direction, resulting in continuous envelope extreme value zones along the line. The intersection of the two types of zonal regions yields a discrete set of points in the overlapping area, serving as the candidate point set for the interface. To construct a continuous and smooth layer boundary from the left mileage to the right mileage, dynamic programming is performed on the candidate point set: each mileage step allows the selection of only candidate points whose depth difference does not exceed one unit; if multiple candidate points exist in the same mileage step, they are selected and retained in the following priority order: higher envelope strength, larger water content gradient transition, and smaller depth difference from the previous mileage step point. To avoid abrupt curvature changes, the angle change of the broken line formed by three adjacent mileage steps is constrained to no more than 5 degrees per 0.02 meters; if this is exceeded, an interpolation point in the depth direction is inserted in the middle step, or the selection reverts to the second-best candidate point.

[0059] For sections with fractures, if the fracture length does not exceed 0.10 meters, it is completed by linear interpolation based on the depths of the two endpoints; if the length exceeds this, the section is marked as low confidence and will be re-evaluated in subsequent rounds. After completing the full-segment search, a continuous and smooth optimal layer boundary is obtained, and this boundary replaces the original layer boundary. This strategy limits the candidate range by the overlap of the envelope extreme value zone and the water content gradient extreme value zone, ensuring that the selected path simultaneously satisfies the requirements of concentrated reflection energy and significant water content transition. Dynamic programming guarantees overall optimality and geometric smoothness across the entire segment.

[0060] Example 9: In the convergence determination during iteration, a convergence check is uniformly performed after each round of dual-domain quantile alignment interface iterations. The first check is the maximum displacement detection of the layer boundary: calculate the depth difference between the current layer boundary and the previous layer boundary at each mileage position, take the maximum absolute value, and if the maximum value is less than the displacement threshold of 0.02 meters determined by the depth step, then the check passes; otherwise, it fails.

[0061] The second test is the difference map change amplitude test: calculate the maximum absolute value change amplitude on the envelope difference map and the phase difference map respectively. The change amplitude is defined as the maximum absolute value of the difference between the current difference map and the previous difference map. Compare the change amplitude of each with the corresponding tenths place value of their respective distributions. If both are less than the corresponding tenths place value, the test is passed.

[0062] Convergence is determined only when both conditions are met simultaneously, and the result is output. If either condition fails, the next iteration begins. To avoid misjudgments due to occasional fluctuations, a verification round is performed after the initial condition is met. If both conditions are met in the verification round, convergence is finally confirmed. If both conditions are not met for 10 consecutive rounds, the non-converged mileage interval is located, with the interval length in 0.50-meter units. The intralayer water content field and layer boundary of this interval are locally reinitialized, and the cycle for this interval restarts. The existing results are retained for the remaining intervals. This approach ensures strict stopping criteria while limiting the iteration cost in the worst-case scenario.

[0063] Example 10: After iterative convergence, the results are stitched together. The layer boundaries of adjacent data blocks are overlapped with the intralayer water content field according to mileage order, with an overlap length of 0.20 meters to 0.40 meters. In the overlap area, the layer boundaries are first aligned at the endpoints: using the depths of the left and right endpoints of the overlap area as anchor points, if the depth difference between the two layer boundaries at the anchor points does not exceed 0.02 meters, a direct linear transition is performed; if it exceeds 0.02 meters, the difference is evenly distributed within the overlap area in three proportions: one-quarter, one-half, and three-quarters, to avoid abrupt changes. For the intralayer water content field, a linear transition is performed in the overlap area according to the mileage direction, linearly interpolating the values ​​of the previous and subsequent data blocks according to their mileage positions. To suppress stitching marks, after the transition is completed, a median smoothing with a sliding window is performed on the entire intralayer water content field, with a window length of 0.20 meters and a step size of 0.02 meters. To maintain the transition characteristics of the layer boundaries, one depth unit is retained on each side of the layer boundary during smoothing and is not included in the window statistics. For potential voids and overlapping coverage areas, a priority rule is applied: if the same cell has two valid values, the arithmetic mean of the two is taken; if a cell is empty and has valid values ​​on both sides within 0.10 meters of each other along the mileage direction, it is filled by linear interpolation on both sides; if the length exceeds this, it is marked as low confidence and noted in the results table. After stitching, a consistency check is performed on the entire segment, including whether the curvature of the layer boundary is continuous, whether the layer thickness variation is smooth within 1 meter, and whether there are isolated gradient peaks exceeding 0.10 in the water content field within the layer. If any non-compliance is found, the median smoothing window is reduced to 0.10 meters in the corresponding interval and the transition is redone until all checks pass. The final output includes a continuous profile, continuous layer boundary vectors, and a layered water content table, ensuring seamless connection of adjacent data blocks in mileage and depth.

[0064] The present invention has been described in detail above. Specific examples have been used to illustrate the principles and implementation methods of the invention. The descriptions of the embodiments above are merely for the purpose of helping to understand the method and core ideas of the present invention. It should be noted that those skilled in the art can make various improvements and modifications to the present invention without departing from its principles, and these improvements and modifications also fall within the protection scope of the claims of the present invention.

Claims

1. A layered inversion method for formation water content based on shallow seismic full waveform, characterized in that, Includes the following steps: Step 1: Under the same survey line conditions, simultaneously acquire the full waveform, complex impedance, and induced polarization spectrum of the seismic wave; Step 2: Construct a layered initial model that includes the layer boundary and the initial value of the water content within the layer; and based on the layered initial model, perform a two-domain quantile alignment interface maintenance iteration until convergence to obtain the final layer boundary and the water content field within the layer. The dual-domain stratification alignment interface maintains iteration, including updating the intralayer water content field and layer boundaries based on the full waveform of the seismic wave. Based on the complex impedance and induced polarization spectrum, the intralayer water content field is updated by quantile alignment of spectral fingerprint and envelope fingerprint; and total variation regularization is performed on the updated intralayer water content field to maintain the transition characteristics of the layer boundary.

2. The method according to claim 1, characterized in that, The steps for constructing a layered initial model that includes layer boundaries and initial values ​​of water content within the layers include: obtaining candidate reflection zones by peak-valley tracking based on the envelope peak, phase plateau, and instantaneous frequency steps of the full seismic wave waveform; obtaining strong response zones by frequency band clustering based on the arc shape and amplitude-phase peak shape of the complex impedance and induced polarization spectrum; performing voxel-level registration between the candidate reflection zones and the strong response zones; extracting connected ridges through connected component analysis; and generating a set of continuous layer boundaries as layer boundaries using shortest path search.

3. The method according to claim 1, characterized in that, Before constructing the layered initial model, a data standardization step is also included, which includes: performing DC removal, bandpass filtering, first arrival template subtraction, inter-channel amplitude equalization and time zero-point correction on the full waveform of the seismic wave in sequence to obtain the standardized full waveform; and performing electrode geometry correction, amplitude and phase consistency and frequency resampling on the complex impedance and induced polarization spectrum to obtain the standardized spectrum.

4. The method according to claim 1, characterized in that, The process of updating the intralayer water content field and layer boundary based on the full waveform of seismic waves specifically includes: generating a one-dimensional reflection sequence under the current layer boundary and intralayer water content field, and using zero-phase source wavelets to convolve and synthesize the full waveform; cross-correlating and aligning the synthesized full waveform with the standardized full waveform to obtain the envelope difference map, phase difference map, and instantaneous frequency difference map; constructing a pyramid sequence from the low-frequency sub-band to the high-frequency sub-band, and matching and tracking according to the sub-band order to locate the deviation area in the difference map, and locally updating the intralayer water content field in the corresponding intralayer region.

5. The method according to claim 4, characterized in that, After the process of locally updating the water content field within the corresponding layer region, the process also includes: calculating the quantile distribution of the envelope time shift within the deviation region, pairing it with the quantile distribution of the phase difference to generate a time-shift-phase pairing table, adjusting the water content field within the layer based on the time-shift-phase pairing table, and straightening the time-to-depth conversion relationship based on the phase anchoring results to update the layer boundary.

6. The method according to claim 1, characterized in that, Based on the complex impedance and induced polarization spectrum, the process of updating the intralayer water cut field through quantile alignment of spectral fingerprint and envelope fingerprint specifically includes: dividing the normalized spectrum into several frequency bands according to the logarithmic frequency; performing circular arc fitting on the Nyquist trajectory within each frequency band and removing the relaxation component; within each intralayer region, statistically analyzing the arc radius obtained from the circular arc fitting and the area enclosed by the trajectory as the spectral fingerprint, and statistically analyzing the empirical quantile curve of the spectral fingerprint; quantile alignment of the empirical quantile curve of the spectral fingerprint with the quantile curve of the envelope intensity to form a monotonic mapping table of spectral fingerprint-envelope fingerprint; and based on the monotonic mapping table of spectral fingerprint-envelope fingerprint, converting the spectral fingerprint of each voxel into a water cut increment and updating the intralayer water cut field.

7. The method according to claim 1, characterized in that, The process of performing total variation regularization on the updated intralayer water cut field specifically includes: calculating the first-order difference map in the mileage direction and the depth direction within each intralayer region; using the upper quartile of each first-order difference map as the shrinkage threshold, setting the difference below the threshold to zero, and reducing the amplitude of the difference above the threshold to obtain the shrunken difference map; performing integral reconstruction on the shrunken difference map, and using the water cut value at the boundary of the intralayer region as the reconstruction anchor point to restore the intralayer water cut field after total variation regularization.

8. The method according to claim 7, characterized in that, After the process of restoring the intralayer water content field after total variation regularization, the process also includes: extracting the intersection points of the envelope extreme value zone and the water content gradient extreme value zone on the reconstructed water content field to form a set of interface candidate points; performing dynamic programming on the set of interface candidate points and constraining the mileage difference and curvature between adjacent points to obtain a continuous and smooth optimal layer boundary, and replacing the original layer boundary.

9. The method according to claim 1, characterized in that, The convergence criterion for iteration until convergence is: the maximum displacement of the layer boundary between two cycles is less than the displacement threshold determined by the mesh resolution, and the maximum change amplitude between the envelope difference map and the phase difference map is less than the amplitude threshold determined by the difference map distribution.

10. The method according to claim 1, characterized in that, After iterative convergence, the results are further stitched together, including: overlapping and stitching the layer boundaries of adjacent data blocks with the intralayer water content field in mileage order; aligning the overlapping areas by endpoints and performing a linear transition; and performing median smoothing of a sliding window on the intralayer water content field to form a continuous profile.

Images