基于快速步進(jìn)法的三維反射波追蹤

姜莉

摘 要： 引入了一種改進(jìn)的快速步進(jìn)算法，對(duì)三維復(fù)雜介質(zhì)中反射波進(jìn)行波前追蹤，并利用QT和OpenGL相結(jié)合的方法對(duì)速度場(chǎng)以及射線路徑進(jìn)行可視化模擬。通過模型仿真，驗(yàn)證了該方法的效率、精度以及可靠性。

關(guān)鍵詞： 快速步進(jìn)法； 反射波； 射線追蹤； 波前追蹤

中圖分類號(hào)： TN911.7?34； P314.3 文獻(xiàn)標(biāo)識(shí)碼： A 文章編號(hào)： 1004?373X（2014）11?0132?03

Abstract： In the paper， an improved fast stepping algorithm is introduced to make wavefront tracing of the reflected wave in 3D heterogeneous media， and perform the visualization simulation of velocity field and ray path by using QT and OpenGL. The efficiency， accuracy and reliability of the newalgorithm were verified by model simulation.

Keywords： fast stepping method； reflected wave； ray tracing； wavefront tracing

費(fèi)馬原理指出，地震波沿著所需時(shí)間最短的路徑傳播，在高頻近似條件下，地震波場(chǎng)的主能量沿射線軌跡傳播。地震射線追蹤的目的就是尋找震源到各個(gè)接收點(diǎn)的地震波傳播的最短傳播路徑，從而研究地震波的傳播規(guī)律，對(duì)地層探測(cè)及層析成像具有重要的意義。針對(duì)傳統(tǒng)射線追蹤方法僅能求出地震波的初至旅行時(shí)的缺點(diǎn)，本文在快速步進(jìn)算法的基礎(chǔ)上給出一種多級(jí)的方法，對(duì)反射波進(jìn)行射線追蹤。該方法將波前到達(dá)的每一層分割為一個(gè)獨(dú)立計(jì)算的區(qū)域，把波前到達(dá)的每層界面重新作為波源重新計(jì)算，從而逐層推進(jìn)式地追蹤反射波或透射波的波前。

1 算法實(shí)現(xiàn)

快速步進(jìn)算法的主要思想就是通過系統(tǒng)地用上風(fēng)格式構(gòu)建旅行時(shí)場(chǎng)來解程函方程。程函方程的上風(fēng)差分格式實(shí)際上意味著信息傳播是從走時(shí)的小值到大值。因此，快速步進(jìn)算法重點(diǎn)在通過從[T（x，y，z）]的小值向外擴(kuò)展來構(gòu)建方程的解[1]。

1.1 旅行時(shí)場(chǎng)方程

在3D旅行時(shí)計(jì)算時(shí)，程函方程為：

1.2 模型建立

把目標(biāo)區(qū)域離散為網(wǎng)格，將沒被界面穿過的網(wǎng)格稱作規(guī)則網(wǎng)格，每個(gè)網(wǎng)格頂點(diǎn)稱作規(guī)則結(jié)點(diǎn)，如圖1中大實(shí)心點(diǎn)所示。在規(guī)則網(wǎng)格部分，每一個(gè)網(wǎng)格都是具有八個(gè)網(wǎng)格頂點(diǎn)的區(qū)域，同時(shí)每個(gè)規(guī)則結(jié)點(diǎn)屬于八個(gè)網(wǎng)格。如果規(guī)則網(wǎng)格被界面切割，稱作不規(guī)則網(wǎng)格，網(wǎng)格邊與界面相交，產(chǎn)生的交點(diǎn)稱作不規(guī)則網(wǎng)格結(jié)點(diǎn)，如圖1中小實(shí)心點(diǎn)所示，每個(gè)交點(diǎn)屬于四個(gè)網(wǎng)格。如果要更新某一點(diǎn)的走時(shí)，認(rèn)為每個(gè)網(wǎng)格相互獨(dú)立。

1.3 算法實(shí)現(xiàn)

如圖2所示，為了加快計(jì)算，把目標(biāo)區(qū)域的所有網(wǎng)格結(jié)點(diǎn)劃分為三個(gè)子集，首先把具有初始旅行時(shí)值的點(diǎn)記為確認(rèn)值點(diǎn)（實(shí)心黑點(diǎn)），這些點(diǎn)的旅行時(shí)在計(jì)算中不再改變；然后把與確認(rèn)值點(diǎn)相鄰一個(gè)網(wǎng)格的點(diǎn)記為近點(diǎn)（實(shí)心綠點(diǎn)），近點(diǎn)需要計(jì)算一個(gè)旅行時(shí)的值，但它不一定是最小值，是可改變的；其他點(diǎn)記為遠(yuǎn)點(diǎn)（空心點(diǎn)），遠(yuǎn)點(diǎn)的旅行時(shí)值還沒有計(jì)算，需要逐次推進(jìn)計(jì)算。

計(jì)算最小旅行時(shí)值的過程如下：首先在近點(diǎn)中搜索旅行時(shí)最小的點(diǎn)，把這個(gè)點(diǎn)從近點(diǎn)子集中去除，加入確認(rèn)值點(diǎn)子集，然后把此點(diǎn)的所有非確認(rèn)值的臨近點(diǎn)歸入近點(diǎn)子集，如果臨近點(diǎn)屬于遠(yuǎn)點(diǎn)子集，則從中去除；接著重新計(jì)算該點(diǎn)臨近點(diǎn)的旅行時(shí)值，如果比原值大，則保留原值，如果小，則取而代之。重復(fù)前面的步驟，直到確認(rèn)值點(diǎn)集包含整個(gè)目標(biāo)區(qū)域。本算法利用堆排序技術(shù)，可有效地在近點(diǎn)集合中搜索最小元素，使算法加快速度。

1.4 反射波前追蹤

對(duì)反射波射線追蹤時(shí)如果直接利用快速步進(jìn)法計(jì)算時(shí)，會(huì)存在以下問題：使用快速步進(jìn)方法的有限差分解必須服從絕對(duì)初至?xí)r間。對(duì)于反射波，波前至少穿過計(jì)算區(qū)域兩次，因此只能通過多級(jí)方法來解決。多級(jí)方法的原理如圖3所示。

第一步，波前傳播到反射界面，然后將該區(qū)域進(jìn)行復(fù)制。

第二步，將快速步進(jìn)方法用于復(fù)制的區(qū)域。第二步開始時(shí)，窄帶圍繞著界面，這樣做并不破壞因果性，因?yàn)槌踔练瓷洳ㄇ爸坏竭_(dá)每個(gè)界面結(jié)點(diǎn)一次。

1.5 走時(shí)更新

1.5.1 規(guī)則網(wǎng)格

對(duì)于規(guī)則網(wǎng)格點(diǎn)的走時(shí)更新，用上風(fēng)格式近似表示式（1）的梯度函數(shù)，可得到如下的表達(dá)式：

1.5.2 不規(guī)則網(wǎng)格

如果被更新結(jié)點(diǎn)處于不規(guī)則網(wǎng)格中，使用相鄰的不規(guī)則點(diǎn)的走時(shí)來更新某一結(jié)點(diǎn)時(shí)，假設(shè)波前近似為平面。當(dāng)使用不規(guī)則網(wǎng)格中的確認(rèn)值點(diǎn)更新時(shí)，此網(wǎng)格單元中的確認(rèn)值點(diǎn)數(shù)可能是一個(gè)，兩個(gè)，三個(gè)或大于三個(gè)。如果由三個(gè)確認(rèn)值點(diǎn)的信息來更新此點(diǎn)走時(shí)，這三點(diǎn)可以惟一確定一個(gè)波陣面，準(zhǔn)確計(jì)算該點(diǎn)走時(shí)。如果僅有一個(gè)或兩個(gè)點(diǎn)是確認(rèn)值點(diǎn)，必須增加額外的約束信息，找到波陣面的解。當(dāng)確認(rèn)值點(diǎn)多于三個(gè)時(shí)，將產(chǎn)生冗余信息。下面分別考慮這四種情況。

（1） 由三個(gè)確認(rèn)值點(diǎn)更新

