A method for broadband seismic constrained inversion imaging of seafloor sulfide ore body structure

Abstract

The application discloses a seabed sulfide ore body structure wideband seismic constrained inversion imaging method and belongs to the field of marine seismic detection. The application aims at the current seabed sulfide ore body structure detection problem, is based on wideband near-bottom vertical cable seismic detection data, uses a compressed sensing reconstruction method to overcome the data missing problem caused by the influence of the instrument itself, other carried equipment and environmental conditions such as sea waves, improves the subsequent processing accuracy, then constructs an initial constrained smooth velocity model based on rock sample physical property measurement, electrical detection, a plurality of station CTD and other multi-source data, further proposes a multi-step frequency damping multi-scale full waveform inversion velocity modeling method, and uses wave field illumination to correct the model in the iterative inversion, realizes high-resolution imaging of the seabed sulfide ore body from the outline to the details. The application proposes an optimization scheme of seabed sulfide seismic imaging from the angles of data accuracy, initial model construction and inversion method technology, and can effectively improve the accuracy of sulfide ore body structure imaging.

Description

Technical Field

[0001] This invention relates to the field of marine seismic detection, specifically a broadband seismic-constrained inversion imaging method for the structure of seafloor sulfide ore bodies. Background Technology

[0002] Hydrothermal circulation systems on the seabed contain important seabed resources, including hydrothermal sulfides rich in copper, zinc, and gold, representing a significant frontier in international marine science. Precise imaging of the structure of these seabed sulfide ore bodies is a crucial problem that urgently needs to be solved and is the most important foundation for resource assessment.

[0003] Currently, direct research on orebody structure internationally largely relies on deep drilling, which places high demands on the number and distribution of drilling operations, significantly limiting costs and operational methods. Due to the characteristics of small distribution areas, deep waters, and complex structures of submarine hydrothermal sulfide orebodies, their geophysical response profiles are only a few hundred meters or even tens of meters wide. Near-bottom geophysical methods based on deep-sea towed platforms or various submersibles (AUVs, ROVs, etc.) play a significant role in sulfide orebody exploration. However, in sulfide orebody exploration, the resolution of the methods themselves is limited. Near-bottom magnetic methods can only roughly delineate the alteration zone range, and near-bottom electrical methods can only roughly delineate the lateral range of the orebody (Zhu et al., JGR, 2020, ZL201910589392.2, ZL202010801649.9), but they cannot accurately obtain orebody depth information, nor can they accurately image the orebody structure.

[0004] Near-bottom seismic surveys are a novel method for studying the structure of sulfide ore bodies in recent years (ZL201821239127.9). By towing the seismic source and / or receiver array in near-seabed water or fixing it vertically to the seabed, and using a medium- to high-frequency seismic source, the detection resolution can be effectively improved. Existing seismic survey results for sulfide ore bodies indicate that vertical seismic surveys may be an effective way to obtain the internal structure of sulfide ore bodies. However, the vertical seismic velocity modeling method based on equivalent offset (EOM) (Asakawa and McIntyre, SEG Annual Meeting, 2015) is a Kirchhoff time migration method. Insufficient velocity modeling accuracy results in incorrect alignment of the phase axis in the ore body imaging results (producing "horizontal bar" artifacts), which restricts the interpretation of the ore body structure.

[0005] Theoretically, full wavefield inversion (FWI) is an ideal solution for velocity modeling and shows promise for near-bottom seismic imaging of seafloor sulfides. However, limitations such as the small distribution area of ​​seafloor sulfides, the complex structure of ore bodies, and unclear physical properties make direct use of FWI for inversion imaging of sulfide ore body structures impractical. Only broadband seismic imaging containing both high-frequency and low-frequency information is suitable for imaging sulfide ore body structures. Therefore, developing a broadband seismic inversion imaging method for seafloor sulfide ore body structures has significant practical value, will strongly promote the implementation of my country's major task of seafloor sulfide resource assessment, and will also be of great importance to improving seafloor exploration capabilities. Summary of the Invention

[0006] To address the shortcomings mentioned above, this invention provides a broadband seismic-constrained inversion imaging method for the structure of seafloor sulfide ore bodies.

[0007] The technical solution of the present invention is as follows:

