CN105652319A - Estimation method of near-surface stratum Q value of complex mediums - Google Patents

Abstract

The invention relates to an estimation method of a near-surface stratum Q value of complex mediums. The estimation method comprises the steps of: 1) collecting an actual measurement Rayleigh surface wave single shot seismic record of a target area, processing the actual measurement Rayleigh surface wave single shot seismic record with a multichannel surface wave analysis method, and obtaining that a near-surface stratum has N kinds of different mediums; 2) establishing an initial vector of the near-surface stratum Q value, and synthesizing a Rayleigh surface wave single shot seismic record with a seismic wave forward modeling tool; 3) carrying out mode separation on the actual measurement Rayleigh surface wave single shot seismic record and the synthesized Rayleigh surface wave single shot seismic record by means of high resolution Radon transformation, obtaining the actual measurement Rayleigh surface wave single shot seismic records and the synthesized Rayleigh surface wave single shot seismic records of a fundamental mode and various high-order modes, calculating initial target function values of the fundamental mode and the various high-order modes, carrying out weighted addition on the initial target function values of the fundamental mode and the various high-order modes, and obtaining a total initial target function value; and 4) utilizing an simulated annealing method to invert the near-surface stratum Q value.

Description

A kind of method of estimation of complex dielectrics near surface stratum Q-value

Technical field

The present invention relates to a kind of method of seismic prospecting, particularly relate to the method for estimation of a kind of complex dielectrics near surface stratum Q-value.

Background technology

Method of seismic prospecting is the highly important method of class finding underground oil and gas reservoir, but in oil-gas exploration, owing to surface conditions has multiformity, there are the complex situations such as open structure or low velocity layer region near surface, the energy causing seismic wave is significantly absorbed or is decayed, and reduces the resolution of seismic data. The Q-value near surface stratum represents that seismic wave stores the ratio of energy and consume energy, reflects seismic wave and often propagates the degree of energy loss during the distance of a wavelength. Estimate near surface stratum Q-value, it is possible to for inverse Q filtering, thus compensate the absorption effect on stratum, correction wavelet stretches mutually, improve lineups seriality and improve seismic data resolution. Rayleigh (Rayleigh) face ripple can provide the information such as shear wave velocity and the Q-value of near surface, Xia et al. proposed the method utilizing Rayleigh face ripple to estimate near surface stratum Q-value in 2002, but the method based on underground be layered medium it is assumed that the situation that near surface stratum is complex dielectrics cannot be adapted to.

Summary of the invention

For the problems referred to above, it is an object of the invention to provide the method for estimation of a kind of complex dielectrics near surface stratum Q-value, based on a kind of new object function of Rayleigh face wave mode separation structure, simulated annealing inverting is adopted to estimate near surface stratum Q-value, need not suppose that underground medium is layered medium, solve Q-value inestimable problem near surface stratum in complex dielectrics situation.

For achieving the above object, the present invention takes techniques below scheme: the method for estimation of a kind of complex dielectrics near surface stratum Q-value, and it comprises the following steps:

1) carry out seismic experiment in target area, gathered the actual measurement Rayleigh face ripple list big gun earthquake record of target area by multiple cymoscopes, and record the distance between each cymoscope and focus and seismic wave and propagate the time used by each cymoscope from focus;Adopt wave analysis method in multiple tracks face to process actual measurement Rayleigh face ripple list big gun earthquake record, it is thus achieved that the velocity of longitudinal wave near surface stratum, shear wave velocity and Media density, show that near surface stratum has the medium that N kind is different;

2) set the initial value of a Q-value to every kind of medium, the initial vector setting up near surface stratum Q-value isWherein,Represent the initial value of the Q-value of jth kind medium, 1��j��N; And set the initial temperature T for simulated annealing₀; Initial vector Q according near surface stratum Q-value⁰, use seismic forward modeling simulation tool to simulate the synthesis Rayleigh face ripple list big gun earthquake record being made up of multiple tracks face ripple record;

3) use high-resolution radon transform respectively actual measurement Rayleigh face ripple list big gun earthquake record and synthesis Rayleigh face ripple list big gun earthquake record to be carried out modal cutoff, obtain actual measurement Rayleigh face ripple list big gun earthquake record and the synthesis Rayleigh face ripple list big gun earthquake record of base order mode and each higher order mode; Actual measurement Rayleigh face ripple list big gun earthquake record according to base order mode and each higher order mode and synthesis Rayleigh face ripple list big gun earthquake record, calculate the initial target functional value of base order mode and each higher order mode respectively; By the initial target functional value weighting summation of base order mode and each higher order mode, obtain total initial target functional value;

4) use simulated annealing inverting, obtain final near surface stratum Q-value.

Described step 3) in use high-resolution radon transform that actual measurement Rayleigh face ripple list big gun earthquake record or synthesis Rayleigh face ripple list big gun earthquake record are carried out modal cutoff, specifically include following steps:

1. use high-resolution radon transform that from time and space territory, actual measurement Rayleigh face ripple list big gun earthquake record or synthesis Rayleigh face ripple list big gun earthquake record are transformed to frequency slowness domain, obtain actual measurement Rayleigh face ripple list big gun earthquake record or the frequency-slowness section of synthesis Rayleigh face ripple list big gun earthquake record, actual measurement Rayleigh face ripple list big gun earthquake record or the synthesis Rayleigh face ripple list big gun earthquake record of base order mode and each higher order mode are respectively provided with scope clearly in frequency-slowness section, and not overlapping each other;

2. in frequency slowness section, the scope of the actual measurement Rayleigh face ripple list big gun earthquake record or synthesis Rayleigh face ripple list big gun earthquake record that represent base order mode and each higher order mode is determined respectively, the value of actual measurement Rayleigh face ripple list big gun earthquake record or synthesis Rayleigh face ripple list big gun earthquake record by removing base order mode or a certain higher order mode other order mode extraneous in frequency slowness section, obtain the frequency slowness section only comprising the actual measurement Rayleigh face ripple list big gun earthquake record of base order mode or a certain higher order mode or synthesis Rayleigh face ripple list big gun earthquake record respectively,

3. adopt and draw east inverse transformation that from frequency-slowness domain, the actual measurement Rayleigh face ripple list big gun earthquake record of base order mode and each higher order mode or synthesis Rayleigh face ripple list big gun earthquake record are switched back to time and space territory, obtain actual measurement Rayleigh face ripple list big gun earthquake record or the synthesis Rayleigh face ripple list big gun earthquake record of base order mode and each higher order mode separated.

Described step 3) in calculate target function value according to the actual measurement Rayleigh face ripple list big gun earthquake record of base order mode and each higher order mode and synthesis Rayleigh face ripple list big gun earthquake record, specifically include following steps:

A, in the actual measurement Rayleigh face ripple list big gun earthquake record and synthesis Rayleigh face ripple list big gun earthquake record of base order mode and each higher order mode the spread speed of a selected Rayleigh face ripple respectively, the actual measurement Rayleigh face ripple list big gun earthquake record and synthesis Rayleigh face ripple list big gun earthquake record of base order mode and each higher order mode select lineups respectively clearly, and picks up amplitude respectively along respective lineups in the actual measurement Rayleigh face ripple list big gun earthquake record and synthesis Rayleigh face ripple list big gun earthquake record of base order mode and each higher order mode;

B, amplitude according to the actual measurement Rayleigh face ripple list big gun earthquake record of base order mode and each higher order mode and synthesis Rayleigh face ripple list big gun earthquake record, calculate the target function value of base order mode and each higher order mode respectively;

C, by the target function value weighting summation of base order mode and each higher order mode, obtain total target function value.

In described step b, the computing formula of target function value is:

$E\left(Q\right)=\underset{t=1}{\overset{n}{\Σ}}\left|log\left(\frac{A^{obs}\left(T\right(H_t))}{A^{syn}\left(T\right(H_t))}\right)\right|$

In formula, E (Q) represents the target function value that Q-value is corresponding; N represents the number of channels of the face ripple record comprised in actual measurement Rayleigh face ripple list big gun earthquake record or synthesis Rayleigh face ripple list big gun earthquake record, t represents t road actual measurement Rayleigh face ripple list big gun earthquake record or the number of channels of synthesis Rayleigh face ripple list big gun earthquake record, 1��t��n; H_tRepresent the distance between cymoscope and the focus of actual measurement Rayleigh face, record t road ripple list big gun earthquake tracer signal; T (H_t) represent that seismic wave propagates the time used by cymoscope recording actual measurement Rayleigh face, t road ripple list big gun earthquake tracer signal from focus; A^obs(T(H_t)) and A^syn(T(H_t)) represent T (H in actual measurement Rayleigh face ripple list big gun earthquake record and synthesis Rayleigh face ripple list big gun earthquake record respectively_t) moment pickup amplitude;

Or the computing formula of target function value is:

$E(Q,W_t)=\underset{t=1}{\overset{n}{\Σ}}\left|log\left(\frac{{\mathrm{\ϵ}}^{obs}\left(T\right(H_t\left)\right)}{{\mathrm{\ϵ}}^{syn}\left(T\right(H_t\left)\right)}\right)\right|$

In formula, W_tRepresent with T (H_t) centered by the time window that takes; ��^obs(T(H_t)) and ��^syn(T(H_t)) represent that actual measurement Rayleigh face ripple list big gun earthquake record and synthesis Rayleigh face ripple list big gun earthquake record are at time window W respectively_tInterior amplitude energy.

Described step 4) middle use simulated annealing inverting near surface stratum Q-value, specifically include following steps:

The first step, set up the iterative vectorized of near surface stratum Q-valueWherein,Represent the Q-value of jth kind medium, 1��j��N, i >=1 during ith iteration; Make i=1;

Second step, make Qⁱ=Q^i-1, iterative target function E (Qⁱ)=E (Q^i-1), iteration temperature T_i=�� T_i-1; Wherein, �� is constant, takes ��=0.99;

3rd step, generate the comparison vector of near surface stratum Q-valueWherein,The comparison Q-value of the jth kind medium generated when representing ith iteration, 1��j��N, i < 1; Comparison vector according near surface stratum Q-valueUse seismic forward modeling simulation tool simulation comparison Rayleigh face ripple list big gun earthquake record; Use high-resolution radon transform respectively actual measurement Rayleigh face ripple list big gun earthquake record and comparison Rayleigh face ripple list big gun earthquake record to be carried out modal cutoff, obtain actual measurement Rayleigh face ripple list big gun earthquake record and the comparison Rayleigh face ripple list big gun earthquake record of base order mode and each higher order mode; Actual measurement Rayleigh face ripple list big gun earthquake record according to base order mode and each higher order mode and comparison Rayleigh face ripple list big gun earthquake record, the amplitude surveying Rayleigh face ripple list big gun earthquake record and comparison Rayleigh face ripple list big gun earthquake record of pickup base order mode and each higher order mode, calculates the comparison target function value of base order mode and each higher order mode respectively; By the comparison target function value weighting summation of base order mode and each higher order mode, obtain total comparison target function value

4th step, calculating $\mathrm{\&Delta;}E=E\left({\stackrel{\&OverBar;}{Q}}^i\right)-E\left(Q^i\right);$

If �� E��0, accept the comparison vector generated, take

If �� E > 0, calculate�� is intermediate computations parameter;

If �� (�� E) > R, accept the comparison vector generated, takeR is equally distributed random number between (0,1);

Otherwise do not accept the comparison vector generated, Q^uIt is worth constant;

If the comparison vector that the 5th step generates for continuous 20 timesDo not accepted, i.e. continuous 20 QⁱBeing worth constant, terminate inverting, output final result is Q^est=Qⁱ;Otherwise, i is made_new=i+1, returns second step.

Described 3rd step adopts random method generate comparison vector, or adopt the Ingber Fast Simulated Annealing proposed to generate comparison vector; The formula that quickly generates of described comparison vector is:

${\stackrel{\&OverBar;}{Q}}^i=Q^i+y(Q_{max}-Q_{min})$

Y=T_isgn(u-0.5)[(1+1/T_i)^|2u-1|-1]

In formula, y is intermediate computations parameter; Q_maxAnd Q_minThe respectively upper limit value and lower limit value of the Q-value span of each medium, Q_min=10, Q_max=200; U is random number, is evenly distributed on [0,1] interval; Sgn is for taking sign function.

Described step 2) in seismic forward modeling simulation tool be finite difference calculus or FInite Element.

Described step 4) in the algorithm that realizes of simulated annealing that adopts be Metropolis algorithm.

Due to the fact that and take above technical scheme, it has the advantage that 1, the method for estimation of the complex dielectrics near surface stratum Q-value of the present invention, need not supposing that underground medium is layered medium, solve Q-value inestimable problem near surface stratum in complex dielectrics situation, practicality is stronger. 2, the method for estimation of the complex dielectrics near surface stratum Q-value of the present invention, it is possible to be widely used in the estimation of various complex dielectrics near surface stratum Q-value, there is important practical significance.

Accompanying drawing explanation

Fig. 1 is the flow chart of the inventive method;

Fig. 2 is Rayleigh face ripple list big gun earthquake record;

Fig. 3 is the result schematic diagram that base order mode separates and amplitude information extracts of face ripple;

Fig. 4 is the result schematic diagram that higher order mode separates and amplitude information extracts of face ripple.

Detailed description of the invention

Below in conjunction with drawings and Examples, the present invention is described in detail.

As it is shown in figure 1, the method for estimation of a kind of complex dielectrics near surface stratum Q-value of the present invention, comprise the following steps:

1) seismic experiment is carried out in target area, gathered the actual measurement Rayleigh face ripple list big gun earthquake record (be called for short effective surface ripple record) of target area by cymoscope, and record the distance between cymoscope and focus and seismic wave and propagate the relevant parameters such as time used by cymoscope from focus. Wherein it is possible to arrange multiple cymoscope, effective surface ripple record is made up of multiple tracks face ripple record, and every one face ripple record is the time signal of a cymoscope record. Wave analysis method in multiple tracks face is adopted to process effective surface ripple record, it is thus achieved that the parameters such as the density of the velocity of longitudinal wave near surface stratum, shear wave velocity and medium, such that it is able to show that near surface stratum has the medium that N kind is different.

Seismic wave includes the polytypes such as face ripple, shear wave and compressional wave, and near surface stratum media is also different to the Q-value of dissimilar seismic wave, and the present invention to estimate is the medium Q-value to shear wave and compressional wave; Present invention assumes that and every kind of medium all has the Q-value to compressional wave and shear wave equal, be therefore collectively referred to as the Q-value of medium.

The present invention has n cymoscope (namely having face, n road ripple record); T represents the number of channels of face ripple record, 1��t��n; H_tRepresent the distance between cymoscope and the focus of record face, t road ripple tracer signal; T (H_t) represent that seismic wave propagates the time used by the cymoscope recording face, t road ripple tracer signal from focus.

2) initial value of a Q-value is set to every kind of medium, thus the initial vector setting up near surface stratum Q-value isWherein,Represent the initial value of the Q-value of jth kind medium, 1��j��N; And set the initial temperature T for simulated annealing₀. Initial vector Q according near surface stratum Q-value⁰, use seismic forward modeling simulation tool to simulate the Rayleigh face ripple list big gun earthquake record (being called for short synthesis face ripple record) being made up of multiple tracks face ripple record.

The Q-value of practical application medium has a general span. The minima being currently known the Q-value estimating the medium obtained is about 30, it is believed that Q-value will not less than 10; Meanwhile, the Q-value more than 200 is only small for the impact of seismic wave, such as when some required precisions are not high, Q=200 and Q=150 on the impact of seismic wave it is believed that difference is little. Therefore, inventive method assumes that the Q-value of each medium is all [Q_min, Q_max] integer value in interval, and in this interval, set the initial value of the Q-value of each medium; Wherein, Q_min=10, Q_max=200. Solving that simulated annealing is tried to achieve is unrelated with the initial value of the Medium Q value chosen, and the initial value of the Q-value of each medium can be arbitrarily designated in above-mentioned interval. Initial temperature T₀It is the parameter in simulated annealing, initial temperature T₀More big, it is thus achieved that the probability of high-quality solution is more big, but the calculating time of cost is by increase; Initial temperature T₀Determination need through repeatedly debugging obtain.

3) high-resolution radon transform is used respectively effective surface ripple record and synthesis face ripple record to be carried out modal cutoff, obtain effective surface ripple record and the synthesis face ripple record of base order mode, and the effective surface ripple record of each higher order mode and synthesis face ripple record, here, higher order mode can have the first high-order, the second high-order even pattern of higher order; Effective surface ripple record according to base order mode and each higher order mode and synthesis face ripple record, calculate the initial target functional value of base order mode and the initial target functional value of each higher order mode respectively; By the initial target functional value weighting summation of base order mode and each higher order mode, obtain total initial target functional value E (Q⁰)��

Face ripple record is a fluctuation signal, has amplitude and phase place, and lineups are the lines of the extreme value (being commonly called as crest or trough) that face, each road ripple vibration phase of face ripple record is identical. As shown in Figure 2, effective surface ripple record and synthesis face ripple record all include the record of the face ripple of various modes, the face ripple record of different mode has different spread speeds, mutually overlap in time and space territory, cannot differentiate, it is difficult to determine the spread speed of Rayleigh face ripple in effective surface ripple record and synthesis face ripple record, causes can not find suitable lineups to pick up the amplitude of Rayleigh face ripple in effective surface ripple record and synthesis face ripple record. Utilize the propagation characteristic of Rayleigh face ripple, use high-resolution radon transform (Luo, 2009) Rayleigh face ripple is carried out modal cutoff, after separation, the various pattern faces ripple record of Rayleigh face ripple there are lineups clearly, therefore the spread speed of a Rayleigh face ripple can be determined in the face ripple record of each pattern, choose lineups clearly, pick up the amplitude of this pattern Rayleigh face ripple along these lineups, finally the face ripple record of each pattern can calculate and obtain a target function value. By the target function value weighting summation of all pattern faces ripple record, thus obtaining the target function value that Rayleigh face ripple is total.

As shown in Figure 3, Figure 4, the method using high-resolution radon transform that effective surface ripple record or synthesis face ripple record are carried out modal cutoff, comprise the following steps:

1. high-resolution radon transform (Luo is used, 2009) effective surface ripple record or synthesis face ripple record are transformed to frequency slowness domain from time and space territory, obtain effective surface ripple record or the frequency-slowness section of synthesis face ripple record, effective surface ripple record or the synthesis face ripple record of each order mode are respectively provided with scope clearly in frequency-slowness section, and not overlapping each other;Wherein, frequency-slowness domain is mathematical description, for the parameter space of variable in representation formula, is an abstract conception, and frequency-slowness section refers to the coordinate in length and breadth during the display of calculated result.

2. in frequency slowness section, the scope of the effective surface ripple record or synthesis face ripple record that represent base order mode and each higher order mode is determined respectively, the value of effective surface ripple record or synthesis face ripple record by removing base order mode or a certain higher order mode other order mode extraneous in frequency slowness section, obtains the frequency slowness section only comprising the effective surface ripple record of base order mode or a certain higher order mode or synthesis face ripple record respectively.

3. adopt and draw east inverse transformation that from frequency-slowness domain, the effective surface ripple record of each order mode or synthesis face ripple record are switched back to time and space territory, obtain effective surface ripple record or the synthesis face ripple record of base order mode and each higher order mode separated.

The method that effective surface ripple record according to base order mode and each higher order mode and synthesis face ripple record calculate target function value, comprises the following steps:

A, in the effective surface ripple record and synthesis face ripple record of base order mode and each higher order mode the spread speed of a selected Rayleigh face ripple respectively, the effective surface ripple record and synthesis face ripple record of base order mode and each higher order mode select lineups respectively clearly, and picks up amplitude respectively along respective lineups in the effective surface ripple record and synthesis face ripple record of base order mode and each higher order mode.

B, amplitude according to the effective surface ripple record of base order mode and each higher order mode and synthesis face ripple record, be utilized respectively formula (2) and calculate the target function value of base order mode and each higher order mode.

$E_k\left(Q\right)={\mathrm{\&Sigma;}}_{t=1}^n\left|log\left(\frac{A_k^{obs}\left(T\right(H_t\left)\right)}{A_k^{syn}\left(T\right(H_t\left)\right)}\right)\right|---\left(1\right)$

In formula, k represents the exponent number of the pattern of face ripple record, k=1,1��k��m during base order mode; M represents the face ripple record with m order mode; E_k(Q) target function value that kth order mode face ripple record is corresponding is represented;WithRepresent the effective surface ripple record of kth order mode respectively and synthesize T (H in the ripple record of face_t) moment pickup amplitude.

In practical application, face ripple records affected by many factors, only by the amplitude calculating target function value of a point improper. Therefore, it is possible to use the amplitude energy in a time window replaces the amplitude of a point. Time window is a time range, is sized to the crest at lineups place or the time span of trough, and by the effective surface ripple record of kth order mode and synthesis face ripple record, the amplitude energy in this time window is designated as respectivelyWith

With T (H_t) centered by take time window W_t, define two analytic signalsWithAt time window W_tInterior inner product is:

$\begin{array}{c}<\tilde{x},\tilde{y}{>}_{W_t}={\&Integral;}_{-W_t/2}^{W_t/2}\tilde{x}\left(t\right)\stackrel{\&OverBar;}{\tilde{y}}\left(t\right)dt\\ \tilde{x}\left(t\right)=x\left(t\right)+iH\left(x\right)\left(t\right)\end{array}---\left(2\right)$

In formula,It isConjugation, H (x) is the Hilbert transform of x (t), then time window W_tInterior amplitude energy can be expressed as $\mathrm{\&epsiv;}={(\tilde{A},A)}_{W_t}.$

Then the object function of kth order mode can be written as:

$E_k(Q,W_t)={\mathrm{\&Sigma;}}_{t=1}^n\left|log\left(\frac{{\mathrm{\&epsiv;}}_k^{obs}\left(T\right(H_t\left)\right)}{{\mathrm{\&epsiv;}}_k^{syn}\left(T\right(H_t\left)\right)}\right)\right|---\left(3\right)$

C, the target function value of base order mode and each higher order mode is weighted be added, obtaining total target function value is:

$E\left(Q\right)={\mathrm{\&Sigma;}}_{k=1}^m{q}_k{E}_k---\left(4\right)$

In formula, q_kFor the weight coefficient of the object function of kth order mode, the weight coefficient of the object function of different order mode has multiple value, it is possible to according to practical situation sweetly disposition; In the present embodiment, the weight coefficient of the object function of base order mode and each higher order mode all takes 1.

4) use simulated annealing inverting, obtain final near surface stratum Q-value, specifically include following steps:

The first step, set up the iterative vectorized of near surface stratum Q-valueWherein,Represent the Q-value of jth kind medium, 1��j��N, i >=1 during ith iteration;Make i=1.

Second step, make Qⁱ=Q^i-1, iterative target function E (Qⁱ)=E (Q^i-1), iteration temperature T_i=�� T_i-1. Wherein, �� is constant, iteration temperature T_iShould slowly reduce, therefore, generally take ��=0.99.

3rd step, one near surface stratum Q-value of stochastic generation comparison vectorWherein,Represent the comparison Q-value of the jth kind medium of stochastic generation, 1��j��N, i >=1 during ith iteration. Comparison vector according near surface stratum Q-valueUse seismic forward modeling simulation tool simulation comparison Rayleigh face ripple list big gun earthquake record (being called for short than opposite ripple record); Use high-resolution radon transform respectively to effective surface ripple record with carry out modal cutoff than opposite ripple record, obtain effective surface ripple record and the ratio opposite ripple record of base order mode and each higher order mode; Effective surface ripple record according to base order mode and each higher order mode and ratio opposite ripple record, pick up base order mode and the effective surface ripple record of each higher order mode and the amplitude than opposite ripple record respectively, calculate the comparison target function value of base order mode and each higher order mode; By the comparison target function value weighting summation of base order mode and each higher order mode, obtain total comparison target function value

4th step, calculating $\mathrm{\&Delta;}E=E\left({\stackrel{\&OverBar;}{Q}}^i\right)-E\left(Q^i\right).$

If �� E��0, illustrate that the comparison vector of stochastic generation can accept, take

If �� E > 0, calculate�� is intermediate computations parameter;

If �� (�� E) > R, illustrate that the comparison vector of stochastic generation can accept, takeR is equally distributed random number between (0,1);

Otherwise do not accept the comparison vector of stochastic generation, QⁱIt is worth constant.

If the comparison vector of the 5th continuous 20 stochastic generation of stepDo not accepted, i.e. continuous 20 QⁱBeing worth constant, terminate inverting, output final result is Q^est=Qⁱ; Otherwise, i is made_new=i+1, returns second step.

In above-described embodiment, described step 2) in seismic forward modeling simulation tool can be finite difference calculus, FInite Element or other method.

In above-described embodiment, described step 4) in the algorithm that realizes of simulated annealing that adopts be Metropolis (in Mai Teluobo) algorithm; In order to accelerate described step 4) in the iterative process of simulated annealing, it is possible to adopt Ingber (Engelberg) a kind of Fast Simulated Annealing proposed that comparison vector is updated. Fast Simulated Annealing uses and depends on being distributed in around current iteration vector like Cauchy (Cauchy) and producing new comparison vector of temperature, and hunting zone tapers into the reduction of temperature.

The quick more new formula of vector is by ratio:

$\begin{array}{c}{\stackrel{\&OverBar;}{Q}}^i=Q^i+y(Q_{\max }-Q_{\min })\\ y=T_i\mathrm{sgn}(u-0.5)\&lsqb;{(1+1/T_i)}^{|2u-1|}-1\&rsqb;\end{array}---\left(5\right)$

In formula, y is intermediate computations parameter, and u is random number, is evenly distributed on [0,1] interval; Sgn is for taking sign function.

The various embodiments described above are merely to illustrate the present invention; wherein each parts structure, position and connected mode etc. thereof are set all can be varied from; every equivalents carried out on the basis of technical solution of the present invention and improvement, all should not get rid of outside protection scope of the present invention.

Claims (10)

1. a method of estimation for complex dielectrics near surface stratum Q-value, it comprises the following steps:

1) carry out seismic experiment in target area, gathered the actual measurement Rayleigh face ripple list big gun earthquake record of target area by multiple cymoscopes, and record the distance between each cymoscope and focus and seismic wave and propagate the time used by each cymoscope from focus; Adopt wave analysis method in multiple tracks face to process actual measurement Rayleigh face ripple list big gun earthquake record, it is thus achieved that the velocity of longitudinal wave near surface stratum, shear wave velocity and Media density, show that near surface stratum has the medium that N kind is different;

2) set the initial value of a Q-value to every kind of medium, the initial vector setting up near surface stratum Q-value isWherein,Represent the initial value of the Q-value of jth kind medium, 1��j��N; And set the initial temperature T for simulated annealing₀; Initial vector Q according near surface stratum Q-value⁰, use seismic forward modeling simulation tool to simulate the synthesis Rayleigh face ripple list big gun earthquake record being made up of multiple tracks face ripple record;

3) use high-resolution radon transform respectively actual measurement Rayleigh face ripple list big gun earthquake record and synthesis Rayleigh face ripple list big gun earthquake record to be carried out modal cutoff, obtain actual measurement Rayleigh face ripple list big gun earthquake record and the synthesis Rayleigh face ripple list big gun earthquake record of base order mode and each higher order mode; Actual measurement Rayleigh face ripple list big gun earthquake record according to base order mode and each higher order mode and synthesis Rayleigh face ripple list big gun earthquake record, calculate the initial target functional value of base order mode and each higher order mode respectively; By the initial target functional value weighting summation of base order mode and each higher order mode, obtain total initial target functional value;

4) use simulated annealing inverting, obtain final near surface stratum Q-value.

