CN104216011A - Reverse time migration method of stable qP wave in TTI (tilted transversely isotropic) media - Google Patents

Abstract

The invention relates to the technical field of processing of seismic exploration data and relates to a reverse time migration method of stable qP wave in TTI (tilted transversely isotropic) media. The method includes the steps: a qP wave equation with regularization items is deducted on the basis of an existing stable qP wave equation in TTI media; adaptive regularization parameters are selected: Sigma=max{sqrt( <2>8( )),sqrt( <2>Phi( )), with 8 referring to a tilt angle of a formation and Phi referring to an azimuth angle of the formation. On the basis of the existing stable qP propagator, regularization items are added to the qP wave equation, a more stable qP wave equation is obtained without excess calculation and storage, TTI media RTM (reverse time migration) is more widely applicable to the model, and the algorithm is more practical; through tilt angle and azimuth angle information of the formation, magnitudes of regularization coefficients are automatically determined without excess manual intervention, and a regularization equation is more steady in processing actual data.

Description

A kind of qP ripple reverse-time migration method of TTI media stabilize

Technical field:

The present invention relates to seismic prospecting data processing technology field, be specifically related to a kind of qP ripple reverse-time migration method of TTI media stabilize.

Background technology:

In TTI medium, stable propagation qP ripple is that can TTI RTM is current practical key, is also study hotspot problem.QP ripple propagation operator research stable in current TTI medium mainly contains following two kinds of thinkings: (1) is from equations for elastic waves.Directly from anisotropic medium equations for elastic waves, do acoustic approximation and obtain the onomatopoeia wave equation corresponding with equations for elastic waves, the equation obtaining has like this retained the most of feature of former equations for elastic waves system, such as in the situation that meeting certain boundary conditions, the conservativeness that keeps elastic potential energy and kinetic energy summation, stability when this can strengthen equation solution.But simultaneously the qP wave equation corresponding with elastic wave more complicated still, solves efficiency lower, the original intention of this and acoustic approximation is run counter to a bit.On the onomatopoeia wave equation basis corresponding with elastic wave, can ignore angle derivative term, obtain approximate equation, but only comprise first order differential operator in this equation, cannot solve this equation with common single order central difference schemes is stable, this is comparatively unfavorable for reverse-time migration; (2) from dispersion relation.QP wave equation in TTI medium can also be derived from the dispersion relation of coupling.But the qP wave equation of deriving from dispersion relation is excessive to the former equations for elastic waves system reform, cannot keep the good feature of elastic wave system.Improper such as choosing due to intermediate variable, even boundary value condition meets in some equation communication process being derived by dispersion relation, potential energy and kinetic energy summation conservation that also cannot keeping system, this can make numerical solution extremely unstable.Stable in order to make by dispersion relation derived equation numerical solution, can do some to acoustic approximation and compromise, introduce limited shear wave velocity, obtain so-called limited shear wave equation.Limited shear wave equation solution is more stable, but also has the problem of himself, such as remaining shear wave energy is stronger than acoustic approximation equation, solves calculated amount also larger.In addition, limited shear wave equation can strengthen solving stability, but can not guarantee absolute.

Summary of the invention:

A kind of qP ripple reverse-time migration method that the object of this invention is to provide TTI media stabilize, it is on existing stable qP propagation operator basis, qP wave equation is added to regularization term, do not increasing in too much calculated amount and storage capacity situation, obtain more stable qP wave equation, make TTI medium RTM wider to the adaptability of model, increased the practicality of algorithm; Utilize stratigraphic dip and azimuth information, automatically determine the size of regularization coefficient, without artificial, too much intervene, make regularization equation more sane when processing real data.

In order to solve the existing problem of background technology, the present invention by the following technical solutions: qP wave equation in A, band regularization term TTI medium;

(a) by equations for elastic waves, derive qP wave equation

For the equations for elastic waves of doing in TTI medium after acoustic approximation, be:

$\mathrm{\&rho;}\frac{\&PartialD;v_x}{\&PartialD;t}=d1\left({\mathrm{\&sigma;}}_{11}\right)+d3({\mathrm{\&sigma;}}_{33}-{\mathrm{\&sigma;}}_{11})$

(1)

Wherein differentiating operator is defined as:

For the equation after acoustic approximation in TTI medium (1), for counting yield, ignore the wherein derivative term of relevant angle, (1) becomes so:

$\mathrm{\&rho;}\frac{\&PartialD;v_x}{\&PartialD;t}=d1\left({\mathrm{\&sigma;}}_{11}\right)$

$\mathrm{\&rho;}\frac{{\&PartialD;v}_y}{\&PartialD;t}=d2\left({\mathrm{\&sigma;}}_{11}\right)$ $\left(3\right)$

$\mathrm{\&rho;}\frac{\&PartialD;v_z}{\&PartialD;t}=d3\left({\mathrm{\&sigma;}}_{33}\right)$

${\&PartialD;}_t\left[\begin{array}{c}{\mathrm{\&sigma;}}_{11}\\ {\mathrm{\&sigma;}}_{33}\end{array}\right]=\left[\begin{array}{cc}{C}_{11}& C_{13}\\ C_{13}& C_{33}\end{array}\right]\left[\begin{array}{c}{d}^{T}1\left(v_x\right)+d^{T}2\left(v_y\right)\\ d^{T}3\left(v_z\right)\end{array}\right]$

Under the hypothesis that is 1 in density, its corresponding second-order equation can be written as:

$\frac{{\&PartialD;}^2{\mathrm{\&sigma;}}_H}{\&PartialD;t^2}=v_p^2[(1+2\mathrm{\&epsiv;})(d{1}^{T}d1\left({\mathrm{\&sigma;}}_H\right)+d{2}^{T}d2\left({\mathrm{\&sigma;}}_H\right))+\sqrt{(1+2\mathrm{\&delta;})}d{3}^{T}d3\left({\mathrm{\&sigma;}}_H\right)]$ $\left(4\right)$

$\frac{{\&PartialD;}^2{\mathrm{\&sigma;}}_v}{\&PartialD;t^2}=v_p^2[\sqrt{(1+2\mathrm{\&delta;})}(d{1}^{T}d1\left({\mathrm{\&sigma;}}_H\right)+d{2}^{T}d2\left({\mathrm{\&sigma;}}_H\right))+d{3}^{T}d3\left({\mathrm{\&sigma;}}_H\right)]$

(b) by dispersion relation, derive qP wave equation;

Flecter (2009) proposes to utilize following dispersion relation to derive corresponding limited shear wave equation:

${\mathrm{\&omega;}}^4=[(v_{\mathrm{px}}^2+v_{\mathrm{sz}}^2)(k_x^{\&prime;2}+k_y^{\&prime;2})+(v_{\mathrm{pz}}^2+v_{\mathrm{sz}}^2)k_z^{\&prime;2}]{\mathrm{\&omega;}}^2-v_{\mathrm{px}}^2{v}_{\mathrm{sz}}^2{(k_x^{\&prime;2}+k_y^{\&prime;2})}^2-$

$v_{\mathrm{pz}}^2{v}_{\mathrm{sz}}^2{k}_z^{\&prime;4}+[v_{\mathrm{pz}}^2(v_{\mathrm{pn}}^2-v_{\mathrm{px}}^2)-v_{\mathrm{sz}}^2(v_{\mathrm{pn}}^2+v_{\mathrm{pz}}^2)](k_x^{\&prime;2}+k_y^{\&prime;2})k_z^{\&prime;2}$

Introduce after suitable intermediate wave field variable, the limited shear wave equation of its correspondence is:

$\frac{{\&PartialD;}^{2}p}{\&PartialD;t^2}=[v_p^2(1+2\mathrm{\&epsiv;})H_2+v_{\mathrm{sz}}^2{H}_1]p+[v_p^2-v_{\mathrm{sz}}^2]H_{1}q$

$\frac{{\&PartialD;}^{2}q}{\&PartialD;t^2}=[(1+2\mathrm{\&delta;})v_p^2-v_{\mathrm{sz}}^2]H_{2}p+(v_p^2{H}_1+v_{\mathrm{sz}}^2{H}_2)q$

Wherein,

$H_1={\sin }^2\mathrm{\&theta;}{\cos }^2\mathrm{\&phi;}\frac{{\&PartialD;}^2}{\&PartialD;x^2}+{\sin }^2\mathrm{\&theta;}{\sin }^2\mathrm{\&phi;}\frac{{\&PartialD;}^2}{\&PartialD;y^2}+{\cos }^2\mathrm{\&theta;}\frac{{\&PartialD;}^2}{\&PartialD;z^2}+$

${\sin }^2\mathrm{\&theta;}\sin 2\mathrm{\&phi;}\frac{{\&PartialD;}^2}{\&PartialD;x\&PartialD;y}+\sin 2\mathrm{\&theta;}\sin \mathrm{\&phi;}\frac{{\&PartialD;}^2}{\&PartialD;y\&PartialD;z}+\sin 2\mathrm{\&theta;}\cos \mathrm{\&phi;}\frac{{\&PartialD;}^2}{\&PartialD;x\&PartialD;z}$

$H_2=\frac{{\&PartialD;}^2}{\&PartialD;x^2}+\frac{{\&PartialD;}^2}{\&PartialD;y^2}+\frac{{\&PartialD;}^2}{\&PartialD;z^2}-H_1$

(c) with the stable qP wave equation of regularization term

In isotropy ACOUSTIC WAVE EQUATION solution procedure, Liu (2008) has proposed a kind of optimization method of anti-frequency dispersion, by revising dispersion relation corresponding to ACOUSTIC WAVE EQUATION, obtains anti-dispersion equation, and its basic thought is as follows:

For ACOUSTIC WAVE EQUATION in isotropic medium:

$\frac{1}{c_0^2}\frac{{\&PartialD;}^{2}p}{{\&PartialD;}^{2}t}={\&dtri;}^{2}p---\left(7\right)$

Plane wave solves form:

$p=a{e}^{i(\mathrm{\&omega;t}\&PlusMinus;\stackrel{\&OverBar;}{k}x)}---\left(8\right)$

For being with ACOUSTIC WAVE EQUATION in regularization term isotropic medium:

$\frac{1}{c_0^2}\frac{{\&PartialD;}^{2}p}{{\&PartialD;}^{2}t}={\&dtri;}^{2}p+\frac{\mathrm{\&sigma;}}{v}\frac{\&PartialD;}{\&PartialD;t}{\&dtri;}^{2}p---\left(9\right)$

Its corresponding plane wave solution is:

$p=e^{-\frac{\mathrm{\&sigma;}}{2}{k}^2{c}_{0}t}{e}^{i\left[\sqrt{1-\frac{{\mathrm{\&sigma;}}^2{k}^2}{4}}\right]\mathrm{\&omega;t}\&PlusMinus;\stackrel{\&RightArrow;}{k}x]}---\left(10\right)$