[0008] This invention provides a broadband seismic-constrained inversion imaging method for the structure of seafloor sulfide ore bodies, comprising the following steps:

[0009] Step 1: Conduct physical property tests on sulfide ore body samples to obtain sound velocity, density, and porosity within a limited range of the sulfide ore body;

[0010] Step 2: Perform tilt state estimation, location information fusion, direct wave suppression, data denoising, amplitude recovery, and data deconvolution processing on vertical cable seismic data;

[0011] Step 3: Perform compressed sensing high-precision reconstruction processing on the vertical cable seismic data;

[0012] Step 4: Initial velocity modeling constrained by multi-source data;

[0013] Based on the acoustic velocity characteristics of the surrounding rock, substrate, and core obtained from the physical property tests in step 1, a smooth velocity model V is constructed. a The smooth velocity model includes the water portion above the seabed and the solid portion below the seabed;

[0014] Step 5: Multi-step frequency damping multi-scale full waveform inversion;

[0015] Full waveform inversion involves finding the minimum value of an objective function through optimization methods. The velocity that satisfies the minimum objective function value is the final inversion result. Leveraging the multi-directional and wide-bandwidth characteristics of near-bottom earthquakes, a step-by-step frequency-damped wavefield full waveform inversion method is used to perform frequency domain full waveform inversion from low to high frequencies, achieving a gradual reconstruction from outline to detail and establishing an accurate velocity model V. b .

[0016] As a preferred embodiment of the present invention, step 1 is as follows:

[0017] Based on the approximate distribution range of sulfide ore bodies determined by existing data, bottom sediment samples and drilling core samples were selected at different locations and depths on the sulfide ore bodies. Physical property tests were conducted under normal and high pressure environments to obtain physical property characteristics of the sulfide ore bodies within a limited range, including sound velocity, density, and porosity.

[0018] As a preferred embodiment of the present invention, step 2 is as follows:

[0019] Using ultra-short baseline data, the seismic source and underwater location of the vertical cable for the survey line operation are determined, and the tilt state of the vertical cable is estimated based on attitude instrument data to obtain the location of each receiving point; the near-bottom seismic data are processed by location information fusion, direct wave suppression, data denoising, amplitude recovery, and data deconvolution.

[0020] In a preferred embodiment of the present invention, in step 3:

[0021] Compressed sensing sparse constraint reconstruction method based on iterative thresholding is used to perform compressed reconstruction processing on vertical cable seismic data. Since data may be missing in both the source and receiver directions, compressed sensing reconstruction is required in two aspects to recover missing shot and trace data:

[0022] a) Seismic acquisition data consists of hundreds or thousands of shot data. For some missing or invalid shot data, the common receiver gather data is reconstructed first to obtain the rule common receiver gather data.

[0023] b) Convert the data reconstructed in the previous step into common gun set data, and reconstruct each common gun set data to obtain the regular common gun set data.

[0024] As a preferred embodiment of the present invention, the compressed sensing reconstruction method based on iterative threshold in step 3 specifically includes:

[0025] According to compressed sensing theory, the process of missing data can be represented as follows:

[0026] b = Rf

[0027] Where f represents complete seismic data with a regular grid, assuming a total of M traces; b represents observed seismic data with missing traces, assuming a total of N traces and MN traces missing; R is the missing sampling matrix, assuming a matrix dimension of N×M; for the seismic data, its waveform is sparsely represented using Curvelet or Seislet transform, assuming the transform matrix is ​​S, then:

[0028] b = Ax, with A = RS*

[0029] Where x is the sparse representation of data f in the S domain, with dimension m, and A = RS * For measurement matrix, superscript * Let represent the conjugate of a complex number; thus, the problem of reconstructing the data can be expressed as:

[0030]

[0031] The wavy lines above the variables represent estimates of the variables, and λ is the Lagrange multiplier.

[0032] Since the missing data is random, and R is a random sampling matrix, the sparsity requirement (kDlog2(m / N)≤N, where k is the sparsity and D is a constant) is met. The sparse constraint inversion method is then used to solve the above equation to obtain the final data reconstruction result.

[0033] Sparse constraint inversion is performed using an iterative thresholding method, with the following iterative formula:

[0034]

[0035] Where i is the iteration number. Represents threshold operation, γ i θ is the iteration step size, and θ is the threshold used.

[0036] As a preferred embodiment of the present invention, the modeling method for the solid part and the water part of the model in step 4 is as follows:

[0037] 1) The solid part of the smooth velocity model has a background velocity of the surrounding rock. The ore body range is based on the outline of the sulfide ore body delineated by other existing data (such as electrical resistivity tomography). The core velocity is used for massive sulfides. The velocity of vein-like sulfides is determined based on the relationship between the two layers of the existing chronically spreading mid-ocean ridge sulfide zone and the measured velocity of massive sulfides in the target area.

[0038] 2) In the smooth velocity model, the water body portion is constructed by establishing a water body sound velocity profile using CTD data obtained from CTD detection operations at several stations within the detection area. The method used is as follows: sound velocity data at a certain depth obtained from CTD at several stations are used as control points, and sound velocity data for the entire detection profile is established through data fitting and interpolation.

[0039] As a preferred embodiment of the present invention, the method for implementing multi-step frequency damping multi-scale full waveform inversion in step 5 is as follows:

[0040] a) Based on the general understanding of the use of seismic sources in seismic detection, the range of frequency components of seismic data is determined by performing spectral analysis on the seismic data, and the number of inversion steps to be performed is determined according to the number of frequency distribution bands of the seismic data.

[0041] b) Determine the starting frequency f of each frequency distribution band based on the frequency characteristics of the seismic data. b and termination frequency f e To ensure f b ~f e The range encompasses all effective frequency components of this frequency distribution band;

[0042] c) Perform frequency-domain wavefield damping processing on the seismic data to obtain the damped wavefield u for each frequency distribution band. α (f), the calculation method is as follows:

[0043]

[0044] Where u(f) is the frequency domain wave field, γ d (f) is the frequency damping operator, where f is the frequency and α is the damping factor related to the inversion scale, used to control the inversion scale. 0 < α < 1 indicates that the wave field suppresses high frequencies according to the exponential curve, and α = 0 indicates that the wave field is undamped.

[0045] d) Perform full waveform inversion step by step, with each step corresponding to the frequency distribution band of one seismic data, and the result of the previous inversion is used as the initial velocity model for the next inversion.

[0046] As a preferred embodiment of the present invention, in step 5, during the multi-step frequency-damped multi-scale full waveform inversion, wave field illumination is added to achieve high-resolution imaging of deep information; the total wave field illumination matrix is ​​calculated once in each iteration, and the calculation formula is as follows:

[0047]

[0048] Among them, IM i Let V be the wavefield illumination matrix for the i-th iteration. i The velocity model generated in the i-th iteration, For model V i The forward wave field of the j-th source at position r is For model V i The inverse propagation wave field at position r of the k-th receiver, N s N represents the total number of earthquake focal points. g This represents the total number of receiving points. The following model modification methods with lighting information are used in the full waveform inversion:

[0049]

[0050] Where, ΔV i Let V be the model modification obtained through the i-th iteration. i+1 This is the velocity model generated for the (i+1)th iteration. Iteration stops when the number of iterations reaches a predetermined value or a preset iteration termination condition is met, and the final velocity model V is obtained. b The preset termination condition is:

[0051]

[0052] Where e(i) is the least squares error of the objective function after i iterations, and η is a positive number not less than zero.

[0053] Compared with the prior art, the beneficial effects of the present invention include at least the following:

[0054] (1) In the inversion, the present invention makes full use of mixed source data containing low-frequency and high-frequency information. This wide-band data is very beneficial to improving the reliability and accuracy of the inversion.

[0055] (2) This invention obtains seabed physical property parameter information by measuring the physical properties of rock samples, and makes full use of existing multi-source data such as electrical detection, core data, and CTD data from several stations to construct a more reliable initial velocity model for full waveform inversion, which greatly improves the effectiveness and accuracy of velocity inversion;

[0056] (3) Based on the analysis of the frequency components and structure of the data, a multi-step frequency damping multi-scale full wavefield inversion method was constructed, and wavefield illumination was used to correct the model during the iteration, so as to realize high-resolution imaging of the seabed sulfide ore body from outline to detail. Attached Figure Description

[0057] Figure 1 This is a schematic diagram of a near-bottom earthquake observation method, which uses a low-frequency sea surface source and a high-frequency deep-sea source to excite the earthquake and acquire broadband data through a deep-sea vertical cable.

[0058] Figure 2 Through Figure 1 A diagram showing broadband vertical cable seismic data acquired by the observation system shown.

[0059] Figure 3 This is a schematic diagram of the overall process for broadband vertical cable seismic data processing according to the present invention;

[0060] Figure 4 This is a diagram showing the high-progress reconstruction results of compressed sensing data according to the present invention. The missing gather data has been reconstructed with relatively accurate accuracy.

[0061] Figure 5This is a diagram illustrating the effect of the frequency domain damping factor used in the multi-step frequency damping multi-scale waveform inversion of this invention.

[0062] Figure 6 This is a display of the broadband seismic-constrained inversion imaging results of deep-sea sulfide ore bodies based on this invention; Figure (a) shows the actual model, and Figure (b) shows the initial velocity model V established based on multi-source data. a Figure (c) shows the velocity model (f) established by the first step of full waveform inversion of the air gun source data. b =0Hz, f e =300Hz, α =0.080, 0.008, 0.003), Figure (d) shows the velocity model V established by the second step of deep-sea source data iterative inversion based on Figure (c). b (f b =200Hz, f e =800Hz, α=0.080, 0.008, 0.003). Detailed Implementation

[0063] The present invention will be further described and illustrated below with reference to specific embodiments. The embodiments described are merely examples of the content of this disclosure and do not limit the scope of the invention. The technical features of each embodiment in the present invention can be combined accordingly, provided that there is no mutual conflict.

[0064] like Figure 1 As shown, deep-sea seafloor sulfide exploration utilizes a broadband acoustic source composed of a surface air gun seismic source and a deep-sea towed transducer source to generate acoustic waves. The former primarily excites low-frequency seismic waves below 150Hz, while the latter primarily excites acoustic waves in the 250–800Hz range. A multi-channel vertical seismic cable, positioned close to the target body, receives the acoustic signals. This observation method offers a wide operating bandwidth, improving vertical resolution; the proximity of the excitation and reception points to the seafloor also compresses the radius of the first Fresnel zone, thereby improving lateral resolution and achieving high-resolution exploration.

[0065] like Figure 2 As shown, this embodiment uses seismic data excited by a broadband source and received by a near-seabed vertical cable for subsequent analysis and processing. This data has multi-directional and broadband characteristics, which is of great significance for achieving accurate seabed imaging.

[0066] like Figure 3 As shown in this embodiment, the broadband seismic-constrained inversion imaging method for seafloor sulfide ore body structure includes the following steps:

[0067] Step 1: Conduct physical property tests on sulfide ore body samples to obtain sound velocity, density, and porosity within a limited range of the sulfide ore body;

[0068] Based on the approximate distribution range of sulfide ore bodies determined by existing data, bottom sediment samples and drilling core samples were selected at different locations and depths on the sulfide ore bodies. Physical property tests were conducted under normal and high pressure environments to obtain physical property characteristics of the sulfide ore bodies within a limited range, including sound velocity, density, and porosity.

[0069] Step 2: Perform tilt state estimation, location information fusion, direct wave suppression, data denoising, amplitude recovery, and data deconvolution processing on vertical cable seismic data;

[0070] Ultra-short baseline data was used to determine the seismic source and underwater location of the vertical cable for the survey line operation. Based on attitude sensor data, the tilt state of the vertical cable was estimated to obtain the location of each receiving point. Near-bottom seismic data underwent location information fusion, direct wave suppression, data denoising, amplitude recovery, and data deconvolution processing.

[0071] Step 3: Perform compressed sensing high-precision reconstruction processing on the vertical cable seismic data;

[0072] In actual data acquisition, due to the influence of the instrument itself, other onboard equipment, and environmental conditions such as sea waves, some shot data from air guns and deep-sea acoustic sources are unusable due to excessive noise. Since data loss has a certain degree of randomness, using compressed sensing-based methods for high-precision data reconstruction can largely overcome the problem of data loss and meet the accuracy requirements of subsequent processing. Because data loss may occur in both the source and receiver directions, compressed sensing reconstruction is needed in both directions to recover missing shot and trace data:

[0073] a) Seismic acquisition data consists of hundreds or thousands of shot data. For some missing or invalid shot data, the common receiver gather data is reconstructed first to obtain the rule common receiver gather data.

[0074] b) Convert the data reconstructed in the previous step into common gun set data, and reconstruct each common gun set data to obtain the regular common gun set data.

[0075] like Figure 4 As shown, the top image contains vertical cable seismic data with missing traces. The missing shot and trace data can be recovered in two steps. The bottom image shows the reconstructed regular data, which can effectively improve the accuracy of subsequent inversion processing.

[0076] In step 3, the compressed sensing reconstruction of vertical cable seismic data is performed using a compressed sensing sparse constraint reconstruction method based on iterative thresholds.

[0077] According to compressed sensing theory, the process of data missing can be represented as:

[0078] b = Rf

[0079] Where f represents complete seismic data with a regular grid (assuming a total of M traces), b represents observed seismic data with missing traces (assuming a total of N traces, with MN traces missing), and R is the missing sampling matrix (assuming a matrix dimension of N×M). For seismic data, its waveform can be sparsely represented using Curvelet or Seislet transforms. Let the transform matrix be S, then:

[0080] b = Ax, with A = RS *

[0081] Where x is the sparse representation of data f in the S domain (with dimension m), and A = RS * For the measurement matrix, the superscript * denotes the conjugate of complex numbers. Thus, the data reconstruction problem can be represented as:

[0082]

[0083] The wavy line above the variable represents the estimate of the variable, and λ is the Lagrange multiplier.

[0084] Since the missing data is random, and R is a random sampling matrix, the above equation can be solved using the sparsity constraint inversion method to obtain the final data reconstruction result, provided that the sparsity requirement (kD log2(m / N)≤N, where k is the sparsity and D is a constant) is met.

[0085] Sparse constraint inversion is performed using an iterative thresholding method, with the following iterative formula:

[0086]

[0087] Where i is the iteration number. Represents threshold operation, γ i θ is the iteration step size, and θ is the threshold used.

[0088] Step 4: Initial velocity modeling constrained by multi-source data;

[0089] Based on the acoustic velocity characteristics of the surrounding rock, substrate, and core obtained from the physical property tests in step 1, and combined with the detection results from other existing methods, a smooth velocity model V is constructed. a The smooth velocity model includes both the water portion above the seabed and the solid portion below the seabed. The modeling methods for these two parts are as follows:

[0090] 1) The solid part of the smooth velocity model has a background velocity of the surrounding rock. The ore body range is based on the outline of the sulfide ore body delineated by other existing data (such as electrical resistivity tomography). The core velocity is used for massive sulfides. The velocity of vein-like sulfides is determined based on the relationship between the two layers of the existing chronically spreading mid-ocean ridge sulfide zone and the measured velocity of massive sulfides in the target area.

[0091] 2) In the smooth velocity model, the water body portion is constructed by establishing a water body sound velocity profile using CTD data obtained from CTD detection operations at several stations within the detection area. The method used is as follows: sound velocity data at a certain depth obtained from CTD at several stations are used as control points, and sound velocity data for the entire detection profile is established through data fitting and interpolation.

[0092] The initial velocity model established under multi-source data constraints is a constraint and approximation of the actual velocity model, which can effectively improve the accuracy of subsequent inversion.

[0093] Step 5: Multi-step frequency damping multi-scale full waveform inversion;

[0094] Full waveform inversion involves finding the minimum value of an objective function through optimization methods; the velocity satisfying this minimum value is the final inversion result. Leveraging the multi-directional and broadband characteristics of near-bottom earthquakes, a damped wavefield full waveform inversion method is used step-by-step to perform frequency domain full waveform inversion from low to high frequencies, achieving a gradual reconstruction from outline to detail and establishing an accurate velocity model V. b .

[0095] For the equation of sound waves:

[0096]

[0097] Where V is velocity, t is time, and r is the underground position. s Let u(r,t) be the location of the earthquake source, and let S(t)δ(rr) be the wave field at point r. s Let be the source function. Its objective function for full waveform inversion can be expressed as:

