Two-dimensional sea surface ghost wave water body imaging measurement method, system, terminal and flow measuring device

Abstract

The application belongs to the technical field of image processing, and discloses a two-dimensional sea surface ghost wave water body imaging measurement method, system, terminal and flow measuring equipment. Frequency domain positive transformation data of seismic data is obtained; for any frequency component and ray parameter, the delay time of virtual reflection and primary reflection is calculated to obtain the first downgoing wave field frequency domain Radon transformation result; the average delay time and average amplitude difference coefficient of the virtual reflection and the primary reflection corresponding to the ray parameter are obtained and calculated again, and the iteration is alternately and continuously performed until the convergence condition is met, so that the final frequency domain Radon transformation result corresponding to the virtual reflection is obtained, and then frequency domain Radon inverse transformation and time domain Fourier inverse transformation are performed to obtain the final virtual reflection prediction data. The application can perform full waveform inversion on useful information such as virtual reflection signals received by the ADCP to obtain property information of the surface seawater, and accurately image the water body.

Description

Technical Field

[0001] This invention belongs to the field of image processing technology, and in particular relates to a two-dimensional sea surface ghost wave water imaging measurement method, system, terminal and flow measurement equipment. Background Technology

[0002] Over the past 30 years, ADCPs have evolved into a series of products with various forms and applications. Based on transducer array structure, they can be categorized into conventional arrays consisting of four circular transducers and phased arrays composed of nearly a thousand small transducers. Based on application, they can be classified as self-capacitive ADCPs, direct-reading / mobile ADCPs (mobile ADCPs all utilize phased array technology), and low-profile ADCPs (or simply LowADCPs).

[0003] Based on signal detection technology, the research and development of ADCP can be divided into three stages. The first stage is characterized by phase-locked loop frequency trackers, which have low spatial resolution and accuracy in current measurement. The second stage uses complex correlation techniques to calculate the Doppler frequency shift of the signal, such as the Narrowband Acoustic Doppler Current Profiler (NBADCP) from RD Corporation. Its advantages include simultaneous measurement of multiple current layers, small size, simple equipment, good stability, and long-term operation. Its disadvantages are limitations imposed by frequency and pulse width, restricting the profile thickness and resulting in large variance in single measurements. The third stage employs broadband acoustic Doppler current profiler (BBADCP) technology, characterized by broadband pulse-pair coherence, emitting a pseudo-random coded pulse composed of N symbol pairs. The profile measurement depth depends on the overall pulse length and frequency, while the measurement layer thickness depends on the symbol length. Because the width of the symbol is much smaller than the width of the pulse, each symbol corresponds to a current measurement, thus greatly reducing the stratum thickness measured by BBADCP and decreasing the variance of a single measurement to 1N times that of NBADCP. Therefore, BBADCP is widely welcomed by the oceanographic community.

[0004] ADCP is a novel flow measurement device developed using the Doppler effect of sound waves. It can measure both relative bottom velocity and relative current velocity. However, there is currently no existing technology for imaging the water body between the seismic source and the sea surface based on ADCP technology.

[0005] Based on the above analysis, the existing technology has the following problems and shortcomings: there is room for further optimization in the processing and inversion resolution and accuracy of virtual reflection signals. The existing technology cannot invert useful information such as virtual reflection signals received by ADCP, resulting in inaccurate imaging of the surface seawater.

[0006] The difficulty in solving the above problems and defects lies in establishing a connection between the ADCP information system and conventional virtual reflection, and accurately inverting the image using virtual reflection information.

[0007] The significance of solving the above problems and defects is that it extracts and utilizes the virtual reflection information from the conventional shielding and filtering, thereby improving the signal-to-noise ratio of the signal and providing richer signal information for high-precision analysis. Summary of the Invention

[0008] To overcome the problems existing in related technologies, the present invention discloses a two-dimensional sea surface ghost wave water body imaging measurement method, system, terminal, and current measurement equipment. Specifically, it relates to a water flow profile imaging method in marine exploration.

[0009] The technical solution is as follows: A two-dimensional sea surface ghost wave water body imaging measurement method, comprising:

[0010] The virtual reflection signal rich in water body information received by the current measurement device ADCP is subjected to full waveform inversion to obtain the property information of the surface seawater and to accurately image the water body.

[0011] In one embodiment, the two-dimensional sea surface ghost wave water body imaging measurement method specifically includes:

[0012] (1) For the earthquake record S(x) i Performing a Fourier transform on the seismic data (t) in the time direction yields the vector S(x, ω). Using a high-precision Radon transform, the frequency domain Radon forward transform data U of the seismic data is obtained. s (p k ,ω);

[0013] (2) For any frequency component ω and ray parameter p k Given the amplitude difference coefficient R between virtual reflection and primary reflection. k = -1, and calculate the delay time Δτ between the virtual reflection and the first reflection. k Construct the matrix G(x) i p k The first Radon transform result of the downlink wave field in the frequency domain is obtained by solving for ω).

[0014] (3) The ray parameter p is obtained by alternating iterative inversion of the objective function. k The corresponding virtual reflection and the average delay time of a single reflection and the coefficient of difference in average amplitude

[0015] (4) Reuse step (3) to obtain and The downlink wave field is obtained using step (2).

[0016] (5) Calculate the average delay time of the virtual reflection and the first reflection again. and the coefficient of difference in average amplitude By iterating repeatedly until the convergence condition is met, the final Radon transform result U in the frequency domain corresponding to the virtual reflection is obtained. Down (p, ω);

[0017] (6) For U Down Performing an inverse Radon transform in the frequency domain and an inverse Fourier transform in the time domain on (p, ω) yields the final virtual reflection prediction data u. Down (x, t).

[0018] In one embodiment, step (3) yields the ray parameter p. k The corresponding virtual reflection and the average delay time of a single reflection and the coefficient of difference in average amplitude include:

[0019] During data acquisition, the depth of the i-th detector is z. i The offset distance is x i p represents the ray parameter matched to the detector;

[0020] Among them, the k-th ray parameter θ k V is the angle of incidence of the seismic wave as it propagates to sea level, and the angle of incidence for virtual reflection; w Let be the seawater velocity; the seismic wave field received at this detector is the sum of the ascending and descending wave fields, then...

[0021] u S (p k ,τ)=u Up (p k ,τ)+u Down (p k ,τ)

[0022]

[0023] Among them, u S For the total wave field, u Up For the upward wave field, u Down R represents the downlink wavefield; τ is the time the seismic wave travels at offset 0; Δτ is the time difference between the uplink and downlink wavefields; and R... k It is the ratio coefficient between the amplitude of the ghost wave and the amplitude of the first reflected wave;

[0024] Based on the geometric relationship of seismic wave propagation, the time difference Δτ between the downgoing and upgoing wave fields is:

[0025]

[0026] Furthermore, since the total received wave field is represented by Nm horizontal slowness degrees p:

[0027]

[0028] The expressions for the total wave field and the downlink wave field corresponding to the virtual reflection in the frequency domain are:

[0029]

[0030] S(x i , where ω) is the Fourier transform result of the seismic data in the time direction. If:

[0031]

[0032] The least squares solution for the virtual reflection wave field is:

[0033] U Down (p, ω) = (G H (x, p, ω)G(x, p, ω)) -1 G(x, p, ω) H U s (x, ω)

[0034] Ignoring shot-receiver offset error, the unknown parameters include the time difference Δτ between the downlink and primary wavefields, the energy ratio coefficient R between virtual and primary reflections, and the downlink wavefield U to be determined. Down (p, ω); To obtain the above three parameters, an alternating iterative inversion method is used, fixing the time difference Δτ between the downlink and primary wave fields, the energy ratio coefficient R between virtual and primary reflections, and the downlink wave field U to be determined. Down From the two parameters in (p, ω), one parameter is inverted, and then the inverted parameter is used as a constraint value to invert the other parameter. This process is repeated to obtain the optimal solution.

[0035] In one embodiment, an objective function for the inversion needs to be determined before starting the alternating iterative inversion process; the relationship between the ascending wave field and the descending wave field is expressed in two forms:

[0036] u S (p k ,τ)=u Up (p k ,τ)+u Down (p k ,τ)

[0037]

[0038] Performing a Fourier transform yields:

[0039] U Up (p k ,ω)=U s (pk ,ω)-U Down (p k ,ω)

[0040]

[0041] Among them, U s (p k ω) is the frequency domain positive transform of the total wave field; then the target of alternating inversion obtains the downlink wave field through the above. The formula yields the same first-order wave field, i.e.

[0042]

[0043] In one embodiment, the two-dimensional sea surface ghost wave water body imaging measurement method further includes:

[0044] Virtual reflections incorporate the seismic response of underground structures. The detector is artificially moved to a mirror position above the sea level relative to the real detector, and the velocity on the sea surface is equivalent to the velocity of seawater. The virtual reflection is treated as a new primary reflection and migrated using conventional migration imaging methods. At the same time, the virtual reflection imaging result is matched and correlated with the primary reflection imaging result to obtain a joint imaging result of virtual reflection and primary reflection.

[0045] Another object of the present invention is to provide a two-dimensional sea surface ghost wave imaging measurement system for implementing the aforementioned two-dimensional sea surface ghost wave water body imaging measurement method, the two-dimensional sea surface ghost wave water body imaging measurement system comprising:

[0046] The earthquake data frequency domain forward transform data acquisition module is used to process earthquake records S(x) i Performing a Fourier transform on the seismic data (t) in the time direction yields the vector S(x, ω). Using a high-precision Radon transform, the frequency domain Radon forward transform data U of the seismic data is obtained. s (p k ,ω);

[0047] The first downlink wave field frequency domain transformation result acquisition module is used for arbitrary frequency component ω and ray parameter p. k Given the amplitude difference coefficient R between virtual reflection and primary reflection. k = -1, and calculate the delay time Δτ between the virtual reflection and the first reflection. k Construct the matrix G(x) i p k The first Radon transform result of the downlink wave field in the frequency domain is obtained by solving for ω).

[0048] The prediction parameter acquisition module is used to obtain the ray parameters p by iteratively inverting the objective function.k The corresponding virtual reflection and the average delay time of a single reflection and the coefficient of difference in average amplitude

[0049] The downlink wave field acquisition module is used to reuse the data obtained in step (3). and The downlink wave field is obtained using step (2).

