基于快速多極間接邊界元法的高頻地震波散射模擬方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明涉及地震波散射模擬領(lǐng)域，具體是一種基于快速多極間接邊界元法的高頻 地震波散射模擬方法。
【背景技術(shù)】
[0002] 近年來，隨著計(jì)算機(jī)硬件水平飛速發(fā)展，大規(guī)模科學(xué)計(jì)算對(duì)于現(xiàn)代科技研發(fā)和大 型工程設(shè)計(jì)愈發(fā)重要。其中不規(guī)則地形對(duì)地震波的散射模擬是很多學(xué)科領(lǐng)域的研究熱點(diǎn)， 可廣泛應(yīng)用于復(fù)雜地形地震動(dòng)科學(xué)確定、地球物理勘探、巖土層參數(shù)反演等。理論分析方法 整體上可分為解析法和數(shù)值法，其中數(shù)值法包括域離散型的有限元法、有限差分法和邊界 型的邊界單元法等。同域離散型方法相比，邊界元法（BEM)具有高精度、降低問題求解維 度、便于處理高梯度應(yīng)力變化的優(yōu)點(diǎn)，同時(shí)特別便于處理無限域波動(dòng)問題（無需引入無反 射邊界條件），另外還克服了普通有限元法高頻彌散的難題。但該方法在實(shí)際應(yīng)用中的主要 缺陷在于待求解方程系數(shù)矩陣為非對(duì)稱稠密滿陣，這樣對(duì)于目前越發(fā)流行的大規(guī)模科學(xué)運(yùn) 算，其計(jì)算效率將大大降低，在三維問題上該弱點(diǎn)將更為明顯。
[0003] 針對(duì)傳統(tǒng)BEM的稠密矩陣特性，近些年來，很多學(xué)者開展了快速BEM的研究工作， 包括快速多極子展開（FMM)、小波方法、ACA技術(shù)、PFFT方法等。其中快速多極子方法發(fā)展 較早，目前理論構(gòu)架已臻于成熟。其基本原理是將多個(gè)源點(diǎn)和場(chǎng)點(diǎn)之間的相互作用運(yùn)算由 逐粒計(jì)算，通過多極展開及局部展開（聚合、轉(zhuǎn)移），轉(zhuǎn)換為粒組與粒組之間的計(jì)算，能夠極 大地降低大規(guī)模多點(diǎn)相互作用運(yùn)算所需的計(jì)算量和存儲(chǔ)量。目前該方法已被廣泛應(yīng)用到熱 學(xué)、電磁學(xué)、固體力學(xué)、流體動(dòng)力學(xué)和聲學(xué)等領(lǐng)域的大型復(fù)雜問題分析。需指出的是，對(duì)于固 體中彈性波散射問題，由于涉及到壓縮波和剪切波間的波型轉(zhuǎn)換（在不連續(xù)面附近可能還 衍生面波），波動(dòng)特征將更為復(fù)雜，因而FMM-BEM方法在該領(lǐng)域的發(fā)展較為緩慢。已有的快 速多極直接邊界元法（FMM-DBEM)由于其數(shù)學(xué)公式推導(dǎo)復(fù)雜，實(shí)現(xiàn)難度大以及收斂效果不 佳等原因一直處于待完善狀態(tài)。

【發(fā)明內(nèi)容】

