基于Python語言的瑞利面波頻散反演實行方案

吳為治，婁利，王鵬，王賓

(安徽省勘查技術(shù)院，安徽 合肥 230041)

0 引言

Python語言作為一種解釋型的高級語言，其語法簡單、代碼可讀性較高且易擴展移植，近年來迅速發(fā)展成為科學、工程研究領(lǐng)域內(nèi)流行的編程語言。Python語言還具有功能強大且開源標準函數(shù)庫以及第三方函數(shù)庫，眾多科學計算軟件都提供有Python的安裝調(diào)用接口，Github網(wǎng)站也提供Python軟件包的開源以方便開發(fā)者搭建展開研究工作的平臺。正是這樣，調(diào)用高效完善的函數(shù)庫，既能滿足專業(yè)化研究的需求，又能進一步提高開發(fā)效率。彭一波等[1]提出利用Python 語言地震學數(shù)據(jù)處理軟件包obspy進行背景噪聲研究的工作方案，實現(xiàn)了兩個地震臺站之間噪聲互相關(guān)格林函數(shù)波頻散曲線的提取，并進行背景噪聲源的時空性研究。這是利用Python語言建立地球物理研究平臺一次成功的案例，然而其并未進行進一步頻散反演。我們知道面波頻散曲線含有豐富的地下層狀結(jié)構(gòu)信息,尤其對橫波速度結(jié)構(gòu)特別敏感，所以通過對頻散曲線的反演可以得到地下層狀介質(zhì)的物理性質(zhì)[2]。因此基于Python語言的頻散反演方法的研究很有必要，可以為今后繼續(xù)研究背景噪聲成像、面波勘探提供開源的支持。

1 研究基礎(chǔ)

1.1 頻散正演算法

頻散曲線的反演基于正演算法，面波頻散的正演計算方法眾多，如Knopoff算法[3]，廣義反射—透射系數(shù)法[4]等，均能快速準確地求取面波頻散。Python語言第三方函數(shù)庫提供有disba軟件包(https://pypi.org/project/disba)用于計算頻散，該軟件包使用Dunkin矩陣或者快速增量矩陣算法計算瑞利波相和群速度頻散曲線；使用Thomson-Haskell方法計算勒夫波相和群速度頻散曲線。GitHub上則有提供pysurf96軟件包(https://github.com/miili/pysurf96)用來計算頻散，其算法基于圣路易斯大學的CPS(computer programs in seismology)中surf96程序[5]，雖然需要Fortran編譯環(huán)境，但穩(wěn)定性更好。這里我們選擇用pysurf96程序來測試(附錄示例1)，選取一個簡單的兩層模型(表1)來正演瑞利面波的頻散曲線(圖1)。

表1 兩層層狀介質(zhì)模型參數(shù)

圖1 基于pysurf96計算的兩層層狀介質(zhì)模型的瑞利面波頻散曲線Fig.1 Rayleigh wave dispersion curve of two layered media model based on pysurf96

從圖1中可以看出基階瑞利面波相速度在低頻部分接近模型半空間橫波速度的0.9倍左右，高頻部分接近模型上層橫波速度的0.9倍左右；二階面波相速度在高頻部分趨近模型半空間橫波速度；三階面波相速度則在低頻接近模型上層的橫波速度。這與夏江海在《高頻面波方法》中提到的現(xiàn)象一致[2]。同時，我們在測試不同速度模型中發(fā)現(xiàn)當?shù)貙幽Ｐ椭邪肟臻g的波速過低時，會無法得到某些頻率的瑞利面波相速度，如使用軟件包中surf96函數(shù)計算頻散會報錯，使用surfdisp96函數(shù)則需要先定義返回值為復數(shù)。這是因為頻散方程在這些頻率的根是復數(shù)，這里可以采用在原半空間的下方加上一層速度更高的層作為新半空間來計算頻散。但這樣會帶來一定的誤差[6]，最好是取頻散方程復數(shù)根解的實部作為相速度來解決這一問題[7]。

1.2 頻散反演理論和算法

頻散曲線的反演是通過某種方式調(diào)節(jié)模型參數(shù)達到對實際頻散曲線的最好擬合。這里我們構(gòu)造目標函數(shù):

Φ(m)=‖d-d(m)‖2,

(1)