[0050] The module for obtaining the frequency domain transformation result corresponding to the virtual reflection is used to recalculate the average delay time between the virtual reflection and the first reflection. and the coefficient of difference in average amplitude By iterating repeatedly until the convergence condition is met, the final Radon transform result U in the frequency domain corresponding to the virtual reflection is obtained. Down (p, ω);

[0051] The final virtual reflection prediction data acquisition module is used to obtain U Down Performing an inverse Radon transform in the frequency domain and an inverse Fourier transform in the time domain on (p, ω) yields the final virtual reflection prediction data u. Down (x, t).

[0052] Another objective of the present invention is to provide an information data processing terminal, the information data processing terminal including a memory and a processor, the memory storing a computer program, and when the computer program is executed by the processor, causing the processor to perform the two-dimensional sea surface ghost wave water body imaging measurement method.

[0053] Another object of the present invention is to provide a computer-readable storage medium storing a computer program, which, when executed by a processor, causes the processor to perform the two-dimensional sea surface ghost wave water body imaging measurement method.

[0054] Another object of the present invention is to provide a flow measurement device ADCP, which performs the two-dimensional sea surface ghost wave water body imaging measurement method.

[0055] Another object of the present invention is to provide a computer program product stored on a computer-readable medium, including a computer-readable program that, when executed on an electronic device, provides a user input interface for implementing the two-dimensional sea surface ghost wave water body imaging measurement method.

[0056] Combining all the above technical solutions, the advantages and positive effects of this invention are as follows: It innovatively utilizes virtual reflection information to enrich the effective signal components of exploration, which can not only corroborate the results of previous signal processing, but also improve the signal-to-noise ratio and provide rich signal information.

[0057] Virtual reflection signals are rich in water body information. This invention can perform full waveform inversion on useful information such as virtual reflection signals received by ADCP to obtain the property information of surface seawater and achieve accurate imaging of water bodies.

[0058] It should be understood that the above general description and the following detailed description are exemplary and explanatory only, and are not intended to limit the disclosure of the present invention. Attached Figure Description

[0059] The accompanying drawings, which are incorporated in and form a part of this specification, illustrate embodiments consistent with this disclosure and, together with the description, serve to explain the principles of this disclosure.

[0060] Figure 1 This is a flowchart of the two-dimensional sea surface ghost wave water body imaging measurement method provided in the embodiments of the present invention.

[0061] Figure 2 This is a JANUS structure diagram of ADCP transmitting acoustic waves provided in an embodiment of the present invention.

[0062] Figure 3 This is a schematic diagram of an ADCP with a JANUS structure provided in an embodiment of the present invention. Detailed Implementation

[0063] To make the above-mentioned objects, features, and advantages of the present invention more apparent and understandable, specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Many specific details are set forth in the following description to provide a thorough understanding of the present invention. However, the present invention can be practiced in many other ways different from those described herein, and those skilled in the art can make similar modifications without departing from the spirit of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0064] like Figure 1 As shown, the present invention provides a two-dimensional sea surface ghost wave water body imaging measurement method, including:

[0065] S101, for earthquake records S(x) i Performing a Fourier transform on the seismic data (t) in the time direction yields the vector S(x, ω). Using a high-precision Radon transform, the frequency domain Radon forward transform data U of the seismic data is obtained. s (p k ,ω);

[0066] S102, for any frequency component ω and ray parameter p k Given the amplitude difference coefficient R between virtual reflection and primary reflection. k = -1, and calculate the delay time Δτ between the virtual reflection and the first reflection. k Construct the matrix G(x) i pk The first Radon transform result of the downlink wave field in the frequency domain is obtained by solving for ω).

[0067] S103, by alternately iterating and inverting the objective function, the ray parameter p is obtained. k The corresponding virtual reflection and the average delay time of a single reflection and the coefficient of difference in average amplitude

[0068] S104, Repeat step S103 to obtain a new... and Then, the second down-wave field frequency domain Radon transform result is obtained using step S102.

[0069] S105, recalculate the average delay time of the virtual reflection and the primary reflection. and the coefficient of difference in average amplitude By iterating repeatedly until the convergence condition is met, the final Radon transform result U in the frequency domain corresponding to the virtual reflection is obtained. Down (p, ω);

[0070] S106, the Radon transform result U in the frequency domain corresponding to the virtual reflection. Down Performing an inverse Radon transform in the frequency domain and an inverse Fourier transform in the time domain on (p, ω) yields the final virtual reflection prediction data u. Down (x, t).

[0071] Another object of the present invention is to provide a two-dimensional sea surface ghost wave imaging measurement system for implementing the aforementioned two-dimensional sea surface ghost wave water body imaging measurement method, the two-dimensional sea surface ghost wave water body imaging measurement system comprising:

[0072] The earthquake data frequency domain forward transform data acquisition module is used to process earthquake records S(x) i Performing a Fourier transform on the seismic data (t) in the time direction yields the vector S(x, ω). Using a high-precision Radon transform, the frequency domain Radon forward transform data U of the seismic data is obtained. s (p k ,ω);

[0073] The first downlink wave field frequency domain transformation result acquisition module is used for arbitrary frequency component ω and ray parameter p. k Given the amplitude difference coefficient R between virtual reflection and primary reflection. k = -1, and calculate the delay time Δτ between the virtual reflection and the first reflection. k Construct the matrix G(x) i p kThe first Radon transform result of the downlink wave field in the frequency domain is obtained by solving for ω).