[0098]

[0099] Where, d obs This represents the inverse propagation wave field of the observation data, where u is the synthetic seismic record calculated based on the acoustic wave equation, and t, s, and g represent time, source, and receiver location, respectively. The objective function in the frequency domain can be expressed as:

[0100]

[0101] Full waveform inversion is a method of finding the objective function E through certain optimization techniques. t or E ωThe minimum value of the objective function, and the velocity that satisfies the minimum value of the objective function, is the final inversion result. For the seismic data in this embodiment, which uses a broadband sound source composed of a sea surface air gun source and a deep-sea towed transducer source, the multi-step frequency damping multi-scale full waveform inversion process is as follows:

[0102] a) Utilizing the characteristics of near-bottom earthquakes—multi-directional (multi-angle joint excitation from sea surface and deep-sea sources) and wide bandwidth (low-frequency band of airguns, mid-to-high-frequency band of deep-sea sources)—spectral analysis of the data reveals that the effective frequency band of sea surface airgun source data is approximately 3–150 Hz, while the effective frequency band of deep-sea transducer sound sources is approximately 250–800 Hz. Therefore, based on this understanding, a two-step full waveform inversion process is determined.

[0103] b) Based on the understanding from step a), determine the starting frequency f of the first frequency band. b =0Hz, termination frequency f e =300Hz, determine the starting frequency f of the second frequency band. b =200Hz, termination frequency f e =800Hz;

[0104] c) Perform frequency domain wavefield damping processing on the seismic data to obtain the damped wavefield for each frequency distribution band. The calculation method is as follows:

[0105]

[0106] Where u(f) is the frequency domain wave field, γ d (f) is the frequency damping operator, where f is the frequency, and α is a damping factor related to the inversion scale, used to control the inversion scale. 0 < α < 1 indicates that the wave field suppresses high frequencies according to an exponential curve, and α = 0 indicates that the wave field is undamped. Damping operator γ d like Figure 5 As shown.

[0107] d) Perform full waveform inversion in two steps, with each step corresponding to the frequency distribution band of one seismic data, and the result of the former inversion is used as the initial velocity model for the next inversion.

[0108] Step 1: For frequency distribution band 1, damping processing is applied to the original seismic data at α = 0.003, 0.008, and 0.080 to obtain the corresponding damped wavefield data, which is then used for inversion processing. The initial velocity for the inversion of the damped wavefield data with α = 0.003 is the result of the initial velocity modeling constrained by the multi-source data in Step 4; in the full waveform inversion processing of the three sets of data, the result of the former inversion is used as the input velocity model for the latter inversion. Damping operator γ d The curve is as follows Figure 5 As shown in a.

[0109] Step 2: For frequency distribution band 2, damping processing is applied to the original seismic data at α = 0.003, 0.008, and 0.080 to obtain the corresponding damped wavefield data, which is then used for inversion processing. The initial velocity for the inversion of the damped wavefield data with α = 0.003 is the result of the full waveform inversion in the previous step; in the full waveform inversion processing of the three sets of data, the result of the former inversion is used as the input velocity model for the latter inversion. Damping operator γ d The curve is as follows Figure 5 As shown in b.

[0110] In the two-step inversion described above, the first step mainly targets the low-frequency data of the air gun source, while the second step targets the mid-to-high-frequency data of the transducer, primarily focusing on the starting frequency f. b and termination frequency f e By adjusting α, the desired effect can be achieved through multi-step damping iterative inversion.

[0111] Preferably, in a specific embodiment of the present invention, in step 5, the multi-step frequency-damped multi-scale full waveform inversion, wavefield illumination is added during the waveform inversion process to achieve high-resolution imaging of deep information. The total wavefield illumination matrix is ​​calculated once per iteration, using the following formula:

[0112]

[0113] Among them, IM i Let V be the wavefield illumination matrix for the i-th iteration. i The velocity model generated in the i-th iteration, For model V i The forward wave field of the j-th source at position r is For model V i The inverse propagation wave field at position r of the k-th receiver, N s N represents the total number of earthquake focal points. g This represents the total number of receiving points. The following model modification methods with lighting information are used in the full waveform inversion:

[0114]

[0115] Where, ΔV i Let V be the model modification obtained through the i-th iteration. i+1 This is the velocity model generated for the (i+1)th iteration. Iteration stops when the number of iterations reaches a predetermined value or a preset iteration termination condition is met, and the final velocity model V is obtained. b The preset termination condition is:

[0116]

[0117] Where e(i) is the least squares error of the objective function after i iterations, and η is a positive number not less than zero.

[0118] like Figure 6 As shown, Figure a is the actual velocity model, and Figure b is the initial velocity model V established based on multi-source data. a Figure c shows the velocity model (f) established by the first step of full waveform inversion of the air gun source data. b =0Hz, f e =300Hz, α = 0.080, 0.008, 0.003), Figure d shows the velocity model established by the second step of deep-sea source data iterative inversion based on Figure c (f b =200Hz, f e =800Hz, α = 0.080, 0.008, 0.003), which is the final inversion result V. b It is evident that it has high resolution and accuracy.

[0119] The above-described embodiments are merely illustrative of several implementations of the present invention, and while the descriptions are specific and detailed, they should not be construed as limiting the scope of the present invention. Those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these modifications and improvements all fall within the scope of protection of the present invention.

Claims

1. A method of broadband seismic constrained inversion imaging of a seafloor sulfide ore body structure, characterized by, Includes the following steps: Step 1: Conduct physical property tests on sulfide ore body samples to obtain sound velocity, density, and porosity within a limited range of the sulfide ore body; Step 2: Perform tilt state estimation, location information fusion, direct wave suppression, data denoising, amplitude recovery, and data deconvolution processing on vertical cable seismic data; Step 3: Perform compressed sensing high-precision reconstruction processing on the vertical cable seismic data; Step 4: Model the initial velocity of the ore body using multi-source data constraints; Based on the acoustic velocity characteristics of the surrounding rock, substrate, and core obtained from the physical property tests in step 1, a smooth velocity model is constructed. The smooth velocity model includes the water portion above the seabed and the solid portion below the seabed; Step 5: Multi-step frequency damping multi-scale full waveform inversion; Full waveform inversion is a process of finding the minimum value of the objective function through optimization methods. The speed at which the minimum value of the objective function is satisfied is the final inversion result. Taking advantage of the multidirectional and broadband characteristics of near-bottom earthquakes, a damped wavefield full waveform inversion method is used stepwise to perform frequency domain full waveform inversion from low to high frequencies, achieving a gradual reconstruction from outline to detail and establishing an accurate velocity model. ; In step 5, the multi-step frequency-damped multi-scale full waveform inversion, the implementation method is as follows: a) Based on the general understanding of the use of seismic sources in seismic detection, the range of frequency components of seismic data is determined by performing spectral analysis on the seismic data, and the number of inversion steps to be performed is determined according to the number of frequency distribution bands of the seismic data. b) Determine the starting frequency of each frequency distribution band based on the frequency characteristics of the seismic data. and termination frequency To ensure ~ The range encompasses all effective frequency components of this frequency distribution band; c) performing frequency domain wave field damping processing on the seismic data to obtain a damped wave field for each frequency distribution band The calculation method is as follows: ； in, For frequency domain wave fields, For frequency damping operators, For frequency, This is a damping factor related to the inversion scale, used to control the inversion scale. This indicates that the wave field suppresses high frequencies according to an exponential curve. This indicates that the wave field is undamped; d) Perform full waveform inversion step by step, with each step corresponding to the frequency distribution band of one seismic data, and the result of the previous inversion is used as the initial velocity model for the next inversion.

2. The seabed sulfide ore body structure broadband seismic constrained inversion imaging method according to claim 1, characterized in that, Step 1 is as follows: Based on the distribution range of sulfide ore bodies determined by existing data, bottom sediment samples and drilling core samples were selected at different locations and depths on the sulfide ore bodies. Physical property tests were carried out under normal pressure and high pressure environments to obtain physical property characteristics of the sulfide ore bodies within a limited range, including sound velocity, density, and porosity.