即實際頻散相速度d與速度結(jié)構(gòu)模型m正演得到頻散相速度d(m)的L2范數(shù),來定義理論計算數(shù)據(jù)對實際數(shù)據(jù)的擬合程度。反演的過程實際上就是采用算法搜索調(diào)整m,使得目標函數(shù)Φ(m)為極小的過程,若m有多個,則取其目標函數(shù)值最小的那一個作為反演的最終結(jié)果。最初頻散曲線反演方法主要利用半波長法、拐點法、漸進線法等。隨著計算能力的進步，算法的開發(fā)，越來越多反演方法應(yīng)用到頻散反演中，尤其是啟發(fā)式算法的出現(xiàn)，讓反演效率和速度提升到一個新高度：凡友華等[8]使用的BFGS擬牛頓法引入到頻散曲線的反演當中,提出了基于基階模的擬牛頓法反演,在理論上進行了反演模擬;石耀霖等[9]以勒夫波為例討論從面波頻散反演地球內(nèi)部構(gòu)造的遺傳算法;翟佳羽等[10]在蟻群算法引入集中度的概念并對層狀介質(zhì)模型以及實測數(shù)據(jù)的頻散曲線進行反演,驗證了方法的有效性。這些頻散反演的方法大致可以分為以微分為基礎(chǔ)的方法、枚舉算法和隨機算法。前面提到的BFGS擬牛頓法即利用了目標函數(shù)的一階導數(shù)信息，應(yīng)用于頻散反演這種非線性化反演時,依賴初始模型的選擇,容易出現(xiàn)不收斂的現(xiàn)象。遺傳算法屬于全局搜索方法，容易陷入局部極小值中,且反演用時較長、效率較低，但可進行優(yōu)化。模擬退火法則是一種偽隨機算法，通過用溫度代表概率進行演算來搜尋最優(yōu)解，更依賴初始模型和參數(shù)選擇。

郭飛在Github上提供的scikit-opt(https://github.com/guofei9987/scikit-opt)軟件包封裝了多種啟發(fā)式算法函數(shù)，并且在軟件包的網(wǎng)絡(luò)在線說明文檔(https://scikit-opt.github.io/scikit-opt/#/zh/README)中給出了算法特性、詳細的參數(shù)說明和使用案例。通過對提供的算法函數(shù)分析，其中可用于頻散反演的有：遺傳算法、差分進化算法、粒子群算法和模擬退火法。

2 算法設(shè)計和模型測試

首先需要實現(xiàn)頻散正演的子程序，輸入速度模型和需要計算的周期數(shù)數(shù)組，調(diào)用pysurf96程序包函數(shù)進行正演計算：如使用的是surf96函數(shù)，則要在模型下方增加一個更高速的半空間進行計算；如使用的是surfdisp96函數(shù)，則要取返回值的實數(shù)部分。然后實現(xiàn)輸入不同模型時計算目標函數(shù)的子程序，最后在主程序?qū)崿F(xiàn)文件數(shù)據(jù)的讀寫、參數(shù)的設(shè)定和反演算法的運算。其中，Python的numpy擴展庫提供豐富的數(shù)組計算函數(shù)，為我們計算提供了方便，代碼也顯得更加簡潔。例如在計算目標函數(shù)時需要進行L2范數(shù)運算，python語言可以用numpy.linalg.norm()方便快捷算出，省去了復雜數(shù)組計算的代碼(附錄示例2)。

a—差分進化算法;b—遺傳算法;c—粒子群算法;d—模擬退火算法a—differential evolution algorithm;b—genetic algorithm;c—particle swarm optimization;d—simulated annealing algorithm圖2 不同啟發(fā)式算法反演的結(jié)果Fig.2 Results of different heuristic algorithms

可以看出深度越淺，橫波速度結(jié)構(gòu)反演的精度越高，例如圖2a差分進化算法左圖第一、二、三層反演的效果都很好，且無論是高速層還是低速層，都很好地反演出來。這是因為瑞利面波高頻對淺層敏感；最后一層即半空間層的橫波速度也都反演的較好。圖2a右圖則是反演得到模型的頻散與實際頻散曲線的差別，可以看出擬合很好。至于中間兩層誤差較大，我們可以利用disba軟件計算瑞利面波的相速度敏感核進行探討，圖3是頻率為2 Hz的瑞利面波相速度敏感核。從圖中紅線和黃線可以看出，縱波速度和密度的改變總體都對相速度影響很?。痪G線的變化說明橫波速度在0～20 m淺層變化對相速度影響大，30～50 m范圍影響則較小，半空間層橫波速度變化則又對相速度起一定影響。這從另一個方向解釋了上述反演的結(jié)果，也可以看出各層反演得到速度結(jié)果的可靠性。

圖3 頻率為2 Hz瑞利面波相速度敏感核Fig.3 Phase sensitivity kernel of Rayleigh wave on 2 Hz

3 算法穩(wěn)定性和效率

