A method for inverting seismic passive source waveforms based on geological information constraints

Abstract

This invention discloses a method for seismic passive source waveform inversion based on geological information constraints, belonging to the field of seismic passive source waveform inversion. The method includes: acquiring a passive source waveform inversion dataset; collecting geological information of the study area and obtaining a geological model based on the geological information; performing a cross product between the gradient of the normalized geological model and the gradient of the normalized velocity model to be solved, obtaining a geological information constraint vector; using the square of the geological information constraint vector as a geological constraint term; performing passive source waveform inversion based on the passive source waveform inversion dataset and the geological constraint term to solve the velocity model; performing inversion iterations until a termination condition is met to obtain the optimal velocity model, thus completing the seismic passive source waveform inversion. This invention overcomes the shortcomings of existing inversion methods in terms of insufficient depth resolution, obtains a velocity model more consistent with the geological model, and improves the interpretation accuracy of the inversion results.

Description

Technical Field

[0001] This invention belongs to the field of earthquake passive source waveform inversion, and particularly relates to an earthquake passive source waveform inversion method based on geological information constraints. Background Technology

[0002] Seismic passive source imaging technology originated in the mid-20th century, initially used for seismic wave propagation research. With advancements in array observation techniques and inversion algorithms, the introduction of environmental noise imaging technology in the early 21st century significantly expanded its application scope. This method offers a wide range of data sources, high imaging resolution, requires no artificial seismic source, is low-cost, and highly applicable. It is suitable for seismically active areas and regions where artificial seismic sources are difficult to implement, and can simultaneously study structures from the surface to deep crust. With the widespread application of dense array methods and the gradual improvement in the accuracy of inversion methods, this method can achieve detailed characterization of subsurface structures, especially excelling in shallow geology and complex tectonic zones.

[0003] Seismic passive source imaging technology uses waveform data from seismic events or environmental noise acquired by seismic instruments. This data contains rich information about the Earth's interior. Simulation and inversion algorithms are used to compare the observed waveforms with seismic wave propagation models to infer information such as the velocity or density of the subsurface medium. Inversion methods based on seismic data include travel-time analysis, full waveform inversion, reflection wave inversion, sparse inversion, geometric inversion, machine learning inversion, and source inversion, among others. Waveform inversion methods, in particular, can utilize all the information of the waveform (including amplitude and phase) to infer subsurface medium parameters, thus capturing complex geological structures more accurately, improving the resolution of imaging results, and making them more suitable for research on deep Earth exploration. In recent years, due to the continuous improvement of computing power, this method has been widely applied in the field of deep Earth interior exploration and has achieved some results. However, the method itself is computationally expensive, easily affected by initial model and data noise, and prone to getting trapped in local minima, resulting in poor convergence and limited inversion accuracy. Meanwhile, passive source imaging single-method inversion is limited by the method itself, with its resolution decreasing with increasing depth and insufficient resolution for subsurface impedance interfaces. Consequently, the inversion results obtained based on this single method have poor correspondence with geological data. To overcome these limitations, it is urgent to explore an effective method to improve the accuracy of the inversion.

[0004] Several invention patents exist for waveform inversion methods based on seismic data. One existing full-waveform inversion method, based on seismic record integration, solves the problem of waveform inversion failing due to local extrema caused by the lack of low-frequency seismic data. Another existing full-waveform inversion method addresses the issue of multiple solutions in waveform inversion algorithms. A further existing full-waveform inversion method and system provides accurate full-waveform inversion even when low-frequency information from land data is lacking. Finally, a more accurate step-by-step initial model is provided, resulting in more reliable velocity and density inversion results, among other improvements.

[0005] The aforementioned methods, by improving the accuracy of the initial model and employing distributed inversion techniques, address situations where single methods are highly dependent on the initial model and where changes to the inversion calculation can easily lead to local minima and fail to provide a globally optimal solution. For passive source waveform inversion methods, the method's resolution decreases with increasing depth, and it also has certain shortcomings in interpreting and integrating with geological data. Therefore, a method that uses known geological data to establish a model constraining waveform inversion leverages the clear spatial structure and physical interpretability of the geological model. Compared to simple velocity inversion, this method achieves multi-parameter collaborative inversion, improving the lateral and longitudinal resolution of subsurface medium parameters and properties. This makes the inversion results more consistent with geological laws and more conducive to the interpretation and analysis of geological structures in complex regions. Summary of the Invention

[0006] To address the aforementioned shortcomings in existing technologies, this invention provides a seismic passive source waveform inversion method based on geological information constraints. This method overcomes the deficiency of depth resolution in existing inversion methods, obtains a velocity model that is more consistent with the geological model, and improves the interpretation accuracy of the inversion results.

[0007] To achieve the aforementioned objectives, the technical solution adopted by this invention is: a seismic passive source waveform inversion method based on geological information constraints, comprising:

[0008] Continuous seismic waveform data from various seismic stations within the study area were collected, and the continuous seismic waveform data were preprocessed to obtain the seismic data from each seismic station.

[0009] The seismic stations in the study area were paired up to obtain several sets of seismic station pairs;

[0010] Based on the earthquake data corresponding to each seismic station pair, the cross-correlation waveform data of each seismic station pair is obtained, and the cross-correlation waveform data of each seismic station pair is filtered based on the signal-to-noise ratio threshold to obtain the measured cross-correlation waveform dataset.

[0011] Based on the measured cross-correlation waveform dataset, frequency band filtering is performed to obtain the passive source waveform inversion dataset;