2. the method for estimation of a kind of complex dielectrics near surface stratum as claimed in claim 1 Q-value, it is characterized in that, described step 3) in use high-resolution radon transform that actual measurement Rayleigh face ripple list big gun earthquake record or synthesis Rayleigh face ripple list big gun earthquake record are carried out modal cutoff, specifically include following steps:

1. use high-resolution radon transform that from time and space territory, actual measurement Rayleigh face ripple list big gun earthquake record or synthesis Rayleigh face ripple list big gun earthquake record are transformed to frequency slowness domain, obtain actual measurement Rayleigh face ripple list big gun earthquake record or the frequency-slowness section of synthesis Rayleigh face ripple list big gun earthquake record, actual measurement Rayleigh face ripple list big gun earthquake record or the synthesis Rayleigh face ripple list big gun earthquake record of base order mode and each higher order mode are respectively provided with scope clearly in frequency-slowness section, and not overlapping each other;

2. in frequency slowness section, the scope of the actual measurement Rayleigh face ripple list big gun earthquake record or synthesis Rayleigh face ripple list big gun earthquake record that represent base order mode and each higher order mode is determined respectively, the value of actual measurement Rayleigh face ripple list big gun earthquake record or synthesis Rayleigh face ripple list big gun earthquake record by removing base order mode or a certain higher order mode other order mode extraneous in frequency slowness section, obtain the frequency slowness section only comprising the actual measurement Rayleigh face ripple list big gun earthquake record of base order mode or a certain higher order mode or synthesis Rayleigh face ripple list big gun earthquake record respectively,

3. adopt and draw east inverse transformation that from frequency-slowness domain, the actual measurement Rayleigh face ripple list big gun earthquake record of base order mode and each higher order mode or synthesis Rayleigh face ripple list big gun earthquake record are switched back to time and space territory, obtain actual measurement Rayleigh face ripple list big gun earthquake record or the synthesis Rayleigh face ripple list big gun earthquake record of base order mode and each higher order mode separated.

3. the method for estimation of a kind of complex dielectrics near surface stratum as claimed in claim 1 or 2 Q-value, it is characterized in that, described step 3) in calculate target function value according to the actual measurement Rayleigh face ripple list big gun earthquake record of base order mode and each higher order mode and synthesis Rayleigh face ripple list big gun earthquake record, specifically include following steps:

A, in the actual measurement Rayleigh face ripple list big gun earthquake record and synthesis Rayleigh face ripple list big gun earthquake record of base order mode and each higher order mode the spread speed of a selected Rayleigh face ripple respectively, the actual measurement Rayleigh face ripple list big gun earthquake record and synthesis Rayleigh face ripple list big gun earthquake record of base order mode and each higher order mode select lineups respectively clearly, and picks up amplitude respectively along respective lineups in the actual measurement Rayleigh face ripple list big gun earthquake record and synthesis Rayleigh face ripple list big gun earthquake record of base order mode and each higher order mode;

B, amplitude according to the actual measurement Rayleigh face ripple list big gun earthquake record of base order mode and each higher order mode and synthesis Rayleigh face ripple list big gun earthquake record, calculate the target function value of base order mode and each higher order mode respectively;

C, by the target function value weighting summation of base order mode and each higher order mode, obtain total target function value.

4. the method for estimation of a kind of complex dielectrics near surface stratum as claimed in claim 3 Q-value, it is characterised in that in described step b, the computing formula of target function value is:

$E\left(Q\right)=\underset{t=1}{\overset{n}{\&Sigma;}}\left|log\left(\frac{A^{obs}\left(T\right(H_t\left)\right)}{A^{syn}\left(T\right(H_t\left)\right)}\right)\right|$

In formula, E (Q) represents the target function value that Q-value is corresponding; N represents the number of channels of the face ripple record comprised in actual measurement Rayleigh face ripple list big gun earthquake record or synthesis Rayleigh face ripple list big gun earthquake record, t represents t road actual measurement Rayleigh face ripple list big gun earthquake record or the number of channels of synthesis Rayleigh face ripple list big gun earthquake record, 1��t��n; H_tRepresent the distance between cymoscope and the focus of actual measurement Rayleigh face, record t road ripple list big gun earthquake tracer signal; T (H_t) represent that seismic wave propagates the time used by cymoscope recording actual measurement Rayleigh face, t road ripple list big gun earthquake tracer signal from focus; A^obs(T(H_t)) and A^syn(T(H_t)) represent T (H in actual measurement Rayleigh face ripple list big gun earthquake record and synthesis Rayleigh face ripple list big gun earthquake record respectively_t) moment pickup amplitude;

Or the computing formula of target function value is:

$E(Q,W_t)=\underset{t=1}{\overset{n}{\&Sigma;}}\left|log\left(\frac{{\mathrm{\&epsiv;}}^{obs}\left(T\right(H_t\left)\right)}{{\mathrm{\&epsiv;}}^{syn}\left(T\right(H_t\left)\right)}\right)\right|$

