外文翻譯--波動方程反射走時(shí)反演

1、<p>　　2248單詞，4098漢字</p><p>　　出處：Zhang S, Schuster G, Luo Y. Wave-equation reflection traveltime inversion[C]//2011 SEG Annual Meeting. Society of Exploration Geophysicists, 2011.</p><p><

2、;b>　　畢 業(yè) 論 文</b></p><p><b>　　文獻(xiàn)翻譯</b></p><p>　　完成日期 2014 年 5 月 28 日 </p><p>　　Wave-equation Reflection Traveltime Inversion</p><p>　　Sanz

3、ong Zhang?, Gerard Schuster, King Abdullah University of Science and Technology, and Yi Luo, Saudi Aramco</p><p><b>　　SUMMARY</b></p><p>　　The main difficulty with iterative waveform

4、 inversion using a gradient optimization method is that it tends to get stuck in local minima associated within the waveform misfit function. This is because the waveform misfit function is highly nonlinear with respect

5、to changes in the velocity model. To reduce this nonlinearity, we present a reflection traveltime tomography method based on the wave equation which enjoys a more quasi-linear relationship between the model and the data.

6、 A local crosscor</p><p>　　INTRODUCTION</p><p>　　Prestack depth migration of 3D seismic data is the industry standard for computing detailed estimates of the earth’s reflectivity distribution. H

7、owever, an accurate velocity model is a precondition for imaging complex geological structures. To estimate this velocity model, there are three primary inversion methods: migration velocity analysis (MVA), traveltime in

8、version, and full waveform inversion. For migration velocity analysis (Symes and Kern, 1994; Sava and Biondi, 2004; Shen and Calandra, </p><p>　　A more detailed analysis shows that traveltime inversion is co

9、nstrained by a high-frequency approximation, and so it fails to invert for the earth’s velocity variations having nearly the same wavelength or less than that of the source wavelet. Consequently, the resolution of the ve

10、locity model constructed from the traveltimes is much less than that of full waveform inversion. The merit is that the traveltime misfit function (normed squared error between observed and calculated traveltimes) is qu&l

11、t;/p><p>　　To exploit the strengths and ameliorate the weaknesses of both ray-based traveltime tomography and full waveform inversion, wave-equation-based traveltime inversion was developed to invert the veloci

12、ty model (Luo and Schuster, 1991a and 1991b; Zhou et al, 1995; Zhang and Wang, 2009; Leeuwen and Mulder, 2010). This kind of inversion methods inverts traveltime using the gradient calculated from the wave equation. It i

13、s not constrained by a high-frequency approximation and traveltime picking is not</p><p>　　This paper is organized into three sections. The first section describes the basic theory of image-domain wave-equat

14、ion reflection traveltime inversion. The second section shows numerical examples to verify the effectiveness of this method. The last section draws some conclusions.</p><p><b>　　THEORY</b></p&

15、gt;<p>　　The key step in WRTI is to transform the reflection data into that recorded by a virtual transmission experiment. This transmission data can then be inverted by WTI (Luo and Schuster, 1991a).</p>&

16、lt;p>　　1). Assume an initial velocity model.</p><p>　　2). Migrate the recorded upgoing reflection data to get the image points at x.</p><p>　　3). Forward propagate the source at xs to x to ge

17、t the downgoing direct wave ps(x, t) as shown in Figure 1(a). Now we have the virtual source wavelet ps(x, t) where the virtual source isat x, which will be used to update the velocity on the receiver side.</p>&l

18、t;p>　　4). Backpropagate the observed reflection data from xg to the image point x and get the upgoing reflection wave pg(x, t) as shown in Figure 1(b). Now we have the virtual reflection data at x which can be used to

19、 find the the traveltime difference between the downgoing direct wave ps(x, t) and the upgoing reflection wave pg(x, t).</p><p>　　5). Crosscorrelate the downgoing direct event in ps(x, t) and the upgoing ref

20、lection event in pg(x, t) to find the time shift ?? between them as shown in Figure 1(c).</p><p>　　6). Update the velocity model by smearing ???along the weighted wavepath between xs and x and between x and

21、xg as shown in Figure 1(d). This step is actually the application of WTI to the virtual transmission data.</p><p>　　7). Repeat steps (3)-(6) for all source and image points.</p><p>　　8). Go back

22、 to step (2) until the norm of the traveltime residual satisfies the specified minimum. </p><p>　　In summary, WRTI can be decomposed into two steps. The first step is to redatum the geophones from the free s

