曲線坐標(biāo)系程函方程的求解方法研究

劉一峰，蘭海強(qiáng)

1中國(guó)科學(xué)院地質(zhì)與地球物理研究所，巖石圈演化國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100029

2中海石油（中國(guó)）有限公司北京研究中心，北京 100027

3中國(guó)科學(xué)院研究生院，北京 100049

曲線坐標(biāo)系程函方程的求解方法研究

劉一峰1，2，3，蘭海強(qiáng)1，3

1中國(guó)科學(xué)院地質(zhì)與地球物理研究所，巖石圈演化國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100029

2中海石油（中國(guó)）有限公司北京研究中心，北京 100027

3中國(guó)科學(xué)院研究生院，北京 100049

笛卡爾坐標(biāo)系中經(jīng)典的程函方程在靜校正、疊前偏移、走時(shí)反演、地震定位、層析成像等許多地球物理工作都有應(yīng)用，然而用其計(jì)算起伏地表的地震波走時(shí)時(shí)卻比較困難.我們通過把曲線坐標(biāo)系中的矩形網(wǎng)格映射到笛卡爾坐標(biāo)系的貼體網(wǎng)格推導(dǎo)出了曲線坐標(biāo)中的程函方程，此時(shí)，曲線坐標(biāo)系的程函方程呈現(xiàn)為各向異性的程函方程（盡管在笛卡爾坐標(biāo)系中介質(zhì)是各向同同性的）.然后嘗試用求解各向同性程函方程的快速推進(jìn)法和Lax－Friedrichs快速掃描算法來分別求解該方程.數(shù)值試驗(yàn)表明未加考慮各向異性程函方程與各向同性程函方程的差別而把求解各向同性程函方程的快速推進(jìn)法直接拓展到曲線坐標(biāo)中的程函方程的做法是錯(cuò)誤的，而Lax－Friedrichs快速掃描算法總能穩(wěn)定地求解曲線坐標(biāo)系的程函方程，進(jìn)而有效地處理了地表起伏的情況，得到穩(wěn)定準(zhǔn)確的計(jì)算結(jié)果.

程函方程，起伏地表，各向異性，快速推進(jìn)法，Lax－Friedrichs快速掃描算法

1 引 言

地震波走時(shí)的計(jì)算在地球物理學(xué)中起著極其重要的作用[1－6].例如，確定震中位置時(shí)，較古老的方法是，以所有地震臺(tái)觀測(cè)到的走時(shí)的剩余值最小為準(zhǔn)則來求解；Kirchhoff偏移需要計(jì)算Green函數(shù)，而Green函數(shù)描述了觀測(cè)點(diǎn)與地下介質(zhì)內(nèi)部點(diǎn)之間的走時(shí)[7－8]；走時(shí)反演和層析成像更是直接利用地表觀測(cè)到的走時(shí)數(shù)據(jù)來直接反演地下的速度分布，進(jìn)而構(gòu)建地下的巖性結(jié)構(gòu)[9－11].傳統(tǒng)的計(jì)算走時(shí)的方法為射線方法，其原理是沿程函方程的特征方向，即射線方向求解走時(shí)，經(jīng)插值計(jì)算后得到地下規(guī)則網(wǎng)格點(diǎn)上的地震波走時(shí).但是，對(duì)于復(fù)雜構(gòu)造的地質(zhì)模型，這種古典射線方法將遇到焦散面及陰影區(qū)等問題[8，12－13].Vidale[14－15]基 于 擴(kuò) 張 波 前 的 思 想 開 創(chuàng) 性地提出一種用有限差分方法來近似程函方程，但該方法采用擴(kuò)張矩形來追蹤波前，在一定的速度分布情況下，計(jì)算出的走時(shí)并不是最小，且在處理強(qiáng)速度界面時(shí)會(huì)出現(xiàn)不穩(wěn)定現(xiàn)象.Qin等[16]改進(jìn)了Vidale的方法，盡可能沿?cái)U(kuò)張的波前面來計(jì)算走時(shí)，方法和Vidale基本相同但他考慮到了因果關(guān)系，首先尋求上一個(gè)近似波前面的走時(shí)全局極小點(diǎn)，然后向外擴(kuò)張；但該方法計(jì)算機(jī)實(shí)現(xiàn)較困難，大部分時(shí)間要用于尋找全局極值，效率不高.Van Trier等[17]將迎風(fēng)有限差分法引入解程函方程，大大的提高了差分格式的穩(wěn)定性.Sethian等[18－21]提出了一種稱之為快速推進(jìn)的方法（Fast marching method，F(xiàn)MM），該方法利用迎風(fēng)差分格式求解局部程函方程，采用窄帶延拓重建走時(shí)波前，利用堆選排技術(shù)保存走時(shí)，將最小走時(shí)放在堆的頂部.該方法顯著縮短了尋找極小值的時(shí)間，計(jì)算量由波前擴(kuò)展法的O（N3）減少到O（N·logN）（logN由堆的排序算法產(chǎn)生），其中N是節(jié)點(diǎn)數(shù).近年來，很多學(xué)者對(duì)快速推進(jìn)法進(jìn)行了推廣和應(yīng)用[22－29].最近，Zhao等[30]提出了一種稱之為快速掃描的方法（Fast sweeping method，F(xiàn)SM）用于求解一階雙曲型偏微分方程，并指出該方法的計(jì)算量為O（N），且證明了該算法的單調(diào)性和穩(wěn)定性.實(shí)際上快速掃描法和快速推進(jìn)法都是求解相同的離散方程，且都是基于因果關(guān)系沿程函方程的特征方向求解，但前者更強(qiáng)健且對(duì)高階方程更靈活，效率也一般較后者高[31－32].

