CN116466390A - 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 PDF

Info

Publication number
CN116466390A
CN116466390A CN202310129707.1A CN202310129707A CN116466390A CN 116466390 A CN116466390 A CN 116466390A CN 202310129707 A CN202310129707 A CN 202310129707A CN 116466390 A CN116466390 A CN 116466390A
Authority
CN
China
Prior art keywords
seismic
event
model
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.)
Granted
Application number
CN202310129707.1A
Other languages
Chinese (zh)
Other versions
CN116466390B (en
Inventor
蒋星达
杨华勇
李超
杨得厚
李跃金
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Southern Marine Science and Engineering Guangdong Laboratory Guangzhou
Original Assignee
Southern Marine Science and Engineering Guangdong Laboratory Guangzhou
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Southern Marine Science and Engineering Guangdong Laboratory Guangzhou filed Critical Southern Marine Science and Engineering Guangdong Laboratory Guangzhou
Priority to CN202310129707.1A priority Critical patent/CN116466390B/en
Publication of CN116466390A publication Critical patent/CN116466390A/en
Application granted granted Critical
Publication of CN116466390B publication Critical patent/CN116466390B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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 invention 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 invention aims to solve the problem of low positioning precision of real-time monitoring of induced seismic events in reservoirs.

Description

Real-time monitoring and positioning method for earthquake induced by large reservoir
Technical Field
The invention 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 invention 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 invention 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 invention has the beneficial effects that: the invention 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 invention, 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 that 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 invention;
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 invention;
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 invention 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 one of ordinary skill in the art without undue burden from the present disclosure, are within the scope of the present disclosure.
Examples:
it should be noted that the terms "first," "second," and the like in the description and the claims of the present invention 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 invention 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 invention, 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 invention 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 invention 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 invention 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 invention,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, the liquid crystal display device comprises a liquid crystal display device,
"birth": randomly generating a new horizon based on the original horizon, the boundary thereofThe face depth obeys a 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 an updated modelType a priori information and likelihood functions. 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->Representing forward playThe qSH wave times of the obtained detectors r to event i and event j are calculated. 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 invention,representing the travel time difference of event i and event j to trace k. Δτ 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 invention. 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 invention, and are intended to enable those skilled in the art to understand the content of the present invention and implement the same, and are not intended to limit the scope of the present invention. All equivalent changes or modifications made in accordance with the essence of the present invention are intended to be included within the scope of the present invention.

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 The seismic traces are monitored for different DASs,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:
where d represents the observed data, m is the model parameters, p (m) is the model prior information, p (d|m) is the likelihood function, p (m|d) is the model posterior probability, and p (d) is the overall probability of the observed data in the model space, which 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 a layerProbability of number of bits; 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, cauchy distribution and the like 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 dissimilarity parameter [ alpha ] in iterative process 0 ,β 0 ,ε,γ,δ]Includes four choices of birth, death, movement and change, among which,
"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 a horizon of the existing speed model based on the original horizon, deleting the horizon, and selecting probabilities 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 a variable [ 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 information and 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 invention,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.
CN202310129707.1A 2023-02-17 2023-02-17 Real-time monitoring and positioning method for earthquake induced by large reservoir Active CN116466390B (en)

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 true CN116466390A (en) 2023-07-21
CN116466390B 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 (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20170074997A1 (en) * 2015-09-16 2017-03-16 Schlumberger Technology Corporation Bayseian microseismic source inversion
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
US20230099919A1 (en) * 2020-03-27 2023-03-30 Chevron U.S.A. Inc. System and method for stochastic full waveform inversion

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20170074997A1 (en) * 2015-09-16 2017-03-16 Schlumberger Technology Corporation Bayseian microseismic source inversion
CN109085642A (en) * 2017-06-14 2018-12-25 中国石油化工股份有限公司 A kind of anisotropic medium micro-seismic event localization method
US20230099919A1 (en) * 2020-03-27 2023-03-30 Chevron U.S.A. Inc. System and method for stochastic full waveform inversion
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
CN116466390B (en) 2023-11-03

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
Coccia et al. Application of Refraction Microtremor (ReMi) technique for determination of 1-D shear wave velocity in a landslide area
US20080247269A1 (en) Analysis of Uncertainty of Hypocenter Location Using the Combination of a VSP and a Subsurface Array
CN109884710B (en) Micro-logging tomography method aiming at excitation well depth design
CN103109207A (en) Method for detection of subsurface seismic events in vertically transversely isotropic media
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)
CN105719433A (en) In-hole seismic wave based advanced prediction method
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
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
Zujun et al. Preliminary analysis on the source properties and seismogenic structure of the 2017 M s 7.0 Jiuzhaigou earthquake.
CN113589375B (en) VSP layer speed inversion method based on calculation during constraint travel of inclined layer
CN111352160B (en) Automatic repositioning device and method for ocean bottom seismograph
CN116009111A (en) High-density electrical method and micro-motion combined near-surface water-rich region exploration method
CN110780345A (en) Three-dimensional velocity analysis method for tunnel advanced seismic exploration seismic data
CN116500673B (en) Artificial island micro-seismic monitoring method, equipment and medium based on distributed optical fiber acoustic wave sensing
CN111310361A (en) Drilling guidance method, system, equipment and storage medium based on earthquake while drilling
CN111126793A (en) Landslide risk assessment method based on ultralow frequency electromagnetic wave
CN109901221A (en) A kind of seismic data anisotropy modeling method based on NMO velocity parameter
Koedel et al. Project report on seismic tomography data interpretation and conceptual model for integrating DAS into borehole seismic tomography surveying
Guo et al. Using Wide-Angle Reflection Wave Technology to Detect Seismic Inversion Data of Complex Geological Structure Zone for a Future Smart World
Bakulin et al. Advances in near-surface characterization and deep imaging with smart DAS upholes

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