23、urface to the image points. The second step is to redatum the source to the image point. Hence, two virtual transmission experiments are formed and used to update the velocity model. The potential benefit is that reflect

24、ion traveltime inversion might enjoy robust convergence properties and not require the tedious picking of reflection traveltimes.</p><p>　　(a) Forward extrapolate source</p><p>　　(b) Backward ex

25、trapolate geophones</p><p>　　(c) Crosscorrelate (d) Searing Δτalong wavepath</p><p>　　Figure 1: (a). The forward extrapolation of the source field. (b). The backward extrapolati

26、on of the geophone field. (c). The crosscorrelation of the downgoing direct wave and the upgoing reflection wave. (d). The misfit gradient is propertial to the ???weighted wavepath functions between the source and the im

27、age points, and the image point and the geophones.</p><p>　　Connective function</p><p>　　The following analysis assumes that the propagation of seismic waves honors the 2-D acoustic wave equatio

28、n. Let p(xr, t|xs)obs be the pressure at time t observed at the receiver location xr due to a source at xs. The source is always assumed to be initiated at zero time. For a given velocity model, p(xr, t|xs)cal denotes th

29、e calculated seismogram that honors the 2D acoustic wave equation. The crosscorrelation function between the forward wavefield and the backward wavefield can be used to determ</p><p><b>　　(1)</b>

30、</p><p>　　where ps(x, t +??) is the forward wavefield initiated by the source at xs</p><p><b>　　(2)</b></p><p>　　Here w(t) is the source wavelet, and g(x, t|xs,0) is the

31、 Green’s function. pg(x, t) is the backward wavefield by the time-reversed propagation of the observed data p(xg, t|xs)obs</p><p><b>　　(3)</b></p><p>　　and ??is the time lag of the c

32、rosscorrelation function. When ??= 0, equation (2) is the conventional correlation imaging condition. The nonzero time lag indicates the inaccuracy of the velocity model. The extremum of f (x,??) should satisfy</p>

33、<p><b>　　(4)</b></p><p><b>　　or</b></p><p><b>　　(5)</b></p><p>　　where T is the estimated maximum time lag between the forward modeled wav

34、e from the source and the backward propagated wave from the receivers. Note ???= 0 indicates that the correct velocity model has been found which generates a downgoing direct wave and upgoing reflection wave arriving at

35、the same time. The derivative of f (x,??) with respect to ? should be zero at ???unless its maximum or minimum is at an end point T or ?T:</p><p><b>　　(6)</b></p><p>　　where ˙ ps(x,

36、t +??) represents the time derivative of the calculated downgoing wave.</p><p>　　Misfit function</p><p>　　The inverse problem is defined as finding a velocity model that minimizes the following

37、misfit function:</p><p><b>　　(7)</b></p><p>　　Here x is the image point, and s is the source position. The reflection traveltime inversion is computed by finding c(x′) that minimizin

38、g the sum of the squared traveltime residuals. For simplicity, a steepest descent non-linear optimization method is used to describe the iterative minimization of equation (7), with the understanding that a preconditione

39、d conjugate gradient method is used in practice. To update the velocity model, the steepest descent method gives</p><p><b>　　(8)</b></p><p>　　where gk(x′) is the steepest descent dir

40、ection for the misfit function S, x′ represents any location in the velocity model, ?k is the step length, and k denotes the kth iteration.</p><p>　　Gradient function</p><p>　　Taking the Frech ?

41、 et derivative of S with respect to velocity perturbations yields the misfit gradient</p><p><b>　　(9)</b></p><p>　　where g (x′) represents the traveltime misfit gradient. Using (6) a

42、nd the rule for an implicit function derivative, we get</p><p><b>　　(10)</b></p><p><b>　　Where</b></p><p><b>　　(11)</b></p><p><b

43、>　　And</b></p><p><b>　　(12)</b></p><p>　　where E is the constant. Under the Born approximation, we can rewrite the misfit gradient (10) as</p><p><b>　　(13

44、)</b></p><p>　　Equation (13) indicates that the gradient function of WRTI inversion consists of two gradient functions of WTI for two virtual transmission experiments. One virtual seismic experiment is

45、 where the geophones are redatumed to the image point, and the source is on the free surface. The other one is where the source are redatumed to the image point, and the geophones are still on the free surface. The veloc