[0074] The prediction parameter acquisition module is used to obtain the ray parameters p by iteratively inverting the objective function. k The corresponding virtual reflection and the average delay time of a single reflection and the coefficient of difference in average amplitude

[0075] The downlink wave field acquisition module is used to reuse the data obtained in step (3). and The downlink wave field is obtained using step (2).

[0076] The module for obtaining the frequency domain transformation result corresponding to the virtual reflection is used to recalculate the average delay time between the virtual reflection and the first reflection. and the coefficient of difference in average amplitude By iterating repeatedly until the convergence condition is met, the final Radon transform result U in the frequency domain corresponding to the virtual reflection is obtained. Down (p, ω);

[0077] The final virtual reflection prediction data acquisition module is used to obtain U Down Performing an inverse Radon transform in the frequency domain and an inverse Fourier transform in the time domain on (p, ω) yields the final virtual reflection prediction data u. Down (x, t).

[0078] The technical solution of the present invention will be further described below with reference to specific embodiments.

[0079] Since this invention proposes precise water imaging based on full waveform inversion of useful information received by an ADCP device, the relevant content is described below to enable those skilled in the art to have a clearer understanding of the invention.

[0080] I. The principle of ghost wave generation

[0081] To increase excitation energy and reduce seawater swell noise, marine seismic exploration often places the seismic source and detector in the seawater. However, seawater impedance is approximately 1, while air impedance is 0, meaning the former is greater than the latter, resulting in a reflection interface with a reflection coefficient of approximately -1. When the seismic source excites the seismic source below sea level, it propagates in both upward and downward directions. When it propagates upward to sea level, a polarity-reversed seismic reflection occurs, known as a virtual reflection or ghost wave. Ghost waves generally have three propagation paths, which can be defined by their generation mechanism as: source virtual reflection, detector virtual reflection, and both source and detector virtual reflection. Among them, source virtual reflection refers to the energy propagating upward after the source is excited being reflected by the seawater-air impedance interface, and then reflected by the seabed to the detector; detector virtual reflection refers to part of the source energy propagating downward, being reflected by the underground interface, and then reflected by the seawater-air wave impedance interface to the detector; source and detector virtual reflection refers to the energy of the source being reflected by the seawater-air impedance interface, propagating downward to the seabed, being reflected again to the seawater-air wave impedance interface, and finally being received by the detector.

[0082] II. ADCP Flow Measurement Principle

[0083] ADCP—Acoustic Doppler Current Profiler—is currently the most effective method for observing upper ocean current profiles. Its basic principle is to measure the Doppler frequency shift of sound waves backscattered by microparticles in seawater, thereby obtaining the velocity of water in different water layers. Its advantages include high accuracy, high resolution, and ease of operation.

[0084] ADCP measures water velocity by measuring the Doppler frequency shift of acoustic backscattering from microparticles in the water. In the instrument coordinate system, the velocity component of the water at each depth along each acoustic beam direction can be calculated using the following equation based on the Doppler frequency (D) measured at that depth:

[0085]

[0086] In the formula F t Where is the transmission frequency, in kHz. C is the speed of sound in water, which can be calculated using the following formula:

[0087] C = 1449.2 + 4.6T - 5.5 * 10 -3 T 3 +2.9*10 -4 T 3 +(1.34-10 -2 T)

[0088] In the formula:

[0089] D = Depth, in meters (m);

[0090] S = Salinity, in units of 10. -3 ;

[0091] T = Temperature, in °C.

[0092] By using the Doppler frequency shift of the water backscattered from four sound beams (at least three), the three-dimensional flow velocities, Vx, Vy, and Vz, can be determined. These can then be converted to Earth coordinates as U (east component), V (north component), and W (vertical component). Since the speed of sound remains essentially constant within a certain depth range (from the surface to several hundred meters), the depth can be determined based on the time difference between sound wave emission and reception. By using continuously emitted sound pulses, determining specific emission time intervals and hysteresis, and estimating the Doppler frequency shift and spectral width, the flow velocity of the water in each layer of the entire water profile can be obtained.

[0093] Generally, ocean current velocity can be considered constant within a certain range. Therefore, ocean current velocity can be defined as a three-dimensional velocity vector in the east-west, north-south, and vertical directions. To measure the velocity vector in these three directions, ADCP typically emits four pulse signals simultaneously in the ocean in the east, west, south, and north directions, each forming an angle φ with the vertical. Figure 2 As shown, this beam structure is called the JANUS structure. Its operational schematic is shown below. Figure 3 As shown.

[0094] For marine earthquake processing, the common approach is to eliminate virtual reflections and then use conventional migration imaging methods to migrate and image the primary reflected wave. However, virtual reflections also contain the seismic response of the water layer. To image them separately, a mirror cable method can be used, in which the detector is moved to a mirror position above the sea level relative to the real detector. In this case, the velocity on the sea surface is equivalent to the velocity of the seawater. Thus, the virtual reflection can be regarded as a new type of primary reflection, and it can be migrated and imaged using conventional migration imaging methods.

[0095] III. ADCP Waveform Determination

[0096] 1. Algorithm Principle

[0097] The traditional FWI algorithm is based on finding the minimum value of the objective function formula (1) constructed using the residual between the simulated data and the measured data. The simulated synthetic record is usually obtained by forward modeling the wave equation using the current velocity model.