[0012] Collect geological information about the study area and derive a geological model based on that information;

[0013] The gradient of the normalized geological model is cross-multiplied with the gradient of the normalized velocity model to be solved to obtain the geological information constraint vector.

[0014] The square of the geological information constraint vector is used as the geological constraint term;

[0015] Based on the passive source waveform inversion dataset and geological constraints, passive source waveform inversion is performed to solve the velocity model; the inversion is iterated until the termination condition is met to obtain the optimal velocity model, thus completing the seismic passive source waveform inversion.

[0016] The beneficial effects of this invention are as follows: This invention fully combines the advantages of passive source waveform inversion method in fine characterization of underground structure and constraint of prior geological information model. By incorporating prior geological information as a constraint into the passive source waveform imaging inversion process through normalized cross gradient constraint, a velocity model more consistent with the geological model is obtained. This effectively overcomes the shortcomings of existing inversion methods in terms of insufficient depth resolution, improves the geological interpretation accuracy of inversion results, and can provide important technical support for high-precision imaging of seismic dense arrays.

[0017] Furthermore, the obtained measured cross-correlation waveform dataset is specifically as follows:

[0018] A1. Based on the seismic data from each seismic station, obtain the seismic data sets for each pair of seismic stations, and set the segment length and signal-to-noise ratio thresholds.

[0019] A2. Obtain the current seismic station pair and divide the seismic data group of the current seismic station pair according to the segment length to obtain several seismic data segment groups; the start and end times of two seismic data segments in the seismic data segment group are the same;

[0020] A3. Perform cross-correlation calculations on each group of seismic data segments to obtain several cross-correlation waveform data;

[0021] A4. Superimpose the cross-correlation waveform data with a signal-to-noise ratio greater than the signal-to-noise ratio threshold to obtain the cross-correlation waveform data of the current seismic station pair;

[0022] A5. Return to step A2 to obtain the cross-correlation waveform data of the next seismic station pair, until the cross-correlation waveform data of all seismic station pairs are obtained;

[0023] A6. Retain the cross-correlation waveform data of seismic station pairs with a signal-to-noise ratio greater than the signal-to-noise ratio threshold to obtain the measured cross-correlation waveform dataset.

[0024] The beneficial effects of the above-mentioned further scheme are as follows: based on the continuous seismic waveform data of the seismic array, the cross-correlation waveform data between the seismic station pairs is obtained through a series of data processing, which prepares for obtaining the data required for inversion.

[0025] Furthermore, obtaining the passive source waveform inversion dataset specifically involves: selecting different filtering windows for each cross-correlation waveform data in the measured cross-correlation waveform dataset to filter different frequency bands, thereby obtaining waveform data of different frequency bands corresponding to each cross-correlation waveform data;

[0026] By integrating waveform data from different frequency bands corresponding to all cross-correlation waveform data in the measured cross-correlation waveform dataset, a passive source waveform inversion dataset is obtained.

[0027] The beneficial effect of the above-mentioned further scheme is that by performing frequency division filtering on the cross-correlation waveform data of the seismic station pairs, waveform datasets of different frequency bands of all seismic station pairs are obtained. The dataset obtained in this step is a passive source waveform inversion dataset.

[0028] Furthermore, the expression for the geological information constraint vector is:

[0029]

[0030]

[0031] in, h is the geological information constraint vector; x h represents the x-component of the geological information constraint vector; y h represents the y-component of the geological information constraint vector; z The z-component of the geological information constraint vector; × is the gradient operator; × is the cross product operator; | is the modulus calculation operator; The gradient of the geological model; m R For geological models; For the velocity model gradient; m S The velocity model to be solved; Gradient vector The modulus; The gradient of the normalized geological model; Gradient vector The modulus; The gradient of the normalized velocity model to be solved; For partial derivative operators; For m R Partial derivatives of the model along the x-axis; For m R The partial derivatives of the model along the y-axis; For mR The partial derivatives of the model along the z-axis; For m S Partial derivatives of the model along the x-axis; For m S The partial derivatives of the model along the y-axis; For m S The partial derivatives of the model along the z-axis; and All are intermediate parameters; Δx is the grid spacing along the x-axis of the rectangular coordinate system; Δy is the grid spacing along the y-axis of the rectangular coordinate system; Δz is the grid spacing along the z-axis of the rectangular coordinate system; m Sr For m S The right mesh of the current central hexahedral mesh of the model; m Sl For m S The left mesh of the current central hexahedral mesh of the model; m Su For m S The upper mesh of the current central hexahedral mesh of the model; m Sb For m S The lower mesh of the current central hexahedral mesh of the model; m Sa For m S The model's current central hexahedral mesh is the front mesh; m Sn For m S The rear mesh of the current central hexahedral mesh of the model; m Rr For m R The right mesh of the current central hexahedral mesh of the model; m Rl For m R The left mesh of the current central hexahedral mesh of the model; m Ru For m R The upper mesh of the current central hexahedral mesh of the model; m Rb For m R The lower mesh of the current central hexahedral mesh of the model; m Ra For m R The model's current central hexahedral mesh is the front mesh; m Rn For m R The mesh following the current central hexahedral mesh of the model.

[0032] The beneficial effects of the above-mentioned further scheme are as follows: the collected geological information is transformed into prior geological information constraints through the normalized cross gradient function, and the mathematical expression of the normalized cross gradient constraint function is obtained under the discrete grid, which prepares for the subsequent construction of prior geological information constraint terms.