46、ity model is updated by smearing the time shifts at the image point along the wavepath between</p><p>　　NUMERICAL EXAMPLES</p><p>　　The first example is associated with a three-layer model. The

47、model in Figure 2(a) is discretized into a mesh with 201x401 gridpoints, with 100 line sources and 401 receivers on the top surface of the model, respectively. A 40-gridpoint wide absorbing sponge zone is added along eac

48、h side, and the grid interval is 20 meters. The source wavelet is a Ricker wavelet with a peak frequency of 10 Hz, and the starting model is shown in Figure 2(b) which is a constant velocity model. The observed seismog&l

49、t;/p><p>　　as shown in Figure 2(d). Figure 2(e) shows the forward modeled wavefield recorded on two reflectors. A time window is used to separate out the downgoing direct wave and the upgoing reflection wave fr

50、om the calculated data and the redatumed data indicated by the dashed lines in Figure 2(d) and Figure 2(e). The first arrival traveltime at reflectors calculated from the eikonal solver is consistent with the center of t

51、he time window. The direct downgoing waves are crosscorrelated with the correspon</p><p>　　CONCLUSION</p><p>　　A new seismic reflection traveltime tomography is presented which reconstructs velo

52、cities from reflection traveltimes computed from solutions to the wave equation. No high-frequency assumption to the data is needed, and traveltime picking and event identification are sometimes unnecessary. The mathemat

53、ical derivation demonstrates that WRTI is roughly equivalent to that of transmission tomography for two virtual transmission experiments. The synthetic data tests illustrate that it converges robus</p><p>　　

54、ACKNOWLEDGMENTS</p><p>　　We thank the 2011 sponsors of Center for Subsurface Imaging and Fluid Modeling (CSIM) at KAUST for their support.</p><p>　　Figure 2: (a). Three-layer true velocity model

55、. (b). The initial constant velocity model. (c). The observed data. The direct wave is removed. (d). The upgoing reflection wave which is obtained by redatuming the observed data from the free surface to the reflectors.

56、(e). The calculated downgoing wave on the reflectors. (f). The inversion result after seven iterations.</p><p>　　Figure 3: (a). The true velocity model with a fault. (b). The initial velocity model. (c). The

57、 inversion result after ten iterations.</p><p>　　REFERENCES</p><p>　　Bishop, T., K. Bube, R. Cutler, R. Langan, P. Love, J. Resnick, R. Shuey, D. Spindler, and H. Wyld, 1985, Tomographic determi

58、nation of velocity and depth in laterally varying media: Geophysics, 50, 903–923, doi:10.1190/1.1441970.</p><p>　　Crase, E., C. Wideman, M. Noble, and A. Tarantola, 1992, Nonlinear elastic waveform inversion

59、 of land seismic reflection data: Journal of Geophysics Research, 97, 4685–4704.</p><p>　　Dines, K., and R. Lytle, 1979, Computerized geophysical tomography: Proceedings of the IEEE, 67, 1065–1073, doi:10.11

60、09/PROC.1979.11390.</p><p>　　Ivansson, S., 1985, A study of methods for tomographic velocity estimation in the presence of lowvelocity zones: Geophysics, 50, 969–988, doi:10.1190/1.1441975.</p><p&

61、gt;　　Levander, A. R., 1988, Fourth-order finite-difference P-SV seismograms: Geophysics, 53, 1425–1436, doi:10.1190/1.1442422.</p><p>　　Lines, L., 1988, Inversion of geophysical data: SEG Geophysics Reprint

62、Series No. 9.</p><p>　　Luo, Y., and G. Schuster, 1991a, Wave equation traveltime inversion: Geophysics, 56, 645–653, doi:10.1190/1.1443081.</p><p>　　Luo, Y., and G. Schuster, 1991b, Wave equatio

63、n inversion of skeletonized geophysical data: Geophysical Journal International, 105, 289–294, doi:10.1111/j.1365-246X.1991.tb06713.x.</p><p>　　Mora, P., 1987, Nonlinear two-dimensional elastic inversion of

64、multioffset seismic data: Geophysics, 52, 1211–1228, doi:10.1190/1.1442384.</p><p>　　Paulsson, B., N. Cook, and T. McEvilly, 1985, Elastic wave velocities and attenuation in an underground granitic repositor

65、y for nuclear waste: Geophysics, 50, 551–570, doi:10.1190/1.1441932.</p><p>　　Pratt, R. G., C. Shin, and G. J. Hicks, 1998, Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion:

66、 Geophysical Journal International, 133, 341–362, doi:10.1046/j.1365- 246X.1998.00498.x.</p><p>　　Sava, P., and B. Biondi, 2004, Wave-equation migration velocity analysis, I: Theory: Geophysical Prospecting,

67、 52, 593–606.</p><p>　　Shen, P., and H. Calandra, 2005, One-way waveform inversion within the framework of adjoint state differential migration: 75th Annual International Meeting, SEG, Expanded Abstracts, 17

68、09–1712.</p><p>　　Symes, W., and M. Kern, 1994, Inversion of reflection seismograms by differential semblance analysis: Algorithm structure and synthetic examples: Geophysical Prospecting, 42, 565–614, doi:1

69、0.1111/j.1365-2478.1994.tb00231.x.</p><p>　　Tarantola, A., 1986, A strategy for nonlinear elastic inversion of seismic reflection data: Geophysics, 51, 1893–1903, doi:10.1190/1.1442046.</p><p>　

70、　Tarantola, A., 1987, Inverse problem theory: Elsevier.</p><p>　　van Leeuwen, T., and W. A. Mulder, 2010, A correlation-based misfit criterion for wave-equation traveltime tomography: Geophysical Journal Int

71、ernational, 182, 1383–1394, doi:10.1111/j.1365- 246X.2010.04681.x.</p><p>　　Zhou, C., W. Cai, Y. Luo, G. Schuster, and S. Hassanzadeh, 1995, Acoustic wave-equation traveltime and waveform inversion of crossh

72、ole seismic data: Geophysics, 60, 765–773, doi:10.1190/1.1443815.</p><p>　　波動方程反射走時(shí)反演</p><p><b>　　摘要</b></p><p>　　迭代梯形反演運(yùn)用一個(gè)梯度優(yōu)化的方法，其主要的困難是該方法往往限制在局部極小值，與此對應(yīng)的問題就是函數(shù)與波形不匹

73、配（即波形失配函數(shù)）。這是因?yàn)橄鄬τ谒俣饶Ｐ偷母淖?，波形失配函?shù)是高度非線性的。為了降低這種非線性，我們提出了一種基于波動方程的反射走時(shí)層析成像的方法，該波動方程在模型和數(shù)據(jù)之間擁有更多的準(zhǔn)線性關(guān)系。在圖像點(diǎn)上面的一個(gè)帶窗口的下行直達(dá)波和上行反射波的局部互相關(guān)會產(chǎn)生滯后的時(shí)間，這個(gè)滯后的時(shí)間可以最大限度的提高相關(guān)性。這種滯后的時(shí)間呈現(xiàn)的反射走時(shí)殘差是陸地模型的背投影，以此作為波動方程旅行時(shí)反演的相同方式來更新速度模型。沒有走時(shí)的采集的需

74、要，沒有高頻的近似的假設(shè)。數(shù)學(xué)的推導(dǎo)和數(shù)值的例子用以部分展現(xiàn)該方法的有效性和魯棒性。</p><p><b>　　引言</b></p><p>　　三維地震資料疊前深度偏移是計(jì)算地球反射率分布的詳細(xì)評估的行業(yè)標(biāo)準(zhǔn)。然而一個(gè)準(zhǔn)確的速度模型是復(fù)雜地質(zhì)構(gòu)造成像的前提。為了顧及這個(gè)速度模型，這里主要有三種反演的方法：偏移速度分析（MVA），走時(shí)反演和全波形反演。對于偏移速度分

75、析（賽姆斯和克恩，1994;薩瓦河和比昂迪，2004; 申和克拉的若，2005），最佳的偏移速度是一種可以在普通圖像上讓反射事件的聚集可以最好的變得平坦。對于走時(shí)反演（戴恩斯和萊特爾，1979;波爾森等，1985; 伊瓦森，1985;比舍普等人，1985;萊恩斯，1988），折射和反射走時(shí)的引進(jìn)被用于反演速度模型的平滑特性，然而全波形反演（塔蘭特拉，1986，1987;莫拉，1987; 克萊斯等人，1992; 周等人，1995;普拉特，

76、1998）反演與陸地模型的細(xì)節(jié)有關(guān)的波形信息。</p><p>　　更詳細(xì)的分析顯示，走時(shí)反演受一個(gè)高頻近似的制約，所以它不能反演陸地速度的變化，其擁有和震源子波幾乎相同或者較小的波長。因此，利用走時(shí)構(gòu)造的速度模型的分辨率比全波形反演的分辨率低得多。走時(shí)的優(yōu)點(diǎn)是失配函數(shù)（觀測和計(jì)算走時(shí)之間賦范誤差）對于速度的擾動是準(zhǔn)線性的，所以，盡管起始模型和真實(shí)模型相差很大也可以實(shí)現(xiàn)一個(gè)有效的速度反演（羅和舒斯特，1991a和

77、1991b; 周等，1995）。盡管對于起始模型和富含噪音的振幅的選擇是非常的敏感的，全波形反演有時(shí)可以重建陸地模型的一個(gè)非常詳細(xì)的估計(jì)。這是因?yàn)闆]有對數(shù)據(jù)的高頻假定，幾乎所有地震事件是嵌在失配函數(shù)。然而，相對于全波形反演的問題就是其失配函數(shù)（觀察到的和合成地震記錄的賦范誤差）可相對于變化速度模型高度非線性。在這種情況下，如果起始模型和真實(shí)模型之間差距較大，一個(gè)梯度優(yōu)化的方法往往限制在局部極小值。</p><p>

78、;　　利用射線走時(shí)成像和全波形反演的優(yōu)點(diǎn)，并改進(jìn)其缺點(diǎn)，從而基于波動方程走時(shí)的反演被用于反演速度模型（羅和舒斯特，1991a，1991b；周等人，1995；張和王2009；李文和穆德2010）。這種反演方法是利用波動方程計(jì)算出的梯度來反演走時(shí)的。該方法沒有受到高頻近似的限制，也沒有進(jìn)行走時(shí)采集的需要。其他還有一些重要的優(yōu)勢，收斂的速度在某種程度上對于原始模型不敏感，而且在噪音存在的情況下，數(shù)據(jù)具有魯棒性。然而，這些反演的方法的設(shè)計(jì)是用于

79、反演傳輸波在地震數(shù)據(jù)中形式，而不是用于反演反射走時(shí)的設(shè)計(jì)。不同于折射波和直達(dá)波，對于反演模型，反射波可以提供更多有關(guān)于更深地下界面的速度信息。然而，如果初始模型和真實(shí)模型相差很大，反射波的全波形反演是困難的。為了克服這些局限性，本文在波動方程傳輸走時(shí)反演（WTI）（羅和舒斯特，1991a和1991b）的基礎(chǔ)之上進(jìn)行了延伸，提出了波動方程反射走時(shí)反演（WRTI）的方法。</p><p>　　本文分為三個(gè)部分。第一部

80、分介紹了在圖像域波動方程反射走時(shí)反演的基本理論。第二部分展現(xiàn)了一些數(shù)值例子來驗(yàn)證這種反演方法的有效性。最后一部分得出了一些結(jié)論。</p><p><b>　　原理</b></p><p>　　WRTI的關(guān)鍵步驟是將反射的地震數(shù)據(jù)轉(zhuǎn)換成為通過一個(gè)虛擬傳輸實(shí)驗(yàn)得出的記錄。這種傳輸數(shù)據(jù)在之后可以運(yùn)用WRTI（羅和舒斯特，1991a）方法來反演。</p><

81、;p>　　假設(shè)有一個(gè)初始的速度模型；</p><p>　　移動記錄的上行反射數(shù)據(jù)來得到圖像上點(diǎn)x；</p><p>　　在xs點(diǎn)的震源向前傳播到達(dá)x點(diǎn)，以獲得下行直達(dá)波ps（x，t），如圖1（a）所示?，F(xiàn)在我們有了虛擬的震源子波ps（x，t），其虛擬的震源點(diǎn)在x點(diǎn)，這將被用于更新在接收器一側(cè)的速度；</p><p>　　反向傳播觀察到的反射數(shù)據(jù)從xg點(diǎn)到達(dá)圖像

82、上的點(diǎn)x，以獲得上行反射波pg（x，t），如圖1（b）所示?，F(xiàn)在我們有了x點(diǎn)的虛擬的反射數(shù)據(jù)，這些數(shù)據(jù)可以被用于找出下行直達(dá)波ps（x，t）和上行反射波pg（x，t）的走時(shí)誤差；</p><p>　　對下行直達(dá)波ps（x，t）和上行反射波pg（x，t）做互相關(guān)，可以找出兩者之間的時(shí)移Δτ，如圖1（c）所示；</p><p>　　沿著xs和x之間和xg和x之間的加權(quán)的傳播路徑，通過去除時(shí)差Δ

83、τ來更新速度模型，如圖1（d）所示。這一步是WTI對于虛擬傳輸數(shù)據(jù)的真實(shí)應(yīng)用；</p><p>　　對于所有的震源和圖像點(diǎn)重復(fù)步驟（3）-（6）；</p><p>　　回到步驟（2），知道走時(shí)殘差滿足要求的最低標(biāo)準(zhǔn)。</p><p>　　總的來說，WRTI可以分為兩個(gè)步驟。第一步就是從自由表面的檢波器到圖像點(diǎn)重建基準(zhǔn)面。第二步是從震源到圖像點(diǎn)重建基準(zhǔn)面。所以，形成了

84、兩個(gè)虛擬傳輸實(shí)驗(yàn)，并用此來更新速度模型。潛在的好處就是反射走時(shí)反演可能具有強(qiáng)大的收斂性，而且不需要繁瑣的反射走時(shí)的的采集。</p><p><b>　　(a)震源向前傳播</b></p><p>　　(b)檢波器向后傳播</p><p>　　(c)互相關(guān)(d)去除時(shí)差Δτ</p><p>　　圖1（a）震

85、源場向前傳播。（b）檢波器場向后傳播。（c）下行直達(dá)波和上行反射波的互相關(guān)。（d）震源和圖像點(diǎn)以及圖像點(diǎn)和檢波器之間的加權(quán)傳播路徑函數(shù)Δτ正比于失配梯度。</p><p><b>　　相關(guān)函數(shù)</b></p><p>　　下面的分析假設(shè)地震波的傳播非常好的滿足二維聲波方程。設(shè) p(xr, t|xs)obs在時(shí)間t的壓力，t為從xs點(diǎn)的震源激發(fā)到xr點(diǎn)的檢波器接收的時(shí)間

86、。震源通常被假設(shè)為在零時(shí)刻開始激發(fā)。對于給定的速度模型， p(xr, t|xs)cal表示計(jì)算出的地震記錄，其和二維聲波方程匹配的很好。向前傳播的波場和向后傳播的波場的互相關(guān)函數(shù)可以用來確定在x點(diǎn)的圖像。</p><p><b>　?。?）</b></p><p>　　ps(x, t +t ) 是由在xs點(diǎn)的震源激發(fā)的向前傳播的波場</p><p&

87、gt;<b>　?。?）</b></p><p>　　在這里 w(t)為震源的子波，g(x, t|xs,0) 為格林函數(shù)。pg(x, t)為反向傳播的波場，該波場是由觀測到的數(shù)據(jù)p(xg, t|xs)反向傳播形成的時(shí)間場</p><p><b>　?。?）</b></p><p>　　τ是互相關(guān)函數(shù)的滯后時(shí)間。當(dāng)τ=0的時(shí)

88、候，方程（2）是現(xiàn)有的互相關(guān)成像條件。非零的滯后時(shí)間表示速度模型的誤差。 f (x,t )的極值應(yīng)該滿足</p><p><b>　?。?）</b></p><p><b>　　或者</b></p><p><b>　　（5）</b></p><p>　　其中T是從震源向前傳播

89、的波和從檢波器向后傳播的波的最大滯后時(shí)間的估計(jì)。注意Δτ=0說明此事的速度模型就是你所尋找的真正的速度模型，在該速度模型情況下，向下傳播的直達(dá)波和向上傳播的反射波在相同的時(shí)間到達(dá)。在Δτ=0時(shí)，f (x,t )的導(dǎo)數(shù)等于零，除非它的最大值或者最小值在點(diǎn)T或者-T結(jié)束：</p><p><b>　?。?）</b></p><p>　　其中ps(x, t +t )表示計(jì)算

90、的下行波的導(dǎo)數(shù)時(shí)間。</p><p><b>　　失配函數(shù)</b></p><p>　　反演的問題就是被定義為發(fā)現(xiàn)了一個(gè)速度模型，該速度模型使得下列的失配函數(shù)最小化：</p><p><b>　?。?）</b></p><p>　　在這里，x是圖像點(diǎn)所在的位置，s是震源所在的位置。反射走時(shí)反演是通過

91、尋找 c(x′)來計(jì)算的， c(x′)是最小平方走時(shí)殘差之和的最小值。為了簡單起見，一個(gè)最速下降的非線性優(yōu)化方法被用于描述使得方程（7）最小化的迭代算法，據(jù)了解，一個(gè)預(yù)條件共軛梯度算法在實(shí)踐中通常被應(yīng)用。</p><p>　　為了更新速度模型，最速下降算法給出了</p><p><b>　　（8）</b></p><p>　　其中γk(x′)是

92、失配函數(shù)S的最速下降的方向， x′ 表示速度模型中的任何位置，αk是步長，k表示第k次的迭代。</p><p><b>　　梯度函數(shù)</b></p><p>　　針對于速度的擾動，利用S作為遞推方程來產(chǎn)生失配梯度</p><p><b>　　（9）</b></p><p>　　其中g(shù) (x′) 表示

93、走時(shí)失配梯度，利用方程（6）和隱函數(shù)的導(dǎo)數(shù)規(guī)則，我們可以得到</p><p><b>　?。?0）</b></p><p>　　其中 （11）</p><p><b>　?。?2）</b></p><p>　　其中E是恒定的。根據(jù)Born近似，我

94、們可以將失配梯度方程（10）重寫為</p><p><b>　?。?3）</b></p><p>　　方程（13）表示W(wǎng)RTI反演的梯度函數(shù)，該梯度函數(shù)由兩個(gè)對于虛擬傳輸實(shí)驗(yàn)的WTI形成的兩個(gè)梯度函數(shù)組成，一個(gè)虛擬的地震實(shí)驗(yàn)就是由檢波器的接收演變到圖像點(diǎn)，而且震源在自由表面。另一個(gè)虛擬的地震實(shí)驗(yàn)就是由震源的激發(fā)演變到圖像點(diǎn)，而且檢波器仍然在自由表面。速度模型運(yùn)用滯后時(shí)

95、間在圖像點(diǎn)上沿著震源和圖像點(diǎn)的傳播路徑以及圖像點(diǎn)和檢波點(diǎn)的傳播路徑的切除來更新速度模型。</p><p><b>　　數(shù)值例子</b></p><p>　　第一個(gè)模型是一個(gè)三層的模型。圖2（a）中的模型是一個(gè)201×401的網(wǎng)格模型，該模型在模型的最上層的表面上有100條震源線和401個(gè)接收點(diǎn)。在每一邊增加了一個(gè)40個(gè)網(wǎng)格寬度的吸收邊界，而且網(wǎng)格間距為20m

96、。震源子波是一個(gè)主頻為10Hz的雷克子波，圖2（b）中顯示的是起始模型，該模型是一個(gè)速度恒定的模型。所觀察到的地震記錄是由對一個(gè)二維聲波方程（該波動方程具有恒定的密度）進(jìn)行四階有限差分的求解產(chǎn)生的。圖2（c）是一個(gè)典型的自由表面的炮集記錄，其中直達(dá)波是被去除的。所觀察到的數(shù)據(jù)是由自由表面演變到發(fā)射面的，如圖2（d）所示。圖2（e）展示的是在兩個(gè)反射界面上的向前傳播的模型的波場記錄。一個(gè)時(shí)間窗口被用來從計(jì)算的數(shù)據(jù)中分離出下行直達(dá)波和上行反

97、射波，演變的數(shù)據(jù)如圖2（d）和圖2（e）的虛線所示。從程涵方程計(jì)算得來的反射界面的首波和時(shí)間窗口的中心一致。下行直達(dá)波和上行發(fā)射波的互相關(guān)可以找到兩者的滯后時(shí)差。圖2（f）顯示的是經(jīng)過7次反演迭代的結(jié)果。它顯示出波動方程反射走時(shí)反演對于速度模型的構(gòu)建是一種有效的方式。接下來，我們在一個(gè)更實(shí)際的存在問題的模型上來測試我們的反演算法。圖3（a）顯示的是一個(gè)存</p><p>　　圖2：（a）三層真實(shí)速度模型。（b）最

98、初的恒定速度模型（c）觀測到的數(shù)據(jù)，直達(dá)波被去除。（d）通過從自由表面到反射面觀測到的數(shù)據(jù)的演化得到的上行反射波。（e）計(jì)算的反射界面的下行波。（f）七次迭代之后的反演結(jié)果</p><p>　　圖3：（a）有問題的真正的速度模型。（b）起始速度模型。（c）十次迭代之后的反演。</p><p><b>　　結(jié)論</b></p><p>　　一個(gè)新