曹 健,陳景波

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

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

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

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

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

從物理的角度而言,自由表面是地表以上的空氣介質(zhì)和地下彈性介質(zhì)之間的一個(gè)接觸界面,即流體/固體邊界。由于地表以上的空氣介質(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è)兩個(gè)具有不同介質(zhì)參數(shù)的各向同性彈性半空間,其分界面位于z=0處,上層介質(zhì)的密度和拉梅常數(shù)分別為ρ-,λ-和μ-,而下層介質(zhì)的密度和拉梅常數(shù)分別為ρ+,λ+和μ+(圖1)。當(dāng)上層介質(zhì)為真空時(shí),即ρ-=0,λ-=0和μ-=0,分界面z=0自然就變成了自由表面。根據(jù)MOCZO等[21]中提出的介質(zhì)參數(shù)平均方式,分界面z=0處的彈性參數(shù)可以表示為:

(1)

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

(2)

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

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

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

(3)

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

(4)

而EA是平均介質(zhì)的彈性矩陣,同時(shí)包含了兩種介質(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會(huì)使得(8)式中0出現(xiàn)在分母的位置上,進(jìn)而使得整個(gè)彈性矩陣(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)一個(gè)有趣的現(xiàn)象。當(dāng)?shù)叵陆橘|(zhì)是非常致密堅(jiān)硬的巖石時(shí),即泊松比接近于0,公式(12)在這一極限條件下變?yōu)?

(13)

這一結(jié)果與MITTET[18]從TI介質(zhì)近似角度提出的自由表面處理方式的結(jié)果相同。而當(dāng)?shù)叵陆橘|(zhì)為松散的沉積物時(shí),即泊松比接近于0.5時(shí),公式(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á)方式能夠很好地整合到原始的彈性動(dòng)力學(xué)方程中,從而得到自由表面處的彈性波控制方程,即:

(15)

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

(16)

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

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

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

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

(17)

和

(18)

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

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

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

(19)

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

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

(20)

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

3 數(shù)值實(shí)驗(yàn)與討論

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

實(shí)驗(yàn)一中介質(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吸收邊界條件來消除人工反射。震源是一個(gè)單點(diǎn)垂直力,施加在自由表面下的第11個(gè)網(wǎng)格層,其震源時(shí)間函數(shù)是主頻為30Hz的雷克子波。在模擬過程中,我們對3種自由表面的處理方式均使用二階交錯(cuò)網(wǎng)格有限差分格式。此外,為了突出3種數(shù)值模擬結(jié)果的精度差異,這里不考慮將計(jì)算網(wǎng)格加密,僅采用二階格式體波頻散分析得出的采樣準(zhǔn)則[23],即最小波長內(nèi)網(wǎng)格點(diǎn)數(shù)為10。

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

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

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

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

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

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

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

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

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

4 結(jié)論

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