為了測試算法穩(wěn)定性和運行效率，這里所有運算過程均使用單CPU核心并且沒有其他計算任務(wù)在系統(tǒng)平臺上運行，系統(tǒng)平臺選擇Ubuntu以避免其他后臺程序干擾。測試使用的處理器為英特爾二代酷睿i5(i5-6200U)主頻2.3 GHz，且并沒有讓GPU參與運算。算法種群規(guī)模、迭代次數(shù)等參數(shù)均使用默認值,且均不加入外部判斷是否繼續(xù)或跳出反演代碼。測試結(jié)果如圖4所示。可以看出差分進化算法擬合度最高，但算法耗時最長；遺傳算法和粒子群算法反演耗時大幅度縮短，但擬合度較差；模擬退火法擬合度較好，耗時由于初始模型的選擇，這里接近差分進化算法?？紤]到默認的種群規(guī)模較大，差分進化算法實際使用可以考慮減少數(shù)值來提升計算效率。遺傳算法和粒子群算法則可以加入合適判定增加反演迭代次數(shù)來提高擬合度的穩(wěn)定性。模擬退火法效率十分依靠初始模型的選擇，很適合已有初始模型的反演，但要注意模擬退火其他參數(shù)選擇也很重要，否則容易出現(xiàn)迭代次數(shù)或者溫度達到后跳出反演。

a—差分進化算法;b—遺傳算法;c—粒子群算法;d—模擬退火算法a—differential evolution algorithm;b—genetic algorithm;c—particle swarm optimization;d—simulated annealing algorithm圖4 不同啟發(fā)式算法的效率和穩(wěn)定性Fig.4 Efficiency and stability of different heuristic algorithms

scikit-opt軟件包還提供多CPU 目標函數(shù)計算和GPU加速功能，能成倍提升反演計算速度，并且可以進行斷點續(xù)行，能有效避免運算過程中部分函數(shù)出錯導致程序崩潰問題。隨著算力的提升，計算機科學的發(fā)展，能更好地進行地球物理學反演問題的研究。Zhang等[11]2021年應(yīng)用大數(shù)據(jù)先行計算好區(qū)域模型的格林函數(shù)并儲存，然后利于人工智能快速識別并反演震源破裂機制參數(shù)，在智能預警地震上取得了突破的進展。如今Python語言已有眾多大數(shù)據(jù)和人工智能算法庫支持，需要我們提供封包好的地球物理學正演和反演方法。引入大數(shù)據(jù)和人工智能，這將是未來地球物理研究重要的方向。

4 實例

前文中提到彭一波[1]在研究中利用obspy軟件包提取了海拉爾盆地的兩個地震臺站間背景噪聲互相關(guān)格林函數(shù)的頻散曲線，這里我們對其頻散曲線進行一個多層厚度5 km的速度模型反演(圖5)?？梢钥闯鰞蓚€臺站之間的地殼以及上地幔平均速度結(jié)構(gòu)：區(qū)域的淺層橫波速度較小,波速在3 km/s附近,這可能是因為沉積層較厚；在15～25 km處存在低速速層，莫霍面深度在40 km附近。這些與李英康[12]、李皎皎[13]等人研究結(jié)論一致?？紤]到背景噪聲提取的長周期面波頻散可靠性較低，上地幔以及深部結(jié)構(gòu)還需結(jié)合地震面波進行反演確定[14]。

圖5 遺傳算法反演海拉爾盆地面瑞利波頻散結(jié)果Fig.5 Inversion of Rayleigh wave dispersion curve in Hailar Basin by Genetic Algorithm

5 結(jié)論及討論

本文基于Python語言scikit-opt軟件包進行瑞利面波頻散曲線反演實行方案的研究，給出了頻散正演算法在實際反演過程中所遇到頻散函數(shù)解虛數(shù)問題的解決方法，搭建了瑞利面波頻散反演地下層狀速度結(jié)構(gòu)的平臺，為今后研究背景噪聲成像、面波勘探提供開源基礎(chǔ)。同時對比了scikit-opt軟件包中不同啟發(fā)式算法的效率和穩(wěn)定性，并利用相速度敏感核探討了反演結(jié)果的可靠性，為使用者在反演算法的選擇上提供一些依據(jù)。最后對彭一波[1]在海拉爾盆地背景噪聲研究提取的頻散曲線進行了反演，得到了該區(qū)域的地殼和上地幔層狀速度模型。

可以看出地球物理正演和反演算法的封包，能為廣大研究者提供一個可靠且高效的平臺，也為今后接入大數(shù)據(jù)和人工智能打下基礎(chǔ)。本文只進行了基階瑞利面波頻散的一維速度模型反演，現(xiàn)在研究趨勢是加入高階頻散進行二維、三維模型的反演，這是將今后繼續(xù)研究的方向。

致謝：感謝R.Hermann、郭飛等人算法的開源共享。

附錄

示例1：基于pysurf96程序包計算瑞利面波基階頻散腳本示例，通過修改mode數(shù)值可以計算高階瑞利面波頻散。

示例2：基于scikit-opt程序包中遺傳算法反演一維層狀速度模型腳本示例。腳本目錄需包含三個文件：模型上下邊界模型l.model和u.model；反演的目標頻散文件m.disp。反演的速度模型結(jié)果寫入本地文件fit.model。