Notice:

$1-\frac{{\mathrm{\&sigma;}}^2{k}^2}{4}\&ap;1---\left(11\right)$

(10) can write:

$p=e^{-\frac{\mathrm{\&sigma;}}{2}{k}^2{c}_{0}t}{e}^{i[\mathrm{\&omega;t}\&PlusMinus;\stackrel{\&RightArrow;}{k}x]}---\left(12\right)$

In solving TTI medium, when qP wave equation, introduce regularization term, high wave number composition useless in unstable solution is done to suitable decay, obtain more stable equation.The limited shear wave equation proposing for Flecter, adds following regularization term:

Succinct in order to describe, (6) are rewritten as:

$\frac{{\&PartialD;}^{2}p}{\&PartialD;t^2}=G_p^{1}p+G_p^{2}q$ $\left(13\right)$

$\frac{{\&PartialD;}^{2}q}{\&PartialD;t^2}=G_q^{1}p+G_q^{2}q$

Add after regularization term, (9) can become:

$\frac{{\&PartialD;}^{2}p}{\&PartialD;t^2}=G_p^{1}p+G_p^{2}q+\mathrm{\&sigma;}\frac{\&PartialD;}{\&PartialD;t}(G_p^{1}p+G_p^{2}q)---\left(14\right)$