[0033] Furthermore, the expression for the geological constraint term is:

[0034]

[0035] Where, Θ h Geological constraints; h is the geological information constraint vector; x h represents the x-component of the geological information constraint vector; y h represents the y-component of the geological information constraint vector; z is the z-component of the geological information constraint vector; T is the transpose.

[0036] The beneficial effect of the above-mentioned further scheme is that the square of the normalized cross ladder function is used to construct the geological information constraint term, which prepares for the subsequent construction of the seismic passive source constraint inversion objective function based on geological information constraints.

[0037] Furthermore, the objective function for the passive source waveform inversion is:

[0038] Θ S =Θ t +ξΘ h

[0039]

[0040] Where, Θ S The objective function for passive source waveform inversion based on geological information constraints is Θ. t ξ represents the data fitting term for passive source waveform inversion; ξ represents the weighting factor for the geological model constraint term; Θ represents... h For geological constraints; d obs N is the passive source waveform inversion dataset; f is the theoretical waveform set obtained by the forward modeling operator; c N represents the total number of filter windows used to filter the cross-correlation waveform data; c is the filter window index; N p p represents the total number of cross-correlated waveform data in the passive waveform inversion dataset; p is the index of the cross-correlated waveform data. f is the data after the p-th cross-correlation waveform data has been filtered by the c-th filter window; pc For in m S Calculated in the model The corresponding theoretical waveform data; σ is the travel time difference between the p-th cross-correlation waveform data and the corresponding theoretical waveform data; pc Let be the covariance of the p-th cross-correlation waveform data after being filtered by the c-th filter window.

[0041] The beneficial effect of the above-mentioned further scheme is that the joint construction of the passive source waveform data term and the geological information constraint term to construct the seismic passive source waveform inversion objective function based on geological information constraints is an important part of the present invention.

[0042] Furthermore, the formula for the inversion iteration is:

[0043]

[0044] δΘ S =δΘ t +ξδΘ h

[0045]

[0046] in, This is the velocity model updated after the k-th iteration; The velocity model before the k-th iteration update; α k The update step size is used to control the magnitude of model updates; H k This is the inverse of the approximate Hessian matrix used to adjust the gradient direction in the k-th iteration; g represents the gradient of the objective function in the k-th iteration. S The gradient of the objective function; Θ S Let m be the objective function for passive source waveform inversion based on geological information constraints; S For velocity model; K σ R is the sensitive kernel function for density σ; P and R σ Both are constants; m σ For a three-dimensional density model; K S K is the sensitive kernel function for S-wave velocity; P For P-wave velocity-sensitive kernel function; m P For a three-dimensional P-wave velocity model; ξ is the weighting factor of the geological model constraint terms; h x h represents the x-component of the geological information constraint vector; y h represents the y-component of the geological information constraint vector; z The component of the geological information constraint vector along the z-axis; For partial derivative operators; Θ t For passive source waveform inversion, the data fitting term is Θ. h For geological constraints; N c N represents the total number of filter windows used to filter the cross-correlation waveform data; c is the filter window index; N p p represents the total number of cross-correlated waveform data in the passive waveform inversion dataset; p is the index of the cross-correlated waveform data. f is the data after the p-th cross-correlation waveform data has been filtered by the c-th filter window; pc For in m S Calculated in the model The corresponding theoretical waveform data; σ is the travel time difference between the p-th cross-correlation waveform data and the corresponding theoretical waveform data; pc Let be the covariance of the p-th cross-correlation waveform data after being filtered by the c-th filter window.

[0047] The beneficial effects of the above-mentioned further scheme are as follows: Using the finite-memory quasi-Newton iterative method to solve the objective function of passive source waveform inversion based on geological information constraints can effectively improve computational efficiency and obtain a high-precision subsurface velocity model. By combining the advantages of the passive source waveform inversion method and the constraints of the prior geological information model, the shortcomings of traditional passive source imaging methods in terms of insufficient resolution with increasing depth can be effectively overcome, improving inversion accuracy and providing important technical support for high-precision detection of dense seismic arrays.

[0048] Furthermore, the termination condition is: reaching the maximum number of iterations, or the data fit difference being lower than the minimum threshold, or the speed model similarity between two adjacent iterations being higher than the minimum threshold.

[0049] The beneficial effects of the above-mentioned further solutions are: guiding the end of the inversion, providing an algorithm exit point, and avoiding entering an infinite loop. Attached Figure Description

[0050] Figure 1 This is a flowchart of passive source waveform inversion based on geological information constraints.

[0051] Figure 2 This is a theoretical geological model diagram.

[0052] Figure 3 This is a graph showing the results of a single-method inversion of a passive source waveform.

[0053] Figure 4 This is a diagram showing the passive source waveform inversion results based on geological model constraints.

[0054] Figure 5 This is a schematic diagram of waveform fitting after filtering the theoretical waveform and the measured cross-correlation waveform through a 10-30s filter window. Detailed Implementation

[0055] The specific embodiments of the present invention are described below to enable those skilled in the art to understand the present invention. However, it should be understood that the present invention is not limited to the scope of the specific embodiments. For those skilled in the art, various changes are obvious as long as they are within the spirit and scope of the present invention as defined and determined by the appended claims. All inventions utilizing the concept of the present invention are protected.

[0056] like Figure 1 As shown, in one embodiment of the present invention, a method for inverting seismic passive source waveforms based on geological information constraints includes:

[0057] Continuous seismic waveform data from various seismic stations within the study area were collected, and the continuous seismic waveform data were preprocessed to obtain the seismic data from each seismic station.

[0058] The seismic stations in the study area were paired up to obtain several sets of seismic station pairs;

[0059] Based on the earthquake data corresponding to each seismic station pair, the cross-correlation waveform data of each seismic station pair is obtained, and the cross-correlation waveform data of each seismic station pair is filtered based on the signal-to-noise ratio threshold to obtain the measured cross-correlation waveform dataset.

[0060] Based on the measured cross-correlation waveform dataset, frequency band filtering is performed to obtain the passive source waveform inversion dataset;

[0061] Collect geological information about the study area and derive a geological model based on that information;

[0062] The gradient of the normalized geological model is cross-multiplied with the gradient of the normalized velocity model to be solved to obtain the geological information constraint vector.

[0063] The square of the geological information constraint vector is used as the geological constraint term;

[0064] Based on the passive source waveform inversion dataset and geological constraints, passive source waveform inversion is performed to solve the velocity model; the inversion is iterated until the termination condition is met to obtain the optimal velocity model, thus completing the seismic passive source waveform inversion.

[0065] The obtained measured cross-correlation waveform dataset is specifically as follows:

[0066] A1. Based on the seismic data from each seismic station, obtain the seismic data sets for each pair of seismic stations, and set the segment length and signal-to-noise ratio thresholds.

[0067] A2. Obtain the current seismic station pair and divide the seismic data group of the current seismic station pair according to the segment length to obtain several seismic data segment groups; the start and end times of two seismic data segments in the seismic data segment group are the same;

[0068] A3. Perform cross-correlation calculations on each group of seismic data segments to obtain several cross-correlation waveform data;

[0069] A4. Superimpose the cross-correlation waveform data with a signal-to-noise ratio greater than the signal-to-noise ratio threshold to obtain the cross-correlation waveform data of the current seismic station pair;

[0070] A5. Return to step A2 to obtain the cross-correlation waveform data of the next seismic station pair, until the cross-correlation waveform data of all seismic station pairs are obtained;

[0071] A6. Retain the cross-correlation waveform data of seismic station pairs with a signal-to-noise ratio greater than the signal-to-noise ratio threshold to obtain the measured cross-correlation waveform dataset.

[0072] The process of obtaining the passive source waveform inversion dataset is as follows: for each cross-correlation waveform data in the measured cross-correlation waveform dataset, different filter windows are selected to filter different frequency bands to obtain waveform data of different frequency bands corresponding to each cross-correlation waveform data;

[0073] By integrating waveform data from different frequency bands corresponding to all cross-correlation waveform data in the measured cross-correlation waveform dataset, a passive source waveform inversion dataset is obtained.

[0074] In this embodiment, obtaining passive source waveform inversion data requires five steps:

[0075] Step 1: Preprocessing of continuous seismic data from a single seismic station. Continuous seismic waveform data are collected from each seismic station within the study area. Preprocessing is performed on the individual seismic data, including downsampling, mean removal, linear trend removal, one-bit time-domain normalization, spectral whitening, and bandpass filtering, to obtain the individual seismic data.

[0076] Step 2: Pair up the seismic stations in the study area, select two seismic stations to form a pair, and pair up all the seismic stations in the area. Hereinafter, the paired seismic stations are referred to as seismic station pairs.

