專(zhuān)利名稱(chēng)：二維復(fù)雜結(jié)構(gòu)三角網(wǎng)射線(xiàn)追蹤全局方法
技術(shù)領(lǐng)域：
本發(fā)明涉及一種二維復(fù)雜結(jié)構(gòu)三角網(wǎng)射線(xiàn)追蹤全局方法，具體是一種二維 復(fù)雜結(jié)構(gòu)的走時(shí)層析成像反演和以Maslov射線(xiàn)理論為基礎(chǔ)的波場(chǎng)計(jì)算的射線(xiàn) 追蹤方法，屬于勘探地球物理技術(shù)領(lǐng)域。
背景技術(shù)：
射線(xiàn)理論和射線(xiàn)方法是認(rèn)識(shí)波場(chǎng)傳播規(guī)律的重要途徑，是研究地下復(fù)雜構(gòu) 造和不均勻介質(zhì)中的波場(chǎng)傳播問(wèn)題的重要手段，也是走時(shí)層析成像反演的基 礎(chǔ)。
無(wú)論是正演模擬還是層析反演，模型參數(shù)化是基本的問(wèn)題。由于工程地質(zhì) 探測(cè)現(xiàn)場(chǎng)條件的復(fù)雜性，所探測(cè)研究對(duì)象的外部幾何形狀大部分情況下是不規(guī) 則的形狀，如"工"字型結(jié)構(gòu)；內(nèi)部地質(zhì)異常體在空間結(jié)構(gòu)上更加復(fù)雜，有的呈 不規(guī)則柱狀體(如溶洞、陷落柱等)，有的呈條帶狀(如斷層、裂縫等)、有 的呈層狀(如軟弱夾層等)。這就模型參數(shù)化提出更高的要求，以滿(mǎn)足復(fù)雜結(jié) 構(gòu)體的高分辨率探測(cè)。目前模型參數(shù)化一般都采用矩形網(wǎng)或三角網(wǎng)，其中矩形 網(wǎng)剖分是最普遍的形式，而三角網(wǎng)剖分也主要是在矩形網(wǎng)基礎(chǔ)上簡(jiǎn)單的三角 化。這些簡(jiǎn)單的矩形網(wǎng)、三角網(wǎng)剖分沒(méi)有充分考慮區(qū)域復(fù)雜邊界，也未將射線(xiàn) 密度、介質(zhì)特性同網(wǎng)格大小有機(jī)結(jié)合起來(lái)，難以滿(mǎn)足復(fù)雜結(jié)構(gòu)的射線(xiàn)追蹤正演 和高分辨率層析成像反演的要求。必須采用不規(guī)則網(wǎng)，如三角形網(wǎng)、不規(guī)則四 邊形網(wǎng)、或三角形網(wǎng)和不規(guī)則四邊形網(wǎng)的混合網(wǎng)，對(duì)模型參數(shù)化。研究表明， 相對(duì)于矩形網(wǎng)絡(luò)參數(shù)化，三角網(wǎng)格參數(shù)化具有如下優(yōu)點(diǎn)①模型剖分的靈活性
強(qiáng)；②對(duì)速度間斷面的描述精度高，速度模型和界面模型具有一致性；③n^演
模擬存儲(chǔ)量小、計(jì)算時(shí)間少；④正演模擬難度降低、誤差減少、效率提高； 用于層析反演時(shí)存儲(chǔ)量小、計(jì)算時(shí)間少、方程性態(tài)好、求解容易。
射線(xiàn)追蹤的理論基礎(chǔ)是，在高頻近似條件下，彈性波場(chǎng)的主能量沿射線(xiàn)軌 跡傳播。傳統(tǒng)的射線(xiàn)追蹤方法，通常意義上包括初值問(wèn)題的試射法(Shootingmethod)和邊值問(wèn)題的彎曲法(Bending method)。 Vidale(l988)在提出程函方程的
有限差分法時(shí)，指出試射法和彎曲法的主要問(wèn)題在于①難于處理介質(zhì)中較強(qiáng) 的速度變化；②難于求出多值走時(shí)中的全局最小走時(shí)；③計(jì)算效率較低； 陰 影區(qū)內(nèi)射線(xiàn)覆蓋密度不足。
對(duì)于三角網(wǎng)模型的射線(xiàn)追蹤問(wèn)題，有關(guān)文獻(xiàn)，作了探討，但其給出的射線(xiàn) 追蹤方法基本歸類(lèi)于通常意義上的初值問(wèn)題的試射法和邊值問(wèn)題的試射法和 彎曲法，不可避免地存在其固有缺點(diǎn)。
對(duì)于矩形網(wǎng)模型的射線(xiàn)追蹤問(wèn)題，有關(guān)文獻(xiàn)對(duì)常見(jiàn)的射線(xiàn)追蹤方法做了比 較詳細(xì)的總結(jié)。為了克服傳統(tǒng)試射法和彎曲法的這缺點(diǎn)，近年來(lái)，許多地球物 理研究者在這方面進(jìn)行了大量的工作，提出了一些精度較高、效率較高而且實(shí) 用的計(jì)算初至旅行時(shí)的波前追蹤方法。研究進(jìn)展主要體現(xiàn)在①在傳統(tǒng)的試射 法及彎曲法的基礎(chǔ)上的改進(jìn)，如各類(lèi)波前重建方法
(Vinje，1992;Sun,1992;Lambare et al,,1996)，除多值走時(shí)外，還較好地解決了計(jì) 算效率及陰影區(qū)覆蓋不足的問(wèn)題；②對(duì)最小走時(shí)算法的改進(jìn)，使之可適應(yīng)多值 走時(shí)計(jì)算，如慢度匹配法(Symes,1998)，可認(rèn)為是最短路徑方法的推廣；③傳 統(tǒng)方法與最小走時(shí)算法的結(jié)合，如HWT方法(Sava， Fomel， 1998)，則是通過(guò) 波前傳播計(jì)算射線(xiàn)路徑。這些方法中，同時(shí)考慮空間所有離散點(diǎn)上的走時(shí)和射 線(xiàn)路徑的全局計(jì)算方法，由于其較高的精度和計(jì)算效率引起了廣泛的關(guān)注。然 而，這些主要針對(duì)規(guī)則網(wǎng)射線(xiàn)追蹤方法，應(yīng)用到三角網(wǎng)中會(huì)使問(wèn)題復(fù)雜化，甚 至難以實(shí)現(xiàn)。
三角網(wǎng)的優(yōu)勢(shì)使得模型三角網(wǎng)參數(shù)化是未來(lái)地球物理模型參數(shù)化的重要趨 勢(shì)，但其波場(chǎng)射線(xiàn)追蹤問(wèn)題卻缺乏深入的研究，因此，三角網(wǎng)模型及其波場(chǎng)射 線(xiàn)追蹤問(wèn)題是現(xiàn)階段地球物理探測(cè)技術(shù)提高正演與反演效率、效果必須要解決 的重要課題。