[x，][y，][z]分別為[x，][y]和[z]三個(gè)方向的單位矢量。由以上的計(jì)算公式可以看出，每個(gè)網(wǎng)格的慢度矢量在計(jì)算旅行時(shí)的過程中，就可求出，不用特別計(jì)算，從而可節(jié)省計(jì)算量。根據(jù)射線理論，射線與慢度矢量的方向一致，射線在模型中某一點(diǎn)的傳播方向，用該點(diǎn)周圍網(wǎng)格點(diǎn)的慢度矢量插值求出，在本文中采用雙線性插值。

2 實(shí)驗(yàn)及仿真

本文采用層狀介質(zhì)模型進(jìn)行實(shí)驗(yàn)仿真，模型為一個(gè)兩層的三維介質(zhì)模型。如圖4所示，上層用單一顏色表示其傳播速度為常數(shù)，下層的速度隨深度變化，顏色越深表示傳播速度越快。目標(biāo)區(qū)域長(zhǎng)寬深分別為2 000 m，2 000 m和1 000 m，在[x，][y]和[z]方向分別劃分為81，81和41個(gè)網(wǎng)格，網(wǎng)格間距在各個(gè)方向上都是25 m，炮點(diǎn)位于地面上（0，1 000，1 000）坐標(biāo)處。圖4是模型速度場(chǎng)的剖面圖，拖動(dòng)右邊滑塊可以控制各方向剖面的移動(dòng)。圖5是射線及界面模擬圖，拖動(dòng)滑塊可以對(duì)圖像進(jìn)行旋轉(zhuǎn)。從圖中可以看出，模擬結(jié)果與理論結(jié)果一致，快速步進(jìn)法確實(shí)是一種快速穩(wěn)定，適用于實(shí)際的方法。

3 結(jié) 論

本文采用Fortran和QT類庫(kù)以及OpenGL相結(jié)合，的方法很好地實(shí)現(xiàn)了快速步進(jìn)法射線追蹤，不僅利用前人已有的成果進(jìn)一步降低開發(fā)難度，提高開發(fā)效率。同時(shí)易于向Unix移植，以適應(yīng)一些用戶的特殊需要。由模型實(shí)驗(yàn)看，所得結(jié)果比較客觀地反映出了實(shí)際情況。所以，此法是一種很好的三維正演方法。

參考文獻(xiàn)

[1] KIM Seongjai. 3?Deikonalsolvers： first?arrival traveltimes [J]. Geophysical， 2002， 67（4）： 1225?1231.

[2] 許琨，吳律，王妙月.改進(jìn)Moser法射線追蹤[J].地球物理學(xué)進(jìn)展，1998（4）：60?64.

[3] 張?jiān)奇櫇h明，師學(xué)明.基于MOSER曲射線追蹤的SIRT聲波層析成像[J].工程地球物理學(xué)報(bào)，2005（3）：167?169.

[4] 張霖斌，劉迎曦，趙振峰，等.有限差分法射線追蹤[J].石油地球物理勘探，1993，42（6）：198?200.

[5] 張美根，賈豫葛，王妙月，等.界面二次源波前擴(kuò)展法全局最小走時(shí)射線追蹤技術(shù)[J].地球物理學(xué)報(bào)，2006，49（4）：1169?1175.

[6] 李艷.地震勘探中二維波動(dòng)方程的小波數(shù)值正演模擬[D].哈爾濱：哈爾濱工業(yè)大學(xué)，2002.

[7] 周竹生，劉喜亮，熊孝雨.彈性介質(zhì)中瑞雷面波有限差分法正演模擬[J].地球物理學(xué)報(bào)，2007，50（2）：567?573.

[8] 劉恩儒，岳建華，劉彥.具有離散裂縫空間分布的二維固體中地震波傳播的有限差分模擬[J].地球物理學(xué)報(bào)，2006，49（1）：180?188.

[9] 朱亞平，楊慧珠.OpenGL技術(shù)在地震數(shù)據(jù)可視化中的應(yīng)用[J].石油地球物理勘探，2000，35（4）：403?414.

[10] WOO Mason， NEIDER Jaekie， DAVIS Tom， et al. OpenGL編程權(quán)威指南[M].北京：中國(guó)電力出版社，2001.

[11] 王博樺，陳振一，王關(guān)橋，等.地震信號(hào)檢測(cè)新機(jī)理研究[J].現(xiàn)代電子技術(shù)，2013，36（21）：111?112.