$\frac{{\&PartialD;}^{2}q}{\&PartialD;t^2}=G_q^{1}p+G_q^{2}q+\mathrm{\&sigma;}\frac{\&PartialD;}{\&PartialD;t}(G_q^{1}p+G_q^{2}q)$

Wherein σ is regularization coefficient, can verify, solves (14) than the labile element that solves high wave number in (13) wave field still less, thereby can make TTI medium RTM applicability when processing real data stronger.For equation (4), also can adopt similar thought to build it with the form of regularization term.

B, adaptive regularization parameter are selected: from (14), can find out, choosing of regularization parameter has a great impact final wave field result of calculation, in order to make algorithm more sane when processing real data, utilize the adaptive calculating regularization of following mode coefficient:

$\mathrm{\&sigma;}=\max \{\mathrm{sqrt}\left({\&dtri;}^2\mathrm{\&theta;}\left(\stackrel{\&RightArrow;}{x}\right)\right),\mathrm{sqrt}\left({\&dtri;}^2\mathrm{\&phi;}\left(\stackrel{\&RightArrow;}{x}\right)\right)---\left(15\right)$

Wherein, θ, φ are respectively inclination angle and the position angle on stratum.

The present invention has following beneficial effect: it is on existing stable qP propagation operator basis, qP wave equation is added to regularization term, do not increasing in too much calculated amount and storage capacity situation, obtain more stable qP wave equation, make TTI medium RTM wider to the adaptability of model, increased the practicality of algorithm; Utilize stratigraphic dip and azimuth information, automatically determine the size of regularization coefficient, without artificial, too much intervene, make regularization equation more sane when processing real data.

Accompanying drawing explanation:

Fig. 1 is single big gun TTI medium qP ripple RTM calculation flow chart that the present invention proposes.

Embodiment:

This embodiment is taked qP wave equation in following technical scheme: A, band regularization term TTI medium;

(a) by equations for elastic waves, derive qP wave equation

For the equations for elastic waves of doing in TTI medium after acoustic approximation, be:

$\mathrm{\&rho;}\frac{{\&PartialD;v}_x}{\&PartialD;t}=d1\left({\mathrm{\&sigma;}}_{11}\right)+d3({\mathrm{\&sigma;}}_{33}-{\mathrm{\&sigma;}}_{11})$

(1)

Wherein differentiating operator is defined as:

For the equation after acoustic approximation in TTI medium (1), for counting yield, ignore the wherein derivative term of relevant angle, (1) becomes so:

$\mathrm{\&rho;}\frac{\&PartialD;v_x}{\&PartialD;t}=d1\left({\mathrm{\&sigma;}}_{11}\right)$

$\mathrm{\&rho;}\frac{{\&PartialD;v}_y}{\&PartialD;t}=d2\left({\mathrm{\&sigma;}}_{11}\right)$ $\left(3\right)$

$\mathrm{\&rho;}\frac{\&PartialD;v_z}{\&PartialD;t}=d3\left({\mathrm{\&sigma;}}_{33}\right)$

${\&PartialD;}_t\left[\begin{array}{c}{\mathrm{\&sigma;}}_{11}\\ {\mathrm{\&sigma;}}_{33}\end{array}\right]=\left[\begin{array}{cc}{C}_{11}& C_{13}\\ C_{13}& C_{33}\end{array}\right]\left[\begin{array}{c}{d}^{T}1\left(v_x\right)+d^{T}2\left(v_y\right)\\ d^{T}3\left(v_z\right)\end{array}\right]$

Under the hypothesis that is 1 in density, its corresponding second-order equation can be written as:

$\frac{{\&PartialD;}^2{\mathrm{\&sigma;}}_H}{\&PartialD;t^2}=v_p^2[(1+2\mathrm{\&epsiv;})(d{1}^{T}d1\left({\mathrm{\&sigma;}}_H\right)+d{2}^{T}d2\left({\mathrm{\&sigma;}}_H\right))+\sqrt{(1+2\mathrm{\&delta;})}d{3}^{T}d3\left({\mathrm{\&sigma;}}_H\right)]$

$\frac{{\&PartialD;}^2{\mathrm{\&sigma;}}_v}{\&PartialD;t^2}=v_p^2[\sqrt{(1+2\mathrm{\&delta;})}(d{1}^{T}d1\left({\mathrm{\&sigma;}}_H\right)+d{2}^{T}d2\left({\mathrm{\&sigma;}}_H\right))+d{3}^{T}d3\left({\mathrm{\&sigma;}}_H\right)]$

(b) by dispersion relation, derive qP wave equation;

Flecter (2009) proposes to utilize following dispersion relation to derive corresponding limited shear wave equation:

${\mathrm{\&omega;}}^4=[(v_{\mathrm{px}}^2+v_{\mathrm{sz}}^2)(k_x^{\&prime;2}+k_y^{\&prime;2})+(v_{\mathrm{pz}}^2+v_{\mathrm{sz}}^2)k_z^{\&prime;2}]{\mathrm{\&omega;}}^2-v_{\mathrm{px}}^2{v}_{\mathrm{sz}}^2{(k_x^{\&prime;2}+k_y^{\&prime;2})}^2-$

$v_{\mathrm{pz}}^2{v}_{\mathrm{sz}}^2{k}_z^{\&prime;4}+[v_{\mathrm{pz}}^2(v_{\mathrm{pn}}^2-v_{\mathrm{px}}^2)-v_{\mathrm{sz}}^2(v_{\mathrm{pn}}^2+v_{\mathrm{pz}}^2)](k_x^{\&prime;2}+k_y^{\&prime;2})k_z^{\&prime;2}$

Introduce after suitable intermediate wave field variable, the limited shear wave equation of its correspondence is:

$\frac{{\&PartialD;}^{2}p}{\&PartialD;t^2}=[v_p^2(1+2\mathrm{\&epsiv;})H_2+v_{\mathrm{sz}}^2{H}_1]p+[v_p^2-v_{\mathrm{sz}}^2]H_{1}q$

$\frac{{\&PartialD;}^{2}q}{\&PartialD;t^2}=[(1+2\mathrm{\&delta;})v_p^2-v_{\mathrm{sz}}^2]H_{2}p+(v_p^2{H}_1+v_{\mathrm{sz}}^2{H}_2)q$

Wherein,

$H_1={\sin }^2\mathrm{\&theta;}{\cos }^2\mathrm{\&phi;}\frac{{\&PartialD;}^2}{\&PartialD;x^2}+{\sin }^2\mathrm{\&theta;}{\sin }^2\mathrm{\&phi;}\frac{{\&PartialD;}^2}{\&PartialD;y^2}+{\cos }^2\mathrm{\&theta;}\frac{{\&PartialD;}^2}{\&PartialD;z^2}+$

${\sin }^2\mathrm{\&theta;}\sin 2\mathrm{\&phi;}\frac{{\&PartialD;}^2}{\&PartialD;x\&PartialD;y}+\sin 2\mathrm{\&theta;}\sin \mathrm{\&phi;}\frac{{\&PartialD;}^2}{\&PartialD;y\&PartialD;z}+\sin 2\mathrm{\&theta;}\cos \mathrm{\&phi;}\frac{{\&PartialD;}^2}{\&PartialD;x\&PartialD;z}$

$H_2=\frac{{\&PartialD;}^2}{\&PartialD;x^2}+\frac{{\&PartialD;}^2}{\&PartialD;y^2}+\frac{{\&PartialD;}^2}{\&PartialD;z^2}-H_1$

(c) with the stable qP wave equation of regularization term

In isotropy ACOUSTIC WAVE EQUATION solution procedure, Liu (2008) has proposed a kind of optimization method of anti-frequency dispersion, by revising dispersion relation corresponding to ACOUSTIC WAVE EQUATION, obtains anti-dispersion equation, and its basic thought is as follows:

For ACOUSTIC WAVE EQUATION in isotropic medium:

$\frac{1}{c_0^2}\frac{{\&PartialD;}^{2}p}{{\&PartialD;}^{2}t}={\&dtri;}^{2}p---\left(7\right)$

Plane wave solves form:

$p=a{e}^{i(\mathrm{\&omega;t}\&PlusMinus;\stackrel{\&OverBar;}{k}x)}---\left(8\right)$

For being with ACOUSTIC WAVE EQUATION in regularization term isotropic medium:

$\frac{1}{c_0^2}\frac{{\&PartialD;}^{2}p}{{\&PartialD;}^{2}t}={\&dtri;}^{2}p+\frac{\mathrm{\&sigma;}}{v}\frac{\&PartialD;}{\&PartialD;t}{\&dtri;}^{2}p---\left(9\right)$

Its corresponding plane wave solution is:

$p=e^{-\frac{\mathrm{\&sigma;}}{2}{k}^2{c}_{0}t}{e}^{i\left[\sqrt{1-\frac{{\mathrm{\&sigma;}}^2{k}^2}{4}}\right]\mathrm{\&omega;t}\&PlusMinus;\stackrel{\&RightArrow;}{k}x]}---\left(10\right)$