[0098]

[0099] The process of finding the extremum is accomplished by constructing a gradient function formula (2) from the objective function. Currently, commonly used methods include the conjugate gradient method, the steepest descent method, and the Gauss-Newton method. In contrast to global algorithms, the biggest difference between various local extremum algorithms lies in the computer time cost.

[0100]

[0101] 2. Implementation process

[0102] 2.1 Wavelet Calculation

[0103] As can be seen from the basic principle of ADCP waveform acquisition, the residual between the simulated synthetic record and the measured record drives the inversion algorithm. Therefore, obtaining an accurate wavelet for forward modeling is the foundation and a crucial step in ADCP waveform acquisition technology. Common wavelet acquisition algorithms include: simulating the wavelet using air gun array parameters, extracting the wavelet from the direct wave of ADCP data, and using far-field wavelets synthesized from near-field records (NFH). This method extracts the wavelet from the information received by the ADCP device.

[0104] 2.2 ADCP Data Preparation

[0105] ADCP data preparation refers to the basic processing performed on the data used for ADCP waveform acquisition; it can also be called ADCP data preprocessing. A potential problem with ADCP waveform acquisition techniques is period skipping, which means the optimization process gets stuck in local extrema.

[0106] Effective low-frequency signals are crucial for the success of inversion. Therefore, obtaining ADCP data rich in low-frequency components is fundamental. For the acquired data, a low-cutoff filter is first applied to eliminate low-frequency DC components and strong-energy surge noise. Multiples are also a major interference signal in general oceanographic data; whether to suppress data multiples is a frequently discussed aspect, requiring comprehensive consideration in conjunction with other parameters, most importantly, the selection of the source wavelet. Therefore, in the data preparation process, only some essential basic work needs to be done, significantly reducing the workload compared to conventional production processing. This is why some in the industry suggest that inversion can shorten the production cycle.

[0107] 2.3 Initial Model Establishment

[0108] ADCP waveform acquisition strictly requires a smooth but relatively accurate velocity model. A good initial model is crucial for updating ADCP waveforms. If the simulated data obtained from forward modeling using the initial model has too large an error compared to the measured data, the inversion process can easily get trapped in local extrema. Therefore, conventional velocity modeling methods are needed to appropriately optimize the velocity model before using it as the initial model for inversion, and a careful check and analysis should be performed before using it in the ADCP waveform acquisition. The initial model can be obtained through well data, time-domain processing, and travel-time inversion methods.

[0109] The method of this invention uses a time-domain ADCP waveform acquisition algorithm and employs a strategy of gradually increasing the inversion frequency in frequency bands. It begins the inversion update with the lowest possible effective signal energy to ensure stable convergence of the inversion results, and then gradually increases the inversion frequency. Therefore, the starting frequency of the inversion is an important parameter.

[0110] Other key parameters include the maximum inversion frequency. The ADCP waveform is obtained through an iterative optimization process, and each iteration contains two wave equation algorithms.

[0111] 2.4 Virtual Reflection Extraction Method

[0112] The time difference and waveform characteristics between virtual reflections and primary reflections differ only slightly in the time domain profile, making direct calculation and comparison difficult to distinguish, thus hindering the prediction of virtual reflection seismic records. Therefore, filtering methods for multiple wave elimination often transform the time domain profile to the FK or Radon domain for processing, obtaining better-quality characteristics, and using the time difference relationship between primary waves and virtual reflections to separate them. Simultaneously, the estimation of the time difference between primary and virtual reflections and the virtual reflection reflection coefficient also affects the accuracy of virtual reflection prediction, adapting to the complexity brought about by seafloor tectonic undulations and velocity variations in virtual reflection prediction and elimination.

[0113] 2.4.1 High-precision Radon transform based on the Majorization-Minimization framework

[0114] The Radon transform is a commonly used method for multiple wave elimination in seismic processing. Any equation or image can be reconstructed using its own projection. This theory was first applied to CT scans, calculating the integral projection value of an object at a certain angle by setting the X-ray path, and then performing an inverse transform on these values ​​to recover the original image. Assuming the position of the ray is r = xsin θ + ycos θ, the total attenuation of the ray projection at that angle is:

[0115]

[0116] Where f(x, y) is the two-dimensional data to be reconstructed.

[0117] 2.4.2 Analysis of Radon Transform and its Key Issues

[0118] The Radon transform in the time domain is a linear integral over a seismic profile following a specific ray trajectory, and its expression is:

[0119]

[0120] Where m(τ, p) represents the transformed data of time-domain data d(t, x) after Radon forward transform, t = φ(τ, p, x) is the integration path, and t, τ, p, x represent time, offset, two-way time intercept, and curvature parameter, respectively. The transformation formula for recovering d(t, x) from Radon domain data m(τ, p) is:

[0121]

[0122] Discretize the above formula along the offset x i (i = 1, ..., N) x To sum the earthquake amplitudes, the specific time-domain formula for the Radon transform is:

[0123]

[0124] Let the number of ray parameters p be N. M The corresponding inverse transform form is:

[0125]

[0126] Representing the Radon forward transform as a matrix, we can obtain...

[0127] m = L T d (3-6)

[0128] To find m, we typically first define the inverse transform operator L.

[0129] d=Lm (3-7)

