From e461a26fe9b898acb2cbec7ebd7caa29cfae9002 Mon Sep 17 00:00:00 2001 From: sfomel Date: Thu, 4 Apr 2024 12:51:04 -0500 Subject: [PATCH] sage --- book/tccs/owe/SConstruct | 3 + book/tccs/owe/Sage/Fig/ostolt.pdf | Bin 0 -> 11798 bytes book/tccs/owe/Sage/ostolt.sage | 8 + book/tccs/owe/owe.bib | 74 +++++ book/tccs/owe/paper.tex | 447 ++++++++++++++++++++++++++++++ 5 files changed, 532 insertions(+) create mode 100644 book/tccs/owe/SConstruct create mode 100644 book/tccs/owe/Sage/Fig/ostolt.pdf create mode 100644 book/tccs/owe/Sage/ostolt.sage create mode 100644 book/tccs/owe/owe.bib create mode 100644 book/tccs/owe/paper.tex diff --git a/book/tccs/owe/SConstruct b/book/tccs/owe/SConstruct new file mode 100644 index 0000000000..f2ec5346ac --- /dev/null +++ b/book/tccs/owe/SConstruct @@ -0,0 +1,3 @@ +from rsf.tex import * + +End(use='amsmath') diff --git a/book/tccs/owe/Sage/Fig/ostolt.pdf b/book/tccs/owe/Sage/Fig/ostolt.pdf new file mode 100644 index 0000000000000000000000000000000000000000..3e97c22ce019449259bd4bc908d4986e59cd7502 GIT binary patch literal 11798 zcmb_?2{@Er^gkgQOV%P{kR;3OVg4lbxx4a5U&(3RhKyQN0|Yj)dO^XnT9n;23%bgr$*# z8kA$g8@Smp6DqLE}te zzcd5A{-H{Js+Tj(1&&%;Ufb0J2oFbSdjO4SQYqe!R7eItnlF{?2@5=Z&)UFGM3L*> zxegFNTna1mLrOLCytEZVSq(4cDD-oQc$JwpX-MR&52;N|35RzxE?YJ1w+so|D9QObo* z!oJw&>`xwlh}&t`q0rbN*I1F0IyckXqGhwi^j=6*W2M~oBBM9c1z-0Lo*GR*TyCWr zdZGmtywE)pP{Uxfx2^o}Y-8iqj`zVqL9_4L`-}#UE?l{qwp;JiLk&&CE$L{}(N3#@ zTctMWnv8pmiRTR!f;k6tMRpja%`htKopA`LXpYp8T$9*6k860QT6&ta6`Lg1G(IF- zR&7IhEdG-6X7;MXfk!hFNrTlDInJ#Z+7{;6Ow%`92euy<)eq@3KC1VTM6NoY(?`_t zmk~%HQZg50xISf&?k5*lAMHVAITIs~+X$Sai08jU^&1^|A4d9WD0z{YK7y^R|8ph( z-V(|Bb9SYB37WY*5~@~QOhIX+YBMcN7P59mnlR>ZB%YByjqB_>11?duZJQ1wxdXTO zktLNwweLufYI*PRh^k5Ii|S#rN(4Xg-1?xOpsZzfV|5dZm0*q*l3m-msj#(`-KvYp zB`W96$t2u+2J5bP+#?gYJ}&F7q%%_|p$7SRp-8LmNS44A-SJbqWL2`XPp=u#ijg{X zWA&>zj7Yo!-10Nytq+V`k8ak*t-T~ryW?|Rr^O+)qQL0H6M36Qn6|}MJ=wy^Hwiki zge0&LWL2{)Y>dgo9yh+x&<&qEPq4Ussqd`=ubI8mt-}SxjAfts>`ypRI2CIBIW&^H z;C9=NS8?us+Pu;uOl~!pwqYiNHr7{h>or2;99qU>toNzU@2I&9EqN!v`NSpz|MNN8 zBExwFr_HX(03zW*d!#0{28{6ASG>6xOr|iL-!8LecuRxmbv#F;b%REUrsE!okaBB3 zpW>miX=(4nm3L%Y0>im;4C2b2 z3^y@R#bMJ<8^0r&hUD$t`NfqYxV9rB>INd>*o<_VI|>WLoQ$on?cSt`iQqanB*inG zUn0bB%N%vhI!?Lt#;EdfJ#3IYQeYdv}<^VZz4G8`yY=Z^i^W z$CbmqKbw%A_3syb_i4YV?vq3IaXTe9<{Z;n-dF@w_yow4z9p zmBra%K5)M3hQ{kliB$9p#xvCs zEm?^dED;=Ig=+Hs{>~KsZOXDf>_Pcp%|p3O;T6$X8IjKO1FjDaUcQhx@cmh|Wui*s zkX1#Qb@7YokSAY0!`$!J1nW(G?5*lMRW`VA=556xOQZ5BHN3A#mvU*gtv;k{F>{jPUSvawL-1^>v5>PhQ;4qJ z@OM-B*^r>?>?hsaZUwPdZYjQhC1~q`&Dv9#={D2xbs;7Wb&8^=@V-*ImdlP z#*XsJ)}#)Z?2(#Yhs#)GX)?aS~|M@8IAc_+E_02oa4}iKmQygE|+`bu1 zNwT|M9Tc)i#@v-^mVM~O2N<-Ql`FO*?LRd|jonzvLMl|;SSCO9;@jIvi)eH2rhq+T z)T%)}F*nPvP64-B_6Fd%`U1Q4=DKWaSI@0?xe(wtu<2IL*yS+Xh}5>nBE^*O+nw)U zZBS5Y{J@);Kg=%GY)W(1doxZp-NK*~`M~^g2%_zc!JaWM5QtUR4TTQ`uDhFbi_Kjw%= z`8DUb=g;=V5!NCGr;v{Ew>I{ClJ$6$!M4tjCpG?V=4q905teFJi7@*p{>e*Ogrvka zLy?KCJ6=B!(mtNGnKC86>zmuEH!!oOMD=a!8YH4#<=!PU@n7`}wbP4zk)6ndFJduV z%QLKAXR@VTef8-1sdd91aTu2wwsb`iF8yI;?X@K&&AC)yG?K3q+3-9LE0q~SZ2B&C_G`@C|f%jr`FRWKP4JH$9qox{@ z`h$Q6z&SO(mBA(j=G z{~eky0Z=0Pzj;zQ_JH#~l8tFW{{b8Ks_NNY=3<}bJ2vwSu!^m26|^Vcz;!-oy`TQ# zt>V30t>m|C7Pqw8jn9_xU$C^@aY6r~QKr&iQ<})I57)bn<{pAUcmLkWrtgevcZ#yF zz~c%a8G}^(2PTvByqHiNq)F8%1@8MZt!Kw`x0GhQKm1Nw>QkIZ>pCsBv{HR5GtuMN zEy}T()^c7BFjUsJ#4pjzAxjcLL#lg9Fbf=rXt%qsBK`~`JP4oVK z)*Z1Dl=Ss|-Rf&T?R@BB9h16Wi%EUt%gx6J*NnV-?^bu6H9zB-l8n=}L7J(8*ZL#k zqIF*`jrHWbao1>^x+W_%q4Oph<10Pl@~HP!-G}4z?~gI-C0=n6AMd>C)vRD&ZnjO! z!O_sS(DkUf;iNJD%ez(gIJm!hXpv<6PP`u2TlI|HJJ zsY#XGk#e9&u0o)uh%s$bNw}QNAhFi72;htBu$@cr@d_9op`a8`LS6^$ER>!&w?4h z*+O>D*|&U#0S5-Rks+;^exVv zu@4?>Y2q5XJM8?yT-?RBJ60f62EoGcZ$9 z@|siaJ8MU{Q>U4JP35r3bJ!kt{e%q=PVQH{@D z9|j{H!V5^I@Ml45V`xJc@hAAUZ9I+HPFAEO+}A>_KK+SyQH6CgarPlB(odK=WT7J+ zds!l30FxkdcI&}&f$QXk>+pM|lT<#!yG%`E>o1&3dA!~8aSCjk({pAy-n*TSdcF1b99lu7DaXnD*PMPy!cH7PkQ!;*b1DpE{NW$J{upW=7nknGkn ztDvO%)`g5d2{V%^(@P2t8n!XQjpupBRL;BeS{2{>Oq9#X^s8VMi@L*%InPK5)=Xsc&T z!<_Y;*7hy=!X`(FDtKI*yNSIfY?QK}jF znE@k78bMtKN)ipCuw1l~yzH%m)dOB(Guxk-j(xA6oZ5JG?bA7%(<^8Yx}00CK^)<~ zZBXlko-Z?;O>^>WxMb+uN9E+5RMg!CJm)hg#xpV<=g9rxwL`|tAK>Z9M7`qiX; zI>0I=cn231_NLHJ#G~`EXO*C|_Nk%~uE|)qk+j*4nqheB7_0QL;M4+HpEDlcp6emG zn7?cAWqdrUHpmh7=&9MnN8X|C?ts29d+!50I<4GwC&UFeh6kLnS@_UAeQ#l3v{K_O zp(`SzANFT72h{F$ND}a4JlZC6sn1=(x_1&WHyVpHL*jAfbvtu`uW`2w!)SW*&%bXtCil9HMh@{ zIjOB_Gb6q&P}wd9+xBe*Syuqbgntn=z6=7%ng=i1Gly*crfQP%@O>pGx0Kf%>2?*H zjcYPxyf9Uy)n_^6qHXtuO{5{&ZFU$ZoEgd?zpfZyKKAvZ>GYmg3;D zM*A8h{h;*c*$T4){6yjY1rQ`Y;RNE5Mo)&bGdDF}S2l{vIzn4;m&#+|$UxRhx2}sg z5l_;$AIhL~N$oN!Qr#qZyHuCRoo}qjY~-TPfjOP4*O~gMf~tHwkXZDJ<;?NAb}LV; zrqEO6=$yPE0mB)DU-G3?Au3m0p;^0O-Q#e9jP|w$#fr>W9gt|V2`Q70~F$=}tau!V4 z%3=7raqI3`L>~3c)@y;tGZNdcpLEIFTPhMbD8s7prY*HBNZf#JCSEqHH$jgLcWA?ijF0Eq|K3VwW)`OXaeAeb2 z*`HQnMM~EhJPzV_QT5t1h@Z`fvJ z?%I08wCkp8&BMx8>(_Of$9#^>#eQ17YtCrECS}YZ;%X9VMsQcTQbfG&l$}3qW|qj<;{4W%M*HZ7NmvV^cAQf`Tyby(L05oBsDFst1#rd|LnI;V%1{Q!QR<_}3pg%451`KCr>|Ie<>1kzfZaHa;eK`Wf6(hKWy$J+(%it zY2PLVnqz-_g6t`g?)x|3ln>7d3-xzYt5NFh+47d@0ky3SynYYSs8uD}^{vX;_t)0i zu<<&ki&fl*OV}n(TZ~&JPc+_EYELY^>9kiR`+n>YTg^4b-J#!0Ys2&&PTRfwn7&nt z!gu6L*nos7Ug#WPi&0JM(M9_Wm$xrE z)70arypY=EKE+?7CZ$-HY0q^b_r=B?CxRRZf;OTzg%ds$_s{T@B{>B*J<=$Z zatS?HXy)2>aUWmxMjXq4%)sfKa_+n5@}1vcliIO2lKk#)Rvv6LYgx!ED1S1|a&N!F z{0dU704~s|f6z3}gvbRj-q57LyxOLCM4(ajmZ$TTPvf6kh(551zL_GP`NcH(l6a4q zYEY)a>C(2A1~+V9+D&=C&^smNh@zzXu}(1&4_q_uFrHLqB(*%hWI|dm+Vtvz*_LZ< z%&Kv>z1Fr2}{wKo`4sZkoj&Ot{ zs34aDTy%h=kQW?<0!rBjM*xkXKsSFlA^@C(z}I>)u;p3B+z!N16+eW-;yLSM~*6 z!G9jOmTdujOa2#q1Li!)&VE4Rf6~0V9|cM-A)*1>&?Z9<@@NcD1&QpVL!~PQ2}4-^ zXoA!gkW}MJQTKB8phC7{Mx%OKKn(+OAYCFL5ElJ=UdsraA-K1h=YLw$K&vD91ctZ&aaBVAQgeh z;gO&X)zEtkC?IJ_5sZiCBoc_AUZg=aT^Nia45}d}Afcdw2L=a?#R0(==Ky`dy7V>% z3lu@17n~xIUK2s=#sfMYg$4pHu8Rfh;(#qJvq~Bi+ zz~JH(Xdp542I2~E{=*D@Y}bz(a0~mv2%RgeA_@|l&OT%Zf3OP4L^md2P|&EKvp`HC zkudtZr3r|5BK${z3}Mj*AUpW8&~**$0@1z;2uDK)MnI#&!6t*& zl9NU+1u0KB*mHUZFJR1zEy$&#m!|&}34jBn5C1J7(9ii`?xmnWKNf?>eNZ%^&*=*c ze6h6{TY$9mIiPp}B&GN92Oeax6#xfc1D$pN4)%!Nq6ZnwqJ+O9zz>fJ#WihUc8gQ} zdjI?RM;-Nh`1_fwh#;wh6%Z7lKgbs&)Y1QIHxQDRy8Y(0e*+P}*8Od>ei47iz1k@G zU!X$X)78r#q+KqG|7*bSNBl*6vY4NV&~tT!Vl}W+NH_F+MSl>E=mISLtPP^~FK|Nw zIfoD;*#e-0kB2wS!_@&EfKebQVC3K~G@6g!HUz@+=PLzoUuPK@0ze2ye+o!p{C$*< zqZ8bLOmPRREX@FV(x-Rz_R@rMKeuRZLnFbZB=|$2aiArS#7QHO(lWp4EYO$g1OpkE zC>ZcVOaFiy#{pAtg8!sJ=mq)*_gbQ%kaz+Rc^M6h16Hw|M*O=Ui)g^=Ddnj**#{)>kG zn?FSyNDNs%mH_!(=-*F!Mk*3P%w9$V=m2=FKl*{t4i2L~X*6H5s|VGWp1*D88U&C! sa7?D&-jD;J`$}CeCvOmO=?}o0_|eF|H2NopLK27sn3R;3kv8mq06SZ2*Z=?k literal 0 HcmV?d00001 diff --git a/book/tccs/owe/Sage/ostolt.sage b/book/tccs/owe/Sage/ostolt.sage new file mode 100644 index 0000000000..a06615083f --- /dev/null +++ b/book/tccs/owe/Sage/ostolt.sage @@ -0,0 +1,8 @@ +k = var('k') +def ost(t, w): + return plot(w*cos(t) - k*sin(t), k, -2, 2) +def env(w): + return plot(sqrt(w^2 + k^2), k, -2, 2, thickness=5) +lines = [ost(t*pi/20,1) for t in range(-10, 11)] +plot = sum(lines)+env(1) +plot.save(frame=True, axes=False, axes_labels=['Wavenumber (rad)', 'Frequency (rad)'], filename='junk_sage.pdf') diff --git a/book/tccs/owe/owe.bib b/book/tccs/owe/owe.bib new file mode 100644 index 0000000000..4351336055 --- /dev/null +++ b/book/tccs/owe/owe.bib @@ -0,0 +1,74 @@ +@Article{pnas, + author = {S[ergey] Fo\-mel and J[ames] A[] Sethian}, + title = {Fast phase-space computation of multiple + arrivals}, + journal = {Proc. of the Nat. Acad. of Sciences}, + year = 2002, + volume = 99, + pages = {7329-7334} +} + +@Article{jcp, + author = {S[ergey] Fo\-mel and J[ames] A[] Sethian}, + title = {Fast Algorithms for Computing +Multiple-Arrival Multiple-Source Non-Viscosity +Solutions of Static {Hamilton-Jacobi} Equations}, + journal = {Journal of Computational Physics, submitted}, + year = 2003 +} + +@Book{maslov, + author = {V[] P[] Maslov and M[] V[] Fedoriuk}, + title = {Semi-classical approximations in +quantum mechanics}, + publisher = {Reidel}, + year = 1981, + address = {Dordrecht} +} + +@Book{norm, + author = {N[] Bleisten and J[] K[] Cohen and J[] W[] Stockwell}, + title = {Mathematics of multidimensional seismic imaging, migration, and inversion}, + publisher = {Springer}, + year = 2001 +} + +@Book{iei, + author = {J[on] F[] Claerbout}, + title = {Imaging the {Earth's} Interior}, + publisher = {Blackwell}, + year = 1985 +} + +@Book{santo, + author = {J[] A[] DeSanto}, + title = {Scalar wave theory: {Green's} functions and applications}, + publisher = {Springer-Verlag}, + year = 1992 +} + +@Article{curvelet, + author = {E[mmanuel] J[] Cand\'{e}s and L[aurent] Demanet}, + title = {The curvelet representation of wave propagators is optimally sparse}, + journal = {Communications on Pure and Applied Mathematics}, + year = 2005, + volume = 58, + pages = {1472-1528} +} + +@inproceedings{seislet, + author = {S[] Fomel}, + booktitle = {76rd Ann. Internat. Mtg.}, + pages = {2847-2850}, + publisher = {Soc. of Expl. Geophys.}, + title = {Towards the seislet transform}, + year = {2006} +} + +@Book{sethian, + author = {J[] A[] Sethian}, + title = {Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science}, + publisher = {Cambridge Univ. Press}, + year = 1999 +} + diff --git a/book/tccs/owe/paper.tex b/book/tccs/owe/paper.tex new file mode 100644 index 0000000000..2ee7f6a00f --- /dev/null +++ b/book/tccs/owe/paper.tex @@ -0,0 +1,447 @@ +\title{Seismic imaging using the oriented wave equation} +\author{Sergey Fomel} + +\maketitle + +\begin{abstract} +Seismic waves propagate not only in space and time but also in +direction. The mathematical description of waves propagating in phase +space (space-time-direction) is provided by the oriented wave +equation. I derive this equation and apply it for seismic images using +synthetic and field data examples. The advantages of imaging in the +phase-space domain include the ease of separating different wave +components, the ease of forming angle gathers, and the ease of +handling seismic anisotropy. The disadvantage is the computational +expense of extra dimensions. +\end{abstract} + +\section{Introduction} + +Imaging in the angle domain has become a recurrent theme in the recent +research on seismic imaging. Accounting for the illumination (wave scattering) +angles provides an illumination compensation in Kirchhoff migration +\cite[]{SEG-1999-13581361,SEG-1999-13621365,SEG-2002-11881191,SEG-2003-09210924}. +Formulating the migration operator directly in the angle coordinates leads to +a powerful imaging method with the ability to handle multiple-path +arrivals~\cite[]{SEG-1998-1538,SEG-2002-11961199,GEO68-01-02320254}. It is also +possible to extract angle-domain common-image gathers from the output of +wave-equation (do\-wnward extrapolation) imaging +\cite[]{SEG-1999-08240827,GEO67-03-08830889,SEG-2002-13601363,GEO68-03-10651074}. + +Studying the problem of multiple-arrival traveltime computation, +\cite{pnas,jcp} realized that the regularity of the phase space, +previously explored in the theoretical studies of the asymptotic wave +propagation \cite[]{maslov}, can be used for constructing an efficient numerical +algorithm. The algorithm is based on a finite-difference discretization of +the set of specially constructed static partial differential equations +(``escape equations'') for the single-valued phase-space traveltime and +location functions. The algorithm of \cite{pnas} serves as a fast +alternative to tracing multiple rays from subsurface imaging point locations +for common-angle Kirchhoff-style imaging. However, its connection with the +wave-equation methods is not immediately apparent. + +In this paper, I extend the theory of the phase-space computations to handle +the wave extrapolation problem. Applying a simple Huygens-principle +construction, I derive a linear partial differential equation that describes +the wave propagation effects in the phase space, where propagating waves are +assigned not only with the time and space coordinates but also with +orientation directions. The oriented wave equation complements the previously +found set of escape equations. It belongs to the transport type and possesses +other attractive theoretical properties. This equation can form the +theoretical basis for angle-domain seismic imaging methods. + +\section{Traveltimes and waves in the physical space} + + +To derive the oriented wave equation, I follow the analogy with the +conventional wave equation operating in the physical space. The first step is +to establish a connection between the traveltime and wave descriptions. In an +isotropic medium, the traveltime $T(\mathbf{x},\mathbf{y})$ from source +$\mathbf{y}$ to a physical point $\mathbf{x}$ is governed by the eikonal +equation +\begin{equation} + \label{eq:eikonal} + \left|\nabla_x T\right|^2 = n^2(\mathbf{x})\;, +\end{equation} +where $n(\mathbf{x})$ is the slowness at $\mathbf{x}$. Let us consider the +problem of extrapolation the wavefield recorded on the surface of the Earth to +some subsurface location. The Huygens principle implies that the continued +wavefield is a superposition of waves from the surface sources. We can express +this superposition formally as +\begin{equation} +\label{eq:cont} +W^{(\pm)}(\mathbf{x},t) = \int\limits_{\partial \mathcal{D}} +A(\mathbf{x},\mathbf{y})\, +\widehat{W}\left(\mathbf{y},t \mp T(\mathbf{x},\mathbf{y})\right)\, +\mathbf{d y}\;, +\end{equation} +where $\mathbf{y}$ is a point on the surface $\partial \mathcal{D}$, $t$ is +time, $\widehat{W}$ is the surface wavefield, $W$ is the wavefield in the +subsurface $\mathbf{x} \in \mathcal{D}$, and $A$ is the appropriate amplitude +function. The $\pm$ sign corresponds to the two possible wave extrapolation +schemes: waves traveling inside or outside $\mathcal{D}$. The sign defines the +time delay of the secondary sources on the surface. + +The integral representation~(\ref{eq:cont}) is commonly associated with the +high-frequency asymptotic solution of the wave equation and involves +additional filtering to balance the spectra of $W$ and $\widehat{W}$ +\cite[]{norm}. In the high-frequency regime, it is equally possible to use it +to derive the wave equation. Start by rewriting equation~(\ref{eq:cont}) in +the frequency domain: +\begin{equation} +\label{eq:contw} +\mathcal{W}^{(\pm)}(\mathbf{x},\omega) = \int\limits_{\partial \mathcal{D}} +A(\mathbf{x},\mathbf{y})\, +\widehat{\mathcal{W}}\left(\mathbf{y},\omega\right)\, +e^{\mp i \omega\,T(\mathbf{x},\mathbf{y})}\, +\mathbf{d y}\;, +\end{equation} +where $\omega$ is the temporal frequency, and $\mathcal{W}$ is the Fourier +transform of $W$. Differentiating both sides of equation~(\ref{eq:contw}) with +respect to the spatial coordinates $\mathbf{x}$ and retaining only the terms +of the highest order in $\omega$, we obtain +\begin{equation} +\label{eq:nabla} +\nabla \mathcal{W}^{(\pm)} = \mp i \omega\, +\int\limits_{\partial \mathcal{D}} +\nabla_{\mathbf{x}} T\,A(\mathbf{x},\mathbf{y})\, +\widehat{\mathcal{W}}\left(\mathbf{y},\omega\right)\, +e^{\mp i \omega\,T(\mathbf{x},\mathbf{y})}\, +\mathbf{d y} + \ldots +\end{equation} +Repeated differentiation produces +\begin{equation} +\label{eq:nabla2} +\nabla^2 \mathcal{W}^{(\pm)} = - \omega^2\, +\int\limits_{\partial \mathcal{D}} +\left|\nabla_{\mathbf{x}} T\right|^2\,A(\mathbf{x},\mathbf{y})\, +\widehat{\mathcal{W}}\left(\mathbf{y},\omega\right)\, +e^{\mp i \omega\,T(\mathbf{x},\mathbf{y})}\, +\mathbf{d y} + \ldots +\end{equation} +According to the eikonal equation~(\ref{eq:eikonal}), the term +$\left|\nabla_{\mathbf{x}} T\right|^2$ does not depend on $\mathbf{y}$ and can +be taken out of the integral sign. The final result is the celebrated wave +equation, commonly used in many different seismic imaging algorithms: +\begin{equation} +\label{eq:nabla} +\nabla^2 \mathcal{W} += - \omega^2\,n^2(\mathbf{x})\,\mathcal{W}(\mathbf{x},\omega) ++ \ldots +\end{equation} +This simplified derivation highlights the formal connection between +traveltimes (wavefronts) and waves in the high-frequency asymptotic +description of the wave propagation. + +\section{Traveltimes and waves in the phase space} +A major problem with the integral wavefield +continuation~(\ref{eq:cont}-\ref{eq:contw}) follows from the fact that, in +complex media, the traveltime $T(\mathbf{x},\mathbf{y})$ can take multiple +values because of the possible multiplicity of ray paths between $\mathbf{x}$ +and $\mathbf{y}$. The remedy is to rewrite the wave continuation integral in +the angle coordinates \cite[]{SEG-1998-1538}, as follows: +\begin{equation} +\label{eq:pscont} + \mathcal{W}^{(\pm)}(\mathbf{x},\omega) = + \int\limits_{|\mathbf{p}|=n(\mathbf{x})} + \widehat{A}(\mathbf{x},\mathbf{p})\, + \widehat{\mathcal{W}}\left(\widehat{\mathbf{y}}(\mathbf{x},\mathbf{p}), + \omega\right)\, +e^{\mp i \omega\,\widehat{T}(\mathbf{x},\mathbf{p})}\,\mathbf{d p}\;. +\end{equation} +Here $\widehat{\mathbf{y}}$ and $\widehat{T}$ are the escape location and the +escape time, uniquely defined by the initial location $\mathbf{x}$ and the +initial direction $\mathbf{p} = \nabla_x T$. They correspond to the location +and traveltime of the ray started at $\mathbf{x}$ and traced in the initial +direction $-\mathbf{p}$ at the ray intersection with the surface $\partial +\mathcal{D}$. By virtue of ray tracing, the escape functions are strictly +single-valued under the assumption of differentiable slowness $n(\mathbf{x})$. +The integral in equation~(\ref{eq:pscont}) is connected with the integral in +equation~(\ref{eq:contw}) via a change of variables. The change of variables +from $\mathbf{p}$ to $\mathbf{y}$ is possible only when the two-way traveltime +function $T(\mathbf{x},\mathbf{y})$ is single-valued \cite[]{GEO52-07-09430964}. + +\cite{pnas} derive the following set of partial differential equations +for describing the wavefront propagation in the phase space: +\begin{equation} +\label{eq:esct} +\mathbf{p} \cdot \nabla_{\mathbf{x}} \widehat{T} + +n(\mathbf{x})\,\nabla n \cdot \nabla_{\mathbf{p}} \widehat{T} + n^2(\mathbf{x}) += 0 +\end{equation} +\begin{eqnarray} +\label{eq:escy} + \mathbf{p} \cdot \nabla_{\mathbf{x}} \widehat{\mathbf{y}} + +n(\mathbf{x})\,\nabla n \cdot \nabla_{\mathbf{p}} \widehat{\mathbf{y}} +& = & 0 \\ +\label{eq:escp} + \mathbf{p} \cdot \nabla_{\mathbf{x}} \widehat{\mathbf{p}} + +n(\mathbf{x})\,\nabla n \cdot \nabla_{\mathbf{p}} \widehat{\mathbf{p}} +& = & 0 +\end{eqnarray} +The escape equations~(\ref{eq:esct}-\ref{eq:escp}) lift the eikonal +equation~(\ref{eq:eikonal}) from the physical space $\mathbf{x}$ to the +higher-dimensional phase space $\{\mathbf{x},\mathbf{p}\}$. They are linear, +transport (convection-type) equations, assuming only single-valued solutions +and suitable for a fast numerical computation. The computation produces +surface arrivals simultaneously from all possible subsurface locations and all +possible initial directions \cite[]{jcp}. + +Figure~\ref{fig:symes5} shows an example output of the multiple-arrival +computation algorithm. The 2-D slowness model has a smooth slow region +surrounded by fast regions. This causes five distinct arrivals for both point +sources and plane-wave sources located at the surface. Figure~\ref{fig:slices} +displays horizontal slices from the actual computational output: the solutions +of the escape equations. The multi-valued traveltime function for each +individual source can be extracted from these solutions by simple +interpolation. + +%% \sideplot{symes5}{width=\columnwidth}{Multiple-arrival traveltimes. Left: +%% wavefronts extracted for a point source at the surface. Right: wavefronts +%% extracted for a plane-wave source at the surface. The slowness model is +%% shown in the background.} \sideplot{slices}{width=\columnwidth}{Horizontal +%% slices out of the phase-space multiple-arrival solution. The slices +%% correspond to 1.8~km depth in Figure~\ref{fig:symes5}. Left: escape +%% location. The curve shows the part of the solution corresponding to the +%% point source result in the left plot of Figure~ \ref{fig:symes5}. Right: +%% escape direction. The curve shows the part of the solution corresponding to +%% the plane source result in the right plot of Figure~ \ref{fig:symes5}. Five +%% branches of the solution are identified by projecting the curves to the +%% horizontal axes.} + +To establish the path from traveltimes to wavefields in the phase space, I +separate the integrand of equation~(\ref{eq:pscont}) and call it an +\emph{oriented wave}: +\begin{equation} +\label{eq:ow} +\mathcal{V}^{(\pm)}(\mathbf{x},\mathbf{p},\omega) = + \widehat{A}(\mathbf{x},\mathbf{p})\, + \widehat{\mathcal{W}}\left(\widehat{\mathbf{y}}(\mathbf{x},\mathbf{p}), + \omega\right)\,e^{\mp i \omega\,\widehat{T}(\mathbf{x},\mathbf{p})} +\end{equation} +Differentiating equation~(\ref{eq:ow}) with respect to both phase-space +components $\mathbf{x}$ and $\mathbf{p}$ and retaining only the terms with the +highest order of frequency $\omega$, we obtain +\begin{eqnarray} + \label{eq:psnabla1} + \nabla_{\mathbf{x}} \mathcal{V}^{(\pm)} & = & + \mp i \omega\,\nabla_{\mathbf{x}} \widehat{T} \mathcal{V}^{(\pm)} + \ldots + \\ + \label{eq:psnabla2} + \nabla_{\mathbf{p}} \mathcal{V}^{(\pm)} & = & + \mp i \omega\,\nabla_{\mathbf{p}} \widehat{T} \mathcal{V}^{(\pm)} + \ldots +\end{eqnarray} +Equations~(\ref{eq:psnabla1}-\ref{eq:psnabla2}) can be combined with the help +of the escape equation~(\ref{eq:esct}) to produce the following equation +formally independent of the traveltime: +\begin{equation} + \label{eq:owew} + \mathbf{p} \cdot \nabla_{\mathbf{x}} \mathcal{V}^{(\pm)} + + n(\mathbf{x})\,\nabla n \cdot \nabla_{\mathbf{p}} \mathcal{V}^{(\pm)} + = \pm i \omega\,n^2(\mathbf{x})\, +\mathcal{V}^{(\pm)}(\mathbf{x},\mathbf{p},\omega) + \ldots +\end{equation} +Equation~(\ref{eq:owew}) is the oriented wave equation: the analog of the wave +equation for propagating oriented waves in the phase space. After the inverse +Fourier transform, it takes the following form in the original time domain: +\begin{equation} + \label{eq:owe} + \boxed{ +\mathbf{p} \cdot \nabla_{\mathbf{x}} V^{(\pm)} + + n(\mathbf{x})\,\nabla n \cdot \nabla_{\mathbf{p}} V^{(\pm)} + = \pm n^2(\mathbf{x})\,\frac{\partial V^{(\pm)}}{\partial t}\;.} +\end{equation} +The oriented wave $\mathcal{V}$ ``knows'' not only its position in space +$\mathbf{x}$ but also its directional orientation $\mathbf{p}$. According to +equation~(\ref{eq:pscont}), the usual propagating wave $\mathcal{W}$ is simply +a superposition of waves oriented in all possible directions or, in other +words, a projection of the oriented waves from the phase space to the +observable physical space. Additional weighting may be required to account +for the amplitude and spectrum differences. The dimensionality analysis and +the analogy with the Weyl angular-spectrum representation \cite[]{santo} suggest +the balanced connection +\begin{equation} + \label{eq:v2w} + \mathcal{W}(\mathbf{x},\omega) \sim i\,\omega\,\int + \mathcal{V}(\mathbf{x},\mathbf{p}_y,\omega)\, + \frac{d\,\mathbf{p}_y} + {\sqrt{n^2(\mathbf{x})-\mathbf{p}_y \cdot \mathbf{p}_y}}\;. +\end{equation} +In equation~(\ref{eq:v2w}), the direction vector $\mathbf{p}$ is +partitioned as $\{\mathbf{p}_y,p_z\}$ with +\[ +p_z=\sqrt{n^2(\mathbf{x})-\mathbf{p}_y \cdot \mathbf{p}_y}\;, +\] +and the integration over a sphere in~(\ref{eq:pscont}) is replaced with +integration over a plane. + +\section{Oriented waves in layered media} + +%\inputdir{Math} + +The connection between physical and oriented waves is particularly clear in +the case when the slowness function changes with only one spatial coordinate. +Let $\mathbf{x}=\{\mathbf{y},z\}$. If the slowness $n(\mathbf{x})$ is a +function of $z$, we can make a change of variables in equation~(\ref{eq:owe}) +from $p_z$ to $q = \sqrt{s^(z)-p_z^2}$ and then restrict $q$ to comply with +the eikonal equation~(\ref{eq:eikonal}). This results in the simplified +equation +\begin{equation} + \label{eq:owevz} +\sqrt{n^2(z)- \mathbf{p}_y \cdot \mathbf{p}_y}\, +\frac{\partial V^{(\pm)}}{\partial z} + +\mathbf{p}_y \cdot \nabla_y \partial V^{(\pm)} += \pm n^2(z)\,\frac{\partial V^{(\pm)}}{\partial t}\;. +\end{equation} +In the frequency-wavenumber domain, the solution to equation~(\ref{eq:owevz}) +for propagating waves in the $z$-direction is given by the simple phase-shift +operator: +\begin{equation} + \label{eq:oweps} + \widetilde{\mathcal{V}}^{(\pm)}(z+\Delta z,\mathbf{k}_y,\mathbf{p}_y,\omega) = + \widetilde{\mathcal{V}}^{(\pm)}(z,\mathbf{k}_y,\mathbf{p}_y,\omega)\, + e^{\,k_z(\mathbf{k}_y,\mathbf{p}_y,\omega)\,\Delta z}\;, +\end{equation} +where +\begin{equation} + \label{eq:owekz} + k_z(\mathbf{k}_y,\mathbf{p}_y,\omega) = + \frac{\pm \omega\,n^2(z) + \mathbf{k}_y \cdot \mathbf{p}_y} + {\sqrt{n^2(z)-\mathbf{p}_y \cdot \mathbf{p}_y}} +\end{equation} +The meaning of equation~(\ref{eq:owekz}) becomes clearer when we notice that +substituting $\mathbf{p}_y= \pm \nabla_y T = \mp \mathbf{k}_y/\omega$ reduces +it to the expression +\begin{equation} + \label{eq:wekz} + k_z(\mathbf{k}_y,\omega) = \pm \omega\,n(z)\, + \sqrt{1-\frac{\mathbf{k}_y \cdot \mathbf{k}_y}{n^2(z)\,\omega^2}}\;, +\end{equation} +and equation~(\ref{eq:oweps}) becomes the familiar phase-shift continuation in +the physical frequency-wavenumber domain \cite[]{GEO43-07-13421351}. Similarly, +the inverse transformation for the constant-slowness case +\begin{equation} + \label{eq:owest} + \pm \omega(\mathbf{k}_y,\mathbf{p}_y) = \frac{k_z\,\sqrt{n^2-\mathbf{p}_y + \cdot \mathbf{p}_y} - \mathbf{k}_y \cdot \mathbf{p}_y}{n^2} +\end{equation} +reduces to the familiar frequency-wavenumber relationship +\cite[]{GEO43-01-00230048} +\begin{equation} + \label{eq:west} + \omega(\mathbf{k}_y) = \pm \mbox{sign}{(k_z)}\frac{\sqrt{k_z^2 + + \mathbf{k}_y \cdot \mathbf{k}_y}}{n}\;. +\end{equation} +Equation~(\ref{eq:owekz}) decomposes these physical-space-domain expressions +into their individual phase-space components. Figure~\ref{fig:ostolt} shows +how the Stolt hyperbola~(\ref{eq:west}) gets composed from individual +phase-space planes~(\ref{eq:owest}). Figure~\ref{fig:omig} displays the +kinematic impulse response of the constant-velocity zero-offset migration to +show how, in the simplest of all possible situations, the energy gets +distributed from the physical space to the phase space. This distribution is +additionally demonstrated in Figure~\ref{fig:ogaz} by creating a decomposition +of the impulse response by partial stacking in the phase-space domain. + +\begin{figure}[htbp] +\centering +\includegraphics[width=\columnwidth]{Math/Fig/ostolt.pdf} +\caption{Stolt's frequency-wavenumber hyperbola + is the envelop of frequency-wavenumber planes from oriented plane waves.} +\label{fig:ostolt} +\end{figure} + +\begin{figure}[htbp] +\centering +\includegraphics[width=0.8\columnwidth]{Math/Fig/omig.pdf} +\caption{Kinematic impulse response of the + constant-velocity zero-offset phase-space migration.} +\label{fig:omig} +\end{figure} + +%\sideplot{ostolt}{width=\columnwidth}{Stolt's frequency-wavenumber hyperbola +% is the envelop of frequency-wavenumber planes from oriented plane waves.} +%\sideplot{omig}{width=0.8\columnwidth}{Kinematic impulse response of the +% constant-velocity zero-offset phase-space migration.} +%% \sideplot{ogaz}{width=0.8\columnwidth,height=1.6\columnwidth}{Decomposition of +%% the constant-velocity zero-offset migration impulse response by partial +%% phase-space stacks. Displayed from top to bottom are stacks from +%% $-60^{\circ}$ to $-30^{\circ}$, $-30^{\circ}$ to $0^{\circ}$ degrees, +%% $0^{\circ}$ to $30^{\circ}$, and $30^{\circ}$ to $60^{\circ}$.} + +Prestack downward continuation of individual phase-space components can be +accomplished by simultaneous continuation of sources and receivers: +\begin{equation} + \label{eq:owedsr} + k_z(\mathbf{k}_s,\mathbf{k}_r,\mathbf{p}_s,\mathbf{p}_r,\omega) = + \frac{\pm \omega\,n^2(z) + \mathbf{k}_s \cdot \mathbf{p}_s} + {\sqrt{n^2(z)-\mathbf{p}_s \cdot \mathbf{p}_s}} + + \frac{\pm \omega\,n^2(z) + \mathbf{k}_s \cdot \mathbf{p}_s} + {\sqrt{n^2(z)-\mathbf{p}_s \cdot \mathbf{p}_s}}\;, +\end{equation} +Equation~(\ref{eq:owedsr}) is related to the double-square-root method +\cite[]{iei}. + +In the 2-D constant-velocity case, it is convenient to express the direction +vectors in terms of angles: +\begin{equation} + \label{eq:angles} +p_r + p_s = 2\,n\,\sin{\alpha}\,\cos{\theta}\;;\quad +p_r - p_s = 2\,n\,\cos{\alpha}\,\sin{\theta}\;, +\end{equation} +where $\alpha$ is the dip angle and $\theta$ is the scattering angle. The +oriented frequency-wavenumber relationship is then +\begin{equation} + \label{eq:owedsr} + k_z(k_m,k_h,\alpha,\gamma,\omega) = + \frac{\pm 2\,\omega\,n^2 + k_m\,\sin{\alpha}\,\cos{\alpha} + + k_h\,\sin{\gamma}\,\cos{\gamma}}{\cos^2{\alpha}-\sin^2{\gamma}}\;, +\end{equation} +where $k_m = k_s + k_r$ and $k_h = k_r - k_s$. The time-space-domain analog of +expression~(\ref{eq:owedsr}) corresponds precisely to the constant-velocity +angle-migration Kirchhoff operators \cite[]{paul}. + +\section{Discussion} + +The oriented wave equation describes wave propagation in the phase space, +where propagating waves are associated not only with time and space +coordinates but also with orientation directions. The ordinary time-and-space +waves are projections from the phase space: a sum of oriented waves from all +possible directions. This concept relates directly to the concept of +angle-domain seismic imaging. It has the following attractive properties: +\begin{itemize} +\item The two directions of wave propagation (up- and down-going waves) are + explicitly separated by the choice of sign in equation~(\ref{eq:owe}). +\item Equation~(\ref{eq:owe}) is the first-order transport + (con\-ve\-cti\-on-type) equation, convenient for both mathematical and + numerical analysis. +\item In the case of lateral velocity variations, equation~(\ref{eq:owe}) + requires the slowness field $n(\mathbf{x})$ to be continuously + differentiable. This is a reasonable requirement in the seismic imaging + problem, because it assures the absence of multiple reflections in the + propagated wavefield. +\item It is possible to extend the oriented wave concept in a straightforward + way to accommodate the orientation-related effects such as anisotropy and + vector wavefield polarization. +\end{itemize} + +The practical advantage of the oriented wave approach follows from the fact +that recorded seismic wavefields allow for a sparse representation in the +phase-space domain. The sparseness is easily understood in the case of layered +or constant-velocity media. In the general case, it should be possible to +preserve sparseness in the downward wavefield continuation by decomposing the +wavefield into spatially constrained oriented beams \cite[]{wg}. The oriented +wave equation allows Gaussian beams \cite[]{GEO66-04-12401250} or other sparse +phase-space data representations to be downward continued without +computationally troublesome ray tracing. Exploring this opportunity remains an +open direction for future research. + +\bibliographystyle{seg} +\bibliography{SEP2,SEG,owe} + +%%% Local Variables: +%%% mode: latex +%%% TeX-master: t +%%% TeX-master: t +%%% TeX-master: t +%%% End: +