發(fā)明內(nèi)容
針對(duì)上述存在問(wèn)題，本發(fā)明的目的在于提供一種可對(duì)二維復(fù)雜結(jié)構(gòu)三角網(wǎng) 進(jìn)行射線(xiàn)追蹤的全局算法。 其技術(shù)解決方案如下
二維復(fù)雜結(jié)構(gòu)三角網(wǎng)射線(xiàn)追蹤全局方法是以波行面代替波前面，波行面的擴(kuò)展、波的傳遞通過(guò)一個(gè)三角單元與其相鄰三角單元的"振動(dòng)"傳遞進(jìn)行，遇回 轉(zhuǎn)波，會(huì)自行向后傳遞，無(wú)需另行"收縮"，各節(jié)點(diǎn)的次級(jí)源確定和最小走時(shí)計(jì) 算在波行面擴(kuò)展過(guò)程中實(shí)現(xiàn)，其計(jì)算方法按以下步驟進(jìn)行
a) 將二維復(fù)雜結(jié)構(gòu)區(qū)域，用慢度分塊均勻的三角形模型將介質(zhì)參數(shù)化， 在三角形單元的邊界上設(shè)置節(jié)點(diǎn)，
b) 設(shè)源節(jié)點(diǎn)的走時(shí)為零，除源節(jié)點(diǎn)外的其他節(jié)點(diǎn)的走時(shí)為無(wú)窮大，將源 節(jié)點(diǎn)所在的三角形納入波行面三角單元數(shù)組或優(yōu)先隊(duì)列，并使源節(jié)點(diǎn)所在的三 角單元處于激活狀態(tài)，作出存在于波行面數(shù)組的標(biāo)志，
c) 在波行面三角單元數(shù)組中査找處于激活狀態(tài)的三角單元中的最小走時(shí)
節(jié)點(diǎn)，確定這一節(jié)點(diǎn)所在三角單元，將此三角單元作為當(dāng)前三角單元，如果這 一節(jié)點(diǎn)是多個(gè)三角單元的公共節(jié)點(diǎn)，選其中一個(gè)三角單元即可，其它三角單元 會(huì)自行遍歷，將這個(gè)最小走時(shí)的當(dāng)前三角單元及其相鄰三角單元納入波行面三 角單元數(shù)組，并做出存在于波行面數(shù)組的標(biāo)志，以便于不多次存入波行面三角 單元數(shù)組中，
d) 從最小走時(shí)三角單元的最小走時(shí)節(jié)點(diǎn)開(kāi)始，計(jì)算、比較當(dāng)前最小走時(shí) 三角單元每個(gè)節(jié)點(diǎn)的旅行時(shí)間，并確定各節(jié)點(diǎn)的次級(jí)源，
e) 再計(jì)算、比較最小走時(shí)三角單元的相鄰三角單元節(jié)點(diǎn)的旅行時(shí)間及確定 各節(jié)點(diǎn)的次級(jí)源，當(dāng)相鄰三角單元節(jié)點(diǎn)走時(shí)有變化時(shí)，將這一相鄰三角單元作 出激活標(biāo)志，
f) 再一次計(jì)算、比較當(dāng)前最小走時(shí)三角單元節(jié)點(diǎn)旅行時(shí)間，當(dāng)當(dāng)前最小走 時(shí)三角單元各節(jié)點(diǎn)走時(shí)無(wú)變化時(shí)，對(duì)當(dāng)前最小走時(shí)三角單元作出"休眠"標(biāo)志， 回到第c)步，當(dāng)當(dāng)前最小走時(shí)三角單元各節(jié)點(diǎn)走時(shí)有變化時(shí)，對(duì)當(dāng)前最小走 時(shí)三角單元作出"激活"標(biāo)志，回到第c)步，當(dāng)所有波行面的節(jié)點(diǎn)走時(shí)都不再 變化時(shí)且波行面中所有的三角單元均處于"休眠"狀態(tài)，波行面擴(kuò)展完畢。
g) 從接收點(diǎn)開(kāi)始，根據(jù)各節(jié)點(diǎn)走時(shí)及次級(jí)源信息，進(jìn)行向源檢索，拾取 從接收點(diǎn)到源點(diǎn)的射線(xiàn)路徑。
所述的波形面是指二維平面內(nèi)，由波源點(diǎn)及某一時(shí)刻波到達(dá)的每個(gè)節(jié)點(diǎn)構(gòu) 成的平面。
所述的源節(jié)點(diǎn)是指激發(fā)源點(diǎn)。所述的相鄰三角單元是指與一個(gè)三角單元頂點(diǎn)或邊相關(guān)的所有三角單元。
由于采用了以上技術(shù)方案，本發(fā)明針對(duì)二維復(fù)雜結(jié)構(gòu)的射線(xiàn)追蹤問(wèn)題，提 出了以波行面代替波前面的三角網(wǎng)射線(xiàn)追蹤全局方法，數(shù)值模擬表明該算法具 有很高的精度和可靠性，克服了傳統(tǒng)的基于規(guī)則網(wǎng)的射線(xiàn)追蹤方法適應(yīng)性差的 弱點(diǎn)，該方法使規(guī)則網(wǎng)成為此方法中的一個(gè)特例，適用于任意多邊形網(wǎng)，如四 邊形網(wǎng)、四邊形與三角形混合網(wǎng)。對(duì)于三維問(wèn)題，稍加改造也同樣適用。
與一些將波前刻畫(huà)成一條曲線(xiàn)不同，本方法的特點(diǎn)是將所遍歷的三角單 元作為一個(gè)波陣平面， 一方面波行面的擴(kuò)展、波的傳遞通過(guò)最小走時(shí)三角單元 與其相鄰三角單元"振動(dòng)"傳遞進(jìn)行，另一方面當(dāng)遇回轉(zhuǎn)波時(shí)，自行向后傳遞， 無(wú)需另行"收縮"，亦即通過(guò)三角單元的相鄰三角單元將先前已經(jīng)作出"休眠"標(biāo) 志的三角單元"激活"，再次成為最小走時(shí)三角單元。這樣就有效的解決了波的 擴(kuò)展和回傳的問(wèn)題。
二維復(fù)雜區(qū)域模型三角化及二維復(fù)雜結(jié)構(gòu)三角網(wǎng)射線(xiàn)追蹤全局方法，將能 促進(jìn)以射線(xiàn)理論為基礎(chǔ)的波場(chǎng)正演研究工作，并使其正演效率極大提高。