前面提到的計(jì)算地震波走時(shí)的方法基本都是基于水平地表的.然而，在我國(guó)隨著東部平原地帶大規(guī)模油氣勘探工作已接近尾聲，目前的油氣勘探重點(diǎn)已轉(zhuǎn)移到了西部、西南部地區(qū).與東部地區(qū)不同，西部地區(qū)劇烈的地形起伏給地震勘探工作提出了嚴(yán)峻的挑戰(zhàn)[33]，傳統(tǒng)的基于平緩構(gòu)造勘探的地震數(shù)據(jù)采集、處理和解釋方法在這類復(fù)雜地表地區(qū)已不再適用[34－35].研究復(fù)雜地表?xiàng)l件下地震波走時(shí)的計(jì)算問題，對(duì)于在這些地區(qū)進(jìn)行反射地震偏移成像、走時(shí)層析反演均有非常重要的意義.到目前為至，絕大多數(shù)的可以處理起伏地表的地震波走時(shí)計(jì)算方法都是基于非規(guī)則 網(wǎng) 格 的[19，23，25，36－41].然 而，目 前 許 多 地 震波場(chǎng)數(shù)值模擬方法都是基于規(guī)則網(wǎng)格的[34，42－48]，這使在這些網(wǎng)格上的走時(shí)計(jì)算變得愈加重要，因?yàn)榛谝陨系恼萁Y(jié)果進(jìn)行Kirchhoff偏移時(shí)需要這些走時(shí)信息.因此，基于規(guī)則網(wǎng)格剖分的起伏地表下的走時(shí)計(jì)算具有重要意義.最近，Sun等[49]提出了一種在規(guī)則網(wǎng)格上求解起伏地表下地震波走時(shí)的方法，該方法的關(guān)鍵是在起伏地表附近采用非均一的網(wǎng)格間距且通過引入地表點(diǎn)、地表上點(diǎn)、界面點(diǎn)、界面下點(diǎn)等概念使其適于用FMM求解.Lan ＆Zhang[50－51]發(fā)展了一種求解起伏地表下地震波走時(shí)的方法.首先，他們借助于貼體網(wǎng)格和坐標(biāo)變換，把笛卡爾坐標(biāo)系中的程函方程變換到曲線坐標(biāo)系中，推導(dǎo)出了曲線坐標(biāo)系的程函方程.此時(shí)，曲線坐標(biāo)系的程函方程呈現(xiàn)為各向異性的程函方程（盡管在笛卡爾坐標(biāo)系中介質(zhì)是各向同同性的）.在各向同性介質(zhì)中，射線速度方向（群速度）與走時(shí)梯度的方向（相速度）一致，此時(shí)，在傳播的走時(shí)場(chǎng)中，走時(shí)的梯度方向指示了能量的傳播方向[50]，因此可以用快速推進(jìn)法或快速掃描法來求解各向同性介質(zhì)的程函方程.然而，在各向異性介質(zhì)中，這些關(guān)系將不復(fù)存在[50].因此，曲線坐標(biāo)系程函方程的求解是一個(gè)比較復(fù)雜的問題，本文將探討求解曲線坐標(biāo)系程函方程的方法.

2 曲線坐標(biāo)系的程函方程及其嘗試求解

2.1 笛卡爾坐標(biāo)和曲線坐標(biāo)的轉(zhuǎn)換

對(duì)于復(fù)雜地表的介質(zhì)，離散網(wǎng)格邊界需與起伏的地表吻合以避免人為的產(chǎn)生誤差.這種網(wǎng)格被稱作 貼 體 網(wǎng) 格[51－52]，且 廣 泛 應(yīng) 用 于 數(shù) 值 模 擬中[34，42－47，5354].貼體網(wǎng) 格 可 以 通 過 由 計(jì) 算 空 間 到 物理空間的坐標(biāo)變換來獲得[51－52]（圖1）.通過坐標(biāo)變換曲線坐標(biāo)系的q，r映射到物理空間的笛卡爾坐標(biāo)系，這兩個(gè)坐標(biāo)系中垂直軸都取向上方向?yàn)檎?

圖1 計(jì)算域和物理域映射示意圖[46]Fig.1 Mapping between computational and physical space in two dimensions[46]

貼體網(wǎng)格生成之后，笛卡爾坐標(biāo)系中的網(wǎng)格點(diǎn)與曲線坐標(biāo)系中的網(wǎng)格點(diǎn)就一一對(duì)應(yīng)，即

由鏈鎖規(guī)則，有

這里qx表示的意義也類似.這些導(dǎo)數(shù)叫做度量導(dǎo)數(shù)，把式（2a，b）式分別代入（3a，b），經(jīng)過化簡(jiǎn)，可得

J是變換的雅克比，它可以表示為J＝xqzr－xrzq，詳細(xì)的推導(dǎo)過程請(qǐng)參見[46].

2.2 曲線坐標(biāo)系中的程函方程

笛卡爾坐標(biāo)系各向同性介質(zhì)中的程函方程為

這一方程給出了射線經(jīng)過慢度為s（x，z）的介質(zhì)中的點(diǎn) （x，z）時(shí)的走時(shí)T（x，z）.

利用關(guān)系式（2a，b），方程（5）可變換為

利用方程（4）來替換上式中的度量導(dǎo)數(shù)qx，qz，rx，rz，則程函方程（6）可進(jìn)一步變換為

此方程即為曲線坐標(biāo)系中的程函方程，是Hamilton－Jacobi方程中常見的一種.上述方程可記為下列更為緊湊的形式

2.3 嘗試求解曲線坐標(biāo)系程函方程的數(shù)值方法

2.3.1 （不考慮各向異性和各向同性程函方程的差別）直接拓展求解各向同性程函方程的快速推進(jìn)法來求解曲線坐標(biāo)系的程函方程

如果不考慮各向異性的程函方程（7）式與各向同性程函方程的差別，則可以把程函方程（8）式離散為下列的形式

其中A、B、C的含義與（7）式相同，h為曲線坐標(biāo)系中的網(wǎng)格間距，Tqmin＝ min（Ti－1，j，Ti＋1，j），Trmin＝min（Ti，j－1，Ti，j＋1）且