Notice:

$1-\frac{{\mathrm{\&sigma;}}^2{k}^2}{4}\&ap;1---\left(11\right)$

(10) can write:

$p=e^{-\frac{\mathrm{\&sigma;}}{2}{k}^2{c}_{0}t}{e}^{i[\mathrm{\&omega;t}\&PlusMinus;\stackrel{\&RightArrow;}{k}x]}---\left(12\right)$

In solving TTI medium, when qP wave equation, introduce regularization term, high wave number composition useless in unstable solution is done to suitable decay, obtain more stable equation.The limited shear wave equation proposing for Flecter, adds following regularization term:

Succinct in order to describe, (6) are rewritten as:

$\frac{{\&PartialD;}^{2}p}{\&PartialD;t^2}=G_p^{1}p+G_p^{2}q$ $\left(13\right)$

$\frac{{\&PartialD;}^{2}q}{\&PartialD;t^2}=G_q^{1}p+G_q^{2}q$

Add after regularization term, (9) can become:

$\frac{{\&PartialD;}^{2}p}{\&PartialD;t^2}=G_p^{1}p+G_p^{2}q+\mathrm{\&sigma;}\frac{\&PartialD;}{\&PartialD;t}(G_p^{1}p+G_p^{2}q)$ $\left(14\right)$