[0130] Constructing the inversion objective function

[0131]

[0132] To ensure that m can recover the original data, the least squares solution of m can be obtained using the principle of optimization.

[0133] m=(LL T ) -1 L T d (3-9)

[0134] It can be seen that the dimension of the two-dimensional seismic data d in the time-domain transformation is N. x ×N t The dimension of m depends on the discrete ray parameters, so the size of m is N. M ×N t Therefore, L T The size should be (N) M ×N t )×(N x ×N tAccording to the matrix expression of the formula, the LL that satisfies the matrix multiplication law T The size is (N) M ×N t )×(N x ×N t ).

[0135] To recover the original data, the ray parameters need to be sufficiently rich, which requires M to be large enough so that LL... T The scale is too large, which will lead to LL T It requires a huge amount of storage space, and L T d also needs (N) M ×N t ) times (N) x ×N t The dot product and summation of 3D vectors are computationally expensive. Therefore, the Radon domain transform is often computed in the frequency domain, provided that the integral path used satisfies the time-invariant characteristic.

[0136] Converting the tx data to the frequency domain allows for the performance of a forward Fourier transform in the time domain for each channel.

[0137] Therefore, the frequency domain expression of the Radon forward transform is:

[0138]

[0139] The corresponding inverse transform form is:

[0140]

[0141] Where ω = 2πf is the angular frequency, and f is the frequency.

[0142] Its matrix representation is as follows:

[0143]

[0144] For any given frequency, L H (ω) has a magnitude of N x ×N t .

[0145] At this point, the objective function becomes

[0146]

[0147] Then, according to the least squares solution method:

[0148] M = (LL H ) -1 L H d (3-14)

[0149] The Radon forward transform can be performed at frequency division points in the frequency domain. Since this solution process is a multi-solution problem, a damping term u is usually applied to ensure the stability of the solution.

[0150] M = (LL H +μI) -1 L H d (3-15)

[0151] Based on the above analysis, the calculation process of Radon's frequency domain forward transform can be obtained as follows:

[0152] Step 1: Perform a time-domain Fourier transform on the TX seismic data;

[0153] Step 2: Solve for M using the formula for each frequency;

[0154] Step 3: Apply the inverse Fourier transform in the time domain to the obtained M to obtain the Radon forward transform τ-p data.

[0155] It is clear from the above calculation process that the form of the ray path, the range of values ​​for the ray parameter p, the magnitude of the damping term, and the choice of the inversion solution algorithm all affect the accuracy and resolution of the Radon transform application, which is crucial for subsequent processing.

[0156] Radon transform has wide applications in seismic exploration, and numerous improved algorithms have emerged. In ray path selection, linear, parabolic, and hyperbolic ray paths are primarily used, specifically for practical applications such as plane wave decomposition, velocity superposition, and multiple wave elimination after normal time difference correction.

[0157] 2.5 Calculation Method of Key Parameters for Radon Domain Virtual Reflection Prediction Based on Alternating Inversion

[0158] 2.5.1 Basic Principles

[0159] Although the Radon domain processing algorithm can predict virtual reflections well, the prediction process requires converting the original seismic data into Radon domain data. However, the actual data has complex morphological characteristics. In the calculation process of Radon domain data, improper selection of ray parameters and inversion constraint parameters (damping factor, sparsity constraint weight) can lead to deviations in the transformation results, thus causing the virtual reflection prediction results to deviate from the expected effect.

[0160] Without loss of generality, we assume that variable depth cable acquisition is used, where the burial depth of the i-th detector is z. i The offset distance is x i p represents the matching ray parameter.

[0161] in, θ kV is the angle of incidence of the seismic wave as it propagates to sea level, and the angle of incidence for virtual reflection; w Let be the velocity of the seawater. It is also known that the seismic wave field received at this detector is the sum of the ascending and descending wave fields.

[0162]

[0163] Among them, u S For the total wave field, u Up For the upward wave field, u Down R represents the downlink wavefield; τ is the time the seismic wave travels at offset 0; Δτ is the time difference between the uplink and downlink wavefields; and R... k It is the ratio coefficient between the amplitude of the ghost wave and the amplitude of the first reflected wave;

[0164] Based on the geometric relationship of seismic wave propagation, the time difference Δτ between the descending wave field and the ascending wave field can be determined as:

[0165]

[0166] Furthermore, since the total received wave field can be represented by Nm horizontal slowness degrees p:

[0167]

[0168] The expressions for the total wave field and the downlink wave field corresponding to the virtual reflection in the frequency domain are:

[0169]

[0170] S(x i , where ω) is the Fourier transform result of the seismic data in the time direction. Assume:

[0171]

[0172] The least squares solution for the virtual reflection wave field is:

[0173] U Down (p, ω) = (G H (x, p, ω)G(x, p, ω)) -1 G(x, p, ω) H U s (x, ω) (3-21)

[0174] Ignoring shot-receiver offset error, the unknown parameters include the time difference Δτ between the downlink and primary wavefields, the energy ratio coefficient R between virtual and primary reflections, and the downlink wavefield U to be determined. Down (p, ω). To obtain the above three parameters, an alternating iterative inversion method is needed, that is, fix two of the parameters, invert one parameter, and then use the inverted parameter as a constraint value to invert the other parameter, and so on to obtain the optimal solution.