[0077] Step 3: Segment the seismic data corresponding to the seismic station pairs according to the set time intervals. Select the seismic waveform data of the seismic station pair and segment the seismic data of the two seismic stations according to the set time intervals (for example, if a seismic station's observation period is 10 days and the set time interval is 1 hour, then the seismic data of one seismic station will be segmented into 10 days × 24 hours / day = 240 segments). It is particularly important to ensure that the absolute start and end times of each seismic data segment are completely consistent (for example, if station 1 in the seismic station pair segments every hour starting from 00:00, then station 2 must also segment every hour starting from 00:00). Perform seismic data segmentation processing on all seismic station pairs.

[0078] Step 4: Obtaining the cross-correlation waveform data of seismic station pairs. The data from the seismic station pairs segmented in Step 3 are cross-correlated for seismic data within the same time period. (For example, after segmenting the data into 1-hour intervals, data within the same time period are selected for cross-correlation calculation for each seismic station pair, resulting in 240 cross-correlation waveform segments.) The cross-correlation waveforms obtained from each time period are then superimposed by selecting data with high signal-to-noise ratios. The superimposed data is the cross-correlation waveform data for that seismic station pair. The same cross-correlation and superposition calculations are performed on the seismic data of all seismic station pairs to obtain the cross-correlation waveform data for all seismic station pairs. Finally, the signal-to-noise ratio (SNR) is calculated for all cross-correlation waveform data, and only the measured cross-correlation waveform data that meets the SNR requirements are retained. In the following formula, d is used as the SNR value. obs Represents the measured cross-correlation waveform data, where N is in the following formula. p The total number of measured cross-correlation waveform data to meet the signal-to-noise ratio requirements, This represents the p-th measured cross-correlation waveform data.

[0079] Step 5, frequency-division waveform acquisition (the frequency-division waveform data is the inversion data used for passive source waveform inversion). The main waveform component in the measured cross-correlation data obtained in Step 4 is a surface wave. Given the dispersion characteristics of surface waves (i.e., the propagation speed of surface waves in different frequency bands is different), different filter windows need to be selected to filter the measured cross-correlation waveform to obtain the frequency-division waveform data. In the following formula, N... c represents the number of selected filter windows, where c is the c-th filter window. The waveform data after filtering by the c-th filter window for the p-th measured cross-correlation waveform data. Figure 5 The waveforms of cross-correlation data from different seismic stations are shown after passing through a 10s to 20s filter window.

[0080] The expression for the geological information constraint vector is:

[0081]

[0082]

[0083] in, h is the geological information constraint vector; x h represents the x-component of the geological information constraint vector; y h represents the y-component of the geological information constraint vector; z The component of the geological information constraint vector along the z-axis; × is the gradient operator; × is the cross product operator; | is the modulus calculation operator; The gradient of the geological model; m R For geological models; For the velocity model gradient; m S The velocity model to be solved; Gradient vector The modulus; The gradient for a normalized geological model; Gradient vector The modulus; The gradient of the normalized velocity model to be solved; For partial derivative operators; For m R Partial derivatives of the model along the x-axis; For m R The partial derivatives of the model along the y-axis; For m R The partial derivatives of the model along the z-axis; For m S Partial derivatives of the model along the x-axis; For m S The partial derivatives of the model along the y-axis; For m S The partial derivatives of the model along the z-axis; and All are intermediate parameters; Δx is the grid spacing along the x-axis of the rectangular coordinate system; Δy is the grid spacing along the y-axis of the rectangular coordinate system; Δz is the grid spacing along the z-axis of the rectangular coordinate system; m Sr For m S The right mesh of the current central hexahedral mesh of the model; m Sl For m S The left mesh of the current central hexahedral mesh of the model; m Su For m S The upper mesh of the current central hexahedral mesh of the model; m Sb For m S The lower mesh of the current central hexahedral mesh of the model; m Sa For m S The model's current central hexahedral mesh is the front mesh; m Sn For m S The rear mesh of the current central hexahedral mesh of the model; m Rr For m R The right mesh of the current central hexahedral mesh of the model; m Rl For m R The left mesh of the current central hexahedral mesh of the model; m Ru For m R The upper mesh of the current central hexahedral mesh of the model; m Rb For m R The lower mesh of the current central hexahedral mesh of the model; m Ra For mR The model's current central hexahedral mesh is the front mesh; m Rn For m R The mesh following the current central hexahedral mesh of the model.

[0084] The expression for the geological constraint term is:

[0085]

[0086] Where, Θ h Geological constraints; h is the geological information constraint vector; x h represents the x-component of the geological information constraint vector; y h represents the y-component of the geological information constraint vector; z is the z-component of the geological information constraint vector; T is the transpose.

[0087] In this embodiment, the geological information includes fault spatial distribution, subsurface rock properties, basin structure, etc. The collected information is used to model, for example... Figure 2 A theoretical geological model is presented, and the mathematical form of this geological model is denoted as m. R Let m be the model to be solved. S (m S This is also the unknown quantity to be solved in the passive source waveform inversion of this invention patent). This invention patent constructs m R and m S The mathematical relationship is based on the "normalized geological model m". R The gradient and the normalized velocity model m to be solved S The matrix obtained by cross-product calculation of the gradient is used as the geological information constraint vector, which is the geological model m. R With the velocity model to be solved m S The structural coupling relationship is represented by vectors. This represents the structural coupling relationship. In three-dimensional space, a vector... To write along a three-component rectangular coordinate system (x, y, z):

[0088]

[0089] and Specifically, it can be expressed as follows:

[0090]

[0091] set up The three components (h) x h y h z ),,h x for The component along the x-axis, h y for The component along the y-axis, h z for The components along the z-axis, the mathematical form of these three quantities can be rewritten as:

[0092]

[0093] In a 3D discretized mesh, Δx, Δy, and Δz represent the mesh spacing along the Cartesian coordinate system (x, y, z axes), respectively. Let the character 'o' specifically refer to the current central hexahedral mesh. The mesh to the right of the central mesh 'o' is denoted as 'r', the mesh to the left of the central mesh 'o' as 'l', the mesh above the central mesh 'o' as 'u', the mesh below the central mesh 'o' as 'b', the mesh in front of the central mesh 'o' as 'a', and the mesh after the central mesh 'o' as 'n'. Then, in m... S In the model, the current central hexahedral mesh is m So The right, left, top, bottom, front, and back grids of the current central hexahedral grid are denoted as m. Sr m Sl m Su m Sb m Sa m Sn The discretized formula can then be rewritten as:

[0094]

[0095]

[0096] Similarly, in m R In the model, let the current central hexahedral mesh be... The right, left, top, bottom, front, and back grids of the current central hexahedral grid are denoted as m. Rr m Rl m Ru m Rb m Ra m Rn The discretized formula can then be rewritten as:

[0097]

[0098] Pick The square of is used as the geological information constraint term of this invention (denoted as Θ). h Its mathematical form can be written as

[0099]

[0100] T is the transpose operator. In this invention, the geological information model m... R It is a known quantity, m S It is the unknown quantity to be solved. Then, under the discretized grid, Θ h There is only m in the middle S Unknown quantity.

[0101] The objective function for passive source waveform inversion is:

[0102] Θ S =Θ t +ξΘ h

[0103]