[0004] 本發(fā)明針對(duì)不規(guī)則地形對(duì)高頻地震波散射問題，提出了一種基于快速多極間接邊 界元法的高頻地震波散射模擬方法。本方法采用快速多極間接邊界元方法（FMM-IBEM)4B 夠快速求解其他方法（解析法，有限差分法，有限元法以及常規(guī)邊界元法）無法完成的高頻 地震波散射問題，能在一臺(tái)普通個(gè)人電腦上完成上百萬自由度問題的三維散射問題，為高 頻或大尺度地震波散射模擬提供了一種有效方法，提高計(jì)算速度的同時(shí)大大降低電腦資源 占用，為實(shí)際工程的快速模擬帶來顯著的經(jīng)濟(jì)效益。
[0005] 為實(shí)現(xiàn)上述目的，本發(fā)明提供如下技術(shù)方案：
[0006] -種基于快速多極間接邊界元法的高頻地震波散射模擬方法，具體步驟為：
[0007] (1)建立邊界元模型
[0008] 考慮平面波入射，基于單層位勢(shì)理論，在邊界表面上施加三個(gè)正交方向的虛擬荷 載以構(gòu)造散射波場(chǎng)；在建立邊界元模型中，首先將邊界表面離散為三角形或四邊形，得到離 散單元；采用斜面圓盤均布荷載近似覆蓋在離散單元上；
[0009] ⑵FMM-IBEM實(shí)施過程
[0010] (2. 1)邊界離散并生成樹結(jié)構(gòu)
[0011] 將邊界離散后，采用第0層正方體包含整個(gè)計(jì)算域，將第0層正方體等分為四個(gè)第 一層正方體，再將每個(gè)第一層正方體等分為四個(gè)第二層正方體，這樣逐級(jí)劃分下去；當(dāng)最小 正方體的邊長(zhǎng)等于給定的數(shù)值時(shí)，停止劃分；然后刪去不含任何單元點(diǎn)的正方體，即可生成 樹結(jié)構(gòu)；定義第L層的正方體為父細(xì)胞，L彡0 ;而L+1層的正方體為其子結(jié)點(diǎn)，1彡1 ;包 含邊界單元點(diǎn)的最小正方體稱之為葉子結(jié)點(diǎn)；第0層正方體為根結(jié)點(diǎn)；其他正方體為樹結(jié) 占.
[0012] (2. 2)近場(chǎng)計(jì)算
[0013] 對(duì)于葉子結(jié)點(diǎn)，其鄰居以及自身立方體內(nèi)場(chǎng)點(diǎn)和源點(diǎn)間相互作用通過間接邊界元 法IBEM計(jì)算；這里的場(chǎng)點(diǎn)指單元內(nèi)部的積分點(diǎn)，源點(diǎn)指虛擬荷載和流量源作用圓盤面上的 離散點(diǎn)；
[0014] (2. 3)下行遍歷計(jì)算多極展開系數(shù)
[0015] 從葉子結(jié)點(diǎn)開始，將離散單元上虛擬均布荷載的作用傳遞到葉子結(jié)點(diǎn)，形成最初 的多極展開系數(shù)，然后通過多極展開中心點(diǎn)間的轉(zhuǎn)移公式Μ2Μ，逐層向上遞歸操作，直到求 出第二層樹結(jié)點(diǎn)的多極展開系數(shù)。
[0016] (2. 4)上行遍歷逐層計(jì)算局部展開系數(shù)
[0017] 對(duì)第二層的每一個(gè)樹節(jié)點(diǎn)，利用多極展開向局部展開轉(zhuǎn)移公式M2L，將相互作用 列表中的所有節(jié)點(diǎn)的多極展開系數(shù)傳遞到該樹節(jié)點(diǎn)的局部展開系數(shù)；對(duì)于父結(jié)點(diǎn)中所有源 點(diǎn)的影響，通過局部展開中心點(diǎn)間的轉(zhuǎn)移公式L2L得到，最后累加到該樹節(jié)點(diǎn)的局部展開 系數(shù)；如此直到求出葉子結(jié)點(diǎn)的局部展開系數(shù)；
[0018] (2. 5)邊界積分和迭代求解
[0019] 當(dāng)樹結(jié)構(gòu)的多極展開系數(shù)和局部展開系數(shù)計(jì)算完畢后，利用樹結(jié)構(gòu)對(duì)不同單元的 邊界積分求和；對(duì)某個(gè)單元的場(chǎng)點(diǎn)X，根據(jù)距離遠(yuǎn)近，源點(diǎn)類型分為三類：①對(duì)其中的近場(chǎng) 源點(diǎn)，即鄰居及本身，使用常規(guī)方法直接計(jì)算；②對(duì)該場(chǎng)點(diǎn)所在葉子"相互作用列表"中的 源點(diǎn)，采用源點(diǎn)所在葉子結(jié)點(diǎn)的多極展開系數(shù)計(jì)算；③對(duì)于更遠(yuǎn)的源點(diǎn)，即上層相互作用列 表，它們包含的邊界單元積分已通過多極展開系數(shù)的轉(zhuǎn)移、疊加和傳遞，集成到了 X所在的 葉子結(jié)點(diǎn)的局部展開系數(shù)；這樣逐個(gè)單元計(jì)算下去，使用樹結(jié)構(gòu)來代替系數(shù)矩陣與迭代向 量相乘所需的存儲(chǔ)和計(jì)算，即可得到各個(gè)單元上虛擬均布荷載的密度；然后采用GMRES算 法迭代求解直到滿足精度要求；
[0020] (2. 6)總波場(chǎng)計(jì)算
[0021] 得到各個(gè)單元上虛擬均布荷載的密度之后，采用多極展開計(jì)算各個(gè)單元對(duì)所關(guān)心 區(qū)域的位移格林函數(shù)，即得到散射波場(chǎng)；最后將散射波場(chǎng)和自由場(chǎng)疊加，即可得到總體波 場(chǎng)。
[0022] 作為本發(fā)明進(jìn)一步的方案：步驟（2. 5)中GMRES算法迭代求解采用預(yù)處理技術(shù)，即 擬采用近場(chǎng)格林函數(shù)解構(gòu)造一窄帶稀疏矩陣，對(duì)方程組進(jìn)行預(yù)處理，以加速GMRES方程求 解收斂。
[0023] 作為本發(fā)明進(jìn)一步的方案：該方法替換了常規(guī)IBffl的位移格林函數(shù)和應(yīng)力格林 函數(shù)的展開計(jì)算方式；常規(guī)IBEM的位移格林函數(shù)和應(yīng)力格林函數(shù)：
[0026] 上式中i, j = 1，2, 3,對(duì)應(yīng)于X，y, z方向；η;為場(chǎng)點(diǎn)X處邊界單元單位法向量與i 軸正方向夾角余弦；
[0027] 其中位勢(shì)函數(shù)：
[0029] 其中q為給定常數(shù)，后文取為k，h (橫波和縱波波數(shù)），復(fù)數(shù)i2= -1，r = I x-y I，X 和y分別為場(chǎng)點(diǎn)和源點(diǎn)；
[0030] 本方法使用樹結(jié)構(gòu)作為主要存儲(chǔ)和運(yùn)算對(duì)象，對(duì)核函數(shù)進(jìn)行展開和傳遞，借助 GMRES迭代算法，在每一次迭代中以樹結(jié)構(gòu)取代系數(shù)矩陣與迭代量相乘；通過迭代精度控 制得出結(jié)果；三維問題采用八叉樹結(jié)構(gòu)，對(duì)任意幾何形狀的彈性體邊界面，都會(huì)生成與之對(duì) 應(yīng)的、具有特定層次和分布的樹結(jié)構(gòu)；
[0031] 根據(jù)平面波展開理論，勢(shì)函數(shù)（3)可多極展開為：
?"為[-1，1]上高斯積分的橫坐標(biāo)及權(quán)系數(shù)，（ = ]^/|(|，r=|M,滿足且 |5| < · O)為第2類第1階球漢克爾函數(shù)，勒讓德函數(shù)；
[0037] 從而，得到了新的位移格林函數(shù)和應(yīng)力格林函數(shù)展開式：