In formula, W_tRepresent with T (H_t) centered by the time window that takes; ��^obs(T(H_t)) and ��^syn(T(H_t)) represent that actual measurement Rayleigh face ripple list big gun earthquake record and synthesis Rayleigh face ripple list big gun earthquake record are at time window W respectively_tInterior amplitude energy.

5. the method for estimation of a kind of complex dielectrics near surface stratum Q-value as described in claim 1 or 2 or 4, it is characterised in that described step 4) middle use simulated annealing inverting near surface stratum Q-value, specifically include following steps:

The first step, set up the iterative vectorized of near surface stratum Q-valueWherein,Represent the Q-value of jth kind medium, 1��j��N, i >=1 during ith iteration; Make i=1;

Second step, make Qⁱ=Q^i-1, iterative target function E (Qⁱ)=E (Q^i-1), iteration temperature T_i=�� T_i-1; Wherein, �� is constant, takes ��=0.99;

3rd step, generate the comparison vector of near surface stratum Q-valueWherein,The comparison Q-value of the jth kind medium generated when representing ith iteration, 1��j��N, i >=1; Comparison vector according near surface stratum Q-valueUse seismic forward modeling simulation tool simulation comparison Rayleigh face ripple list big gun earthquake record; Use high-resolution radon transform respectively actual measurement Rayleigh face ripple list big gun earthquake record and comparison Rayleigh face ripple list big gun earthquake record to be carried out modal cutoff, obtain actual measurement Rayleigh face ripple list big gun earthquake record and the comparison Rayleigh face ripple list big gun earthquake record of base order mode and each higher order mode; Actual measurement Rayleigh face ripple list big gun earthquake record according to base order mode and each higher order mode and comparison Rayleigh face ripple list big gun earthquake record, the amplitude surveying Rayleigh face ripple list big gun earthquake record and comparison Rayleigh face ripple list big gun earthquake record of pickup base order mode and each higher order mode, calculates the comparison target function value of base order mode and each higher order mode respectively; By the comparison target function value weighting summation of base order mode and each higher order mode, obtain total comparison target function value

4th step, calculating $\mathrm{\&Delta;}E=E\left({\stackrel{\&OverBar;}{Q}}^i\right)-E\left(Q^i\right);$

If �� E��0, accept the comparison vector generated, take

If �� E > 0, calculate $\mathrm{\&rho;}\left(\mathrm{\&Delta;}E\right)=\exp (-\frac{\mathrm{\&Delta;}E}{T_i}),$ �� is intermediate computations parameter;

If �� (�� E) > R, accept the comparison vector generated, takeR is equally distributed random number between (0,1);

Otherwise do not accept the comparison vector generated, QⁱIt is worth constant;

If the comparison vector that the 5th step generates for continuous 20 timesDo not accepted, i.e. continuous 20 QⁱBeing worth constant, terminate inverting, output final result is Q^est=Qⁱ; Otherwise, i is made_new=i+1, returns second step.

6. the method for estimation of a kind of complex dielectrics near surface stratum as claimed in claim 5 Q-value, it is characterised in that adopt random method to generate comparison vector in described 3rd step, or adopt the Ingber Fast Simulated Annealing proposed to generate comparison vector; The formula that quickly generates of described comparison vector is:

${\stackrel{\&OverBar;}{Q}}^i=Q^i+y(Q_{max}-Q_{min})$

Y=T_isgn(u-0.5)[(1+1/T_i)^|2u-1|-1]

In formula, y is intermediate computations parameter; Q_maxAnd Q_minThe respectively upper limit value and lower limit value of the Q-value span of each medium, Q_min=10, Q_max=200; U is random number, is evenly distributed on [0,1] interval; Sgn is for taking sign function.

7. the method for estimation of a kind of complex dielectrics near surface stratum Q-value as described in claim 1 or 2 or 4 or 6, it is characterised in that described step 2) in seismic forward modeling simulation tool be finite difference calculus or FInite Element.

8. the method for estimation of a kind of complex dielectrics near surface stratum as claimed in claim 5 Q-value, it is characterised in that described step 2) in seismic forward modeling simulation tool be finite difference calculus or FInite Element.

9. the method for estimation of a kind of complex dielectrics near surface stratum Q-value as described in claim 1 or 2 or 4 or 6 or 8, it is characterised in that described step 4) in the algorithm that realizes of simulated annealing that adopts be Metropolis algorithm.

10. the method for estimation of a kind of complex dielectrics near surface stratum as claimed in claim 7 Q-value, it is characterised in that described step 4) in the algorithm that realizes of simulated annealing that adopts be Metropolis algorithm.