Method for joint inversion of near-seismic body wave three-dimensional quality factor
By acquiring the displacement spectrum of seismic body waves in the frequency domain and iteratively updating the parameters, the reliability problem of three-dimensional seismic quality factor imaging is solved, achieving more efficient and accurate quality factor inversion and reducing coupling and uncertainty.
Patent Information
- Application Number
- CN202511210035.2
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-08-27
- Publication Date
- 2025-12-05
AI Technical Summary
In existing technologies, the two-step method for extracting seismic mean attenuation travel time has strong uncertainty, and the multi-event-multi-station inversion method cannot consider the consistency of multiple similar rays, resulting in low reliability of seismic three-dimensional quality factor imaging.
By obtaining the target seismic body wave displacement spectrum from the earthquake to the station in the frequency domain, and combining the initial three-dimensional quality factor, source, receiver location and average attenuation travel time formula, the theoretical spectral amplitude is calculated. The parameter perturbation is solved by the residual vector formula and first-order smoothing constraint, and the initial parameters are iteratively updated to obtain the target three-dimensional quality factor.
This avoids the problem of inconsistent average attenuation factors for the same rays, reduces the coupling between the quality factor and the seismic corner frequency, and improves the reliability and inversion accuracy of the quality factor structure and source parameters.
Smart Images
Figure CN121069474A_ABST
Abstract
Description
TECHNICAL FIELD
[0001] The present application relates to the technical field of seismic data processing, and relates to a method for jointly inverting three-dimensional quality factors. BACKGROUND
[0002] The seismic quality factor is a very important medium physical parameter, and has important indicative significance for studying the seismogenic mechanism of earthquakes and determining fault fracture zones. In addition, jointly interpreting the seismic velocity structure and the quality factor of earthquakes, i.e., the Q value structure, can help us more clearly understand the physical state of the underground medium and better understand the fault activity process.
[0003] In the related art, the two-step method needs to extract the average attenuation travel time on each ray from the seismic wave data. Since the seismic average attenuation travel time has strong coupling with the seismic corner frequency, the uncertainty of the average attenuation travel time extracted by the two-step method is strong, which reduces the reliability of the quality factor imaging. The multi-event-multi-station inversion method cannot consider the consistency constraint of the average attenuation travel times obtained by multiple similar rays, so the average attenuation travel times on the similar rays have great differences.
[0004] Therefore, how to efficiently and accurately obtain reliable seismic three-dimensional quality factors has become a problem to be solved. SUMMARY
[0005] Therefore, the embodiments of the present application provide a method for jointly inverting three-dimensional quality factors, which at least solves the problem that the related art cannot efficiently and accurately obtain reliable seismic three-dimensional quality factors.
[0006] The embodiments of the present application provide a method for jointly inverting three-dimensional quality factors, which comprises: obtaining a target seismic body wave displacement spectrum of a seismic wave to a station in a frequency domain; obtaining a first average attenuation travel time through an initial three-dimensional quality factor, a source, a receiver position and a first average attenuation travel time formula; calculating a theoretical spectrum amplitude at a frequency in a natural logarithm domain through the first average attenuation travel time, an initial corner frequency, an initial site effect, an initial low-frequency amplitude, a frequency and an amplitude formula; and converting an amplitude of the target seismic body wave displacement spectrum of the seismic wave to the station into an observed spectrum amplitude in the natural logarithm domain; obtaining a first residual vector formula based on the observed spectrum amplitude and the theoretical spectrum amplitude in the amplitude formula; the first residual vector formula comprises a parameter perturbation amount; grid the seismic wave three-dimensional propagation velocity and the initial three-dimensional quality factor in the first average attenuation travel time formula to obtain a second average attenuation travel time formula; and substitute the second average attenuation travel time formula into the first residual vector formula to obtain a second residual vector formula; add a first-order smoothing constraint and a damping parameter to the second residual vector formula to obtain a third residual vector formula; and solve the third residual formula to obtain a value of a parameter perturbation; superimpose the value of the parameter perturbation to the initial three-dimensional quality factor, the initial corner frequency, the initial low-frequency amplitude, and the initial site effect respectively to obtain a new initial three-dimensional quality factor, a new initial corner frequency, a new initial low-frequency amplitude, and a new initial site effect of a next round of iteration; substitute the new initial three-dimensional quality factor, the new initial corner frequency, the new initial low-frequency amplitude, and the new initial site effect back into the first residual vector formula to obtain a first residual vector, and when a difference between the first residual vector and an amplitude satisfies a preset condition, take the new initial three-dimensional quality factor as a target three-dimensional quality factor.
[0007] The scheme provided by the embodiment of the present application comprises the following steps: obtaining a target seismic body wave displacement spectrum of a seismic wave to a station in a frequency domain; obtaining a first average attenuation travel time through an initial three-dimensional quality factor, a seismic source, a receiver position and a first average attenuation travel time formula; calculating a theoretical spectrum amplitude at a frequency through the first average attenuation travel time, an initial corner frequency, an initial site effect, an initial low-frequency amplitude, a frequency and an amplitude formula; converting an amplitude of the target seismic body wave displacement spectrum of the seismic wave to the station into an observed spectrum amplitude in a natural logarithm domain; obtaining a first residual vector formula based on the observed spectrum amplitude and the theoretical spectrum amplitude in the amplitude formula; the first residual vector formula comprises a parameter perturbation; obtaining a second average attenuation travel time formula after griding the three-dimensional propagation velocity of the seismic wave in the first average attenuation travel time formula and the initial three-dimensional quality factor; obtaining a second residual vector formula by substituting the second average attenuation travel time formula into the first residual vector formula; obtaining a third residual vector formula by adding a first-order smoothing constraint and a damping parameter into the second residual vector formula; obtaining a value of the parameter perturbation by solving the third residual formula; superimposing the value of the parameter perturbation into the initial three-dimensional quality factor, the initial corner frequency, the initial low-frequency amplitude and the initial site effect respectively to obtain a new initial three-dimensional quality factor, a new initial corner frequency, a new initial low-frequency amplitude and a new initial site effect of a next round of iteration; substituting the new initial three-dimensional quality factor, the new initial corner frequency, the new initial low-frequency amplitude and the new initial site effect back into the first residual vector formula to obtain a first residual vector; and taking the new initial three-dimensional quality factor as a target three-dimensional quality factor when a difference between the first residual vector and the amplitude meets a preset condition. In the process, the new technology can avoid the defects that the average attenuation factors of the same ray are inconsistent by jointly inverting the seismic source parameters and the quality factor structure and avoiding the process of extracting the average factor, and the new technology can also reduce the coupling between the quality factor and the seismic corner frequency to some extent, reduce the inversion parameter amount and more reliably obtain the quality factor structure and the seismic source parameters. BRIEF DESCRIPTION OF DRAWINGS
[0008] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the following will briefly introduce the drawings needed in the embodiment description. Obviously, the drawings in the following description are only some embodiments of the present application, and for those skilled in the art, other drawings can also be obtained from these drawings without any creative effort. Figure 1 A flowchart of a joint inversion three-dimensional quality factor method provided by the embodiment of the present application; Figure 2A schematic diagram illustrating the effect of the distribution of seismic events and station geographical locations in a synthetic test provided in an embodiment of the present invention; Figure 3 This diagram illustrates a comparison of the results of the traditional two-step quality factor method and the joint inversion quality factor method for synthesis testing at different depth slices, as provided in this embodiment of the invention. Detailed Implementation
[0009] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. The following embodiments are used to illustrate the present invention, but are not intended to limit the scope of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.
[0010] In the following description, references are made to “some embodiments,” which describe a subset of all possible embodiments. However, it is understood that “some embodiments” may be the same subset or different subsets of all possible embodiments and may be combined with each other without conflict.
[0011] It should be noted that the terms "first, second, and third" used in the embodiments of the present invention are only used to distinguish similar objects and do not represent a specific ordering of objects. It is understood that "first, second, and third" can be interchanged in a specific order or sequence where permitted, so that the embodiments of the present invention described herein can be implemented in an order other than that illustrated or described herein.
[0012] It will be understood by those skilled in the art that, unless otherwise defined, all terms used herein (including technical and scientific terms) have the same meaning as commonly understood by one of ordinary skill in the art to which these embodiments of the invention pertain. It should also be understood that terms such as those defined in general dictionaries should be understood to have the same meaning as in the context of the prior art and should not be interpreted in an idealized or overly formal sense unless specifically defined as herein.
[0013] Figure 1 This is a flowchart illustrating a joint inversion method for three-dimensional quality factors provided in an embodiment of the present invention. The joint inversion method for three-dimensional quality factors provided in this embodiment of the present invention can be executed by an electronic device, such as a computer or a server.
[0014] like Figure 1 As shown, the joint inversion method for three-dimensional quality factors includes: S101. Obtain the target seismic body wave displacement spectrum from the earthquake to the station in the frequency domain.
[0015] In embodiments of the present invention, the target seismic body wave displacement spectrum in the frequency domain is obtained by performing a Fourier transform on the seismic body wave (first-arrival P-wave or first-arrival S-wave) signal in the time domain. This process converts the target seismic wave signal from the time domain to the frequency domain, allowing for the analysis of the energy distribution and phase information of different frequency components. S102. Obtain the first average attenuation time using the initial three-dimensional quality factor, source, receiver location, and the first average attenuation travel time formula.
[0016] In an embodiment of the present invention, the earthquake source is the starting point of the earthquake, the station is the location where the receiver is installed, and the initial three-dimensional quality factor can be randomly selected or obtained experimentally; no limitation is imposed here. The first average attenuation travel time is calculated based on the initial three-dimensional quality factor, the earthquake source, the receiver location, and the first average attenuation travel time formula, which is as follows: In the above formula, The first average attenuation travel time is obtained by integrating 1 / (VQ) along the ray path (i.e., the propagation path from the source to the receiver). `source` is the location of the seismic source, `receiver` is the location of the receiver, `V` is the three-dimensional propagation velocity of the seismic wave, `Q` is the initial three-dimensional quality factor, and `dr` is the infinitesimal element of the ray path from the source to the receiver.
[0017] S103. Calculate the theoretical spectral amplitude at the frequency using the first average attenuation travel time, initial corner frequency, initial site effect, initial low-frequency amplitude, frequency and amplitude formulas; and convert the amplitude of the target seismic body wave displacement spectrum from the earthquake to the station into the observed spectral amplitude in the natural logarithm domain.
[0018] In some embodiments of the present invention, the amplitude formula is as follows: In the above formula, For earthquake, For the station, For the first The earthquake to the first Each station on the frequency The theoretical spectral amplitude at that location, Let i be the corner frequency of the i-th earthquake. , The first average attenuation travel time along the ray path from the i-th earthquake to the j-th station It can be calculated using the first average decay travel time formula. For the first station and frequency related site effect, For the first earthquake to the first station on the low-frequency amplitude , the low-frequency amplitude considers the geometric diffusion factor independent of frequency, and contains the seismic radiation pattern. The low-frequency amplitude in the formula is related to the corner frequency of the source, and here the source model is assumed , , . Where the initial corner frequency, the initial site effect, and the initial low-frequency amplitude can be obtained from experiments during inversion.
[0019] Further, the first average attenuation time, the initial corner frequency, the initial site effect, the initial low-frequency amplitude, and the frequency are substituted into the amplitude formula to obtain the theoretical spectral amplitude at frequency f.
[0020] Further, for the amplitude formula, the partial derivative of the observed spectral amplitude with respect to the corner frequency is , the partial derivative with respect to the low-frequency amplitude is , the partial derivative with respect to the first average attenuation factor is , and the partial derivative with respect to the site response is . Here, it is assumed that the quality factor is independent of frequency.
[0021] S104, based on the observed spectral amplitude and the theoretical spectral amplitude in the amplitude formula, a first residual vector formula is obtained; the first residual vector formula includes parameter perturbation.
[0022] In an embodiment of the present application, based on the Taylor expansion of the observed spectral amplitude with respect to the theoretical spectral amplitude in the amplitude formula at the initial corner frequency, the initial site effect, the initial low-frequency amplitude, and the initial first average attenuation time, and taking the first-order approximation of the parameters of the amplitude formula, a first residual vector formula is obtained, wherein the first residual vector formula is as follows: In the above formula, is the residual vector between the observed spectral amplitude and the theoretical spectral amplitude, , the dimension is , is the total number of body wave amplitude spectrum frequency points, is the index of the earthquake, and the total number of body wave amplitude spectrum frequency points is composed of the number of rates and the number of stations, wherein is the spectral amplitude residual vector of the corresponding nth station, is the first frequencies, and the number of frequency points of each frequency spectrum is , is a sensitivity matrix of the source parameters and the first average attenuation travel time, is the source parameter (low-frequency amplitude and corner frequency ) of the mth earthquake and the first perturbation of the first average attenuation travel time, is a sensitivity matrix of the station site response, is the second perturbation of the site response.
[0023] wherein the spectral amplitude residual vector is the difference between the observed spectral amplitude in the natural logarithm domain and the theoretical spectral amplitude in the natural logarithm domain.
[0024] Exemplarily, it is assumed that there is a group of earthquakes m = 1, 2,..., M, each of which is recorded by n stations, wherein n = 1, 2,..., N, and N is the maximum number of stations. The observed amplitude spectrum of each earthquake to the station has p = 1, 2,..., P frequencies. Therefore, for one earthquake, there are unknowns (the low-frequency amplitude and the corner frequency of the earthquake) on the ray path. It is assumed that the mth earthquake is recorded by n stations, and the total number of observed body wave amplitude spectrum frequency points is .
[0025] S105, the seismic wave three-dimensional propagation speed and the initial three-dimensional quality factor in the first average attenuation travel time formula are gridded to obtain a second average attenuation travel time formula; and the second average attenuation travel time formula is substituted into the first residual vector formula to obtain a second residual vector formula.
[0026] In some embodiments of the present application, the first average attenuation travel time formula is discretized, that is, the seismic wave three-dimensional propagation speed and the initial three-dimensional quality factor in the first average attenuation travel time are gridded to obtain a second average attenuation travel time formula, as follows: In the above formula, is the number of grid points, is the index number of the grid point, and are the seismic wave three-dimensional propagation speed value and the initial three-dimensional quality factor value at the mth grid point, is the length of the ray from the source to the receiver on the ith grid point, is the travel time from the source to the ith grid point, is is the attenuation factor at the i-th grid point. This equation is the second average attenuation traveltime equation.
[0027] Further, the second average attenuation traveltime is used to replace the spectrum fitting inversion part, i.e. the second average attenuation traveltime is used to replace the part in the first residual vector calculation equation to obtain the second residual vector calculation equation. The second residual vector calculation equation is as follows: In the above equation, is the index of the earthquake, is the three-dimensional attenuation structure perturbation, where represents the transpose of the vector, is the perturbation of the attenuation factor at the i-th grid point. is the source parameter perturbation, which contains the low-frequency amplitude perturbation of the i-th earthquake to each station and the corner frequency perturbation of the i-th earthquake. , where is the low-frequency amplitude perturbation of the i-th station. is the sensitivity matrix of the observed spectrum value of the i-th earthquake to all stations with respect to the attenuation structure, which has a dimension of , n is the station index, is the number of frequency points of each spectrum, is the number of grid points, which is consistent with the definition of the first residual vector calculation equation. is the sensitivity matrix of the observed spectrum value of the i-th earthquake to all stations with respect to the source parameters (low-frequency amplitude and corner frequency ), is the sensitivity matrix of the observed spectrum value of the i-th earthquake to all stations with respect to the site effect of each station, which is the same as the first residual equation. where, is as shown in the following equation:
[0028] where, is as shown in the following equation: In the above equation, is the sensitivity matrix of the i-th spectrum to the attenuation structure, which has a dimension of .
[0029] where, is as shown in the following equation: In the above formula, is the partial derivative of the amplitude value of the n-th frequency bin in the m-th spectrum with respect to t*, is the travel time of the m-th earthquake to the n-th station at the L-th grid point.
[0030] wherein is shown in the following formula: wherein, , , denotes transposition, denotes the partial derivative of the spectral amplitude of the n-th frequency of the m-th station spectrum with respect to the corner frequency , denotes the partial derivative of the spectral amplitude of the n-th frequency of the m-th station spectrum with respect to the corner frequency .
[0031] S106, a first-order smoothing constraint and a damping parameter are added to the second residual vector formula to obtain a third residual vector formula, and the value of the parameter perturbation quantity is obtained by solving the third residual formula.
[0032] In the embodiments of the present application, the first-order smoothing constraint generally refers to imposing a limit on the first derivative (i.e., the rate of change) of a function or a signal. The purpose of this constraint is to ensure that the obtained result does not have abrupt changes, but is relatively smooth. The damping parameter mainly appears in methods such as ridge regression (Ridge Regression), Tikhonov regularization, etc., and is used to control the strength of the regularization term. For the observed spectrum data of all earthquakes to all stations, a first-order smoothing constraint and a damping parameter are added to the second residual vector formula to obtain a third residual vector formula, as follows: In the above formula, is a smoothing factor, is a first-order smoothing matrix, is a damping factor, is an identity matrix. Wherein, , , and are known, and the conjugate gradient descent algorithm can be used to solve in the formula. Wherein , , , , wherein Indicates transpose. For the first The spectral amplitude residuals from each earthquake to all stations are shown in the second residual vector formula. Similarly, , , It is consistent with the definition of the second residual vector formula.
[0033] S107. The values of the parameter disturbances are superimposed on the initial three-dimensional quality factor, the initial corner frequency, the initial low-frequency amplitude, and the initial site effect, respectively, to obtain the new initial three-dimensional quality factor, the new initial corner frequency, the new initial low-frequency amplitude, and the new initial site effect for the next iteration.
[0034] S108. Substitute the new initial three-dimensional quality factor, the new initial corner frequency, the new initial low-frequency amplitude, and the new initial site effect back into the first residual vector formula to obtain the first residual vector. When the difference between the first residual vector and the amplitude satisfies the preset condition, the new initial three-dimensional quality factor is used as the target three-dimensional quality factor.
[0035] In an embodiment of the invention, the solved disturbance values are superimposed onto the initial three-dimensional quality factor, initial corner frequency, initial low-frequency amplitude, and initial field effect corresponding to the current round, respectively, to obtain new initial three-dimensional quality factor, new initial corner frequency, new initial low-frequency amplitude, and new initial field effect for the next iteration. These new initial three-dimensional quality factor, new initial corner frequency, new initial low-frequency amplitude, and new initial field effect are then substituted back into the first residual vector formula to obtain the first residual vector. When the difference between the first residual vector and the amplitude satisfies a preset condition, the new initial three-dimensional quality factor is used as the target three-dimensional quality factor. When the difference between the first residual vector and the amplitude does not satisfy the preset condition, steps S102 to S108 are repeated, iterating continuously until the target three-dimensional quality factor is obtained.
[0036] It can be understood that in the embodiments of the present application, the target seismic body wave displacement spectrum of the earthquake to the station in the frequency domain is obtained; the first average attenuation travel time is obtained through the initial three-dimensional quality factor, the source, the receiver position and the first average attenuation travel time formula; the theoretical spectrum amplitude at the frequency is calculated through the first average attenuation travel time, the initial corner frequency, the initial site effect, the initial low-frequency amplitude, the frequency and the amplitude formula; and the amplitude of the target seismic body wave displacement spectrum of the earthquake to the station is converted into the observed spectrum amplitude in the natural logarithm domain; the first residual vector formula is obtained based on the observed spectrum amplitude and the amplitude formula; the first residual vector formula includes the parameter perturbation; the seismic wave three-dimensional propagation velocity in the first average attenuation travel time formula and the initial three-dimensional quality factor are gridded to obtain the second average attenuation travel time formula; and the second average attenuation travel time formula is substituted into the first residual vector formula to obtain the second residual vector formula; The first-order smoothing constraint and the damping parameter are added to the second residual vector formula to obtain the third residual vector formula; and the value of the parameter perturbation is obtained by solving the third residual formula; the value of the parameter perturbation is superimposed into the initial three-dimensional quality factor, the initial corner frequency, the initial low-frequency amplitude and the initial site effect respectively to obtain the new initial three-dimensional quality factor, the new initial corner frequency, the new initial low-frequency amplitude and the new initial site effect of the next round of iteration; the new initial three-dimensional quality factor, the new initial corner frequency, the new initial low-frequency amplitude and the new initial site effect are substituted back into the first residual vector formula to obtain the first residual vector, and when the difference between the first residual vector and the amplitude meets the preset condition, the new initial three-dimensional quality factor is taken as the target three-dimensional quality factor. In this process, by jointly inverting the source parameters and the quality factor structure, the process of extracting the average factor is avoided, the new technology can avoid the defect that the average attenuation factors of the same ray are inconsistent, in addition, the constraint that different rays cross and share the same quality factor is also used, which reduces the coupling between the quality factor and the seismic corner frequency to some extent, reduces the amount of inversion parameters and can more reliably obtain the quality factor structure and the source parameters of the earthquake.
[0037] In the embodiments of the present application, S101 further includes S10, which is described by the following steps.
[0038] S10, obtaining the seismic body wave displacement spectrum of the earthquake to the station in the frequency domain, and obtaining the target seismic body wave displacement spectrum after removing the instrument response in the seismic body wave displacement spectrum.
[0039] In some embodiments of the present application, after removing the instrument response from the seismic body wave displacement spectrum, it usually refers to correcting the original record data to eliminate the influence of the seismograph or other recording equipment on the signal, so as to restore the true ground displacement spectrum, i.e. the target seismic body wave displacement spectrum.
[0040] The following is a comparative example of the joint inversion three-dimensional quality factor method and the traditional two-step quality factor imaging method on synthetic data according to the present application, as follows: Figure 2 The effect schematic diagram of the synthetic test seismic event and station geographical position distribution provided by the embodiment of the present application comprises a plurality of seismic events and a plurality of stations. In Figure 2 , the black circles are seismic events, and the red triangles are stations. The abscissa represents longitude, and the ordinate represents latitude. N and E represent north latitude and east longitude, respectively. It is assumed that each seismic event is recorded by all stations, and a total of 26695 rays are generated. In this experiment, 1405 seismic corner frequencies and 26695 low-frequency amplitudes are randomly generated, the corner frequency range is 3Hz~17Hz, and the low-frequency amplitude is within . The site effect of each station is set to 1. Here, the joint inversion of the P-wave three-dimensional quality factor is taken as an example, in which the three-dimensional P-wave velocity model adopts the local regional P-wave velocity structure (known). In addition, on the basis of a three-dimensional uniform quality factor (Q value size is 400) model, a positive and negative 50% perturbation is performed between adjacent grid points to generate the real P-wave quality factor (Qp) model in the synthetic test, and subsequently, the real P-wave quality factor (Qp) model is used as the basis for forward average attenuation travel time. Figure 3 The result comparison schematic diagram of the traditional two-step quality factor method and the joint inversion quality factor method at different depth slices provided by the synthetic test of the embodiment of the present application is shown. The real quality factor model distribution at depths of 7km and 10km is shown in the left column of Figure 3 . According to the given parameters, the amplitude spectrum of the P-wave is synthesized, and the theoretical average attenuation travel time is obtained, and 5% random noise is added to the synthetic amplitude spectrum to obtain the synthetic observation spectrum. In the case of a given uniform P-wave initial quality factor model of 400, the calculation efficiency and quality factor imaging accuracy of the traditional two-step quality factor imaging method and the joint inversion three-dimensional quality factor method are compared and analyzed. dQp is the perturbation amplitude of Qp.
[0041] The joint inversion three-dimensional quality factor method (new method) of the present application and the traditional two-step quality factor imaging algorithm are compared as shown in Table 1 below: Table 1 Comparison of two methods The model parameter quantity of the new method (the joint inversion three-dimensional quality factor method of the application) is related to the three-dimensional model grid quantity (ncell), the path quantity (npath), the earthquake quantity (neq), the station quantity (nsta) and the frequency point quantity of each frequency spectrum (nfreqs), and the specific model parameter quantity is ncells+npath+neq+nsta*nfreqs. Generally, the model grid quantity of the new method is less than the path quantity, and therefore, compared with the traditional two-step quality factor imaging algorithm, the joint inversion three-dimensional quality factor method of the application has the characteristics of small memory and high inversion efficiency.
[0042] Figure 3 The comparative diagram of the results of the traditional two-step quality factor method and the joint inversion quality factor method in the synthetic test at different depth slices is provided for the embodiments of the application. The left column of diagrams is the quality factor distribution diagram of the real P-wave quality factor model at the depths of 7km and 10km. The middle column of diagrams is the P-wave quality factor recovered by the traditional two-step quality factor method at the depths of 7km and 10km. The right column of diagrams is the P-wave quality factor recovered by the joint inversion quality factor method at the depths of 7km and 10km. As can be seen from the diagrams, the joint inversion quality factor can better recover the real P-wave quality factor in the synthetic test. In addition, from the fitting residual of the average attenuation travel time and the corner frequency in Table 1, the joint inversion quality factor method is better than the traditional two-step quality factor method.
[0043] In summary, the joint inversion quality factor method is superior to the traditional two-step quality factor method in the calculation efficiency and the accuracy.
[0044] The above embodiments are only used for describing the embodiments of the application, but not for limiting the embodiments of the application. The ordinary skilled in the art can make various changes and modifications without departing from the spirit and scope of the embodiments of the application, and all equivalent technical solutions also belong to the scope of the embodiments of the application. The patent protection scope of the embodiments of the application should be defined by the claims.
Claims
1. A joint inversion three-dimensional quality factor method, characterized in that, The method comprises the following steps: obtaining a target seismic body wave displacement spectrum of an earthquake to a station in a frequency domain; obtaining a first average attenuation travel time through an initial three-dimensional quality factor, a source, a receiver position and a first average attenuation travel time formula; calculating a theoretical spectral amplitude at a frequency in a natural logarithm domain through the first average attenuation travel time, an initial corner frequency, an initial site effect, an initial low-frequency amplitude, a frequency and an amplitude formula; and converting an amplitude of the target seismic body wave displacement spectrum of the earthquake to the station into an observed spectral amplitude in the natural logarithm domain; obtaining a first residual vector formula based on the observed spectral amplitude and the theoretical spectral amplitude in the amplitude formula; the first residual vector formula comprises a parameter perturbation; the three-dimensional propagation velocity of the seismic wave in the first average attenuation travel time formula and the initial three-dimensional quality factor are gridded to obtain a second average attenuation travel time formula; the second average attenuation travel time formula is substituted into the first residual vector formula to obtain a second residual vector formula; a first-order smoothing constraint and a damping parameter are added to the second residual vector formula to obtain a third residual vector formula; and the third residual formula is solved to obtain a value of the parameter perturbation; the value of the parameter perturbation is superimposed on the initial three-dimensional quality factor, the initial corner frequency, the initial low-frequency amplitude and the initial site effect respectively to obtain a new initial three-dimensional quality factor, a new initial corner frequency, a new initial low-frequency amplitude and a new initial site effect of the next round of iteration; the new initial three-dimensional quality factor, the new initial corner frequency, the new initial low-frequency amplitude and the new initial site effect are substituted back into the first residual vector formula to obtain a first residual vector; when a difference between the first residual vector and the amplitude satisfies a preset condition, the new initial three-dimensional quality factor is taken as a target three-dimensional quality factor.
2. The method of claim 1, wherein, Before the step of obtaining the target seismic body wave displacement spectrum of the earthquake to the station in the frequency domain, the method further comprises: obtaining a seismic body wave displacement spectrum of the earthquake to the station in the frequency domain, and removing an instrument response in the seismic body wave displacement spectrum to obtain a target seismic body wave displacement spectrum.
3. The method according to claim 1 or 2, characterized in that, The parameter perturbation comprises a source parameter perturbation, an average attenuation factor perturbation, a site response perturbation and a three-dimensional attenuation structure perturbation, wherein the source parameter perturbation comprises a low-frequency amplitude perturbation and a corner frequency perturbation of the earthquake to each station.
4. The method of claim 1, wherein, The step of obtaining the first residual vector formula based on the observed spectrum and the amplitude formula comprises: performing Taylor expansion on the amplitude formula based on the observed spectrum and taking a first-order approximation of parameters in the amplitude formula to obtain a first residual vector formula of a theoretical spectrum, and the first residual vector formula is shown in the following formula: In the above formula, The residual vector between the observed spectral amplitude and the theoretical spectral amplitude. , dimension , This represents the total number of frequency points in the body wave amplitude spectrum. As an index for earthquakes, the total number of frequency points in the body wave amplitude spectrum consists of the number of rates and the number of stations, where This represents the spectral amplitude residual vector of the corresponding nth station. For the first There are 10 frequencies, and the number of frequency points in each spectrum is 10. , The sensitivity matrix includes source parameters and the first average attenuation travel time. For the first The source parameters of the earthquake and the first disturbance of the first mean attenuation travel time, This is the sensitivity matrix for the site response of the station. This is the second disturbance in the site response.
5. The method of claim 1, wherein, the first average attenuation travel time formula is shown in the following formula: In the above formula, is the first average attenuation traveltime, which is obtained by integrating 1 / (VQ) along a raypath. source is the source location, receiver is the receiver location, V is the three-dimensional propagation velocity of the seismic wave, Q is the initial three-dimensional quality factor, and dr is the raypath differential of the source to the receiver location. the three-dimensional propagation velocity of the seismic wave in the first average attenuation travel time formula and the initial three-dimensional quality factor are gridded to obtain the second average attenuation travel time formula, as shown in the following formula: In the above formula, L is the number of grid points, and i is the index number of the grid point. and The first The three-dimensional propagation velocity of seismic waves at each grid point and the initial three-dimensional quality factor. Assign the length of the ray from the seismic source to the receiver location at the i-th grid point. The travel time from the earthquake source to the i-th grid point is... for , which is the attenuation factor at the i-th grid point.
6. The method of claim 1, wherein, the second residual vector formula is shown in the following formula: In the above formula, is the index of the earthquake, is the three-dimensional attenuation structure perturbation quantity, where represents the transpose of the vector, is the perturbation quantity of the attenuation factor at the th grid point. is the source parameter perturbation quantity, which contains the low-frequency amplitude perturbation quantity of the th earthquake to each station and the corner frequency perturbation quantity of the th earthquake , where is the low-frequency amplitude perturbation quantity of the th station . is the sensitivity matrix of the observed spectral amplitude of the th earthquake to the attenuation structure at all stations, which has a dimension of , n is the station index, is the number of frequency points of each spectrum, is the number of grid points, consistent with the definition of the first residual vector calculation formula. is the sensitivity matrix of the observed spectral amplitude of the th earthquake to the source parameter at all stations, is the sensitivity matrix of the observed spectral amplitude of the th earthquake to the site effect of each station at all stations, consistent with the first residual formula; wherein as shown in the following equation: wherein is the sensitivity matrix of the n-th spectral pair attenuation structure, of dimension wherein as shown in the following equation: In the above formula, is the partial derivative of the amplitude value of the mth earthquake to the nth station at the Lth grid point with respect to t*, is the travel time of the mth earthquake to the nth station at the Lth grid point. wherein as shown in the following equation: wherein , , denotes the transpose, denotes the partial derivative of the spectral amplitude of the spectrum of the th station at the th frequency with respect to the corner frequency , denotes the partial derivative of the spectral amplitude of the spectrum of the th station at the th frequency with respect to the corner frequency .
7. The method according to any one of claims 1 to 6, characterized in that, the third residual vector formula is shown in the following formula: where is the observed spectral data for all earthquakes to all stations, the first residual vector formula, and the second residual vector formula, where is the first smoothing constraint and is the damping parameter. is the smoothing factor, is the first smoothing matrix, is the damping factor, is the identity matrix. Where, , , and Given, the solution to the equation may be found using the conjugate gradient descent algorithm. Where , , , where denotes the transpose, is the spectral amplitude residual for the th earthquake to all stations.