附圖1速度模型參數(shù)化三角單元
附圖2三角頂點(diǎn)的鄰域三角頂點(diǎn)
附圖3節(jié)點(diǎn)的鄰域點(diǎn)，其中，(a)為節(jié)點(diǎn)是三角頂點(diǎn)情形，(b)為節(jié)點(diǎn)是三 角邊的插入點(diǎn)情形，(c)為節(jié)點(diǎn)在三角邊的插入點(diǎn)之間情形，(d)為節(jié)點(diǎn)在三角單 元內(nèi)部情形
附圖4 一個(gè)三角單元的相鄰三角單元
附圖5 —個(gè)三角單元的波前鄰域點(diǎn)走時(shí)計(jì)算
附圖6波前不同近似的比較，其中，(a)為離散波前點(diǎn)近似示意圖，(b)為 連續(xù)波前點(diǎn)近似示意圖
附圖7波行面從(a)到(h)的擴(kuò)展過(guò)程 附圖8各節(jié)點(diǎn)次級(jí)源計(jì)算實(shí)例
附圖9節(jié)點(diǎn)次級(jí)源檢索，其中，(a)為節(jié)點(diǎn)E在三角單元內(nèi)情形，(b)為節(jié) 點(diǎn)E在三角頂點(diǎn)或插入點(diǎn)上情形，(c) (d)為節(jié)點(diǎn)E在邊插入點(diǎn)之間
附圖10 二水平層介質(zhì)模型計(jì)算結(jié)果，其中，(a)為三角單元?jiǎng)澐纸Y(jié)果和初至波射線(xiàn)路徑圖，(b)為初至波波陣面圖，(c)為右邊界走時(shí)對(duì)比，(d)為上 邊界走時(shí)對(duì)比
附圖11水平層狀模型初至波射線(xiàn)路徑及波陣面圖，其中，(a)為剖分網(wǎng)格， (b)為初至波射線(xiàn)路徑及波陣面
附圖12含空洞模型初至波射線(xiàn)路徑及波陣面，其中，(a)為剖分網(wǎng)格，(b) 為初至波射線(xiàn)路徑及波陣面
附圖13含斷層模型初至波射線(xiàn)路徑及波陣面，其中，(a)為剖分網(wǎng)格，(b) 為初至波射線(xiàn)路徑及波陣面
具體實(shí)施例方式
下面結(jié)合附圖及具體實(shí)例對(duì)本發(fā)明的計(jì)算方法作進(jìn)一步詳細(xì)說(shuō)明。 見(jiàn)附圖。
為便于更好理解本發(fā)明的技術(shù)方案，首先對(duì)二維復(fù)雜結(jié)構(gòu)三角網(wǎng)射線(xiàn)追蹤 全局方法涉及到的相關(guān)概念進(jìn)行描述。
如附圖1，為二維區(qū)域三角剖分形成的三角單元。規(guī)定三角單元的三個(gè) 頂點(diǎn)按逆時(shí)針?lè)较蚺帕?，如A—B — C，三邊均為有向線(xiàn)段，如AB、 BC、 CA; 每個(gè)單元內(nèi)速度度均勻。在每個(gè)三角單元的3個(gè)邊上設(shè)置相同多的節(jié)點(diǎn)，節(jié)點(diǎn) 距離可以均勻分配，也可以非均勻分配，這里采用均勻分配，且三角單元公共 邊設(shè)置的節(jié)點(diǎn)相同。
二維區(qū)域內(nèi)的所有點(diǎn)稱(chēng)為節(jié)點(diǎn)，包括源點(diǎn)、接收點(diǎn)，三角單元頂點(diǎn)及三角 邊的插入點(diǎn)等，對(duì)應(yīng)地，稱(chēng)之為源節(jié)點(diǎn)、接收節(jié)點(diǎn)、三角單元頂點(diǎn)節(jié)點(diǎn)及三角 邊插入節(jié)點(diǎn)。
與三角單元的一個(gè)頂點(diǎn)相鄰的所有三角單元頂點(diǎn)節(jié)點(diǎn)稱(chēng)為該三角單元頂 點(diǎn)節(jié)點(diǎn)的鄰域三角頂點(diǎn)，其方向定義為從該點(diǎn)指向鄰域三角頂點(diǎn)，并約定鄰域 三角頂點(diǎn)繞該點(diǎn)逆時(shí)針排放。如附圖2， Po的鄰域三角頂點(diǎn)分別為P, 、 P2 、
p3 、 p4 、 p5 、 p6。
稱(chēng)一個(gè)三角形單元內(nèi)部或邊上的任意兩個(gè)相鄰不同節(jié)點(diǎn)為相鄰點(diǎn)，并約定 其指向性。 一個(gè)節(jié)點(diǎn)的所有相鄰點(diǎn)稱(chēng)為這個(gè)節(jié)點(diǎn)的鄰域點(diǎn)，其方向定義為從該 節(jié)點(diǎn)指向鄰域點(diǎn)，并約定鄰域點(diǎn)繞該節(jié)點(diǎn)逆時(shí)針排放，為減少不必要的鄰域點(diǎn)， 同一邊上，只保留與該節(jié)點(diǎn)緊鄰的節(jié)點(diǎn)。對(duì)區(qū)域三角網(wǎng)而言，節(jié)點(diǎn)的鄰域點(diǎn)在不考慮邊界情況下主要有四種情形， 如附圖3，節(jié)點(diǎn)ss(空心圓)的鄰域點(diǎn)rr (所有實(shí)心圓)，即節(jié)點(diǎn)在三角單 元頂點(diǎn)、節(jié)點(diǎn)在三角邊的插入點(diǎn)、節(jié)點(diǎn)在三角邊的插入點(diǎn)之間、節(jié)點(diǎn)在三角單 元內(nèi)部。對(duì)一個(gè)三角單元而言，對(duì)于附圖3中(a)、 (b)、 (c)三種情形， 一個(gè)節(jié) 點(diǎn)在該三角單元中的鄰域點(diǎn)只需考慮在該三角單元中的鄰域點(diǎn)。相應(yīng)的， 一個(gè) 節(jié)點(diǎn)的鄰域點(diǎn)可能包括三角單元頂點(diǎn)、三角邊的插入點(diǎn)、三角邊的插入點(diǎn)之 間的點(diǎn)、三角單元內(nèi)部點(diǎn)這四類(lèi)節(jié)點(diǎn)。
在計(jì)算從一個(gè)己知節(jié)點(diǎn)走時(shí)到其鄰域點(diǎn)的旅行時(shí)過(guò)程中，該已知走時(shí)節(jié)點(diǎn) 稱(chēng)為其鄰域點(diǎn)的波前點(diǎn)。如附圖3中的ss，其鄰域點(diǎn)稱(chēng)為該節(jié)點(diǎn)的波前鄰域點(diǎn)， 如附圖3中的rr。 一個(gè)己知走時(shí)節(jié)點(diǎn)的次級(jí)源即為波由源傳播到該節(jié)點(diǎn)所經(jīng)歷 路徑上的上一個(gè)鄰域節(jié)點(diǎn)。因此，在計(jì)算波前點(diǎn)到其波前鄰域點(diǎn)的旅行時(shí)過(guò)程 中，各波前點(diǎn)可以相繼看做是其波前鄰域點(diǎn)的次級(jí)源，但只有當(dāng)遍歷所有波前 點(diǎn)時(shí)，各節(jié)點(diǎn)的次級(jí)源才是使該節(jié)點(diǎn)走時(shí)最小的上一個(gè)鄰域節(jié)點(diǎn)。
與一個(gè)三角單元頂點(diǎn)或邊相關(guān)的所有三角單元，稱(chēng)為該三角單元的相鄰三 角單元。如附圖4，三角單元AQ的相鄰三角單元包括A, A。12個(gè)三角單元。 在下面將要定義的波行面的擴(kuò)展與回轉(zhuǎn)波的回傳就是通過(guò)相鄰三角單元來(lái)擴(kuò) 展和傳遞的。
在二維平面內(nèi)，某一時(shí)刻由波源點(diǎn)及波到達(dá)的每個(gè)節(jié)點(diǎn)構(gòu)成的平面，稱(chēng)為 波行面。在三角網(wǎng)射線(xiàn)追蹤全局算法中，波行面的表現(xiàn)形式為，波遍歷的三角 單元組成的平面。若區(qū)域剖分為其它形式的多邊形網(wǎng)，如四邊形網(wǎng)、四邊形與 三角形混合網(wǎng)，則波行面的表現(xiàn)形式為，波遍歷的四邊形單元、或四變形單元 與三角形單元混合組成的平面。值得說(shuō)明的是，計(jì)算各節(jié)點(diǎn)的走時(shí)過(guò)程中，在 完成整個(gè)區(qū)域的節(jié)點(diǎn)走時(shí)計(jì)算之前，波行面各節(jié)點(diǎn)的走時(shí)不一定必須是最小走 時(shí)。
由上述概念，就可得到三角網(wǎng)射線(xiàn)追蹤算法中運(yùn)用較頻繁的有向邊、節(jié)點(diǎn)、 三角單元等數(shù)據(jù)結(jié)構(gòu)。 有向邊類(lèi)
class Edge 〃有向邊類(lèi) public:Vertex *dest; double cost; Triangle ^lefttriangle; Triangle *righttriangle;
};
節(jié)點(diǎn)類(lèi)
class Vertex
public:
double xj double y5
〃第二個(gè)點(diǎn)位置
〃邊的代價(jià)，可以是距離，也可以是時(shí)間； 〃邊的左三角形 〃邊的右三角形
〃節(jié)點(diǎn)X坐標(biāo)
〃節(jié)點(diǎn)y坐標(biāo)
dcllist<Edge> adj—tri— vex ;〃該點(diǎn)相鄰三角頂點(diǎn) dcllist<Edge> adj—vex; 〃該點(diǎn)鄰±或點(diǎn)，
Vertex *prev—vex; bool scratchy double density; double travel time;
〃dcllist為雙向循環(huán)鏈表模板
〃該點(diǎn)的次級(jí)源 〃搜索標(biāo)志
〃該點(diǎn)三角單元密度控制大小 〃由源點(diǎn)到該節(jié)點(diǎn)的旅行時(shí)
};
三角單元類(lèi)
class Triangle〃三角形類(lèi) public:
dcllist<Vertex *> tripoly;〃由三角單元頂點(diǎn)及邊插入點(diǎn)構(gòu)成的三角環(huán)
Vertex *a， *b, *c; 〃三角單元三個(gè)頂點(diǎn)
double v; 〃三角單元波速
bool scratch; 〃是否被遍歷標(biāo)志(休眠、激活標(biāo)志)
bool insert; 〃是否進(jìn)入波陣平面標(biāo)志
double min—time; 〃最小走時(shí)
DCLLIST::iterator min—time—itr; 〃最小走時(shí)位置，〃DCLLIST::iterator,三角環(huán)迭代位置 dcllist<Triangle *> adj—tri; 〃相鄰三角單元
};
二維復(fù)雜結(jié)構(gòu)三角網(wǎng)射線(xiàn)追蹤全局方法是利用波行面的擴(kuò)展完成節(jié)點(diǎn)最 小走時(shí)計(jì)算和次級(jí)源的確定，而波行面的擴(kuò)展是通過(guò)一個(gè)單元與其相鄰單元的 "振動(dòng)"傳遞完成。
對(duì)于一個(gè)復(fù)雜區(qū)域，采用三角剖分方法將區(qū)域介質(zhì)參數(shù)化，視每個(gè)三角單 元慢度均勻，并設(shè)置節(jié)點(diǎn)及三角網(wǎng)的一些拓?fù)潢P(guān)系。三角網(wǎng)射線(xiàn)追蹤全局計(jì)算
方法的基本步驟如下
第一步，利用波行面的擴(kuò)展完成節(jié)點(diǎn)最小走時(shí)計(jì)算和次級(jí)源的確定。以波 行面代替波前面，波行面的擴(kuò)展、波的傳遞通過(guò)三角單元與其相鄰三角單元的 "振動(dòng)"傳遞進(jìn)行，遇回轉(zhuǎn)波，會(huì)自行向后傳遞，無(wú)需另行"收縮"。
第二步，利用計(jì)算出的各節(jié)點(diǎn)的旅行時(shí)和次級(jí)源信息拾取射線(xiàn)路徑。即根 據(jù)每個(gè)節(jié)點(diǎn)的最小走時(shí)、次級(jí)源，從接收點(diǎn)出發(fā)，向源拾取射線(xiàn)路徑。
二維復(fù)雜結(jié)構(gòu)三角網(wǎng)射線(xiàn)追蹤全局方法的實(shí)現(xiàn)過(guò)程分以下幾方面。
1設(shè)置各節(jié)點(diǎn)的拓?fù)潢P(guān)系
將二維復(fù)雜結(jié)構(gòu)區(qū)域，用慢度分塊均勻的三角形模型將介質(zhì)參數(shù)化，在三 角形單元的邊界上設(shè)置節(jié)點(diǎn)，同時(shí)設(shè)置各節(jié)點(diǎn)之間的拓?fù)潢P(guān)系，主要包括以下
幾類(lèi)節(jié)點(diǎn)的拓?fù)潢P(guān)系
(1) 按逆時(shí)針?lè)较?，設(shè)置各節(jié)點(diǎn)的鄰域三角頂點(diǎn)，這里的節(jié)點(diǎn)包括源點(diǎn)、 接收點(diǎn)、三角頂點(diǎn)及三角邊的插入點(diǎn)，同時(shí)設(shè)置由節(jié)點(diǎn)至鄰域三角頂點(diǎn)形成的 有向邊的左、右三角形，主要目的是便于査詢(xún)?nèi)切蔚南噜徣切巍?br> (2) 按逆時(shí)針?lè)较?，設(shè)置各節(jié)點(diǎn)的鄰域點(diǎn)，這里的節(jié)點(diǎn)包括源點(diǎn)、接收 點(diǎn)、三角頂點(diǎn)及三角邊的插入點(diǎn)，同時(shí)設(shè)置由節(jié)點(diǎn)至其鄰域點(diǎn)形成的有向邊的 左、右三角形，主要目的是便于計(jì)算由波前點(diǎn)至波前鄰域點(diǎn)的旅行時(shí)。將接收 點(diǎn)和源點(diǎn)等同于三角頂點(diǎn)及三角邊的插入點(diǎn)進(jìn)行設(shè)置的一個(gè)優(yōu)點(diǎn)是使得射線(xiàn) 追蹤不必對(duì)源點(diǎn)和接收點(diǎn)另行計(jì)算，路徑拾取只需給出接收點(diǎn)坐標(biāo)或接收點(diǎn)編 號(hào)即可。
(3) 按逆時(shí)針?lè)较?，設(shè)置各三角單元的相鄰三角單元，根據(jù)三角頂點(diǎn)的 鄰域三角頂點(diǎn)設(shè)置，主要目的是便于三角單元之間波的傳遞和波行面的擴(kuò)展。2波前鄰域點(diǎn)最小走時(shí)及次級(jí)源位置計(jì)算
波前鄰域點(diǎn)最小走時(shí)計(jì)算及次級(jí)源的確定只需考慮在一個(gè)三角單元內(nèi)進(jìn) 行，本步驟藴含于波行面的擴(kuò)展中，波行面擴(kuò)展細(xì)節(jié)見(jiàn)下一部分的描述。
附圖5是一個(gè)三角單元，設(shè)ss點(diǎn)是波前節(jié)點(diǎn)，其波前鄰域點(diǎn)用實(shí)心圓表 示，其中rr節(jié)點(diǎn)是ss的某 一鄰域點(diǎn)，它們的坐標(biāo)分別為(Xss,hs)和(x ，_yrr); (xrs,yrs) 為rr節(jié)點(diǎn)次級(jí)源的坐標(biāo)；^為ss節(jié)點(diǎn)的走時(shí)；^為rr節(jié)點(diǎn)的走時(shí)；v為三角 單元的速度；^為波從ss節(jié)點(diǎn)以本三角單元速度v沿直線(xiàn)傳播到rr節(jié)點(diǎn)所需 的時(shí)間
=ik, - &)2 + )2 ]2 /v ( j )
于是可以求出rr節(jié)點(diǎn)的走時(shí)及相應(yīng)的次級(jí)源節(jié)點(diǎn)坐標(biāo)，即當(dāng)(^+ ^ 時(shí)，有
^ = f + f
V尸 *\sv
'A., =
H、. (2) 在一個(gè)三角單元內(nèi)計(jì)算節(jié)點(diǎn)走時(shí)時(shí)，從這個(gè)三角單元走時(shí)最小的節(jié)點(diǎn)ss 開(kāi)始，然后讓ss遍歷本三角單元的所有節(jié)點(diǎn)，本三角單元的節(jié)點(diǎn)走時(shí)計(jì)算才算 完畢，但各節(jié)點(diǎn)走時(shí)不一定就是最小走時(shí)，只有當(dāng)ss經(jīng)歷所有波前點(diǎn)時(shí)，r 才 變成最小走時(shí)，(AS ,;^)才是真正次級(jí)源坐標(biāo)。
這里的次級(jí)源是用節(jié)點(diǎn)近似，波前鄰域點(diǎn)的走時(shí)是用波前點(diǎn)走時(shí)加上波前 點(diǎn)到波前鄰域點(diǎn)局部走時(shí)的和表示的，相當(dāng)于把波前面看成是一系列離散點(diǎn)次 級(jí)源，如附圖6(a)。實(shí)際上到達(dá)每一節(jié)點(diǎn)的波并不一定正好是從網(wǎng)格上的波前 點(diǎn)發(fā)出，而可能是從網(wǎng)格點(diǎn)之間的點(diǎn)發(fā)出，如附圖6(b)。這樣局部走時(shí)和次級(jí) 源位置就要進(jìn)行修正。
如附圖6(b)，要計(jì)算節(jié)點(diǎn)E的最小走時(shí)和次級(jí)源的位置，假定E的次級(jí)源 在A、 B、 C三點(diǎn)之間，如D點(diǎn)，設(shè)某一點(diǎn)為局部原點(diǎn)(如A點(diǎn))，在A， B， C三 點(diǎn)之間走時(shí)函數(shù)表示為"x)， x為AC方向上局部一維坐標(biāo)。設(shè)前方點(diǎn)E在AC 方向上的局部坐標(biāo)為XE，垂直AC方向的距離為么則經(jīng)AC線(xiàn)上任一點(diǎn)x到 達(dá)鄰域點(diǎn)E的走時(shí)為
^ =,(x) + [(x-x￡)2+4/v (3)v為三角單元內(nèi)的速度，令x在AC之間變化時(shí)，,e取極小，可求出波前 點(diǎn)即次級(jí)源點(diǎn)D:
若A、 B、 C 3點(diǎn)的局部坐標(biāo)為xA、 xb、 xc， 3點(diǎn)的走時(shí)為"，,B， /c，則 用雙曲線(xiàn)近似計(jì)算AC上任一點(diǎn)走時(shí)
容易證明，當(dāng)介質(zhì)均勻時(shí)，"x)是某一點(diǎn)源引起的走時(shí)時(shí)，上式是嚴(yán)格精 確的。當(dāng)介質(zhì)不均勻時(shí)，雙曲線(xiàn)近似也是3點(diǎn)走時(shí)插值方法中精度最高的。(4) 式是一個(gè)極小值求解問(wèn)題，可采用梯度法、黃金搜索等方法求得xd。按邊讓A、 B、 C歷遍除E點(diǎn)所在邊的其它兩邊的節(jié)點(diǎn)就完成了 E點(diǎn)在這個(gè)三角形中的最 小走時(shí)和次級(jí)源位置的計(jì)算。
3波行面擴(kuò)展
在基于矩形網(wǎng)的射線(xiàn)追蹤的全局方法中，常利用一矩形模擬波陣面的擴(kuò)張 與收縮來(lái)計(jì)算節(jié)點(diǎn)最小走時(shí)和次級(jí)源位置，但在三角網(wǎng)中，要擬定某一時(shí)刻的 波陣線(xiàn)將會(huì)非常復(fù)雜，而且采用波陣面，也不利于回轉(zhuǎn)波的傳遞。
本發(fā)明的做法是，利用波行平面，將波行平面的擴(kuò)展視為一個(gè)遞歸過(guò)程， 其算法描述如下
(1) 設(shè)源點(diǎn)的走時(shí)為零，除源點(diǎn)外的其他節(jié)點(diǎn)的走時(shí)為無(wú)窮大。將源點(diǎn) 所在的三角形納入波行面三角單元存儲(chǔ)容器，并使這個(gè)三角單元處于激活狀 態(tài)，作出存在于波行面的標(biāo)志。
(2) 在波行面三角單元存儲(chǔ)容器中查找處于激活狀態(tài)的三角單元中的最
小走時(shí)節(jié)點(diǎn)，確定這一節(jié)點(diǎn)所在三角單元，將此三角單元作為當(dāng)前三角單元；
如果這一節(jié)點(diǎn)是多個(gè)三角單元的公共點(diǎn)，選其中一個(gè)三角單元即可，其它三角 單元會(huì)自行遍歷。將這個(gè)最小走時(shí)的當(dāng)前三角單元及其相鄰三角單元納入波行 面三角單元存儲(chǔ)容器，并做出存在于波行面的標(biāo)志，便于不多次存入波行面三 角單元存儲(chǔ)容器。
=min 《x) + ((x-x五))2(3)利用前述"波前鄰域點(diǎn)最小走時(shí)及次級(jí)源位置計(jì)算"，計(jì)算、比較當(dāng) 前最小走時(shí)三角單元節(jié)點(diǎn)旅行時(shí)(從本三角單元最小走時(shí)節(jié)點(diǎn)開(kāi)始，下同)；再 計(jì)算、比較最小走時(shí)三角單元的相鄰三角單元節(jié)點(diǎn)的旅行時(shí)，當(dāng)相鄰三角單元 節(jié)點(diǎn)走時(shí)有變化時(shí)，將這一相鄰三角單元作出激活標(biāo)志；再一次計(jì)算、比較當(dāng) 前最小走時(shí)三角單元節(jié)點(diǎn)旅行時(shí)，當(dāng)當(dāng)前最小走時(shí)三角單元各節(jié)點(diǎn)走時(shí)無(wú)變化 時(shí)，對(duì)當(dāng)前最小走時(shí)三角單元作出"休眠"標(biāo)志，回到第(2)步，當(dāng)當(dāng)前最小 走時(shí)三角單元各節(jié)點(diǎn)走時(shí)有變化時(shí)，對(duì)當(dāng)前最小走時(shí)三角單元作出"激活"標(biāo) 志，回到第(2)步；以次類(lèi)推，遍歷所有三角單元。當(dāng)所有波行面的節(jié)點(diǎn)走 時(shí)都不再變化時(shí)且所有波行面中的三角單元均處于"休眠"狀態(tài)，波行面擴(kuò)展完 畢。
與一些將波前刻畫(huà)成一條曲線(xiàn)不同，本計(jì)算方法的特點(diǎn)是將所遍歷的三 角單元作為一個(gè)波陣平面，以波行面代替波前面； 一方面通過(guò)最小走時(shí)三角單 元的與其相鄰三角單元間的"振動(dòng)"傳遞擴(kuò)展波行面，另一方面當(dāng)遇回轉(zhuǎn)波時(shí)， 也通過(guò)三角單元的相鄰三角單元之間的傳遞將先前已經(jīng)"休眠"的三角單元"激 活"，而再次成為最小走時(shí)三角單元，無(wú)需另行"收縮"。這樣就有效的解決了波 的擴(kuò)展和回傳的問(wèn)題。
附圖7展示了波行面從(a)到(h)的擴(kuò)展過(guò)程，源點(diǎn)在左上角，粗線(xiàn)三角單元 表示當(dāng)前三角單元。
附圖8為一計(jì)算實(shí)例，源點(diǎn)在左上角，顯示出各節(jié)點(diǎn)的次級(jí)源，如A的次 級(jí)源為B。
4波傳播路徑的向源拾取
任意點(diǎn)到源點(diǎn)的射線(xiàn)路徑通過(guò)節(jié)點(diǎn)最小走時(shí)和次級(jí)源信息向源追索完成。 本發(fā)明將源點(diǎn)和接收點(diǎn)按照節(jié)點(diǎn)統(tǒng)一處理，因而，在向源路徑拾取過(guò)程中不必 對(duì)接收點(diǎn)走時(shí)另行計(jì)算，直接追索即可。附圖9表示節(jié)點(diǎn)E的次級(jí)源F的檢索， 根據(jù)E的位置，有以下幾類(lèi)節(jié)點(diǎn)的次級(jí)源檢索
(1) 當(dāng)節(jié)點(diǎn)E位于三角單元內(nèi)部時(shí)，如附圖9(a)，通過(guò)雙曲線(xiàn)近似和最 小走時(shí)搜索從三角單元三邊上找到次級(jí)源，算出該節(jié)點(diǎn)走時(shí)(這一類(lèi)節(jié)點(diǎn)可在 波行面擴(kuò)展中統(tǒng)一處理，這里只需檢索即可)；
(2) 當(dāng)節(jié)點(diǎn)E是三角頂點(diǎn)或三角邊的插入點(diǎn)時(shí)(9(b))，可直接拾取次級(jí)源；
(3)當(dāng)節(jié)點(diǎn)E位于三角邊插入點(diǎn)A、 C之間時(shí)，檢索E的次級(jí)源F要依 據(jù)A的次級(jí)源B與C的次級(jí)源D。這時(shí)主要會(huì)遇到如附圖9(c)、 (d)兩種情況 (c)情形是A的次級(jí)源B與C的的次級(jí)源D在一個(gè)三角單元內(nèi)；(d)情形就比較 復(fù)雜，由于C點(diǎn)是三角頂點(diǎn)，C的次級(jí)源D可能存在于其它三角單元內(nèi)。對(duì)附 圖9(c)、 (d)兩種情況，也通過(guò)雙曲線(xiàn)近似和最小走時(shí)搜索從三角單元相應(yīng)邊上 尋找次級(jí)源，并算出該節(jié)點(diǎn)走時(shí)。 5算法精度檢驗(yàn)
為了檢驗(yàn)三角網(wǎng)射線(xiàn)追蹤算法的精度和可靠性，對(duì)附圖lO(a)所示模型的三 角網(wǎng)射線(xiàn)追蹤正演模擬結(jié)果同根據(jù)Snell定理計(jì)算的理論結(jié)果進(jìn)行了對(duì)比。源 位于左邊界上端，接收點(diǎn)分別位于模型上邊界和右邊界，接收點(diǎn)距2m。第1 層介質(zhì)速度為2000m/s，第2層介質(zhì)速度為4000m/s，模型水平方向?qū)挾葹?0m， 垂直方向深度分別是第1層10m，第2層40m。在對(duì)模型進(jìn)行網(wǎng)絡(luò)劃分時(shí)，采 用密度為2m的函數(shù)進(jìn)行三角剖分。附圖10(a)為三角網(wǎng)剖分結(jié)果和初至波射線(xiàn) 路徑圖(粗線(xiàn))；附圖10(b)為初至波波陣面圖，等時(shí)線(xiàn)間隔lms。附圖10(c)、 (d) 分別為右邊界和上邊界走時(shí)對(duì)比結(jié)果，圖中實(shí)線(xiàn)為用解析分析公式計(jì)算的理論 初至波走時(shí)，圓點(diǎn)為三角網(wǎng)追蹤算法計(jì)算的初至波走時(shí)。右邊界最大走時(shí)相對(duì) 誤差為0.0009%;上邊界最大走時(shí)相對(duì)誤差僅為0.00028%。由此可以看出，采 用三角網(wǎng)射線(xiàn)追蹤算法具有相當(dāng)高的計(jì)算精度。 具體實(shí)施例
實(shí)施例一 水平層狀均勻介質(zhì)模型
附圖11為水平層狀模型的三角剖分及射線(xiàn)追蹤結(jié)果。源點(diǎn)位于左上角， 接收點(diǎn)分別位于模型上邊界、右邊界及下邊界。模型示于附圖ll(a)，第1層 介質(zhì)速度為2000m/s，第2層介質(zhì)速度為4000m/s，第3層介質(zhì)速度為3000m/s， 第4層介質(zhì)速度為5000m/s，模型水平方向?qū)挾葹?0m，垂直方向4個(gè)層厚度 均為10m，對(duì)模型采用密度為2m的函數(shù)進(jìn)行三角剖分。附圖ll(a)為三角剖分 結(jié)果，附圖ll(b)為采用本算法得到的不同時(shí)刻的波陣面和射線(xiàn)路徑。
實(shí)施例二含空洞模型
附圖12為含空洞模型的三角剖分及射線(xiàn)追蹤結(jié)果。源點(diǎn)位于左上角，接均勻密度控制函數(shù)。附圖12(a)為三 角剖分結(jié)果，附圖12(b)為采用本算法得到的不同時(shí)刻的波陣面和射線(xiàn)路徑。 實(shí)施例三含斷層模型
附圖13為含斷層模型的三角剖分及射線(xiàn)追蹤結(jié)果。源點(diǎn)位于左下角，接 收點(diǎn)分別位于模型外邊界。模型示于附圖13(a)，在區(qū)域內(nèi)設(shè)置了斷層子域Q2 及低速子域Qs，為模擬03到Qs速度較大的梯度變化，增加網(wǎng)格加密子域Q4， 網(wǎng)格剖分采用了不均勻密度控制函數(shù)，各子域的模擬速度Q1: vp=3000m/s， Q2: Vp=2000m/s， Q3: vp=4000m/s， Q4: vp=4000m/s， Q5: vp =3000m/s。附圖13(a) 為三角剖分結(jié)果，附圖13(b)為采用本算法得到的不同時(shí)刻的波陣
權(quán)利要求
1 二維復(fù)雜結(jié)構(gòu)三角網(wǎng)射線(xiàn)追蹤全局方法，其特征在于二維復(fù)雜結(jié)構(gòu)三角網(wǎng)射線(xiàn)追蹤全局方法是以波行面代替波前面，波行面的擴(kuò)展、波的傳遞通過(guò)一個(gè)三角單元與其相鄰三角單元的“振動(dòng)”傳遞進(jìn)行，遇回轉(zhuǎn)波，會(huì)自行向后傳遞，無(wú)需另行“收縮”，各節(jié)點(diǎn)的次級(jí)源確定和最小走時(shí)計(jì)算在波行面擴(kuò)展過(guò)程中實(shí)現(xiàn)，其計(jì)算方法按以下步驟進(jìn)行a)將二維復(fù)雜結(jié)構(gòu)區(qū)域，用慢度分塊均勻的三角形模型將介質(zhì)參數(shù)化，在三角形單元的邊界上設(shè)置節(jié)點(diǎn)，b)設(shè)源節(jié)點(diǎn)的走時(shí)為零，除源節(jié)點(diǎn)外的其他節(jié)點(diǎn)的走時(shí)為無(wú)窮大，將源節(jié)點(diǎn)所在的三角形納入波行面三角單元數(shù)組或優(yōu)先隊(duì)列，并使源節(jié)點(diǎn)所在的三角單元處于激活狀態(tài)，作出存在于波行面數(shù)組的標(biāo)志，c)在波行面三角單元數(shù)組中查找處于激活狀態(tài)的三角單元中的最小走時(shí)節(jié)點(diǎn)，確定這一節(jié)點(diǎn)所在三角單元，將此三角單元作為當(dāng)前三角單元，如果這一節(jié)點(diǎn)是多個(gè)三角單元的公共節(jié)點(diǎn)，選其中一個(gè)三角單元即可，其它三角單元會(huì)自行遍歷，將這個(gè)最小走時(shí)的當(dāng)前三角單元及其相鄰三角單元納入波行面三角單元數(shù)組，并做出存在于波行面數(shù)組的標(biāo)志，以便于不多次存入波行面三角單元數(shù)組中，d)從最小走時(shí)三角單元的最小走時(shí)節(jié)點(diǎn)開(kāi)始，計(jì)算、比較當(dāng)前最小走時(shí)三角單元每個(gè)節(jié)點(diǎn)的旅行時(shí)間，并確定各節(jié)點(diǎn)的次級(jí)源，e)再計(jì)算、比較最小走時(shí)三角單元的相鄰三角單元節(jié)點(diǎn)的旅行時(shí)間及確定各節(jié)點(diǎn)的次級(jí)源，當(dāng)相鄰三角單元節(jié)點(diǎn)走時(shí)有變化時(shí)，將這一相鄰三角單元作出激活標(biāo)志，f)再一次計(jì)算、比較當(dāng)前最小走時(shí)三角單元節(jié)點(diǎn)旅行時(shí)間，當(dāng)當(dāng)前最小走時(shí)三角單元各節(jié)點(diǎn)走時(shí)無(wú)變化時(shí)，對(duì)當(dāng)前最小走時(shí)三角單元作出“休眠”標(biāo)志，回到第c)步，當(dāng)當(dāng)前最小走時(shí)三角單元各節(jié)點(diǎn)走時(shí)有變化時(shí)，對(duì)當(dāng)前最小走時(shí)三角單元作出“激活”標(biāo)志，回到第c)步，當(dāng)所有波行面的節(jié)點(diǎn)走時(shí)都不再變化時(shí)且波行面中所有的三角單元均處于“休眠”狀態(tài)，波行面擴(kuò)展完畢，g)從接收點(diǎn)開(kāi)始，根據(jù)各節(jié)點(diǎn)走時(shí)及次級(jí)源信息，進(jìn)行向源檢索，拾取從接收點(diǎn)到源點(diǎn)的射線(xiàn)路徑。
2 如權(quán)利要求1所述的二維復(fù)雜結(jié)構(gòu)三角網(wǎng)射線(xiàn)追蹤全局方法，其特征在 于所述的波形面是指二維平面內(nèi)，由波源點(diǎn)及某一時(shí)刻波到達(dá)的每個(gè)節(jié)點(diǎn)構(gòu) 成的平面。
3 如權(quán)利要求1所述的二維復(fù)雜結(jié)構(gòu)三角網(wǎng)射線(xiàn)追蹤全局方法，其特征在 于所述的源節(jié)點(diǎn)是指激發(fā)源點(diǎn)。
4 如權(quán)利要求1所述的二維復(fù)雜結(jié)構(gòu)三角網(wǎng)射線(xiàn)追蹤全局方法，其特征在 于所述的相鄰三角單元是指與一個(gè)三角單元頂點(diǎn)或邊相關(guān)的所有三角單元。
全文摘要
本發(fā)明涉及二維復(fù)雜結(jié)構(gòu)的走時(shí)層析成像反演和以Maslov射線(xiàn)理論為基礎(chǔ)的波場(chǎng)計(jì)算的射線(xiàn)追蹤方法。對(duì)于二維復(fù)雜外部幾何邊界及區(qū)域內(nèi)部復(fù)雜空間結(jié)構(gòu)，用慢度分塊均勻的三角形模型將介質(zhì)參數(shù)化，在三角形單元的邊界上設(shè)置節(jié)點(diǎn)。以波行面代替波前面，波行面的擴(kuò)展、波的傳遞通過(guò)三角單元與其相鄰三角單元的“振動(dòng)”傳遞進(jìn)行，遇回轉(zhuǎn)波，會(huì)自行向后傳遞，無(wú)需另行“收縮”，各節(jié)點(diǎn)的次級(jí)源位置和最小走時(shí)計(jì)算在波行面擴(kuò)展過(guò)程中實(shí)現(xiàn)；利用各節(jié)點(diǎn)次級(jí)源走時(shí)與方向信息，通過(guò)最小走時(shí)搜索，拾取從接收點(diǎn)到源點(diǎn)的射線(xiàn)路徑。本計(jì)算方法適用于任意多邊形網(wǎng)。
文檔編號(hào)G01V1/30GK101533102SQ20091006147
公開(kāi)日2009年9月16日 申請(qǐng)日期2009年4月9日 優(yōu)先權(quán)日2009年4月9日
發(fā)明者丁亮清, 于師建, 余才盛, 況碧波, 劉新志, 劉方文, 劉潤(rùn)澤, 周習(xí)軍, 喻維鋼, 智 張, 張建清, 李張明, 李文忠, 熊永紅, 程含發(fā), 陸二男, 華 陳 申請(qǐng)人:長(zhǎng)江工程地球物理勘測(cè)武漢有限公司