[0104] Where, Θ S The objective function for passive source waveform inversion based on geological information constraints is Θ. t ξ represents the data fitting term for passive source waveform inversion; ξ represents the weighting factor for the geological model constraint term; Θ represents... h For geological constraints; d obs N is the passive source waveform inversion dataset; f is the theoretical waveform set obtained by the forward modeling operator; c N represents the total number of filter windows used to filter the cross-correlation waveform data; c is the filter window index; N p p represents the total number of cross-correlated waveform data in the passive waveform inversion dataset; p is the index of the cross-correlated waveform data. f is the data after the p-th cross-correlation waveform data has been filtered by the c-th filter window; pc For in m S Calculated in the model The corresponding theoretical waveform data; σ is the travel time difference between the p-th cross-correlation waveform data and the corresponding theoretical waveform data; pc Let be the covariance of the p-th cross-correlation waveform data after being filtered by the c-th filter window.

[0105] The formula for the inversion iteration is:

[0106]

[0107] δΘ S =δΘ t +ξδΘ h

[0108]

[0109] in, This is the velocity model updated after the k-th iteration; The velocity model before the k-th iteration update; α kThe update step size is used to control the magnitude of model updates; H k This is the inverse of the approximate Hessian matrix used to adjust the gradient direction in the k-th iteration; g represents the gradient of the objective function in the k-th iteration. S The gradient of the objective function; Θ S Let m be the objective function for passive source waveform inversion based on geological information constraints; S For velocity model; K σ R is the sensitive kernel function for density σ; P and R σ Both are constants; m σ For a three-dimensional density model; K S K is the sensitive kernel function for S-wave velocity; P For P-wave velocity-sensitive kernel function; m P For a three-dimensional P-wave velocity model; ξ is the weighting factor of the geological model constraint terms; h x h represents the x-component of the geological information constraint vector; y h represents the y-component of the geological information constraint vector; z The z-component of the geological information constraint vector; For partial derivative operators; Θ t For passive source waveform inversion, the data fitting term is Θ. h For geological constraints; N c N represents the total number of filter windows used to filter the cross-correlation waveform data; c is the filter window index; N p p represents the total number of cross-correlated waveform data in the passive waveform inversion dataset; p is the index of the cross-correlated waveform data. f is the data after the p-th cross-correlation waveform data has been filtered by the c-th filter window; pc For in m S Calculated in the model The corresponding theoretical waveform data; σ is the travel time difference between the p-th cross-correlation waveform data and the corresponding theoretical waveform data; pc Let be the covariance of the p-th cross-correlation waveform data after being filtered by the c-th filter window.

[0110] The termination conditions are: reaching the maximum number of iterations, the data fit difference being lower than the minimum threshold, or the speed model similarity between two adjacent iterations being higher than the minimum threshold.

[0111] In this embodiment, the passive source waveform inversion objective function based on geological information constraints consists of two parts Θ. t and Θ h Θ t Let be the data fitting term for the inversion of passive source waveform data. Then, the expression for the objective function of passive source waveform inversion based on geological information constraints is:

[0112] Θ S =Θ t +ξΘ h

[0113] Where, Θ S The objective function for passive source waveform inversion based on geological information constraints is Θ. t ξ represents the data fitting term for passive source waveform inversion. ξ is the weighting factor for the geological model constraint term; Θ h These are the constraints for the geological model obtained above.

[0114] Where Θ t The mathematical form is:

[0115]

[0116] Θ t The physical meaning of the data fitting term is the sum of the travel time differences between the measured cross-correlation waveform data and the theoretical cross-correlation waveform data involved in the inversion calculation. 2 For the sum of squares operator; m S The unknown quantity to be solved in the passive source waveform inversion is the S-wave velocity model. In the formula, d... obs N is the collective term for all cross-correlated waveform data used in the inversion, f is the collective term for the measured theoretical waveforms corresponding to the forward modeling operator calculations, and N is the collective term for the cross-correlated waveforms used in the inversion. p Let p be the total number of cross-correlation waveform data used in the inversion, and let p be the p-th cross-correlation waveform data among the cross-correlation waveform data used in the inversion. For the p-th cross-correlation waveform data, N c Let c be the total number of filter windows used to filter each cross-correlation waveform data, and c be the c-th filter window selected. σ is the data after the p-th cross-correlation waveform data has been filtered by the c-th filter window. pc Let f be the covariance of the p-th cross-correlated waveform data after being filtered by the c-th filter window. Let f be the theoretical waveform forward modeling operator. pc For in m S The model calculates the theoretical cross-correlation waveform data of the p-th theoretical waveform data and then filters it through the c-th filter window. The time difference is the measured cross-correlation waveform data of the p-th waveform and the theoretical waveform data after being filtered by the c-th filter window. For N p Each cross-correlation waveform data and N c The sum of the travel time differences of each frequency band. Figure 5 The image shows waveform fitting diagrams of cross-correlation waveform data and theoretical cross-correlation waveform data from different seismic stations after being filtered by a 10s to 20s filter window.

[0117] Since the inversion objective function is a strongly nonlinear function, it cannot be solved directly. An iterative method is an effective approach. This invention employs a finite-memory quasi-Newton method for solving the problem. The most crucial step in the iterative solution process is calculating the gradient of the objective function.

[0118] The differential of the objective function can be written as:

[0119] δΘ S =δΘ t +ξδΘ h