$\frac{{\&PartialD;}^{2}q}{\&PartialD;t^2}=G_q^{1}p+G_q^{2}q+\mathrm{\&sigma;}\frac{\&PartialD;}{\&PartialD;t}(G_q^{1}p+G_q^{2}q)$

Wherein σ is regularization coefficient, can verify, solves (14) than the labile element that solves high wave number in (13) wave field still less, thereby can make TTI medium RTM applicability when processing real data stronger.For equation (4), also can adopt similar thought to build it with the form of regularization term.

B, adaptive regularization parameter are selected: from (14), can find out, choosing of regularization parameter has a great impact final wave field result of calculation, in order to make algorithm more sane when processing real data, utilize the adaptive calculating regularization of following mode coefficient:

$\mathrm{\&sigma;}=\max \{\mathrm{sqrt}\left({\&dtri;}^2\mathrm{\&theta;}\left(\stackrel{\&RightArrow;}{x}\right)\right),\mathrm{sqrt}\left({\&dtri;}^2\mathrm{\&phi;}\left(\stackrel{\&RightArrow;}{x}\right)\right)---\left(15\right)$

Wherein, θ, φ are respectively inclination angle and the position angle on stratum.

This embodiment has following beneficial effect: it is on existing stable qP propagation operator basis, qP wave equation is added to regularization term, do not increasing in too much calculated amount and storage capacity situation, obtain more stable qP wave equation, make TTI medium RTM wider to the adaptability of model, increased the practicality of algorithm; Utilize stratigraphic dip and azimuth information, automatically determine the size of regularization coefficient, without artificial, too much intervene, make regularization equation more sane when processing real data.

Claims (1)

1. a qP ripple reverse-time migration method for TTI media stabilize, is characterized in that its main technical content is: (A), qP wave equation in band regularization term TTI medium;

(a) by equations for elastic waves, derive qP wave equation

For the equations for elastic waves of doing in TTI medium after acoustic approximation, be:

$\mathrm{\&rho;}\frac{\&PartialD;v_x}{\&PartialD;t}=d1\left({\mathrm{\&sigma;}}_{11}\right)+d3({\mathrm{\&sigma;}}_{33}-{\mathrm{\&sigma;}}_{11})$

(1)

Wherein differentiating operator is defined as:

For the equation after acoustic approximation in TTI medium (1), for counting yield, ignore the wherein derivative term of relevant angle, (1) becomes so:

$\mathrm{\&rho;}\frac{\&PartialD;v_x}{\&PartialD;t}=d1\left({\mathrm{\&sigma;}}_{11}\right)$

$\mathrm{\&rho;}\frac{{\&PartialD;v}_y}{\&PartialD;t}=d2\left({\mathrm{\&sigma;}}_{11}\right)$ $\left(3\right)$

$\mathrm{\&rho;}\frac{\&PartialD;v_z}{\&PartialD;t}=d3\left({\mathrm{\&sigma;}}_{33}\right)$

${\&PartialD;}_t\left[\begin{array}{c}{\mathrm{\&sigma;}}_{11}\\ {\mathrm{\&sigma;}}_{33}\end{array}\right]=\left[\begin{array}{cc}{C}_{11}& C_{13}\\ C_{13}& C_{33}\end{array}\right]\left[\begin{array}{c}{d}^{T}1\left(v_x\right)+d^{T}2\left(v_y\right)\\ d^{T}3\left(v_z\right)\end{array}\right]$

Under the hypothesis that is 1 in density, its corresponding second-order equation can be written as:

$\frac{{\&PartialD;}^2{\mathrm{\&sigma;}}_H}{\&PartialD;t^2}=v_p^2[(1+2\mathrm{\&epsiv;})(d{1}^{T}d1\left({\mathrm{\&sigma;}}_H\right)+d{2}^{T}d2\left({\mathrm{\&sigma;}}_H\right))+\sqrt{(1+2\mathrm{\&delta;})}d{3}^{T}d3\left({\mathrm{\&sigma;}}_H\right)]$ $\left(4\right)$

$\frac{{\&PartialD;}^2{\mathrm{\&sigma;}}_v}{\&PartialD;t^2}=v_p^2[\sqrt{(1+2\mathrm{\&delta;})}(d{1}^{T}d1\left({\mathrm{\&sigma;}}_H\right)+d{2}^{T}d2\left({\mathrm{\&sigma;}}_H\right))+d{3}^{T}d3\left({\mathrm{\&sigma;}}_H\right)]$

(b) by dispersion relation, derive qP wave equation;

Flecter (2009) proposes to utilize following dispersion relation to derive corresponding limited shear wave equation:

${\mathrm{\&omega;}}^4=[(v_{\mathrm{px}}^2+v_{\mathrm{sz}}^2)(k_x^{\&prime;2}+k_y^{\&prime;2})+(v_{\mathrm{pz}}^2+v_{\mathrm{sz}}^2)k_z^{\&prime;2}]{\mathrm{\&omega;}}^2-v_{\mathrm{px}}^2{v}_{\mathrm{sz}}^2{(k_x^{\&prime;2}+k_y^{\&prime;2})}^2-$ $\left(5\right)$

$v_{\mathrm{pz}}^2{v}_{\mathrm{sz}}^2{k}_z^{\&prime;4}+[v_{\mathrm{pz}}^2(v_{\mathrm{pn}}^2-v_{\mathrm{px}}^2)-v_{\mathrm{sz}}^2(v_{\mathrm{pn}}^2+v_{\mathrm{pz}}^2)](k_x^{\&prime;2}+k_y^{\&prime;2})k_z^{\&prime;2}$

Introduce after suitable intermediate wave field variable, the limited shear wave equation of its correspondence is:

$\frac{{\&PartialD;}^{2}p}{\&PartialD;t^2}=[v_p^2(1+2\mathrm{\&epsiv;})H_2+v_{\mathrm{sz}}^2{H}_1]p+[v_p^2-v_{\mathrm{sz}}^2]H_{1}q$ $\left(6\right)$

$\frac{{\&PartialD;}^{2}q}{\&PartialD;t^2}=[(1+2\mathrm{\&delta;})v_p^2-v_{\mathrm{sz}}^2]H_{2}p+(v_p^2{H}_1+v_{\mathrm{sz}}^2{H}_2)q$

Wherein,

$H_1={\sin }^2\mathrm{\&theta;}{\cos }^2\mathrm{\&phi;}\frac{{\&PartialD;}^2}{\&PartialD;x^2}+{\sin }^2\mathrm{\&theta;}{\sin }^2\mathrm{\&phi;}\frac{{\&PartialD;}^2}{\&PartialD;y^2}+{\cos }^2\mathrm{\&theta;}\frac{{\&PartialD;}^2}{\&PartialD;z^2}+$

${\sin }^2\mathrm{\&theta;}\sin 2\mathrm{\&phi;}\frac{{\&PartialD;}^2}{\&PartialD;x\&PartialD;y}+\sin 2\mathrm{\&theta;}\sin \mathrm{\&phi;}\frac{{\&PartialD;}^2}{\&PartialD;y\&PartialD;z}+\sin 2\mathrm{\&theta;}\cos \mathrm{\&phi;}\frac{{\&PartialD;}^2}{\&PartialD;x\&PartialD;z}$

$H_2=\frac{{\&PartialD;}^2}{\&PartialD;x^2}+\frac{{\&PartialD;}^2}{\&PartialD;y^2}+\frac{{\&PartialD;}^2}{\&PartialD;z^2}-H_1$

(c) with the stable qP wave equation of regularization term

In isotropy ACOUSTIC WAVE EQUATION solution procedure, Liu (2008) has proposed a kind of optimization method of anti-frequency dispersion, by revising dispersion relation corresponding to ACOUSTIC WAVE EQUATION, obtains anti-dispersion equation, and its basic thought is as follows:

For ACOUSTIC WAVE EQUATION in isotropic medium:

$\frac{1}{c_0^2}\frac{{\&PartialD;}^{2}p}{{\&PartialD;}^{2}t}={\&dtri;}^{2}p---\left(7\right)$

Plane wave solves form:

$p=a{e}^{i(\mathrm{\&omega;t}\&PlusMinus;\stackrel{\&RightArrow;}{k}x)---\left(8\right)}$

For being with ACOUSTIC WAVE EQUATION in regularization term isotropic medium:

$\frac{1}{c_0^2}\frac{{\&PartialD;}^{2}p}{{\&PartialD;}^{2}t}={\&dtri;}^{2}p+\frac{\mathrm{\&sigma;}}{v}\frac{\&PartialD;}{\&PartialD;t}{\&dtri;}^{2}p---\left(9\right)$

Its corresponding plane wave solution is:

$p=e^{-\frac{\mathrm{\&sigma;}}{2}{k}^2{c}_{0}t}{e}^{i\left[\sqrt{1-\frac{{\mathrm{\&sigma;}}^2{k}^2}{4}}\right]\mathrm{\&omega;t}\&PlusMinus;\stackrel{\&RightArrow;}{k}x]}---\left(10\right)$