3. The seabed sulfide ore body structure broadband seismic constrained inversion imaging method according to claim 1, characterized in that, Step 2 is as follows: Using ultra-short baseline data, the seismic source and underwater location of the vertical cable for the survey line operation are determined, and the tilt state of the vertical cable is estimated based on attitude instrument data to obtain the location of each receiving point. Near-bottom earthquake data is processed by location information fusion, direct wave suppression, data denoising, amplitude recovery, and data deconvolution.

4. The seabed sulfide ore body structure broadband seismic constrained inversion imaging method according to claim 1, characterized in that, In step 3: A high-precision compressed sensing reconstruction method based on iterative thresholds is used to reconstruct vertical cable seismic data. Since data may be missing in both the source and receiver directions, compressed sensing reconstruction is required in both directions to recover missing shot and trace data. a) Seismic acquisition data consists of hundreds or thousands of shot data. For some missing or invalid shot data, the common receiver gather data is reconstructed first to obtain the rule common receiver gather data. b) Convert the data reconstructed in the previous step into common gun set data, and reconstruct each common gun set data to obtain the regular common gun set data.

5. The seabed sulfide ore body structure broadband seismic constrained inversion imaging method according to claim 4, characterized in that, In step 3, the compressed sensing reconstruction method based on iterative thresholds specifically includes: According to compressed sensing theory, the process of missing data can be represented as follows: ； in, For complete seismic data with a regular grid, let there be a total of Dao data, For seismic observation data containing missing traces, let there be a total of Data missing Dao data, Let the missing sampling matrix have dimension 1. For seismic data, its waveforms are sparsely represented using Curvelet or Seislet transforms. Let the transform matrix be... ,but: ； where, is the data In Sparse representation of the domain, dimension , is the measurement matrix, superscript denotes the conjugate of a complex number; the reconstruction problem of the data is then expressed as: ； where the tilde over the variable represents an estimate of the variable, is a Lagrange multiplier; Because data missing is random, Given a random sampling matrix, while satisfying the sparsity requirement... In this case, the sparse constraint inversion method is used to solve the above equation to obtain the final data reconstruction result. ,in, Sparsity, It is a constant; Sparse constraint inversion is performed using an iterative thresholding method, with the following iterative formula: ； wherein, is the iteration number, denotes a threshold operation, is the iteration step size, is the threshold value used.

6. The seabed sulfide ore body structure broadband seismic constrained inversion imaging method according to claim 1, characterized in that, In step 4, the modeling methods for the solid and water parts of the model are as follows: 1) The solid part of the smooth velocity model has a background velocity of the surrounding rock velocity, the ore body range is based on the outline of the sulfide ore body delineated by existing data, the massive sulfide uses the core velocity, and the velocity of the vein-like sulfide is determined based on the relationship between the two layers of the existing chronically spreading mid-ocean ridge sulfide zone and the measured velocity of the massive sulfide in the target area. 2) In the smooth velocity model, the water body part is established by using CTD data obtained from CTD detection operations at several stations in the detection area to establish the water body sound velocity profile. The method used is as follows: using the sound velocity data at a certain depth obtained from CTD at several stations as control points, the sound velocity data of the entire detection profile is established through data fitting and interpolation.

7. The seabed sulfide ore body structure broadband seismic constrained inversion imaging method according to claim 1, characterized in that, In step 5, the multi-step frequency-damped multi-scale full waveform inversion, wavefield illumination is added during the waveform inversion process to achieve high-resolution imaging of deep information; the total wavefield illumination matrix is ​​calculated once in each iteration, and the calculation formula is as follows: ； in, For the first The wavefield illumination matrix of the next iteration, For the first The velocity model generated in the next iteration For the model The next The epicenter is located at The orthogonal wave field at that location, For the model Next Each detector point is located at The reverse propagation wave field at that location, The total number of earthquake focal points, The total number of receiving points; the following model modification methods with lighting information are used in the full waveform inversion: ； in, To pass the first The model modification amount obtained from the next iteration inversion, model For the first The velocity model is generated in the next iteration; iteration stops when the number of iterations reaches a predetermined value or a preset iteration termination condition is met, and the final velocity model is obtained. The preset termination condition is: ； wherein is the step size, the least square error of the objective function after the is a positive number not less than zero.

Images