[0175] Before beginning the alternating iterative inversion process, it is necessary to first determine an objective function for the inversion. As discussed above, the relationship between the ascending and descending wave fields can be expressed in two forms:

[0176]

[0177] Performing a Fourier transform on it yields:

[0178]

[0179] Among them, U s (p k Let ω) be the frequency domain positive transform of the total wave field. The objective of alternating inversion should be to ensure that the downlink wave field is obtained through inversion. Using the above formula, the obtained primary wave fields are equal, i.e.

[0180]

[0181] Based on the previously analyzed Radon domain virtual reflection prediction method using the Majorization-Minimization optimization framework, when Δτ is fixed, the L1 norm can be used to obtain high-resolution Radon domain transform results. Therefore, the Radon domain virtual reflection prediction process based on alternating inversion is as follows:

[0182] (1) For the earthquake record S(x) i Performing a Fourier transform on the seismic data (t) in the time direction yields the vector S(x, ω). Using a high-precision Radon transform, the frequency domain Radon forward transform data U of the seismic data is obtained. s (p k ,ω);

[0183] (2) For any frequency component ω and ray parameter p k Given the amplitude difference coefficient R between virtual reflection and primary reflection. k = -1, and calculate the delay time Δτ between the virtual reflection and the first reflection. k Thus constructing matrix G(x) i p k The first Radon transform result of the downlink wave field in the frequency domain is obtained by solving for ω).

[0184] (3) The ray parameter p is obtained by alternating iterative inversion of the objective function. k The corresponding virtual reflection and the average delay time of a single reflection and the coefficient of difference in average amplitude

[0185] (4) Reuse step (3) to obtain and The downlink wave field is obtained using step (2).

[0186] (5) Calculate the average delay time of the virtual reflection and the first reflection again. and the coefficient of difference in average amplitude By iterating repeatedly until the convergence condition is met, the final Radon transform result U in the frequency domain corresponding to the virtual reflection is obtained. Down (p, ω);

[0187] (6) For U Down Performing an inverse Radon transform in the frequency domain and an inverse Fourier transform in the time domain on (p, ω) yields the final virtual reflection prediction data u. Down (x, t).

[0188] The energy of the virtual reflection data obtained by the above process is matched with the energy of the seismic data. That is, the virtual reflection data obtained by subtracting the predicted energy from the seismic data is the primary reflection data.

[0189] 2.6 Basic Idea of ​​Virtual Reflection Imaging

[0190] For marine seismic processing, the common approach is to eliminate virtual reflections and then use conventional migration imaging methods to migrate and image the primary reflection wave. Virtual reflections also contain seismic responses from subsurface structures. By artificially moving the detector to a mirror position above the sea level relative to the actual detector, and treating the velocity at the sea surface as equivalent to the seawater velocity, virtual reflections can be considered a new type of primary reflection, which can then be migrated and imaged using conventional migration imaging methods. Simultaneously, the virtual reflection imaging results are matched and correlated with the primary reflection imaging results to obtain a joint imaging result of virtual and primary reflections, further improving the imaging accuracy of subsurface structures.

[0191] Other embodiments of this disclosure will readily occur to those skilled in the art upon consideration of the specification and practice of the disclosure herein. This application is intended to cover any variations, uses, or adaptations of this disclosure that follow the general principles of this disclosure and include common knowledge or customary techniques in the art not disclosed herein. The specification and examples are to be considered exemplary only, and the true scope and spirit of this disclosure are indicated by the appended claims.

[0192] It should be understood that this disclosure is not limited to the precise structures described above and shown in the accompanying drawings, and various modifications and changes can be made without departing from its scope. The scope of this disclosure should be limited by the appended claims.

Claims

1. A two-dimensional sea surface ghost wave water body imaging measurement method, characterized in that, This two-dimensional sea surface ghost wave water imaging measurement method uses an ADCP (Advanced Diffusion Profiler) device to receive water-rich information, extracts the effective information from filtered virtual reflections, performs full waveform inversion to obtain the property information of surface seawater, and processes and images the water signal. Specifically, it includes the following steps: (1) Earthquake records Perform a Fourier transform in the time direction to obtain a vector. By using the Radon transform, the frequency domain Radon forward transform data of the seismic data can be obtained. ; (2) For any frequency component and radiation parameters Given the amplitude difference coefficient between virtual reflection and primary reflection The delay times of the virtual reflection and the first reflection were calculated. Construct a matrix The first Radon transform result of the downlink wave field in the frequency domain was obtained by solving the problem. ; (3) The ray parameters are obtained by iteratively inverting the objective function. The corresponding virtual reflection and the average delay time of a single reflection and the coefficient of difference in average amplitude ; (4) Obtained through step (3) and The downlink wave field is obtained using the calculation in step (2). ; (5) Calculate the average delay time of the virtual reflection and the first reflection again. and the coefficient of difference in average amplitude By iterating alternately until the convergence condition is met, the final Radon transform result in the frequency domain corresponding to the virtual reflection is obtained. ; (6) Performing an inverse Radon transform in the frequency domain and an inverse Fourier transform in the time domain yields the final virtual reflection prediction data. .