Notice:

$1-\frac{{\mathrm{\&sigma;}}^2{k}^2}{4}\&ap;1---\left(11\right)$

(10) can write:

$p=e^{-\frac{\mathrm{\&sigma;}}{2}{k}^2{c}_{0}t}{e}^{i[\mathrm{\&omega;t}\&PlusMinus;\stackrel{\&RightArrow;}{k}x]}---\left(12\right)$

In solving TTI medium, when qP wave equation, introduce regularization term, high wave number composition useless in unstable solution is done to suitable decay, obtain more stable equation.The limited shear wave equation proposing for Flecter, adds following regularization term:

Succinct in order to describe, (6) are rewritten as:

$\frac{{\&PartialD;}^{2}p}{\&PartialD;t^2}=G_p^{1}p+G_p^{2}q$ $\left(13\right)$

$\frac{{\&PartialD;}^{2}q}{\&PartialD;t^2}=G_q^{1}p+G_q^{2}q$

Add after regularization term, (9) can become:

$\frac{{\&PartialD;}^{2}p}{\&PartialD;t^2}=G_p^{1}p+G_p^{2}q+\mathrm{\&sigma;}\frac{\&PartialD;}{\&PartialD;t}(G_p^{1}p+G_p^{2}q)$ $\left(14\right)$

$\frac{{\&PartialD;}^{2}q}{\&PartialD;t^2}=G_q^{1}p+G_q^{2}q+\mathrm{\&sigma;}\frac{\&PartialD;}{\&PartialD;t}(G_q^{1}p+G_q^{2}q)$

Wherein σ is regularization coefficient, can verify, solves (14) than the labile element that solves high wave number in (13) wave field still less, thereby can make TTI medium RTM applicability when processing real data stronger.For equation (4), also can adopt similar thought to build it with the form of regularization term.

(B) adaptive regularization parameter is selected: from (14), can find out, choosing of regularization parameter has a great impact final wave field result of calculation, in order to make algorithm more sane when processing real data, utilize the adaptive calculating regularization of following mode coefficient:

$\mathrm{\&sigma;}=\max \{\mathrm{sqrt}\left({\&dtri;}^2\mathrm{\&theta;}\left(\stackrel{\&RightArrow;}{x}\right)\right),\mathrm{sqrt}\left({\&dtri;}^2\mathrm{\&phi;}\left(\stackrel{\&RightArrow;}{x}\right)\right)---\left(15\right)$

Wherein, θ, φ are respectively inclination angle and the position angle on stratum.