圖2描述了用快速推進(jìn)法求解離散的程函方程的詳細(xì)過程.首先初始化震源點(diǎn)，其走時(shí)在后續(xù)計(jì)算中保持不變.我們把震源點(diǎn)周圍的點(diǎn)的走時(shí)按從小到大的順序存在一個(gè)數(shù)組里，那么它就構(gòu)成了這一時(shí)刻的波前.取出其中的最小值，固定它的值（在后續(xù)迭代中保持不變），所有周圍未被固定的點(diǎn)的值可按下式進(jìn)行迭代：

圖2 快速推進(jìn)法更新步驟[19]Fig.2 Updating procedure for the fast－marching method[19]

其中aa＝A＋B＋C，bb＝2·A·Tqmin＋B·（Tqmin＋Trmin）＋2·C·Trmin，cc＝A·Tqmin2＋B·Tqmin·Trmin＋C·Trmin2－h(huán)2·s2.我 們 把 走 時(shí) 為min（Tqmin，Trmin）的點(diǎn)記為（ii，jj），di，j表示點(diǎn)（i，j）到點(diǎn)（ii，jj）的實(shí)際距離（笛卡爾坐標(biāo)系中兩點(diǎn)間的距離），則，x和z表示笛卡爾坐標(biāo)系中點(diǎn)的坐標(biāo).將每一個(gè)新計(jì)算出來的網(wǎng)格點(diǎn)都放進(jìn)波前數(shù)組里，利用堆選排技術(shù)保存走時(shí)，然后再?gòu)亩褩＠锶〕鲎钚≈?，更新周圍未被固定的點(diǎn)的走時(shí)，依次重復(fù)，直至遍布整個(gè)計(jì)算空間.

2.3.2 Lax－Friedrichs快速掃描法

作為求解Hamilton－Jacobi方程的一種迭代方法，快速掃描法近年來得到了迅速發(fā)展.Tsai等[55]將其應(yīng)用到一類基于均一網(wǎng)格上的Godunov型Hamiltonians的Hamilton－Jacobi方程的求解當(dāng)中.Kao等[56]提出了一種基于Lax－Friedrichs數(shù)值哈密爾頓的（Gauss－Seidel）G－S掃描型算法.它可以求解任意的靜態(tài)Hamilton－Jacobi方程，這種方法簡(jiǎn)單易行、收斂速度快.與基于Godunov數(shù)值哈密爾頓的方法不同的是基于Lax－Friedrichs數(shù)值哈密爾頓的方法需要加一些邊界條件.

2.3.2.1 Lax－Friedrichs哈密爾頓

對(duì)曲 線 坐 標(biāo) 系 中 的 矩 形 區(qū) 域 [qmin，qmax]×[rmin，rmax]進(jìn)行均勻剖分，網(wǎng)格點(diǎn)記為 （qi，rj），i＝0，1，…，m1，m1＋1，j＝0，1，…，m2，m2＋1，這里qi＝ （i－1）hq＋qmin，rj＝ （j－1）hr＋rmin，hq＝

曲線坐標(biāo)系中的程函方程（8）是一種特殊的靜態(tài)Hamilton－Jacobi方程，它可以寫做如下的形式：

其中

為哈密爾頓算子，A，B，C的意義與前述相同；Ω為計(jì)算區(qū)域；Γ 為震源點(diǎn)；s（q，r）是慢度（速度的倒數(shù)）；g（q，r）為震源點(diǎn)的走時(shí)值.

用Lax－Friedrichs格式來離散哈密爾頓算子H（Tq，Tr）：

（1）從左下角向右上角，即i＝1∶m1，j＝1∶m2；

（2）從右下角向左上角，即i＝m1∶1，j＝1∶m2；

（3）從左上角向右下角，即i＝1∶m1，j＝m2∶1；

（4）從右上角向左下角，即i＝m1∶1，j＝m2∶1.

以從左下角向右上角方向掃描，即i＝1∶m1，j＝1∶m2為例，有

要指出的是以上的更新公式不涉及任何非線性求逆，因此本文的求解方法很容易實(shí)現(xiàn).

2.3.2.2 計(jì)算邊界條件

由于利用Lax－Friedrichs數(shù)值哈密爾頓計(jì)算一個(gè)網(wǎng)格點(diǎn)的解時(shí)需要它周圍所有相鄰點(diǎn)的信息，因此必須仔細(xì)地確定那些超出計(jì)算區(qū)域的網(wǎng)格點(diǎn)的值.為此，采用Kao等[56]提出的方法處理計(jì)算邊界，他們結(jié)合外推、最大和最小策略計(jì)算那些超出計(jì)算區(qū)域的網(wǎng)格點(diǎn)的值.二維問題的邊界條件為

最后將Gauss－Seidel型Lax－Friedrichs掃描算法作個(gè)總結(jié).它主要包含初始化、交替掃描和實(shí)施計(jì)算邊界條件這三個(gè)步驟.具體如下：

1.初始化.在邊界Γ上的網(wǎng)格點(diǎn)或鄰近的網(wǎng)格點(diǎn)，賦精確值或插值，并且每次迭代過程中這些值固定不變.其它網(wǎng)格點(diǎn)上，賦較大的值M，并且這些值在迭代過程中被更新，這里M應(yīng)該大于真解的最大值.

2.交替掃描.在第n＋1次迭代時(shí)，除了那些賦了精確值或插值的網(wǎng)格點(diǎn)，按照（11）計(jì)算在所有網(wǎng)格點(diǎn) （qi，rj），1≤i≤m1，1≤j≤m2的值，并且只有當(dāng)它小于以前的值時(shí)才更新值得注意的是，這個(gè)過程需要按交替掃描方向進(jìn)行，即在二維情形需要四個(gè)不同的掃描方向.一般來說，d維情形需要2d次交替掃描.

3.實(shí)施計(jì)算邊界條件.每次掃描之后，按照（16）實(shí)施計(jì)算邊界條件，不同的邊界只需作一些細(xì)微的修改.

4.收斂準(zhǔn)則.如果 ‖Tn＋1－Tn‖L1≤δ，算法收斂并終止，這里δ是一個(gè)給定的收斂閥值.

3 計(jì)算實(shí)例

為了研究本文給出的兩種數(shù)值方法求解曲線坐標(biāo)系程函方程的正確性和適用性，我們首先考慮一個(gè)水平地表模型，然后對(duì)一組地表地形為余弦函數(shù)的起伏地表模型進(jìn)行計(jì)算分析.

為了便于與解析解對(duì)比，我們假設(shè)模型中介質(zhì)都是均勻的，速度均為2000m/s.

3.1 水平地表模型

模型大小為 [－1km，1km ]2，震源坐標(biāo)約為（－0.25km，0.50km）.我們把兩種數(shù)值方法計(jì)算的數(shù)值解與解析解進(jìn)行對(duì)比，結(jié)果見走時(shí)等值線圖3，圖4為兩種方法計(jì)算的結(jié)果的絕對(duì)誤差分布圖.由圖可知，直接拓展求解各向同性程函方程的快速推進(jìn)法與用Lax－Friedrichs快速掃描法所求得的數(shù)值解都與解析解吻合的較好，且前者的精度更高.這也可以從過炮點(diǎn)且分別與縱軸和橫軸平行的兩條測(cè)線上的走時(shí)值中（圖5）觀察到.這是因?yàn)?，?duì)于水平地表的模型，此時(shí)各向異性的程函方程（7）就退化為各向同性的程函方程（5），而（16）也就退化為求解各向同性程函方程的迎風(fēng)差分格式，因此得到更為準(zhǔn)確的解；Lax－Friedrichs快速掃描法不是迎風(fēng)格式的，但是它具有很強(qiáng)的健壯性，對(duì)任意形式的Hamiltonians都可以求解，這可以通過后面的實(shí)例得到驗(yàn)證.

圖3 水平地表模型當(dāng)震源坐標(biāo)為 （－0.25km，0.50km）時(shí)的初至波走時(shí)等值線圖Fig.3 Traveltime contours（in seconds）for the flat surface model with the source located at（－0.25km，0.50km）

圖4 直接拓展求解各向同性程函方程的快速推進(jìn)法（a）和Lax－Friedrichs快速掃描法（b）所得的數(shù)值解的絕對(duì)誤差分布Fig.4 Absolute error maps for the numerical solutions obtained by a simple modification of the fast marching method for the isotropic eikonal equation（a）and Lax－Friedrichs fast sweeping method（b），respectively.The units in the error maps are second（s）.

圖5 用兩種方法計(jì)算的過炮點(diǎn)且分別與縱軸（a）和橫軸（b）平行的兩條測(cè)線上的走時(shí)值及其相應(yīng)的解析解對(duì)比.Fig.5 Comparisons of travel times at the two lines which cross the source and parallel z－and x－coordinates，respectively.

3.2 起伏地表模型

3.2.1 余弦形起伏地表模型

選取一組起伏地表模型進(jìn)行研究.模型表面是由一個(gè)山丘和兩個(gè)凹陷組合成的起伏型地表.它的高程可以用下面的余弦函數(shù)來表示

其中，A為余弦函數(shù)的振幅，A取值不同表示地表地形起伏的劇烈程度不同，分別取A為0.1、0.2來進(jìn)行研究.一般的，復(fù)雜地表下的走時(shí)很難用解析的方法進(jìn)行刻畫，但在地表地形可以用解析式表示的情況下，其解析解還是可以求解的.對(duì)于這兩個(gè)起伏地表的均勻介質(zhì)模型，計(jì)算出其相應(yīng)的解析解（見附錄）并和本文計(jì)算的數(shù)值解對(duì)比來考查本文數(shù)值方法求解曲線坐標(biāo)系程函方程的正確性.

圖6為當(dāng)A取0.2時(shí)模型的結(jié)構(gòu)化網(wǎng)格剖分圖.當(dāng)A＝0.1時(shí)，震源坐標(biāo)分別為（－0.25km，0.50km），（－0.25km，0.88km），把兩種方法計(jì)算的數(shù)值解與解析解做一對(duì)比，結(jié)果見圖7.從圖中可以觀察到，Lax－Friedrichs快速掃描法計(jì)算的結(jié)果與解析解吻合的較好，而直接拓展求解各向同性程函方程的快速推進(jìn)法所得的數(shù)值解與解析解出現(xiàn)了較大偏差，這是因?yàn)榇藭r(shí)曲線坐標(biāo)系中的程函方程（7）為各向異性的，且隨著地表附近貼體網(wǎng)格的扭曲程度逐漸加重方程（7）的各向異性增強(qiáng)[57]（可以從曲線坐標(biāo)系的等值線圖8中觀察到），簡(jiǎn)單的拓展求解各向同性程函方程的快速推進(jìn)法所得的差分格式（12）并不是方程（7）的正確的迎風(fēng)格式的離散形式.而Lax－Friedrichs快速掃描法對(duì)各向異性的程函方程也是適用的.從地表的走時(shí)曲線（圖9）也可以得到類似的結(jié)論.

圖6 地表地形為z（x）＝1.0＋0.2cos（1.5πx）km的模型結(jié)構(gòu)化網(wǎng)格剖分Fig.6 The boundary conforming grids for the model whose surface is described by the function z（x）＝1.0＋0.2cos（1.5πx）km

圖7 振幅為0.1km的余弦型起伏地表模型震源坐標(biāo)分別約為（a）（－0.25km，0.50km），（b）（－0.25km，0.88km）時(shí)的初至波走時(shí)等值線圖黑色實(shí)線表示解析解；紅色虛線表示直接拓展求解各向同性程函方程的快速推進(jìn)法所得的數(shù)值解；綠色虛線表示用Lax－Friedrichs快速掃描法求得的數(shù)值解，走時(shí)等值線的單位為秒.Fig.7 Traveltime contours（in seconds）for the model whose surface is a cosine curve－like feature The source is at（a）（－0.25km，0.50km）and（b）（－0.25km，0.88km），respectively.The black solid line denotes the analytical solution，while the red and green dashed lines denote numerical solutions obtained by a simple modification of the fast marching method for the isotropic eikonal equation and Lax－Friedrichs fast sweeping method，respectively.

圖8 振幅為0.1km的余弦型起伏地表模型震源坐標(biāo)分別約為（a）（－0.25km，0.50km），（b）（－0.25km，0.88km）時(shí)在曲線坐標(biāo)系中的初至波走時(shí)等值線圖Fig.8 Traveltime contours（in seconds）in the curvilinear coordinate for the model whose surface is a cosine curve－like feature

圖9 振幅為0.1（km）的余弦型起伏地表模型地表上走時(shí)的數(shù)值解與解析解對(duì)比震源坐標(biāo)分別約為（a）（－0.25km，0.50km），（b）（－0.25km，0.88km）.黑色實(shí)線表示解析解；紅色虛線表示直接拓展求解各向同性程函方程的快速推進(jìn)法所得的數(shù)值解；綠色虛線表示用Lax－Friedrichs快速掃描法求得的數(shù)值解.Fig.9 Comparison of travel times of the analytical and numerical solutions at the irregular surface for the model whose topography amplitude is 0.1km with the source located at（a）（－0.25km，0.50km）and（b）（－0.25km，0.88km）.The black solid line denotes the analytical solution，while the red and green dashed lines denote numerical solutions obtained by a simple modification of the fast marching method for the isotropic eikonal equation and Lax－Friedrichs fast sweeping method，respectively.

當(dāng)A取0.2時(shí)，地表地形起伏更加劇烈，對(duì)于地表附近的點(diǎn)，方程（7）的各向異性較A＝0.1時(shí)強(qiáng).此時(shí)離散格式（12）已無法求解（根式下出現(xiàn)負(fù)數(shù)），而用而Lax－Friedrichs快速掃描法求解依然可以得到正確穩(wěn)定的解，分別見圖10、11.這進(jìn)一步證實(shí)了我們前面的結(jié)論.

圖10 振幅為0.2（km）的余弦型起伏地表模型當(dāng)震源分別位于（a）（－0.2km，1.1km），（b）（－0.2km，－0.1km）時(shí)走時(shí)等值線圖紅線表示解析解，藍(lán)表示數(shù)值解（Lax－Friedrich快速掃描法求得的結(jié)果），走時(shí)等值線的單位為s.Fig.10 Traveltime contours（in seconds）for the model whose topography amplitude is 0.2km with the source located at（a）（－0.2km，1.0km）and（b）（－0.2km，－0.1km）The red and blue lines denote the analytical and numerical solutions（Computed by the Lax－Friedrichs fast sweeping method），respectively.

圖11 振幅為0.2km的余弦型起伏地表模型地表上走時(shí)的數(shù)值解（Lax－Friedrich快速掃描法求得的結(jié)果）與解析解對(duì)比震源坐標(biāo)分別為 （a）（－0.2km，1.1km），（b）（－0.2km，－0.1km）.黑色虛線和實(shí)線分別表示解析解和數(shù)值解.Fig.11 Comparison of travel times of the analytical and numerical（Computed by the Lax－Friedrichs fast sweeping method）solutions at the irregular surface for the model whose topography amplitude is 0.2km with the source located at（a）（－0.2km，1.1km），（b）（－0.2km，－0.1km），respectively.The dotted and solid lines show，respectively，the analytical and numerical solutions.

3.2.2 劇烈起伏地表下的Marmousi模型

為了進(jìn)一步考查本文的方法，我們選取了一個(gè)地表劇烈起伏的模型，Sun等[49]在模擬起伏地表的走時(shí)也選用了類似的模型.模型大小為1.6km×1.32km，地表在1.1～1.32km之間起伏，表示的地形為兩個(gè)山丘和兩個(gè)凹陷的交替組合.用1280×960的網(wǎng)格來離散此模型，剖分的貼體網(wǎng)格見圖12.震源點(diǎn)的坐標(biāo)為（0.8km，1.1km），把以震源為圓心半徑為0.1km的區(qū)域內(nèi)的點(diǎn)的初值賦為精確值.圖13為劇烈起伏地表下的Marmousi模型的初至波走時(shí)圖.由圖可見，Lax－Friedrichs快速掃描法對(duì)地表地形起伏劇烈，速度縱橫向變化劇烈的介質(zhì)依然具有很好的適用性和穩(wěn)定性.對(duì)于此模型，直接用快速推進(jìn)法依然無法求解.

圖12 劇烈起伏地表模型結(jié)構(gòu)化網(wǎng)格剖分為清楚起見，網(wǎng)格密度減少為原來的1/8.倒立的黑色三角形表示在起伏的地表上布設(shè)的檢波器，水平距離從0到1.6km.Fig.12 The boundary conforming grids in the severe undulated surface model For clarity，the grids are displayed with a reducing density factor of 8.The inverted black triangles denote the receivers located at the irregular surface，with horizontal distances between 0and 1.6km.

圖13 劇烈起伏地表下的Marmousi模型中初至波走時(shí)等值線.震源坐標(biāo)為（0.8km，1.1km）.Fig.13 Travel time contours（in seconds）superimposed on the revised Marmousi model，the source is located at（0.8km，1.1km）.

4 結(jié) 論

復(fù)雜地表?xiàng)l件下的地震波走時(shí)計(jì)算對(duì)于實(shí)現(xiàn)復(fù)雜地表?xiàng)l件下的靜校正、偏移成像、走時(shí)反演、層析成像、地震定位等具有重要意義.借助于貼體網(wǎng)格，把笛卡爾坐標(biāo)系中的程函方程變換到曲線坐標(biāo)系中，然后，嘗試采用直接拓展各向同性程函方程的快速推進(jìn)法和Lax－Friedrichs快速掃描算法求解曲線坐標(biāo)系的程函方程.數(shù)值實(shí)例表明，此時(shí)在曲線坐標(biāo)系中程函方程已變?yōu)楦飨虍愋缘某毯匠?，因此通過直接拓展各向同性程函方程的快速推進(jìn)法來求解該方程是錯(cuò)誤的，而用Lax－Friedrichs快速掃描法求解總能得到穩(wěn)定正確的解，這對(duì)于起伏地表的走時(shí)計(jì)算具有重要意義.

致 謝 本研究工作得到了中國(guó)科學(xué)院地質(zhì)與地球物理研究所張中杰研究員的悉心指導(dǎo)和幫助.三位審稿專家對(duì)本研究工作提供了許多建設(shè)性的意見.筆者在此表示衷心的感謝.

附錄：復(fù)雜地表下均勻介質(zhì)模型中的精確走時(shí)計(jì)算策略

一般的，復(fù)雜地表下的走時(shí)很難用解析的方法進(jìn)行刻畫，但在地表地形可以用解析式表示的情況下，其解析解還是可以求解的.在此，我們就來探討本文模型的解析解的求解.

以地表地形為振幅為0.2km的余弦曲線模型為例來介紹其精確走時(shí)的求解.首先，討論震源在點(diǎn)（－0.2km，1.1km）的情形，對(duì)于其它情形后續(xù)在討論.求解的難點(diǎn)是如何確定從震源到區(qū)域內(nèi)其它點(diǎn)的射線的最短路徑.根據(jù)從震源‘S’到模型內(nèi)其它點(diǎn)的射線路徑的形態(tài)，求解區(qū)域可以被剖分為三部分，如圖A1a顯示的兩個(gè)灰色區(qū)域和白色區(qū)域的部分.對(duì)于白色區(qū)域內(nèi)的點(diǎn)，從震源到該點(diǎn)的最短射線路徑均為直線，因此其走時(shí)很容易求得.對(duì)于兩個(gè)灰色區(qū)域，我們以左上角的灰色區(qū)域?yàn)槔齺斫榻B，右上角的區(qū)域可類似計(jì)算.對(duì)于灰色區(qū)域內(nèi)的一點(diǎn)B（弧下方，線段DE上方），過該點(diǎn)做一條直線與相切，切點(diǎn)為C（在上），那么此時(shí)從震源S到點(diǎn)B的最短射線路徑為SD＋＋CB（易被證明，在此從略），向量SD和CB的長(zhǎng)度可以通過兩點(diǎn)間距離公式求得，而的長(zhǎng)度可以通過弧長(zhǎng)公式計(jì)算得出.進(jìn)而可以求出其相應(yīng)的走時(shí).同理，可以計(jì)算出此區(qū)域內(nèi)其它點(diǎn)的走時(shí).

當(dāng)震源位于點(diǎn) （－0.2km，－0.1km）（圖A1（b））上時(shí)，只有右上角一小部分區(qū)域內(nèi)的點(diǎn)需要按照上述的方法特別處理.

綜上所述，我們就可以求得該模型中的初至波精確走時(shí).對(duì)于其它模型的情況可以類似求解，這里從略.

圖A1 本文中的起伏地表模型走時(shí)計(jì)算策略Fig.A1 The strategy for computing analytical solutions of travel times for the irregular surface models of this paper

（References）

[1] Aki K，Richards P G.Quantitative seismology.2002：Univ Science Books.

[2] Hole J.Nonlinear high－resolution three－dimensional seismic travel time tomography.Journal of Geophysical Research，1992，97（B5）：6553－6562.

[3] Chen Y，Badal J，Hu J.Love and Rayleigh wave tomography of the Qinghai－Tibet Plateau and surrounding areas.Pure and Applied Geophysics，2010，167（10）：1－33.

[4] Zelt C，Smith R.Seismic travel time inversion for 2－D crustal velocity structure.Geophysical Journal International，1992，108（1）：16－34.

[5] Zhang Z J，Wang G J，Teng J W，et al.CDP mapping to obtain the fine structure of the crust and upper mantle from seismic sounding data：an example for the southeastern China.Physics of the Earth and Planetary Interiors，2000，122（1－2）：133－146.

[6] 張中杰，秦義龍，陳赟等.由寬角反射地震資料重建殼幔反射結(jié)構(gòu)的相似性剖面.地球物理學(xué)報(bào)，2004，47（003）：469－474.Zhang Z J，Qin Y L，Chen Y，et al.Reconstruction of the semblance section for the crust and mantle reflection structure by wide－angle reflection seismic data.Chinese J.Geophys.（in Chinese），2004，47（3）：469－474

[7] 朱金明，王麗燕.地震波走時(shí)的有限差分法計(jì)算.地球物理學(xué)報(bào)，1992，35（1）：86－92.Zhu J M，Wang L Y.Finite difference calculation of seismic travel times.Chinese J.Geophys.（in Chinese），1992，35（1）：86－92.

[8] Cerveny V，Molotkov I A，and Psencik I.Ray methods in seismology.1977，Univ.of Karlova Press.

[9] Zhang Z，Badal J，Li Y，et al.Crust－upper mantle seismic velocity structure across Southeastern China.Tectonophysics，2005，395（1－2）：137－157.

[10] Zhang Z，Li Y，Lu D，et al.Velocity and anisotropy structure of the crust in the Dabieshan orogenic belt from wide－angle seismic data.Physics of the Earth and Planetary Interiors，2000，122（1－2）：115－131.

[11] Zhang Z，Wang Y，Chen Y，et al.Crustal structure across Longmenshan fault belt from passive source seismic profiling.Geophys.Res.Lett，2009，36.

[12] Xu T，Xu G，Gao E，et al.Block modeling and segmentally iterative ray tracing in complex 3 Dmedia.Geophysics，2006，71：T41－T51.

[13] Xu T，Zhang Z，Gao E，et al.Segmentally Iterative Ray Tracing in Complex 2D and 3DHeterogeneous Block Models.Bulletin of the Seismological Society of America，100（2）：841－850.

[14] Vidale J E.Finite－difference calculation of travel times.Bulletin of the Seismological Society of America，1988，78（6）：2062－2076.

[15] Vidale J E.Finite－difference calculation of traveltimes in three dimensions.Geophysics，1990，55（5）：521－526.

[16] Qin F，Luo Y，Olsen K B，et al.Finite－difference solution of the eikonal equation along expanding wavefronts.Geophysics，1992，57：478－487.

[17] Van Trier J， Symes W W. Upwind finite－difference calculation of traveltimes.Geophysics，1991，56（6）：812－821.

[18] Sethian J.Evolution，implementation，and application of level set and fast marching methods for advancing fronts.Journal of Computational Physics，2001，169（2）：503－555.

[19] Sethian J A.Fast marching methods.SIAM review，1999，41（2）：199－235.

[20] Sethian J A，Popovici A M.Three dimensional traveltimes computation using the Fast Marching Method.Geophysics，1999，64（2）：516－523.

[21] Sethian J A. A fast marching level set method for monotonically advancing fronts.Proceedings of the National Academy of Sciences of the United States of America，1996，93（4）：1591－1595.

[22] Alkhalifah，T. ＆ Fomel，S.，Implementing the fast marching eikonal solver：spherical versus Cartesian coordinates.Geophys.Prospect.，2001，49：165－178.

[23] Fomel S.A variational formulation of the fast marching eikonal solver.SEP－95：Stanford Exploration Project，1997：127－147.

[24] Kim S.3－D eikonal solvers： First－arrival traveltimes.Geophysics，2002，67（4）：1225－1231.

[25] Rawlinson N，Sambridge M. Wave front evolution in strongly heterogeneous layered media using the fast marching method.Geophysical Journal International，2004，156（3）：631－647.

[26] Rawlinson N， Sambridge M. Multiple reflection and transmission phases in complex layered media using a multistage fast marching method.Geophysics，2004，69（5）：1338.

[27] Rawlinson N，Sambridge M.The fast marching method：an effective tool for tomographic imaging and tracking multiple phases in complex layered media.Exploration Geophysics，2005，36（4）：341－350.

[28] 孫章慶，孫建國(guó)，韓復(fù)興.復(fù)雜地表?xiàng)l件下基于線性插值和窄帶技術(shù)的地震波走時(shí)計(jì)算.地球物理學(xué)報(bào)，2009，52（11）：2846－2853.Sun Z Q，Sun J G，Han F X.Travehimes computation using linear interpoiation and narrow band technique under complex topographical conditions.Chinese J.Geophys.（in Chinese），2009，52（11）：2846－2853

[29] De Kool M，Rawlinson N，Sambridge M.A practical gridbased method for tracking multiple refraction and reflection phases in three－dimensional heterogeneous media.Geophysical Journal International，2006，167（1）：253－270.

[30] Zhao H.A fast sweeping method for eikonal equations.Mathematics of Computation，2005，74（250）：603－628.

[31] Zhang Y T，Zhao H K，Qian J.High order fast sweeping methods for static Hamilton－Jacobi equations.Journal of Scientific Computing，2006，29（1）：25－56.

[32] 蘭海強(qiáng)，張智，徐濤等.地震波走時(shí)場(chǎng)模擬的快速推進(jìn)法和快速掃描法比較研究.地球物理學(xué)進(jìn)展，2012，已接受.Lan H Q，Zhang Z，Xu T，et.al.A comparative study on the fast marching and fast sweeping methods in the calculation of first－arrival traveltime field. Progress in Geophys.（in Chinese），2012，accepted.

[33] Zhang Z J，Yang L Q，Teng J W，Badal J.An overview of the earth crust under China.Earth－Science Reviews，104（1－3）：143－166DOI：10.1016/j.earscirev.2010.10.003，2011.

[34] Lan H，Zhang Z.Three－dimensional wave－field simulation in heterogeneous transversely isotropic medium with irregular free surface.Bulletin of the Seismological Society of America，2011101（3）：1354－1370.

[35] 閻世信，劉懷山，姚雪根.山地地球物理勘探技術(shù).北京：石油工業(yè)出版社，2000.Yan S X，Liu H S，Yao X G.Geophysical Exploration Technology in the Mountainous Area（in Chinese）.Beijing：Petroleum Industry Press，2000

[36] Kao C Y，Osher S，Qian J.Legendre－transform－based fast sweeping methods for static Hamilton－Jacobi equations on triangulated meshes.Journal of Computational Physics，2008，227（24）：10209－10225.

[37] Kimmel R，Sethian J A.Computing geodesic paths on manifolds.Proceedings of the National Academy of Sciences of the United States of America，1998，95（15）：8431－7435.

[38] Lelièvre P G，F(xiàn)arquharson C G，Hurich C A.Computing first－arrival seismic traveltimes on unstructured 3－D tetrahedral grids using the Fast Marching Method.Geophysical Journal International，184（2）：885－896.

[39] Qian J，Zhang Y T，Zhao H K.A fast sweeping method for static convex Hamilton－Jacobi equations. Journal of Scientific Computing，2007，31（1）：237－271.

[40] Qian J，Zhang Y T，Zhao H K.Fast sweeping methods for Eikonal equations on triangular meshes.SIAM Journal on Numerical Analysis，2008，45（1）：83－107.

[41] Sethian J A，Vladimirsky A.Fast methods for the Eikonal and related Hamilton－Jacobi equations on unstructured meshes.Proceedings of the National Academy of Sciences of the United States of America，2000，97（11）：5699.

[42] Appelo D，Petersson N.A stable finite difference method for the elastic wave equation on complex geometries with free surfaces.Communications in Computational Physics，2009，5（1）：84－107.

[43] 孫章慶，孫建國(guó)，張東良.二維起伏地表?xiàng)l件下坐標(biāo)變換法直流電場(chǎng)數(shù)值模擬.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2009，39（3）：528－534.Sun Z Q，Sun J G，Zhang D L.2DDC electric field numerical modeling including surface topography using coordinate transformation method.Journal of Jilin University：Earth Science Edition，2009，39（3）：528－534.

[44] Lan H，Zhang Z.2Dfinite－difference seismic modeling of an open fluid－filled fracture with irregular free surface，in review

[45] Tessmer E，Kosloff D.3－D elastic modeling with surface topography by a Chebychev spectral method.Geophysics，1994，59：464.

[46] 蘭海強(qiáng)，劉佳，白志明.VTI介質(zhì)起伏地表地震波場(chǎng)模擬.地球物理學(xué)報(bào)，2011，54（8）：2072－2084.Lan H Q，Liu J，Bai Z M.Wave－field simulation in VTI media with irregular free surface.Chinese J.Geophys.（in Chinese），2011，54（8）：2072－2084.

[47] 孫建國(guó).復(fù)雜地表?xiàng)l件下地球物理場(chǎng)數(shù)值模擬方法評(píng)述.世界地質(zhì)，2007，26（003）：345－362.Sun J G.Methods for numerical modeling of geophysical fields under complex topographica I conditions：a critical review.Global Geology （in Chinese），2007，26（003）：345－362.

[48] Lan H，Zhang Z.Comparative study of the free－surface boundary condition in two－dimensional finite－difference elastic wave field simulation.Journal of Geophysics and Engineering，2011，8：275－286.

[49] Sun J，Sun Z，Han F.A finite difference scheme for solving the eikonal equation including surface topography.Geophysics，2011，76（4）：T53－T63.

[50] Zhang Z J，Teng J W，Wan Z C.Establishment of parametric equations of seismic－wavefronts in 2－D anisotropic media.Chinese Science Bulletin，1994，39（22）：1890－1894.

[51] Hvid S L.Three dimensional algebraic grid generation[Ph.D.thesis]，Technical University of Denmark，1994.

[52] Thompson J，Warsi Z，Mastin C.Numerical grid generation：foundations and applications.1985：North－h(huán)olland Amsterdam.

[53] Hestholm S，Ruud B.3－D finite－difference elastic wave modeling including surface topography.Geophysics，1998，63：613－622.

[54] Tessmer E，Kosloff D，Behle A.Elastic wave propagation simulation in the presence of surface topography.Geophys.J.Int，1992，108（2）：621－632.

[55] Tsai Y H R，Cheng L T，Osher S，et al.Fast sweeping algorithms for a class of Hamilton－Jacobi equations.SIAM Journal on Numerical Analysis，2004：673－694.

[56] Kao C Y，Osher S，Qian J.Lax－Friedrichs sweeping scheme for static Hamilton－Jacobi equations.Journal of Computational Physics，2004，196（1）：367－391.

[57] 蘭海強(qiáng)，徐濤，白志明.貼體網(wǎng)格各向異性對(duì)坐標(biāo)變換法求解起伏地表下地震初至波走時(shí)的影響研究.地球物理學(xué)報(bào)，已接受.Lan H Q，Xu T，Bai Z M.Study of the effects due to the anisotropic stretching of the surface－fitting grid on the traveltime computation for irregular surface by the coordinate transforming method.Chinese J.Geophys.（in Chinese），in press.

Study on the numerical solutions of the eikonal equation in curvilinear coordinate system

LIU Yi－Feng1，2，3，LAN Hai－Qiang1，3
1 State Key Laboratory of Lithospheric Evolution，Institute of Geology and Geophysics，Chinese Academy of Sciences，Beijing100029，China
2 Beijing Research Center of CNOOC，Beijing100027，China
3 Graduate School of Chinese Academy of Sciences，Beijing100049，China

The classical eikonal equation is commonly used in Cartesian coordinate system for problems that involve static correction，prestack migration，earthquake location and seismic tomography，but is less effective for calculating travel times in an earth model that has an irregular surface.We transform the eikonal equation to a curvilinear coordinate system that conforms with the boundaries of the domain.The eikonal equation in the curvilinear coordinate system displays the mathematical form of an anisotropic eikonal equations（even though the medium is isotropic in the Cartesian coordinate system）.Then，we try to solve the anisotropic eikonal equations with two numerical methods，a simple modification of the fast marching methodand Lax－Friedrichs fast sweeping method，respectively.Through a number of different numerical experiments we find：one may compute wrong solutions by extending fast marching methods designed for isotropic eikonal equations to anisotropic eikonal equations without taking into account the differences between them；the Lax－Friedrichs fast sweeping method is stable and valid for the eikonal equation in the curvilinear coordinate system，which has a significant meaning in the direction of solution of the eikonal equation in the curvilinear coordinate system.

Eikonal equation，Topography，Anisotropy，F(xiàn)ast marching method，Lax－Friedrichs fast sweeping method

10.6038/j.issn.0001－5733.2012.06.023

P631

2012－02－05，2012－03－29收修定稿

國(guó)家自然科學(xué)基金（41074033，40874041，41021063）資助.

劉一峰，男，1964年生，博士，主要從事油氣勘探開發(fā)方法方面的研究.E－mail：liuyifeng＠cnooc.com.cn

＊通訊作者 蘭海強(qiáng)，男，現(xiàn)為中國(guó)科學(xué)院地質(zhì)與地球物理研究所博士研究生，主要從事復(fù)雜介質(zhì)中地震波傳播與成像及地震資料處理的理論與方法等研究.E－mail：lanhq＠m(xù)ail.iggcas.ac.cn

劉一峰，蘭海強(qiáng).曲線坐標(biāo)系程函方程的求解方法研究.地球物理學(xué)報(bào)，2012，55（6）：2014－2026，

10.6038/j.issn.0001－5733.2012.06.023.

Liu Y F，Lan H Q.Study on the numerical solutions of the eikonal equation in curvilinear coordinate system.Chinese J.Geophys.（in Chinese），2012，55（6）：2014－2026，doi：10.6038/j.issn.0001－5733.2012.06.023.

（本文編輯 劉少華）