2. The two-dimensional sea surface ghost wave water body imaging measurement method according to claim 1, characterized in that, In step (3), the ray parameters are obtained. The corresponding virtual reflection and the average delay time of a single reflection and the coefficient of difference in average amplitude include: The signal is received by a detector with independent signal reception and input to the current measurement device ADCP for processing. During acquisition, the depth of the i-th detector is... The offset distance is p represents the ray parameter matched to the detector; Wherein, the k-th ray parameter ; The angle of incidence is the angle at which the seismic wave propagates to the sea level and the angle at which the wave is virtually reflected. Let be the seawater velocity; the seismic wave field received at this detector is the sum of the ascending and descending wave fields, then... ； ； in, For the total wave field, For upward wave field, For downlink wave field; This represents the time it takes for the seismic wave to travel at a distance of 0. The time difference between the downlink wave field and the first wave field. The amplitude difference coefficient between virtual reflection and primary reflection; Based on the geometric relationship of seismic wave propagation, the time difference between the downlink wavefield and the primary wavefield... for: ； The total received wave field is composed of The level of slowness p represents: Then, the expressions for the total wave field and the downlink wave field corresponding to the virtual reflection in the frequency domain are: ； Given the Fourier transform result of the seismic data in the time direction, if: ； The least squares solution for the virtual reflection wave field is: 。 3. The two-dimensional sea surface ghost wave water body imaging measurement method according to claim 2, characterized in that, Calculate the time difference between the downlink wave field and the first wave field. The ratio coefficient R of virtual reflection to primary reflection energy, and the Radon transform result in the frequency domain corresponding to virtual reflection. In this process, an alternating iterative inversion method is used to fix the time difference between the downlink wavefield and the primary wavefield. The energy ratio coefficient R between virtual reflection and primary reflection, and the Radon transform result in the frequency domain corresponding to virtual reflection. The optimal solution is obtained by repeatedly using two parameters, one parameter is inverted, and then the inverted parameter is used as a constraint value to invert the other parameter.

4. The two-dimensional sea surface ghost wave water body imaging measurement method according to claim 2, characterized in that, Before starting the alternating iterative inversion process, an objective function for the inversion needs to be determined; the relationship between the up-row and down-row waves can be expressed in the following forms: ； ； Performing a Fourier transform yields: ； ； in, For the frequency domain positive transform of the total wave field; then the target of alternating inversion obtains the downlink wave field through inversion, via , The first-order wavefield obtained by the formula is equal, as follows: 。 5. The two-dimensional sea surface ghost wave water body imaging measurement method according to claim 1, characterized in that, The two-dimensional sea surface ghost Water body imaging measurement methods also include: Virtual reflections incorporate the seismic response of underground structures. The detector is artificially moved to a mirror position above the sea level relative to the real detector, and the velocity on the sea surface is equivalent to the velocity of seawater. The virtual reflection is treated as a new primary reflection and migrated using conventional migration imaging methods. At the same time, the virtual reflection imaging result is matched and correlated with the primary reflection imaging result to obtain a joint imaging result of virtual reflection and primary reflection.

6. A two-dimensional sea surface ghost wave imaging measurement system for implementing the two-dimensional sea surface ghost wave water body imaging measurement method according to any one of claims 1-5, characterized in that, The two-dimensional sea surface ghost wave water imaging measurement system includes: The seismic data frequency domain forward transform data acquisition module is used to process seismic records. Perform a Fourier transform in the time direction to obtain a vector. By utilizing high-precision Radon transform, frequency domain Radon forward transform data of seismic data are obtained. ; The first downlink wave field frequency domain transformation result acquisition module is used for arbitrary frequency components. and radiation parameters Given the amplitude difference coefficient between virtual reflection and primary reflection The delay times of the virtual reflection and the first reflection were calculated. Construct a matrix The first Radon transform result of the downlink wave field in the frequency domain was obtained by solving the problem. ; The prediction parameter acquisition module is used to obtain ray parameters by iteratively inverting the objective function. The corresponding virtual reflection and the average delay time of a single reflection and the coefficient of difference in average amplitude ; The downlink wave field acquisition module obtains... and Calculate the downlink wave field ; The module for obtaining the frequency domain transformation result corresponding to the virtual reflection is used to recalculate the average delay time between the virtual reflection and the first reflection. and the coefficient of difference in average amplitude By iterating alternately until the convergence condition is met, the final Radon transform result in the frequency domain corresponding to the virtual reflection is obtained. ; The final virtual reflection prediction data acquisition module is used for... Performing an inverse Radon transform in the frequency domain and an inverse Fourier transform in the time domain yields the final virtual reflection prediction data. .

7. An information data processing terminal, characterized in that, The information data processing terminal includes a memory and a processor. The memory stores a computer program. When the computer program is executed by the processor, the processor performs the two-dimensional sea surface ghost wave water body imaging measurement method according to any one of claims 1-5.

8. A computer-readable storage medium storing a computer program, which, when executed by a processor, causes the processor to perform the two-dimensional sea surface ghost wave water body imaging measurement method according to any one of claims 1-5.

9. A flow measurement device ADCP, characterized in that, The flow measurement device ADCP performs the two-dimensional sea surface ghost wave water body imaging measurement method according to any one of claims 1-5.

10. A computer program product stored on a computer-readable medium, comprising a computer-readable program that, when executed on an electronic device, provides a user input interface to implement the two-dimensional sea surface ghost wave water body imaging measurement method according to any one of claims 1-5.

Images