[0120] δ is the differential operator, and the difficulty in solving this formula lies in δΘ. t and δΘ h In formula Θ t In the equation f() represents the density of the medium model required for the theoretical waveform forward modeling (density is represented by σ, m σ (for three-dimensional density model), S-wave velocity (m) S (for the three-dimensional S-wave velocity model) and P-wave velocity (m) P (For a three-dimensional P-wave velocity model). δΘ t It can be written as:

[0121] δΘ t =K σ δln(m σ )+K S δln(m S )+K P δln(m P )

[0122] Where δΘ t For the objective function Θ t Differential of data terms in the middle; m σ For a three-dimensional density model, δm σ For the infinitesimal element of the density model; m S For a three-dimensional S-wave velocity model, δm S For the infinitesimal element of the S-wave velocity model; m P For a three-dimensional P-wave velocity model, δm P Let ln(m) be the infinitesimal element of the P-wave velocity model; ln() is the logarithmic operator, ln(m) σ For model m σ Find the logarithm; δln(m) σ ) is ln(m σ The infinitesimal element of ); ln(m S For model m S Find the logarithm; δln(m) S ) is ln(m S The infinitesimal element of ); ln(m P For model m PFind the logarithm; δln(m) P ) is ln(m P The infinitesimal element of K; σ K is the sensitive kernel function for density σ; S K is the sensitive kernel function for S-wave velocity; P This is the P-wave velocity-sensitive kernel function.

[0123]

[0124] m is usually constructed using empirical formulas. σ m S and m P The relation, through empirical formulas, can be expressed as f() containing three parameters (m S m σ m P Reduction has only the variable m S Introducing a constant R P and R σ To construct m σ m S and m P Relationship, R P and R σ The value needs to be determined based on prior information from different regions. Let m P =R P m S and m σ =R σ m P =R σ R P m S At the same time, the formula and

[0125] Then δΘ t The formula can be written as

[0126]

[0127] δΘ h If Θ is the gradient of the constraint terms in the geological model, then under the discretized grid, Θ h There is only m in the middle S The mathematical form of an unknown quantity can be written as:

[0128]

[0129] but:

[0130]

[0131] The gradient of the objective function is denoted as g. S , but:

[0132]

[0133] The iterative formula, using a finite-memory quasi-Newton method, is simplified as follows:

[0134]

[0135] For the initial model (where the known quantities are used, an initial known S-wave velocity model needs to be provided based on the existing velocity models in the region), According to The gradient of the objective function is calculated, where α0 is the initial given iteration step size, and H0 is the gradient of the objective function. and The calculated Hessian matrix, This is the result of the first iteration. This is for updating the model in the k-th iteration, where the value of k ranges from k = 0, 1, 2, 3, ..., N. max N max This is the maximum number of iterations set. To obtain the k-th model based on the iterative formula A new model obtained through calculation. α k The update step size is used to control the magnitude of model updates; H k It is the inverse of the approximate Hessian matrix used to adjust the gradient direction in the k-th iteration. To make the model Substitute into formula g S The value of the gradient of the objective function is calculated in the kth iteration. The model is continuously updated and iterated until one of the following termination conditions is met: the maximum number of iterations is reached, the data fit difference is lower than the minimum threshold, or the model similarity between two adjacent iterations is higher than the minimum threshold, thus obtaining the optimal speed model. Figure 3 and Figure 4 All of them use an iterative method to obtain m S Speed ​​results, among which Figure 3 This is the result of a single-method inversion of passive source waveform imaging. Figure 4 Display based on Figure 2 The passive source waveform inversion results based on geological model constraints. The figure shows that the inversion results based on geological model constraints have higher accuracy than the single-method inversion results of passive source waveform imaging, indicating that the new technology provided by this invention is reliable, effective, and advanced.

Claims

1. A method for seismic passive source waveform inversion based on geological information constraint, characterized in that, The application relates to a method for seismic passive source waveform inversion. The method comprises the following steps: Collecting continuous seismic waveform data of each seismic station in a research area, and preprocessing each continuous seismic waveform data to obtain seismic data of each seismic station; Pairing each seismic station in the research area two by two to obtain a plurality of seismic station pairs; According to the seismic data corresponding to each seismic station pair, obtaining cross-correlation waveform data of each seismic station pair, and screening the cross-correlation waveform data of each seismic station pair based on a signal-to-noise ratio threshold to obtain a measured cross-correlation waveform data set; According to the measured cross-correlation waveform data set, performing frequency-division filtering processing to obtain a passive source waveform inversion data set; Collecting geological information of the research area, and obtaining a geological model based on the geological information; wherein, is a geological information constraint vector; is a component of the geological information constraint vector in the x-axis; is a component of the geological information constraint vector in the y-axis; is a component of the geological information constraint vector in the z-axis; is a gradient operator; is a cross product operator; is a module length calculation operator; is a gradient of the geological model; is a geological model; is a velocity model gradient; is a velocity model to be solved; is a module value of the gradient vector ; is a normalized gradient of the geological model; is a module value of the gradient vector ; is a normalized gradient of the velocity model to be solved; is a partial derivative operator; is a partial derivative of the model along the x-axis direction; is a partial derivative of the model along the y-axis direction; is a partial derivative of the model along the z-axis direction; is a partial derivative of the model along the x-axis direction; is a partial derivative of the model along the y-axis direction; is a partial derivative of the model along the z-axis direction; is a partial derivative of the model along the x-axis direction; is a partial derivative of the model along the y-axis direction; is a partial derivative of the model along the z-axis direction; is a partial derivative of the model along the x-axis direction; is a partial derivative of the model along the y-axis direction; is a partial derivative of the model along the z-axis direction; and are intermediate variables; is a grid spacing along the x-axis direction of the rectangular coordinate system; is a grid spacing along the y-axis direction of the rectangular coordinate system; is a grid spacing along the z-axis direction of the rectangular coordinate system; is a right grid of a current central hexahedral grid of the model; is a left grid of the current central hexahedral grid of the model; is an upper grid of the current central hexahedral grid of the model; is a lower grid of the current central hexahedral grid of the model; is a front grid of the current central hexahedral grid of the model; is back grid of the model current center hex mesh; is right grid of the model current center hex mesh; is left grid of the model current center hex mesh; is up grid of the model current center hex mesh; is down grid of the model current center hex mesh; is front grid of the model current center hex mesh; is back grid of the model current center hex mesh; Cross-multiplying the gradient of the normalized geological model and the gradient of the to-be-solved velocity model to obtain a geological information constraint vector; the expression of the geological information constraint vector is: Taking the square of the geological information constraint vector as a geological constraint term; in, The objective function for passive source waveform inversion based on geological information constraints; For passive source waveform inversion data fitting term; These are the weighting factors for the constraints in the geological model; Geological constraints; This is a passive source waveform inversion dataset; This is the theoretical waveform set obtained by forward modeling. This represents the total number of filter windows used to filter the cross-correlation waveform data. For the filter window index; This represents the total number of cross-correlated waveform data in the passive waveform inversion dataset. For indexing cross-correlation waveform data; For the first The cross-correlation waveform data is processed by the first... Data filtered by one filter window; In order to be in Calculated in the model The corresponding theoretical waveform data; For the first The time difference between the cross-correlated waveform data and the corresponding theoretical waveform data; For the first The cross-correlation waveform data is processed by the first... The covariance of the data filtered by each filter window.

2. The method according to claim 1, wherein, Performing passive source waveform inversion based on the passive source waveform inversion data set and the geological constraint term to solve the velocity model; performing inversion iteration until a termination condition is met to obtain an optimal velocity model, completing the seismic passive source waveform inversion; the objective function of the passive source waveform inversion is: The measured cross-correlation waveform data set is obtained in the following steps: A1, obtaining seismic data groups of each seismic station pair according to the seismic data of each seismic station, and setting a segment length and a signal-to-noise ratio threshold; A2, obtaining a current seismic station pair, and segmenting the seismic data group of the current seismic station pair based on the segment length to obtain a plurality of seismic data segment groups; the start and end times of two seismic data segments in the seismic data segment group are the same; A3, respectively performing cross-correlation operation on each seismic data segment group to obtain a plurality of cross-correlation waveform data; A4, superimposing the cross-correlation waveform data with a signal-to-noise ratio greater than the signal-to-noise ratio threshold to obtain the cross-correlation waveform data of the current seismic station pair; A5, returning to step A2 to obtain cross-correlation waveform data of the next seismic station pair until the cross-correlation waveform data of all seismic station pairs are obtained; 3. The method according to claim 1, wherein, A6, retaining the cross-correlation waveform data of the seismic station pair with a signal-to-noise ratio greater than the signal-to-noise ratio threshold to obtain the measured cross-correlation waveform data set. The passive source waveform inversion data set is obtained in the following steps: selecting different filtering windows for filtering different frequency bands for each cross-correlation waveform data in the measured cross-correlation waveform data set to obtain different frequency band waveform data corresponding to each cross-correlation waveform data; 4. The method of claim 1, wherein, Integrating the different frequency band waveform data corresponding to all cross-correlation waveform data in the measured cross-correlation waveform data set to obtain the passive source waveform inversion data set. wherein, is a geological constraint term; is a geological information constraint vector; is a component of the geological information constraint vector in the x-axis; is a component of the geological information constraint vector in the y-axis; is a component of the geological information constraint vector in the z-axis; is a transpose.

5. The method of claim 1, wherein, The expression of the geological constraint term is: wherein, is the velocity model updated for the th iteration; is the velocity model updated for the th iteration; is the update step length to control the magnitude of the model update; is the inverse of the approximate Hessian matrix used to adjust the gradient direction for the th iteration; is the inverse of the approximate Hessian matrix used to adjust the gradient direction for the th iteration; is the value of the objective function gradient for the th iteration; is the objective function gradient; is the objective function; is the sensitivity kernel for density and are constants; is the three-dimensional density model; is the S wave velocity sensitivity kernel; is the P wave velocity sensitivity kernel; is the three-dimensional P wave velocity model; is the weight factor for the geological model constraint term; is the x-component of the geological information constraint vector; is the y-component of the geological information constraint vector; is the z-component of the geological information constraint vector; is the partial derivative operator; is the data fitting term for passive source waveform inversion; is the geological constraint term; is the total number of filter windows used to filter the cross-correlation waveform data; is the filter window index; is the total number of cross-correlation waveform data in the passive waveform inversion data set; is the cross-correlation waveform data index; is the data of the th cross-correlation waveform data filtered by the th filter window; is the corresponding theoretical waveform data calculated in the model; is the traveltime difference between the th cross-correlation waveform data and the corresponding theoretical waveform data; is the traveltime difference between the th cross-correlation waveform data and the corresponding theoretical waveform data; is the data of the covariance of the data filtered by the filter window.

6. The method of claim 1, wherein, The formula of the inversion iteration is: The termination condition is that the maximum number of iterations is reached, or the data fitting difference is lower than a minimum threshold, or the similarity of the velocity models of adjacent two times of iteration is higher than a minimum threshold.

Images