[0042] 其中，C；1表示或V為局部展開系數(shù)；上述公式⑶和公式（9)采用樹 結(jié)構(gòu)存儲(chǔ)各局部展開系數(shù)及多極展開系數(shù)。
[0043] 與現(xiàn)有技術(shù)相比，本發(fā)明的有益效果是：
[0044] 本發(fā)明針對(duì)大型實(shí)際復(fù)雜非水平層狀沉積盆地對(duì)地震波的散射模擬存在的不足， 例如塔里木盆地，四川盆地等大型盆地存在的低頻不穩(wěn)定性。本發(fā)明主要解決了場(chǎng)地較為 堅(jiān)硬（剪切波速約為2000m/s~4000m/s)的大小大約在一百公里范圍內(nèi)的實(shí)際盆地的地 震動(dòng)模擬、場(chǎng)地較軟的中小型盆地模擬和高頻地震波散射模擬問題，實(shí)現(xiàn)了在普通個(gè)人電 腦上快速模擬高頻或大型不規(guī)則地形地震動(dòng)問題，使該問題的模擬實(shí)現(xiàn)普遍化和經(jīng)濟(jì)化。
[0045] 本發(fā)明結(jié)合快速多極算法（FMM)和間接邊界元法（IBEM)，給出了基于對(duì)角形式的 平面波展開的快速多極間接邊界元法，該方法數(shù)學(xué)公式推導(dǎo)簡(jiǎn)潔易懂，實(shí)現(xiàn)形式簡(jiǎn)便直觀