CN116466390B - Real-time monitoring and positioning method for earthquake induced by large reservoir - Google Patents
Real-time monitoring and positioning method for earthquake induced by large reservoir Download PDFInfo
- Publication number
- CN116466390B CN116466390B CN202310129707.1A CN202310129707A CN116466390B CN 116466390 B CN116466390 B CN 116466390B CN 202310129707 A CN202310129707 A CN 202310129707A CN 116466390 B CN116466390 B CN 116466390B
- Authority
- CN
- China
- Prior art keywords
- seismic
- model
- event
- time
- positioning
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 88
- 238000012544 monitoring process Methods 0.000 title claims abstract description 49
- 230000002441 reversible effect Effects 0.000 claims abstract description 18
- 239000013307 optical fiber Substances 0.000 claims abstract description 14
- 239000011159 matrix material Substances 0.000 claims abstract description 13
- 238000012937 correction Methods 0.000 claims abstract description 4
- 238000009826 distribution Methods 0.000 claims description 33
- 238000005070 sampling Methods 0.000 claims description 9
- 230000008569 process Effects 0.000 claims description 7
- 230000003287 optical effect Effects 0.000 claims description 6
- 238000004364 calculation method Methods 0.000 claims description 5
- 230000008859 change Effects 0.000 claims description 5
- 238000004891 communication Methods 0.000 claims description 4
- 230000003094 perturbing effect Effects 0.000 claims description 4
- 230000007704 transition Effects 0.000 claims description 4
- RZVAJINKPMORJF-UHFFFAOYSA-N Acetaminophen Chemical compound CC(=O)NC1=CC=C(O)C=C1 RZVAJINKPMORJF-UHFFFAOYSA-N 0.000 claims description 2
- 238000011156 evaluation Methods 0.000 claims description 2
- 230000003203 everyday effect Effects 0.000 claims description 2
- 230000035945 sensitivity Effects 0.000 claims description 2
- 238000009827 uniform distribution Methods 0.000 claims description 2
- 230000004075 alteration Effects 0.000 claims 1
- 238000012217 deletion Methods 0.000 claims 1
- 230000037430 deletion Effects 0.000 claims 1
- 230000009286 beneficial effect Effects 0.000 description 6
- 230000015572 biosynthetic process Effects 0.000 description 4
- 238000010586 diagram Methods 0.000 description 4
- 239000000835 fiber Substances 0.000 description 3
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 3
- 238000012512 characterization method Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 239000003651 drinking water Substances 0.000 description 2
- 235000020188 drinking water Nutrition 0.000 description 2
- 230000009977 dual effect Effects 0.000 description 2
- 230000000763 evoking effect Effects 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000006698 induction Effects 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/282—Application of seismic models, synthetic seismograms
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
- G01H9/00—Measuring mechanical vibrations or ultrasonic, sonic or infrasonic waves by using radiation-sensitive means, e.g. optical means
- G01H9/004—Measuring mechanical vibrations or ultrasonic, sonic or infrasonic waves by using radiation-sensitive means, e.g. optical means using fibre optic sensors
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/288—Event detection in seismic signals, e.g. microseismics
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
- G01V1/362—Effecting static or dynamic corrections; Stacking
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/30—Assessment of water resources
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Remote Sensing (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Geophysics (AREA)
- Business, Economics & Management (AREA)
- Emergency Management (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The application relates to the technical field of reservoir induced earthquake monitoring, and discloses a real-time monitoring and positioning method for large reservoir induced earthquake, which comprises the following steps: an amphibious integrated distributed optical fiber acoustic wave sensing monitoring system is arranged in a multiple reservoir seismic zone, and meter-level space precision continuous observation and seismic signal acquisition are carried out on seismic waves formed by near-surface earthquakes; picking up the relative arrival time of the DAS seismic signals by using a cross-correlation method, and constructing a time difference matrix; correcting the anisotropic VTI speed model by adopting a Bayesian theory and reversible jump Markov chain Monte Carlo algorithm; positioning the natural seismic event by using a double-difference positioning method based on the updated VTI speed model; the method comprises the steps of comprehensive seismic signal acquisition, velocity model correction and double-difference positioning, and the accurate position of the induced earthquake is measured. The application aims to solve the problem of low positioning precision of real-time monitoring of induced seismic events in reservoirs.
Description
Technical Field
The application relates to the technical field of reservoir induced earthquake monitoring, in particular to a real-time monitoring and positioning method for large reservoir induced earthquake.
Background
A large number of reservoirs in China are built on geological fault zones, and the fault zone activities lead to frequent induction of various earthquake events, so that urban drinking water and life and property safety are threatened. Accurate positioning of induced earthquake is helpful to describe a hidden micro fault zone below the reservoir, and dynamic monitoring of reservoir safety is guaranteed. However, the conventional seismograph monitoring method is difficult to arrange the earthquake equipment in the reservoir for a long time, and the earthquake signals cannot be obtained in a short distance. Meanwhile, the earthquake monitoring station has the defect of large interval distance, and the complete earthquake wave field data is difficult to obtain. Various fractures and different lithologies in the formation can also cause the formation velocity to be anisotropic, and the conventional isotropic velocity model can be difficult to accurately locate seismic events. In addition, the seismic event locations obtained with seismic wave travel time are also subject to noise interference. The common influence of observation equipment, speed model errors and earthquake travel time noise causes poor positioning accuracy of reservoir induced earthquakes, and potential micro fault zones are difficult to find. The distributed optical fiber acoustic wave sensing technology (Distributed fiberAcoustic Sensing, DAS) developed in recent years can be distributed in water for a long time, has good environment adaptability, can record a seismic wave field according to meter-level precision, and is beneficial to obtaining the accurate position of a seismic signal through closely monitoring seismic wave inversion. Meanwhile, the anisotropic VTI (Vertical Transverse Isotropy) velocity model is utilized to depict the stratum medium velocity field, so that the real stratum velocity structure can be more approximated. The double-difference positioning method is a relative positioning method, can reduce dependence on a velocity model and the accuracy of the arrival time of the earthquake, and deduces the spatial position of the induced earthquake event by means of the known more accurate main earthquake event position. Therefore, compared with the traditional earthquake positioning method, the combined application of the DAS optical fiber monitoring, the VTI speed model and the double-difference positioning method is more beneficial to obtaining the accurate reservoir induced earthquake event position.
Accordingly, the prior art is still in need of improvement and development.
Disclosure of Invention
Aiming at the defects in the prior art, the application provides a real-time monitoring and positioning method for the induced earthquake of a large reservoir, which aims to solve the problem of low real-time monitoring and positioning precision of the induced earthquake event in the reservoir.
In order to achieve the above purpose, the present application may be performed by the following technical scheme:
a real-time monitoring and positioning method for large reservoir induced earthquake comprises the following steps:
an amphibious integrated distributed optical fiber acoustic wave sensing monitoring system is arranged in a multiple reservoir seismic zone, and meter-level space precision continuous observation and seismic signal acquisition are carried out on seismic waves formed by near-surface earthquakes;
picking up the relative arrival time of the DAS seismic signals by using a cross-correlation method, and constructing a time difference matrix; correcting the anisotropic VTI speed model by adopting a Bayesian theory and reversible jump Markov chain Monte Carlo algorithm;
positioning the natural seismic event by using a double-difference positioning method based on the updated VTI speed model;
the method comprises the steps of comprehensive seismic signal acquisition, velocity model correction and double-difference positioning, and the accurate position of the induced earthquake is measured.
Compared with the prior art, the application has the beneficial effects that: the application provides a novel real-time monitoring and positioning method for large-scale reservoir induced earthquake, which is characterized in that a amphibious distributed optical fiber acoustic wave sensing (Distributed fiber Acoustic Sensing, DAS) monitoring system is arranged in a multiple reservoir seismic zone to continuously observe the meter-level space precision of earthquake waves formed by near-surface earthquake; picking up the relative arrival time of the DAS seismic signals by using a cross-correlation method, and constructing a time difference matrix; correcting the anisotropic VTI speed model by adopting a Bayesian theory and reversible jump Markov chain Monte Carlo (reversible jump Markov Chain Monte Carlo, rjMCMC) algorithm; and positioning the natural seismic event by using a double-difference positioning method based on the updated VTI speed model. The method has the advantages that the continuous wave field of the natural earthquake in the reservoir can be recorded with high spatial resolution, meanwhile, the anisotropic VTI speed model is obtained in an optimized mode, and the accurate reservoir bottom earthquake event position can be obtained by combining the double-difference positioning method with higher positioning precision. The method is beneficial to improving the accuracy of the hidden fault zone characterization below the reservoir and provides a good technical means for reservoir safety assessment.
Drawings
In order to more clearly illustrate the technical solutions of the embodiments of the present application, the following description will briefly explain the drawings needed in the embodiments, and it is obvious that the drawings in the following description are only some embodiments of the present application, and other drawings can be obtained according to these drawings without inventive effort for a person skilled in the art.
FIG. 1 is a flow chart of a method for real-time monitoring and positioning of large reservoir induced earthquakes according to an embodiment of the present application;
FIG. 2 is a schematic diagram of a distributed fiber optic sensing device for monitoring reservoir induced seismic signals according to an embodiment of the present application;
FIG. 3 is a flow chart for correcting a VTI media velocity model using the rjMCMC method;
FIG. 4 is a schematic diagram of a dual differential positioning method for obtaining the accurate position of the reservoir induced earthquake.
Detailed Description
The following description of the embodiments of the present application will be made clearly and completely with reference to the accompanying drawings, in which it is apparent that the embodiments described are only some embodiments of the present application, but not all embodiments. All other embodiments, which can be made by those skilled in the art based on the embodiments of the application without making any inventive effort, are intended to be within the scope of the application.
Examples:
it should be noted that the terms "first," "second," and the like in the description and the claims of the present application and the above figures are used for distinguishing between similar objects and not necessarily for describing a particular sequential or chronological order. It is to be understood that the data so used may be interchanged where appropriate such that the embodiments of the application described herein may be implemented in sequences other than those illustrated or otherwise described herein. Furthermore, the terms "comprises," "comprising," and "having," and any variations thereof, are intended to cover a non-exclusive inclusion, such that a process, method, system, article, or apparatus that comprises a list of steps or elements is not necessarily limited to those steps or elements expressly listed or inherent to such process, method, article, or apparatus.
In the description of the present application, the meaning of "plurality" means at least two, for example, two, three, etc., unless specifically defined otherwise. Furthermore, unless explicitly specified and limited otherwise, the terms "mounted," "connected," and "connected" are to be construed broadly, and may be either fixedly connected, detachably connected, or integrally connected, for example; can be mechanically or electrically connected; can be directly connected or indirectly connected through an intermediate medium, and can be communication between two elements. The specific meaning of the above terms in the present application will be understood in specific cases by those of ordinary skill in the art.
The word "exemplary" is used hereinafter to mean "serving as an example, embodiment, or illustration. Any embodiment described as "exemplary" is not necessarily to be construed as preferred or advantageous over other embodiments.
The earthquake of the reservoir has important influence on the urban drinking water safety and life and property. Conventional seismic monitoring equipment has difficulty in ensuring that seismic signals are monitored in real time in water for a long time. It is also difficult for a land discretely deployed seismic station to obtain full wavefield information for the seismic. The lack of the traditional method for collecting and processing the reservoir seismic signals causes that the hidden fault zone at the bottom of the reservoir is difficult to find, and the healthy operation of the reservoir is threatened at any time. Therefore, the device capable of recording the information of the full wave field of the earthquake and suitable for underwater long-time monitoring is developed correspondingly to develop an accurate real-time positioning algorithm of the earthquake event, and is beneficial to the fine depiction of reservoir stratum. The application utilizes the amphibious integrated distributed optical fiber acoustic wave sensing monitoring system to continuously observe the meter-scale space precision of the wave field formed by the induced earthquake at the bottom of the reservoir. Meanwhile, a cross-correlation method is used for picking up the relative arrival time of the DAS seismic signals, an arrival time difference matrix is constructed, and the influence of noise on a travel time result is reduced. And correcting the reservoir stratum VTI speed model based on the Bayesian theory and the reversible jump Markov chain Monte Carlo algorithm to obtain a more accurate stratum speed field. The updated VTI speed model is used for positioning the natural seismic event by using a double-difference positioning method, so that the influence of speed errors on positioning results is effectively reduced. Based on the method, the application provides a novel reservoir induced earthquake real-time monitoring and accurate positioning method.
Referring to fig. 1, the method for monitoring and positioning the induced earthquake of the large reservoir in real time can comprise the following steps:
step 1: an amphibious integrated distributed optical fiber acoustic wave sensing monitoring system is arranged in a multiple reservoir seismic zone, and meter-scale space precision continuous observation and seismic signal acquisition are carried out on seismic waves formed by near-surface earthquakes.
Specifically, referring to fig. 2, fig. 2 is a schematic diagram of a distributed optical fiber sensing device for monitoring reservoir induced seismic signals. The earthquake positioning method is based on an amphibious integrated distributed optical fiber sound wave sensing monitoring system, and can be used for continuously observing the meter-scale space precision of earthquake waves formed by near-surface earthquakes. The method comprises the following specific steps:
(1) And (3) arranging communication optical cables on land in the areas where the reservoir induces earthquake, and arranging submarine cables at the bottom of the reservoir. Whether land cable or sea cable, it is desirable to maintain good coupling contact with the formation in order to accept changes in formation stress caused by seismic waves.
(2) And the optical cable signal demodulator is placed in a machine room or other safe positions, so that stable power supply of the demodulator is ensured. The demodulator interface is connected in series with the land cable and the sea cable. The strain caused by the seismic waves felt by the land cable and the sea cable is resolved by the signal demodulator.
(3) The demodulator is used for setting the space monitoring strain sensitivity of the optical cable to be about one meter, the time sampling rate is more than 200Hz, and the continuous wave field of the earthquake transmitted to the ground surface is monitored without interval every day.
Step 2: picking up the relative arrival time of the DAS seismic signals by using a cross-correlation method, and constructing a time difference matrix; and correcting the anisotropic VTI speed model by adopting a Bayesian theory and reversible jump Markov chain Monte Carlo algorithm.
Specifically, the method for picking up the relative arrival time of the DAS seismic signals by using the cross-correlation method constructs a time difference matrix, and specifically comprises the following steps:
first, a long-short window ratio method is used to pick up valid seismic events. The long term window (LTA) is an Average value of energy of a longer time signal sampling length, and the short term window (Short Term Average, STA) is an Average value of energy of a shorter time signal sampling length:
wherein S (N) is the signal amplitude obtained by DAS monitoring, N is the number of long-time window signal samples, M is the number of short-time window signal samples, and N > M. Ratio is the Ratio of the energy of the longer time signal to the shorter time signal obtained by recording. When it is greater than a threshold, the monitored signal may be determined to be caused by a seismic event.
And secondly, accurately picking up the relative arrival time of different optical fiber monitoring channels by using a cross-correlation method, and constructing a time difference matrix. The cross-correlation method formula can be expressed as:
wherein x is 1 And x 2 Monitoring seismic traces for different DAS, R x1x2 Is a similarity coefficient. And precisely obtaining similarity coefficients of the seismic events extracted by Ratio among different optical fiber channels by using a cross-correlation formula, and selecting a time difference with the maximum similarity coefficient of different seismic channels to form a matrix.
Referring to fig. 3, fig. 3 is a flow chart for correcting VTI media velocity model using the rjMCMC method. Correcting the anisotropic VTI speed model by adopting a Bayesian theory and reversible jump Markov chain Monte Carlo algorithm, wherein the specific steps comprise:
the seismic waves may be split into qP, qSV and qSH waves in the VTI medium. From five anisotropic parameters [ alpha ] 0 ,β 0 ,ε,γ,δ]And (5) determining. The velocity model includes a number n of horizons and a horizon depth D. Thus, the inverse VTI velocity model parameter expression is determined as: m= [ n, D, α 0 ,β 0 ,ε,γ,δ]. Wherein alpha is 0 And beta 0 Representing the speed of each layer along the symmetry axis P wave and SV wave; epsilon represents the ratio of the velocity of each layer of qP wave in the horizontal and vertical directions; gamma represents the ratio of the velocity of each layer qSH wave in the horizontal and vertical directions; delta represents the rate of change of the qP wave of each layer in the vertical direction.
The Bayesian theory expression is determined as follows:
where d represents the observed data and m is the model parameter. p (m) is prior model information, p (d|m) is a likelihood function, p (m|d) is posterior model probability, and p (d) is overall probability of observed data in model space and is a constant.
And by utilizing a Bayesian inversion algorithm and combining prior information of each model parameter and likelihood functions, a final speed structure can be obtained by solving posterior probability distribution of the model parameters.
Determining a priori probability distribution of model parameters in bayesian inversion can be expressed as:
p(m)＝p(n)p(D|n)p(α 0 |n)p(β 0 |n)p(ε|n)p(γ|n)p(δ|n)
wherein p (n) represents the probability of the number of horizons; p (D|n) represents the probability distribution per layer depth in case the number of layers is n; p (alpha) 0 I n) represents the anisotropy parameter α per layer in case of the number of layers being n 0 Probability distribution of (2); p (beta) 0 I n) represents the anisotropy parameter β for each layer in the case where the number of layers is n 0 Probability distribution of (2); p (ε|n) represents the probability distribution of the anisotropy parameter ε for each layer with the number of layers n; p (γn) represents a probability distribution of the anisotropy parameter γ for each layer in the case where the number of layers is n; p (δ|n) represents the probability distribution of the anisotropy parameter δ per layer in the case of the number of layers n.
All prior information can be designed into uniform distribution, gaussian distribution, cauchy distribution and the like according to known data and used for solving the final posterior probability distribution.
Likelihood functions of model parameters in Bayesian inversion are determined. The likelihood function represents the distribution of noise in the observed data. In most cases, the likelihood function of the noise contribution can be expressed as a gaussian distribution:
in the method, in the process of the application,mean value of noise and sigma standard deviation of noise.
And combining the established likelihood function with prior probability distribution to construct a Bayesian formula posterior probability density distribution formula.
Iterative generation of posterior probability model by reversible jump Markov chain Monte Carlo algorithm, and random updating of velocity horizon number n, horizon depth D or dissimilarity parameter [ alpha ] in iterative process 0 ,β 0 ,ε,γ,δ]The four options of birth, death, movement and change are included. Wherein,,
"birth": randomly generating a new horizon based on the original horizon, wherein the interface depth obeys the uniform probability distribution:l is the total layer number, n is the current velocity model layer number, D new The velocity model horizon depth is newly generated.
"death": and randomly selecting one horizon of the existing speed model based on the original horizon, and deleting the horizon. The selection probability is as follows:D death the velocity model horizon is deleted.
"move": randomly selecting a horizon from the existing velocity models, randomly perturbing the depth of the horizon according to Gaussian probability, wherein the perturbation probability is as follows:D move for newly generated horizon depth σ 1 Is the standard deviation of the depth disturbance.
"change": randomly selecting a horizon from the existing velocity models, perturbing its anisotropic parameters according to gaussian probabilities, the perturbation probabilities being:v is a variable [ alpha ] 0 ,β 0 ,ε,γ,δ]K is an integer from 1 to 5, σ k Representing the standard deviation of five anisotropic parameters.
And judging whether the update speed model is accepted or not by using the acceptance probability. The probability of acceptance can be expressed as:
wherein m is old M is the original model new Is the updated model. p (m) old ) And p (d|m) old ) Is the prior information and likelihood function of the original model. p (m) new ) And p (d|m) new ) Is the updated model prior information and likelihood function. q (m) new |m old ) Generating transition probabilities of new models for original models, q (m old |m new ) The transition probabilities of the original model are generated for the new model.
The calculated acceptance probability accept_ratio is compared with the random number r between 0, 1. If accept_ratio > r, the updated model will be accepted; if accept_ratio < r, then the original model will go to the next cycle.
And performing posterior probability evaluation on model parameters obtained by Bayesian inversion. The specific process comprises the following steps: selecting a model with the maximum layer number occupation probability in the posterior probability as a final speed model; calculating the average value of each layer depth of the final speed model as the final horizon speed; calculating the anisotropy parameter [ alpha ] of each layer 0 ,β 0 ,ε,γ,δ]As the result of the anisotropic parameter inversion of each layer.
By using the steps, the number n of the layers of the VTI speed model, the depth D of the layers and the anisotropic parameters [ alpha ] of each layer can be finally determined 0 ,β 0 ,ε,γ,δ]。
Step 3: and positioning the natural seismic event by using a double-difference positioning method based on the updated VTI speed model.
Specifically, referring to fig. 4, fig. 4 is a schematic diagram of a dual differential positioning method for obtaining accurate positions of reservoir induced earthquakes. The method for positioning the natural seismic event by utilizing the double-difference positioning method based on the updated VTI speed model comprises the following specific steps:
picking up the arrival information of the induced seismic event, and establishing an objective function between the arrival information and the main event:
where ψ is the objective function. nr represents the number of DAS seismic traces, and N1 and N2 represent the number of main and evoked seismic events. t=o+t, T is the arrival time of DAS recordings, O is the start time of the seismic event, and T is the travel time from the source point to the seismic trace.And->The qP wave travel times of observed seismic traces r to event i and event j are represented. />And->The qSV wave travel times of observed traces r through events i and j are shown. />And->The qSH wave travel times of observed traces r through events i and j are shown. O (O) i And O j Is the actual start time of event i and event j. />And->The qP wave travel times of the detectors r to event i and event j obtained by forward calculation are represented. />And->The qSV wave times for the detectors r through event i and event j obtained by forward calculation are represented. />And->The qSH wave times for the detectors r through event i and event j obtained by forward calculation are represented. The obs and cal represent observed and calculated values.
Iteration through the gauss-newton algorithm, the above objective function linearization can be expressed as:
in the method, in the process of the application,representing the travel time difference of event i and event j to trace k. Deltaτ ij Indicating the difference in the starting moments of event i and event j. V is the background VTI velocity model. ΔH ij And DeltaZ ij Representing the relative distance between the horizontal and vertical positions of event i and event j.
According to the above calculation method, the relative position of the evoked seismic event with respect to the main seismic event can be finally obtained. The extensive situation of the underground fault zone can be finally described by a large number of induced seismic event positioning results, and some tiny hidden fault structures are found.
Step 4: the method comprises the steps of comprehensive seismic signal acquisition, velocity model correction and double-difference positioning, and the accurate position of the induced earthquake is measured.
In summary, according to the method, the amphibious distributed optical fiber acoustic wave sensing (Distributed fiber Acoustic Sensing, DAS) monitoring system is arranged in the multiple reservoir earthquake zones, so that meter-level space precision continuous observation is carried out on earthquake waves formed by near-surface earthquakes; picking up the relative arrival time of the DAS seismic signals by using a cross-correlation method, and constructing a time difference matrix; correcting the anisotropic VTI speed model by adopting a Bayesian theory and reversible jump Markov chain Monte Carlo (reversible jump Markov Chain Monte Carlo, rjMCMC) algorithm; and positioning the natural seismic event by using a double-difference positioning method based on the updated VTI speed model. The method fully utilizes the characteristics of good adaptability and high spatial sampling in the water of the DAS system, optimizes the velocity model by adopting continuous seismic wave field data, introduces rjMCMC theory in the process of correcting the velocity model to obtain a variable-dimension VTI (Vertical Transverse Isotropy) anisotropic velocity structure, and finally positions the seismic event by utilizing a double-difference positioning method, accurately depicts the spatial continuity of the seismic event and accurately depicts a hidden micro fault zone in a reservoir. The method has the advantages that the method can record the continuous wave field of the natural earthquake in the reservoir with high spatial resolution, simultaneously optimize and obtain the anisotropic VTI speed model, and combine with the double-difference positioning method with higher positioning precision, the accurate reservoir bottom earthquake event position can be obtained. The method is beneficial to improving the accuracy of the hidden fault zone characterization below the reservoir and provides a good technical means for reservoir safety assessment.
In the description of the present specification, a description referring to terms "one embodiment," "some embodiments," "examples," "specific examples," or "some examples," etc., means that a particular feature, structure, material, or characteristic described in connection with the embodiment or example is included in at least one embodiment or example of the present application. In this specification, schematic representations of the above terms are not necessarily directed to the same embodiment or example. Furthermore, the particular features, structures, materials, or characteristics described may be combined in any suitable manner in any one or more embodiments or examples. Furthermore, the different embodiments or examples described in this specification and the features of the different embodiments or examples may be combined and combined by those skilled in the art without contradiction.
The above embodiments are only for illustrating the technical concept and features of the present application, and are intended to enable those skilled in the art to understand the content of the present application and implement the same, and are not intended to limit the scope of the present application. All equivalent changes or modifications made in accordance with the essence of the present application are intended to be included within the scope of the present application.
Claims (10)
1. A real-time monitoring and positioning method for large reservoir induced earthquake is characterized by comprising the following steps:
an amphibious integrated distributed optical fiber acoustic wave sensing monitoring system is arranged in a multiple reservoir seismic zone, and meter-level space precision continuous observation and seismic signal acquisition are carried out on seismic waves formed by near-surface earthquakes;
picking up the relative arrival time of the DAS seismic signals by using a cross-correlation method, and constructing a time difference matrix; correcting the anisotropic VTI speed model by adopting a Bayesian theory and reversible jump Markov chain Monte Carlo algorithm;
positioning the natural seismic event by using a double-difference positioning method based on the updated VTI speed model;
the method comprises the steps of comprehensive seismic signal acquisition, velocity model correction and double-difference positioning, and the accurate position of the induced earthquake is measured.
2. The method for real-time monitoring and positioning of large reservoir induced earthquake according to claim 1, wherein the method for continuously observing the meter-scale space precision and acquiring the earthquake signal of the earthquake wave formed by the near-surface earthquake comprises the following specific steps: the method comprises the steps of (1) land laying communication optical cables in a reservoir earthquake multiple zone, and laying submarine cables at the bottom of the reservoir; the optical cable signal demodulator is arranged in a machine room and is connected with the land communication cable and the reservoir sea cable in series; analyzing the strain caused by the seismic waves sensed by the sea and land cable by using a signal demodulator; the space monitoring strain sensitivity is set to be about one meter, the time sampling rate is greater than 200Hz, and the continuous wave field of the earthquake transmitted to the surface is monitored without interval every day.
3. The method for real-time monitoring and positioning of large reservoir induced earthquake according to claim 1, wherein the method for picking up the relative arrival time of DAS seismic signals by using a cross-correlation method is constructed into a time difference matrix, and comprises the following specific steps:
and picking up the effective seismic event by using a long-short time window ratio method, wherein the long-short time window is an energy average value of a longer-time signal sampling length, and the short-time window is an energy average value of a shorter-time signal sampling length:
wherein S (N) is the signal amplitude obtained by DAS monitoring, N is the sampling number of long-time window signals, M is the sampling number of short-time window signals, N is more than M, and Ratio is the average energy Ratio of the short-time signals to the long-time signals obtained by recording;
when the signal is larger than the threshold value, the monitored signal can be judged to be caused by the earthquake event;
then picking up accurate relative earthquake arrival time of different optical fiber monitoring channels by using a cross-correlation method, and constructing a time difference matrix, wherein the formula of the cross-correlation method is as follows:
wherein x is 1 And x 2 Monitoring seismic traces for different DAS, R x1x2 And as for the similarity coefficient, accurately obtaining the similarity coefficient of the seismic event extracted by Ratio among different optical channels by using a cross-correlation formula, and selecting the time difference with the maximum similarity coefficient of the seismic channels to form a time difference matrix.
4. The method for real-time monitoring and positioning of large reservoir induced earthquakes according to claim 3, wherein the anisotropic VTI velocity model is corrected by using a markov chain monte carlo algorithm based on bayesian theory and reversible jump, and the specific steps include:
the seismic wave can be split into qP wave, qSV wave and qSH wave in VTI medium, and the VTI medium is composed of five anisotropic parameters [ alpha ] 0 ,β 0 ,ε,γ,δ]Determining; the speed model comprises the number n of horizons and the horizon depth D;
thus, the inversion VTI speed model parameters m are determined as: m= [ n, D, α 0 ,β 0 ,ε,γ,δ]
Wherein alpha is 0 And beta 0 Representing the velocity magnitude of each layer along the symmetry axis qP wave and qSV wave; epsilon represents the ratio of the velocity of each layer of qP wave in the horizontal and vertical directions; gamma represents the ratio of the velocity of each layer qSH wave in the horizontal and vertical directions; delta represents the rate of change of the qP wave of each layer in the vertical direction;
the Bayesian theory expression is determined as follows:
wherein d represents observation data, m is model parameters, p (m) is model priori information, p (d|m) is a likelihood function, p (m|d) is model posterior probability, and p (d) is overall probability of the observation data in a model space and is a constant;
and obtaining a final speed structure by solving posterior probability distribution of model parameters by using a Bayesian inversion algorithm and combining prior information and likelihood functions of the model parameters.
5. The method for real-time monitoring and positioning of large reservoir induced earthquakes according to claim 4, wherein the anisotropic VTI velocity model is corrected by using a markov chain monte carlo algorithm based on bayesian theory and reversible jump, and the specific steps include:
determining the prior probability distribution of model parameters in Bayesian inversion, wherein the prior distribution formula is expressed as follows:
p(m)＝p(n)p(D|n)p(α 0 |n)p(β 0 |n)p(ε|n)p(γ|n)p(δ|n)
wherein p (n) represents the probability of the number of horizons; p (D|n) represents the probability distribution per layer depth in case the number of layers is n; p (alpha) 0 I n) represents the anisotropy parameter α per layer in case of the number of layers being n 0 Probability distribution of (2); p (beta) 0 I n) represents the anisotropy parameter β for each layer in the case where the number of layers is n 0 Probability distribution of (2); p (ε|n) represents the probability distribution of the anisotropy parameter ε for each layer with the number of layers n; p (γn) represents a probability distribution of the anisotropy parameter γ for each layer in the case where the number of layers is n; p (δ|n) represents the probability distribution of the anisotropy parameter δ per layer with the number of layers n;
all prior information is designed into uniform distribution, gaussian distribution and Cauchy distribution and is used for solving the final posterior probability distribution.
6. The method for real-time monitoring and positioning of large reservoir induced earthquakes according to claim 5, wherein the anisotropic VTI velocity model is corrected by using a markov chain monte carlo algorithm based on bayesian theory and reversible jump, and the specific steps include:
determining likelihood functions of model parameters in Bayesian inversion, wherein the likelihood functions represent distribution conditions of noise in observed data, and in most cases, likelihood functions formed by the noise are expressed as Gaussian distribution:
mean value of noise, sigma standard deviation of noise;
and establishing a likelihood function and combining the prior probability distribution to construct a Bayesian formula posterior probability density distribution formula.
7. The method for monitoring and positioning the induced earthquake of the large reservoir in real time according to claim 1, wherein the anisotropic VTI velocity model is corrected by using a markov chain monte carlo algorithm based on the bayesian theory and the reversible jump, and the specific steps include:
iterative generation of posterior probability model by reversible jump Markov chain Monte Carlo algorithm, and random updating of velocity horizon number n, horizon depth D or anisotropic parameter [ alpha ] during iteration 0 ,β 0 ,ε,γ,δ]Including birth, death, movement, and alteration, among others,
"birth": randomly generating a new horizon based on the original horizon, wherein the interface depth obeys the uniform probability distribution:l is the total layer number, n is the current velocity model layer number, D new Generating a velocity model horizon depth for the new generation;
"death": randomly selecting an existing velocity model based on an original horizonA horizon of the model is deleted, and the selection probability is as follows:D death a velocity model horizon for deletion;
"move": randomly selecting a horizon from the existing velocity models, randomly perturbing the depth of the horizon according to Gaussian probability, wherein the perturbation probability is as follows:D move for newly generated horizon depth σ 1 Is the standard deviation of depth disturbance;
"change": randomly selecting a horizon from the existing velocity models, perturbing its anisotropic parameters according to gaussian probabilities, the perturbation probabilities being:v is an anisotropic parameter [ alpha ] 0 ,β 0 ,ε,γ,δ]K is an integer from 1 to 5, σ k Representing the standard deviation of five anisotropic parameters.
8. The method for monitoring and positioning the induced earthquake of the large reservoir in real time according to claim 1, wherein the anisotropic VTI velocity model is corrected by using a markov chain monte carlo algorithm based on the bayesian theory and the reversible jump, and the specific steps include:
judging whether the update speed model is accepted or not by using an acceptance probability, wherein the acceptance probability is expressed as:
wherein m is old M is the original model new For the updated model, p (m old ) And p (d|m) old ) Is the prior information and likelihood function of the original model, p (m new ) And p (d|m) new ) Is updated model prior informationAnd likelihood function, q (m new |m old ) Generating transition probabilities of new models for original models, q (m old |m new ) Generating a transition probability of the original model for the new model;
calculating a random number r between the acceptance probability accept_ratio and [0,1], and if the acceptance_ratio > r, then the updated model will be accepted; if accept_ratio < r, then the original model will go to the next cycle.
9. The method for monitoring and positioning the induced earthquake of the large reservoir in real time according to claim 1, wherein the anisotropic VTI velocity model is corrected by using a markov chain monte carlo algorithm based on the bayesian theory and the reversible jump, and the specific steps include:
the posterior probability evaluation is carried out on model parameters obtained by Bayesian inversion, and the flow comprises the following steps:
(1) Selecting a model with the maximum layer number occupation probability in the posterior probability as a final speed model;
(2) Calculating the average value of each layer depth of the final speed model as the final horizon depth;
(3) Calculating the anisotropy parameter [ alpha ] of each layer 0 ,β 0 ,ε,γ,δ]The average value of (2) is used as the anisotropic parameter inversion result of each layer;
by the steps of finally determining the number n of layers, the depth D of the layers and the anisotropic parameters [ alpha ] of each layer of the VTI speed model 0 ,β 0 ,ε,γ,δ]。
10. The method for real-time monitoring and positioning of large reservoir induced earthquakes according to claim 1, wherein the method for positioning natural earthquake events by using double difference positioning based on the updated VTI velocity model comprises the following specific steps:
picking up the arrival information of the seismic event, and establishing an objective function between the arrival information and the main event:
wherein, ψ is an objective function, nr represents the number of DAS seismic traces, N1 and N2 represent the numbers of main seismic events and induced seismic events, t=O+T, T is the arrival time of DAS records, O is the starting time of seismic events, T is the travel time from a seismic source point to a seismic trace,and->qP wave travel time indicating observed trace r to event i and event j, +.>And->qSV wave travel time, representing observed trace r to event i and event j, +.>And->qSH wave travel time, O, representing observed seismic trace r to event i and event j i And O j Is the actual start time of event i and event j, < +.>And->qP wave travel time representing forward calculated detector r to event i and event j, +.>And->qSV wave-away time representing forward calculated detector r to event i and event j, +.>And->qSH wave times representing forward calculated detector r to event i and event j, obs and cal represent observed and calculated values;
iteration is performed by means of a gauss-newton algorithm, the above linearization of the objective function being expressed as:
in the method, in the process of the application,indicating the travel time difference between event i and event j to trace k, Δτ ij Representing the difference in the starting moments of event i and event j, V is the background VTI speed model, ΔH ij And DeltaZ ij Representing the relative distance between the horizontal position and the vertical position of the event i and the event j;
according to the calculation method, the relative position of the induced seismic event relative to the main seismic event is finally obtained, a large number of induced seismic event positioning results can finally describe the extension condition of the underground fault zone, and some tiny hidden fault structures are found.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310129707.1A CN116466390B (en) | 2023-02-17 | 2023-02-17 | Real-time monitoring and positioning method for earthquake induced by large reservoir |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310129707.1A CN116466390B (en) | 2023-02-17 | 2023-02-17 | Real-time monitoring and positioning method for earthquake induced by large reservoir |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN116466390A CN116466390A (en) | 2023-07-21 |
| CN116466390B true CN116466390B (en) | 2023-11-03 |
Family
ID=87176044
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310129707.1A Active CN116466390B (en) | 2023-02-17 | 2023-02-17 | Real-time monitoring and positioning method for earthquake induced by large reservoir |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116466390B (en) |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109085642A (en) * | 2017-06-14 | 2018-12-25 | 中国石油化工股份有限公司 | A kind of anisotropic medium micro-seismic event localization method |
| CN111722280A (en) * | 2020-06-29 | 2020-09-29 | 重庆大学 | Acoustic emission event Bayes positioning method, system and medium for removing P wave first-motion system observation error |
| CN112904419A (en) * | 2021-01-26 | 2021-06-04 | 南方科技大学 | Microseism imaging method and terminal equipment |
| CN113568037A (en) * | 2021-08-17 | 2021-10-29 | 中油奥博(成都)科技有限公司 | Earthquake and geological disaster monitoring system and method based on optical fiber sensing technology |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US10634803B2 (en) * | 2015-09-16 | 2020-04-28 | Schlumberger Technology Corporation | Bayseian microseismic source inversion |
| US20230099919A1 (en) * | 2020-03-27 | 2023-03-30 | Chevron U.S.A. Inc. | System and method for stochastic full waveform inversion |
-
2023
- 2023-02-17 CN CN202310129707.1A patent/CN116466390B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109085642A (en) * | 2017-06-14 | 2018-12-25 | 中国石油化工股份有限公司 | A kind of anisotropic medium micro-seismic event localization method |
| CN111722280A (en) * | 2020-06-29 | 2020-09-29 | 重庆大学 | Acoustic emission event Bayes positioning method, system and medium for removing P wave first-motion system observation error |
| CN112904419A (en) * | 2021-01-26 | 2021-06-04 | 南方科技大学 | Microseism imaging method and terminal equipment |
| CN113568037A (en) * | 2021-08-17 | 2021-10-29 | 中油奥博(成都)科技有限公司 | Earthquake and geological disaster monitoring system and method based on optical fiber sensing technology |
Also Published As
| Publication number | Publication date |
|---|---|
| CN116466390A (en) | 2023-07-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Xie et al. | Preliminary analysis on the source properties and seismogenic structure of the 2017 M S 7.0 Jiuzhaigou earthquake | |
| Zollo et al. | Bayesian estimation of 2-DP-velocity models from active seismic arrival time data: imaging of the shallow structure of Mt Vesuvius (Southern Italy) | |
| Coccia et al. | Application of Refraction Microtremor (ReMi) technique for determination of 1-D shear wave velocity in a landslide area | |
| WO2008124759A1 (en) | Location of a microseismic event using the combination of a vsp array and a surface array | |
| CN102901985B (en) | A kind of Depth Domain interval velocity modification method being applicable to relief surface | |
| CN105277982A (en) | Shale total organic carbon content earthquake prediction method | |
| Ugalde et al. | Passive seismic monitoring of an experimental CO2 geological storage site in Hontomín (Northern Spain) | |
| CN109884710A (en) | For the micro logging chromatography imaging method of excitation well depth design | |
| CN103913768A (en) | Method and device for modeling superficial layer in earth surface based on seismic wave data | |
| Cercato et al. | Shear‐wave velocity profiling at sites with high stiffness contrasts: a comparison between invasive and non‐invasive methods | |
| CN111352160B (en) | Automatic repositioning device and method for ocean bottom seismograph | |
| CN104749630B (en) | Method for constructing microseism monitoring velocity model | |
| CN113466933B (en) | Depth weighting-based seismic slope tomography method | |
| CN116466390B (en) | Real-time monitoring and positioning method for earthquake induced by large reservoir | |
| CN114861515A (en) | Method, device, equipment and medium for calculating layer speed data volume | |
| Sauvin et al. | Towards joint inversion/interpretation for landslide-prone areas in Norway-integrating geophysics and geotechnique | |
| CN113589375B (en) | VSP layer speed inversion method based on calculation during constraint travel of inclined layer | |
| CN116500673B (en) | Artificial island micro-seismic monitoring method, equipment and medium based on distributed optical fiber acoustic wave sensing | |
| CN110780345A (en) | Three-dimensional velocity analysis method for tunnel advanced seismic exploration seismic data | |
| CN111310361A (en) | Drilling guidance method, system, equipment and storage medium based on earthquake while drilling | |
| Yang et al. | Revealing the shallow soil structure of the Yigong Lake in the Tibetan Plateau using a portable distributed acoustic sensing interrogator | |
| Guo et al. | Using Wide-Angle Reflection Wave Technology to Detect Seismic Inversion Data of Complex Geological Structure Zone for a Future Smart World | |
| Lizhi et al. | The seismic CT method in measuring rock bodies | |
| Koedel et al. | Project report on seismic tomography data interpretation and conceptual model for integrating DAS into borehole seismic tomography surveying | |
| CN118135135A (en) | Landslide hazard prediction method based on seismic noise |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |