彈性波有限差分模擬中自由表面的自適應(yīng)表達(dá)

曹 健,陳景波

(1.中國科學(xué)院地質(zhì)與地球物理研究所,中國科學(xué)院油氣資源研究重點實驗室,北京100029;2.中國科學(xué)院地球科學(xué)研究院,北京100029;3.中國科學(xué)院大學(xué),北京100049)

固體地球表面通常也被稱為自由表面,是一個強(qiáng)阻抗突變界面,它的存在使得其附近的地震波波場變得更為復(fù)雜,例如會產(chǎn)生強(qiáng)振幅面波(瑞利波和洛夫波),在地表處產(chǎn)生反射縱波和轉(zhuǎn)換橫波,以及由于縱橫波干涉造成的場地效應(yīng)等[1]。由于接收到的地震數(shù)據(jù)大多來自地表或近地表,這些復(fù)雜的地震波不可避免地被記錄在地震記錄中,因此需要在地震波正演數(shù)值模擬時考慮自由表面的影響。而自由表面邊界條件則是對這一強(qiáng)阻抗突變界面的數(shù)學(xué)表達(dá)。

對于有限元類方法,例如譜元法[2-3]和間斷有限元法[4],自由表面邊界條件是自然滿足的,但這類方法實現(xiàn)起來相對復(fù)雜,并且對于復(fù)雜模型的計算網(wǎng)格設(shè)計較為困難[5]。相比較而言,有限差分法具有算法簡單、可靈活用于各種復(fù)雜模型的數(shù)值模擬、精確和計算高效等優(yōu)點,目前已被廣泛應(yīng)用于地震波的正演數(shù)值模擬。然而,有限差分法中,自由表面的實現(xiàn)需要進(jìn)行特別處理。為了獲得精確并且穩(wěn)定的自由表面邊界實現(xiàn)方式,在過去的幾十年中研究人員圍繞自由表面邊界條件實現(xiàn)方式開展了大量研究。簡單地講,自由表面邊界實現(xiàn)方式可分為兩大類。第一類是基于彈性波方程的二階位移表達(dá)形式,這種情況下,自由表面邊界條件通過位移的空間導(dǎo)數(shù)來表示,其形式屬于偏微分方程中的Neumann邊界條件。因此,可采用不同種類的差分格式來近似自由表面邊界處的空間導(dǎo)數(shù),例如,中心差分顯式近似格式[6]、單邊顯式近似格式[7-8]、中心差分隱式近似格式[9]以及復(fù)合近似等格式[10-12]。這類自由表面邊界處理方式的主要缺陷在于對高泊松比介質(zhì)進(jìn)行正演模擬時會產(chǎn)生計算不穩(wěn)定現(xiàn)象,并且只能達(dá)到二階的計算精度。第二類實現(xiàn)方式是基于彈性波方程的一階速度-應(yīng)力表達(dá)形式。與第一類方式相比,自由表面邊界條件可以直接采用應(yīng)力為0的形式進(jìn)行表達(dá),即偏微分方向中的Dirichlet邊界條件,這一形式在數(shù)學(xué)上要相對容易處理一些。其中應(yīng)用比較廣泛的有真空格式[13-14]、應(yīng)力鏡像法[15-17]和一些介質(zhì)參數(shù)修正方法[18-19]。它們的主要優(yōu)勢是對于任意的泊松比值介質(zhì)都能提供穩(wěn)定的正演模擬結(jié)果,但存在著某些泊松比值情況下計算結(jié)果不夠準(zhǔn)確,需要進(jìn)行網(wǎng)格加密處理以提高計算精度。

為了解決上述提及的直接針對自由表面邊界條件數(shù)值近似處理帶來的穩(wěn)定性和精度方面的問題,本文采取了CAO等[20]提出的重新建立自由表面表達(dá)方式的想法,從平均介質(zhì)理論[21]出發(fā),并結(jié)合極限思想,從物理角度給出了自由表面的另一種新的自適應(yīng)數(shù)學(xué)表達(dá)方式,從而更易于在有限差分框架下實現(xiàn)。首先基于平均介質(zhì)理論,從物理角度給出自由表面處的密度和本構(gòu)關(guān)系應(yīng)滿足的自適應(yīng)表達(dá)方式;然后具體介紹了該自由表面自適應(yīng)表達(dá)方式在速度-應(yīng)力交錯網(wǎng)格有限差分中的實現(xiàn)方式;最后通過一系列數(shù)值實驗來對比和討論該自由界面自適應(yīng)表達(dá)方式在不同泊松比情況下的數(shù)值模擬效果。

1 自由表面的自適應(yīng)表達(dá)

從物理的角度而言,自由表面是地表以上的空氣介質(zhì)和地下彈性介質(zhì)之間的一個接觸界面,即流體/固體邊界。由于地表以上的空氣介質(zhì),其密度和彈性模量要遠(yuǎn)遠(yuǎn)小于地下巖石介質(zhì)的密度和彈性模量,因此完全可以忽略。根據(jù)這一真空近似[22],自由表面還可以被視為真空/固體邊界。MOCZO等[21]從平均介質(zhì)的角度出發(fā),通過彈性模量調(diào)和平均和密度參數(shù)代數(shù)平均的方式來表示各向同性彈性介質(zhì)內(nèi)部的介質(zhì)參數(shù)不連續(xù)界面,即固體/固體邊界。我們將上述平均介質(zhì)思想引入到自由表面這一強(qiáng)不連續(xù)界面的表示上來,從而推導(dǎo)出自由表面的另一種自適應(yīng)數(shù)學(xué)表達(dá)方式。

假設(shè)兩個具有不同介質(zhì)參數(shù)的各向同性彈性半空間,其分界面位于z=0處,上層介質(zhì)的密度和拉梅常數(shù)分別為ρ-,λ-和μ-,而下層介質(zhì)的密度和拉梅常數(shù)分別為ρ+,λ+和μ+(圖1)。當(dāng)上層介質(zhì)為真空時,即ρ-=0,λ-=0和μ-=0,分界面z=0自然就變成了自由表面。根據(jù)MOCZO等[21]中提出的介質(zhì)參數(shù)平均方式,分界面z=0處的彈性參數(shù)可以表示為:

(1)

式中:κ是體積模量,κ=λ+2μ/3。如果將這一參數(shù)平均方式直接應(yīng)用到自由表面處的介質(zhì)參數(shù)計算上來,即讓ρ-=λ-=μ-=0,可以直接得到:

(2)

這里的密度減半可以被理解為自由表面處的質(zhì)點與內(nèi)部質(zhì)點相比更容易發(fā)生振動[18],但結(jié)果中的拉梅常數(shù)為0顯然沒有意義且不合理。這是因為MOCZO等[21]中最后給出的參數(shù)平均方式是針對介質(zhì)內(nèi)部的固體/固體不連續(xù)界面的,并沒有考慮到自由表面的特殊性。因此,我們需要將自由表面的特殊性合理地引入到平均介質(zhì)的表達(dá)方式中。

圖1 推導(dǎo)自由表面處本構(gòu)關(guān)系時的彈性參數(shù)改變方式

根據(jù)MOCZO等[21]的平均介質(zhì)思想,兩個相連的各向同性彈性介質(zhì)可以用一個平均介質(zhì)來進(jìn)行等效,其滿足的本構(gòu)關(guān)系為:

(3)

式中:τA和εA分別為平均介質(zhì)的應(yīng)力和應(yīng)變張量,也可以理解為兩個介質(zhì)相連界面處的應(yīng)力和應(yīng)變張量,在二維情況下其分量形式為:

(4)

而EA是平均介質(zhì)的彈性矩陣,同時包含了兩種介質(zhì)各自的拉梅常數(shù),即:

(5)

其中:

(6)

接下來通過逐步改變式(5)中上部介質(zhì)的物理參數(shù),使得這一固體/固體邊界的物理性質(zhì)逐步近似于自由表面,進(jìn)而得到可以描述自由表面處應(yīng)力和應(yīng)變張量的本構(gòu)關(guān)系。

首先,使上部介質(zhì)的剪切模量μ=0(如圖1),從而上述的固體/固體邊界轉(zhuǎn)變?yōu)榱黧w/固體邊界,其相應(yīng)的彈性矩陣可以通過將μ-=0代入公式(5)得到,即:

(7)

其中:

(8)

由于空氣介質(zhì)屬于流體介質(zhì)中的一種,所以這一處理方式得到的公式(7)同樣可以用來表示空氣/固體邊界。而自由表面這一概念則是對于地球表面這一空氣/固體邊界在真空近似下的描述,因此一種比較直接的方式是使公式(7)中上部流體介質(zhì)的λ=0,來獲得真空/固體邊界處的彈性矩陣,即自由表面處的本構(gòu)關(guān)系。然而,令λ-=0會使得(8)式中0出現(xiàn)在分母的位置上,進(jìn)而使得整個彈性矩陣(7)無意義。一種數(shù)學(xué)上嚴(yán)格且合理的處理方式是利用極限思想,即讓λ-趨近于0而不是等于0(如圖1),從而可以在極限意義上得到自由表面處的彈性矩陣:

(9)

這里使λ-接近0但不等于0可以很好的表達(dá)出地球表面的物理性質(zhì),即上部空氣介質(zhì)的彈性模量遠(yuǎn)遠(yuǎn)小于地下介質(zhì),而使λ-趨近于0則是體現(xiàn)了自由表面這一概念中的真空近似思想。將(9)式代入到(4)式中,我們可以得到自由表面處應(yīng)力和應(yīng)變張量滿足的本構(gòu)關(guān)系,即:

(10)

此外,引入介質(zhì)泊松比

(11)

并消去λ,(10)式還可以表示為:

(12)

從中我們可以發(fā)現(xiàn)一個有趣的現(xiàn)象。當(dāng)?shù)叵陆橘|(zhì)是非常致密堅硬的巖石時,即泊松比接近于0,公式(12)在這一極限條件下變?yōu)?

(13)

這一結(jié)果與MITTET[18]從TI介質(zhì)近似角度提出的自由表面處理方式的結(jié)果相同。而當(dāng)?shù)叵陆橘|(zhì)為松散的沉積物時,即泊松比接近于0.5時,公式(12)在這一極限條下可以表示為:

(14)

其結(jié)果與XU等[19]提出的聲波/彈性波邊界方法(AEA)相一致。換句話說,MITTET[18]和XU等[19]中的自由表面處理方式分別是新得到的本構(gòu)關(guān)系(公式(10))在小泊松比和大泊松比下的極限情況。這也意味著這一自由表面處的本構(gòu)關(guān)系將適用于更大泊松比范圍內(nèi)的自由表面表達(dá)。

因此,我們可以使用(2)式中的密度修正和(10)式中新推導(dǎo)出的本構(gòu)關(guān)系來表達(dá)自由表面,并且這一表達(dá)方式能夠很好地整合到原始的彈性動力學(xué)方程中,從而得到自由表面處的彈性波控制方程,即:

(15)

和修正后的本構(gòu)關(guān)系:

(16)

其中ux和uz分別是位移向量沿x方向和z方向的分量。由于(16)式是泊松比自適應(yīng)的,所以我們將(15)式和(16)式合稱為自由表面的自適應(yīng)表達(dá)方式。簡單的講,這一表達(dá)方式僅僅是對原始彈性動力學(xué)方程中的密度和本構(gòu)關(guān)系的一個修正,因此在計算過程中是不會增加額外的CPU和內(nèi)存需求。

2 在速度-應(yīng)力交錯網(wǎng)格有限差分中的應(yīng)用

本節(jié)中,我們將給出自由表面自適應(yīng)表達(dá)方式在時間域有限差分?jǐn)?shù)值模擬中的具體實現(xiàn)方式。我們采用顯式的交錯網(wǎng)格有限差分離散格式,主要是考慮到這種差分格式不會出現(xiàn)高泊松比值時數(shù)值模擬精度降低這一問題[20],從而保證自由表面表達(dá)方式的泊松比自適應(yīng)性不會因為數(shù)值格式的離散化而受影響。

以VIRIEUX[23]提出的二階速度-應(yīng)力交錯網(wǎng)格有限差分格式為例,如圖2,i,j分別為x方向和z方向格點坐標(biāo),自由表面位于j=0處,(15)式和(16)式離散后可得:

(17)

和

(18)

式中:速度分量vx=?ux/?t,vz=?uz/?t;Δt是時間步長;Δx和Δz分別為x方向和z方向的網(wǎng)格間距。值得注意的是,由于速度垂直分量vz在離散過程中位于自由表面之下的半網(wǎng)格點上(如圖2),屬于介質(zhì)內(nèi)部點,所以(17)式中的第二個式子不需要對密度進(jìn)行減半修正。此外,(17)式中的應(yīng)力張量組分需要根據(jù)其相對于自由表面的分布情況分為3類進(jìn)行計算。

圖2 自由表面的自適應(yīng)表達(dá)在速度-應(yīng)力交錯網(wǎng)格中的離散方式

第一類,位于自由表面之上(j<0),這種情況下,需要根據(jù)真空近似的思想,將應(yīng)力張量組分全部置0,即:

(19)

第二類,位于自由表面處(j=0),此時應(yīng)力張量組分需要滿足自由表面處的本構(gòu)關(guān)系((10)式),(18)式給出了這一類應(yīng)力張量組分在時間迭代過程中的計算方式。

第三類,位于自由表面之下(j>0),此時應(yīng)力張量組分均位于介質(zhì)內(nèi)部,不再需要考慮自由表面的影響,可以直接根據(jù)這一層介質(zhì)的本構(gòu)關(guān)系進(jìn)行計算,即:

(20)

公式(17)、公式(18)、公式(19)、公式(20)只是用來計算自由表面處的速度波場,而模型內(nèi)部其他位置處的波場仍然采用傳統(tǒng)的速度-應(yīng)力交錯網(wǎng)格有限差分格式[23]計算,以這種方式可以很容易地在地震波正演數(shù)值模擬中引入自由表面的影響,從而實現(xiàn)全波場的模擬計算。類似地,按照以上3類處理方式將很容易得到這一自由表面表達(dá)在高階有限差分格式中的實現(xiàn)。

3 數(shù)值實驗與討論

為了討論自由表面自適應(yīng)表達(dá)方式在彈性波時間域數(shù)值模擬中的計算效果,我們設(shè)計了兩個數(shù)值實驗。實驗一是模擬均勻半空間中的波場,也就是Lamb問題計算。由于Lamb問題存在解析解[24],因此可以將由自由表面自適應(yīng)表達(dá)方式得出的數(shù)值模擬結(jié)果與解析解對比,從而評估這一自由表面自適應(yīng)表達(dá)方式的準(zhǔn)確性。此外,為了測試其在不同泊松比值情況下的效果,我們設(shè)計了4種具有不同泊松比的模型,σ分別為0.200,0.300,0.400和0.489,并將自由表面自適應(yīng)表達(dá)方式計算出的結(jié)果與使用MIT-TET[18]和XU等[19]方法得到的結(jié)果進(jìn)行對比。實驗二采用了比較接近真實地下結(jié)構(gòu)的Overthrust模型,用來測試這一自由表面自適應(yīng)表達(dá)方式在復(fù)雜模型地震波全波場模擬中的可行性。

實驗一中介質(zhì)的縱波速度為5500m/s,密度為1500kg/m3,橫波速度分別為3368,2940,2245,807m/s,分別對應(yīng)于泊松比σ為0.2,0.3,0.4,0.489的情況。模型網(wǎng)格大小為800×200。為了模擬半空間,除頂部的自由表面邊界外,其余邊界處均使用PML吸收邊界條件來消除人工反射。震源是一個單點垂直力,施加在自由表面下的第11個網(wǎng)格層,其震源時間函數(shù)是主頻為30Hz的雷克子波。在模擬過程中,我們對3種自由表面的處理方式均使用二階交錯網(wǎng)格有限差分格式。此外,為了突出3種數(shù)值模擬結(jié)果的精度差異,這里不考慮將計算網(wǎng)格加密,僅采用二階格式體波頻散分析得出的采樣準(zhǔn)則[23],即最小波長內(nèi)網(wǎng)格點數(shù)為10。

圖3和圖4分別給出了不同泊松比值情況下模擬得到的質(zhì)點速度水平分量和垂直分量的地震記錄,其中紅色實線為Lamb問題的解析解結(jié)果,藍(lán)色虛線為MITTET[18]方法結(jié)果,綠色虛線為XU等[19]提出的AEA方法結(jié)果,黑色虛線為自適應(yīng)自由表面表達(dá)方式結(jié)果。對比可以看出,3種方法在體波的模擬上都比較精確,而模擬結(jié)果的差異體現(xiàn)在瑞利波上。當(dāng)泊松比較小時,由MITTET方法[18]計算出的地震記錄與解析解一致性較好,而XU等[19]方法得出的結(jié)果中瑞利波的到達(dá)時間較早,與解析解有明顯的偏差。隨著介質(zhì)泊松比的增加,XU等[19]方法計算出的地震記錄與解析解越來越接近,特別是當(dāng)σ=0.489時,兩者的地震記錄基本上完全重合,但使用MITTET[18]方法模擬得到的瑞利波則顯示出到時延遲和數(shù)值頻散現(xiàn)象,且越來越偏離解析解。這意味著這兩種方法的模擬效果與泊松比有關(guān),二者的適用性受泊松比值范圍的限制。與這兩種方法不同的是,自由表面自適應(yīng)表達(dá)方式計算出的結(jié)果無論是在小泊松比值還是大泊松比值情況下,都顯示出與解析解較好的一致性。為了定量評估3種方法的精度,表1給出了正演模擬地震記錄的L2范數(shù)誤差[14]。可以看出,通過自由表面自適應(yīng)表達(dá)方式計算得到的數(shù)值結(jié)果,其準(zhǔn)確度基本達(dá)到了90%,而在小泊松比值和大泊松比值情況下,其結(jié)果分別與MITTET[18]和XU等[19]方法得到的結(jié)果相近,進(jìn)而證實了第1節(jié)中關(guān)于自由表面自適應(yīng)表達(dá)方式在極端泊松比值情況下的討論結(jié)果。

圖3 不同泊松比值情況下的Lamb問題合成地震記錄(質(zhì)點速度水平分量vx)對比a σ=0.200; b σ=0.300; c σ=0.400; d σ=0.489

圖4 不同泊松比值情況下的Lamb問題合成地震記錄(質(zhì)點速度垂直分量vz)對比a σ=0.200; b σ=0.300; c σ=0.400; d σ=0.489

在實驗二中,我們選取Overthrust模型中的部分結(jié)構(gòu)來進(jìn)行數(shù)值模擬測試,這是考慮到在實際情況中,推覆體構(gòu)造更容易造成地表附近的介質(zhì)泊松比值變化范圍大且不連續(xù)。圖5展示了計算中使用模型的縱波速度、橫波速度、密度和泊松比。在模擬計算中,網(wǎng)格大小為269×1601,網(wǎng)格間距為2.8m。模型頂部邊界采用自適應(yīng)自由表面表達(dá)方式進(jìn)行處理,其余邊界使用PML吸收邊界條件消除人工反射。與實驗一相同,震源為單點垂直力,作用在(x=2990.4m,z=56m)處,其震源時間函數(shù)為主頻15Hz的雷克子波。圖6和圖7分別給出了正演模擬得到的不同時刻速度波場的水平分量和垂直分量,從中可以很好的觀測到由于自由表面引起的反射縱波和轉(zhuǎn)換橫波,以及強(qiáng)振幅瑞利波的產(chǎn)生和傳播。此外,通過與縱波速度和橫波速度的疊加顯示,還可以看出瑞利波在層狀介質(zhì)中的速度頻散和泊松比橫向突變點處的散射。這些都表明自由表面自適應(yīng)表達(dá)方式的引入能夠很好地產(chǎn)生出與自由表面相關(guān)且合乎物理規(guī)律的波現(xiàn)象,從而證實了其在復(fù)雜模型正演模擬中的可行性。

表1 3種自由表面實現(xiàn)方式在不同泊松比情況下的L2誤差比較

注:M為Mittet方法結(jié)果,X為Xu等提出的AEA方法結(jié)果,A為自適應(yīng)自由表面表達(dá)方式結(jié)果。

圖5 Overthrust模型a 縱波速度; b 橫波速度; c 密度; d 泊松比

圖6 波場快照(質(zhì)點速度水平分量vx)a t=0.24s時刻; b t=0.78s時刻; c t=1.05s時刻

圖7 波場快照(質(zhì)點速度垂直分量vz)a t=0.24s時刻,; b t=0.78s時刻; c t=1.05s時刻

4 結(jié)論

本文以平均介質(zhì)理論和極限思想為基礎(chǔ),通過分析自由表面的物理性質(zhì),提出了自由表面邊界的另一種新的自適應(yīng)表達(dá)方式。簡單地講,這一表達(dá)方式是對原始彈性動力學(xué)方程在自由表面處的密度和本構(gòu)關(guān)系的修正,因此在數(shù)值模擬中不會產(chǎn)生額外的計算量和存儲需求。在時間域速度-應(yīng)力交錯網(wǎng)格有限差分中的應(yīng)用證實了這一表達(dá)方式能夠很容易地在有限差分?jǐn)?shù)值模擬方法中實現(xiàn)。數(shù)值實驗和理論分析均表明,這一自由表面自適應(yīng)表達(dá)方式是泊松比自適應(yīng)的,對于不同大小的泊松比值均能給出較為精確的模擬結(jié)果,在空間采樣率為最小波長內(nèi)10個網(wǎng)格點的二階交錯網(wǎng)格差分格式中,其準(zhǔn)確度達(dá)到90%。此外,推覆體模型的正演模擬證實了這一表達(dá)方式在復(fù)雜模型正演中的可行性,進(jìn)而可以應(yīng)用于實際情況中近地表不均勻風(fēng)化層內(nèi)的復(fù)雜波場的模擬和分析。