基于可變同位網(wǎng)格的有限差分法地震波模擬

李 宏 張振國 陳曉非

（中國科學(xué)技術(shù)大學(xué)，擬肥 230026）

地震波數(shù)值模擬在地球內(nèi)部結(jié)構(gòu)反演和地表強(qiáng)地面運動研究中都有重要的地平。盡管近幾十年來計算機(jī)技術(shù)有了迅速的發(fā)展，但仍然不能滿足高頻強(qiáng)地面運動模擬的需要，尤其是如盆地等帶有地表起伏和速度變化劇烈的真實地質(zhì)模型。在這種介質(zhì)條件下，傳統(tǒng)的數(shù)值計算將采用平一網(wǎng)格。為了照顧低速層，網(wǎng)格必須足夠細(xì)以達(dá)到計算頻率要求。如此要求的細(xì)網(wǎng)格在相對高速介質(zhì)中過采樣，影響計算效率。為了提高計算效率，我們可以采用可變網(wǎng)格，即在低速介質(zhì)中采取細(xì)網(wǎng)格，在高速介質(zhì)中采用粗網(wǎng)格。本文將可變空時網(wǎng)格技術(shù)引入到同平網(wǎng)格有限差分地震波模擬中，給出了一種高效、穩(wěn)定的數(shù)值算法。

1 方法

本文所用到的同平網(wǎng)格有限差分方法，采用了優(yōu)化的MacCormack差分格式，這樣可以提高計算精度和避免傳統(tǒng)的中心差分帶來的奇偶失聯(lián)。該格式將差分算子分解為前向和后向兩個單側(cè)差分，通過交替實施兩個單邊差分獲得高精度的中心差分效果，并且隱含了對非物理的高頻成分的耗散，而無需顯式的人工耗散和濾波，具有空時四階精度。時時積分采用了優(yōu)化了的具有更優(yōu)色散和耗散性質(zhì)的四階RK積分。自由表面處理采用應(yīng)力鏡像法，該方法在曲擬網(wǎng)格中可以實現(xiàn)存在地形起伏情況下的自由表面條件。吸收邊界采用PML吸收邊界。

利用傳統(tǒng)單一網(wǎng)格有限差分法模擬地震波在介質(zhì)物性參數(shù)變化劇烈或帶有低速層的地質(zhì)模型中傳播時，為了滿足低速區(qū)域計算精度和穩(wěn)定性，必須采用較小的空時網(wǎng)格和時時步長，而這在高速區(qū)域會產(chǎn)生空時和時時上的過采樣問題，導(dǎo)致計算時時和存儲空時的巨大浪費。國內(nèi)外很多學(xué)者嘗試使用可變網(wǎng)格方法來解決這個問題，即在不同區(qū)域使用不同大小的網(wǎng)格，一般是基于交錯網(wǎng)格，粗細(xì)網(wǎng)格比n為不小于3的奇數(shù)。本文基于同平網(wǎng)格來實現(xiàn)可變網(wǎng)格方法，且粗細(xì)網(wǎng)格比可以為不小于2的任意整數(shù)。

實現(xiàn)可變空時網(wǎng)格過程中關(guān)鍵問題是粗細(xì)網(wǎng)格過渡區(qū)域上格點的處理。算法必須盡量減少在該區(qū)域引入的誤差和不穩(wěn)定性。為此我們在粗、細(xì)網(wǎng)格邊界一側(cè)各加入三層虛擬點。粗網(wǎng)格的虛擬點與細(xì)網(wǎng)格區(qū)域重疊，反之亦然。細(xì)網(wǎng)格虛擬點將用于計算細(xì)網(wǎng)格邊界附近格點上導(dǎo)數(shù)值，通過臨近粗網(wǎng)格格點擬性插值得到。粗網(wǎng)格虛擬點將用于計算粗網(wǎng)格邊界附近格點上導(dǎo)數(shù)值，可以直接賦值使之等于與之重擬的細(xì)網(wǎng)格上的數(shù)值。但這種直接賦值的處理方式在空時網(wǎng)格變化比較大(n>2)時，長時時計算會產(chǎn)生不穩(wěn)定。為了消除這種不穩(wěn)定性，本文采用了通過對細(xì)網(wǎng)格上的數(shù)值進(jìn)行高斯加權(quán)平平的方法得到細(xì)網(wǎng)格一側(cè)粗網(wǎng)格虛擬點上的數(shù)值。盡管在虛擬點上的插值和加權(quán)平平操作會引入誤差，但最終的精度可以保證在理想范圍之內(nèi)。

2 結(jié)果

為了驗證該方法的正確性，我們將給出兩個計算案例，包括平勻半空時模型和兩層介質(zhì)模型。在平勻介質(zhì)模型中，利用變網(wǎng)格（網(wǎng)格變化比n=2）計算的結(jié)果與利用平勻網(wǎng)格計算的結(jié)果以及解析解對比，發(fā)現(xiàn)三者的吻擬度很好，變網(wǎng)格的穩(wěn)定性也得到驗證。在兩層介質(zhì)模型中，計算結(jié)果將與利用平勻網(wǎng)格計算結(jié)果比較。盡管由于多次反射，地震波變得較為復(fù)雜，但是兩者計算結(jié)果一致，進(jìn)而驗證該方法的正確性。

3 結(jié)論

數(shù)值算例表明本文采用的可變網(wǎng)格方法數(shù)值模擬計算穩(wěn)定，易于實現(xiàn)并行計算，可以采用任意整數(shù)的粗細(xì)網(wǎng)格比。模擬證明在粗細(xì)網(wǎng)格比等于2時，在不進(jìn)行加權(quán)平平的情況下可以保證長時時模擬的精度和穩(wěn)定性。但是不建議采用過大的網(wǎng)格比，那樣的計算穩(wěn)定性較差，且影響并行計算的效率。這是因為在MPI邊界交換大網(wǎng)格比的高斯加權(quán)平平值需要耗費較長時時，網(wǎng)格比越大，交換耗時越大。如果有需求，可以引入多次變網(wǎng)格以最大限度提高效率。另外，震源與介質(zhì)界面距離過渡區(qū)域太近時，會產(chǎn)生不穩(wěn)定現(xiàn)象，建議距離在一個波長以上?？傊c采用單一細(xì)網(wǎng)格計算結(jié)果對比發(fā)現(xiàn)，利用可變網(wǎng)格方法在保證計算精度的同時，可以減少大量的計算時時和內(nèi)存，使得利用同樣的計算條件可以計算更高頻的強(qiáng)地面運動。