李 宏 張振國(guó) 陳曉非

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

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

1 方法

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

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

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

2 結(jié)果

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

3 結(jié)論

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