基于lsqr法頻率域波形反演的存儲(chǔ)方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明屬于地震數(shù)據(jù)處理的全波形反演技術(shù)領(lǐng)域，涉及一種基于LSQR法頻率域 波形反演的存儲(chǔ)方法。
【背景技術(shù)】
[0002] 全波形反演方法具有重構(gòu)地下速度結(jié)構(gòu)，揭示復(fù)雜地下構(gòu)造與巖性細(xì)節(jié)信息的潛 力，既可在時(shí)間域也可在頻率域?qū)崿F(xiàn)。頻率域相對(duì)于時(shí)間域反演具有計(jì)算高效、數(shù)據(jù)選擇靈 活等優(yōu)勢(shì)。近十幾年來頻率域全波形反演理論在波場(chǎng)模擬方法、反演頻率選擇策略等方面 取得了極大進(jìn)展。頻率域全波形反演涉及求解大型稀疏線性方程組的問題，常用到的求解 方法是LU分解法，LU分解法需要對(duì)矩陣求逆存在占用內(nèi)存大，計(jì)算量大等不利于存儲(chǔ)和計(jì) 算的問題，這樣會(huì)使得對(duì)計(jì)算機(jī)硬件要求比較高。
[0003] 最小二乘QR(LSQR，Paige，1982)法是一種求解大型稀疏奇異方程的迭代方法，具 有占用內(nèi)存小，計(jì)算量小的優(yōu)點(diǎn)，能很容易地利用矩陣的稀疏性簡化計(jì)算，適合求解大型稀 疏方程問題，對(duì)計(jì)算機(jī)硬件要求很低。但是應(yīng)用到頻率域全波形反演中存在多源重復(fù)求解 方程的問題。

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

[0004] 本發(fā)明的目的是提供一種基于LSQR法頻率域波形反演的存儲(chǔ)方法，解決了現(xiàn)有 技術(shù)中計(jì)算效率差、內(nèi)存占用大、對(duì)計(jì)算機(jī)硬件要求高的缺點(diǎn)。
[0005] 本發(fā)明采用的技術(shù)方案是，一種基于LSQR法頻率域波形反演的存儲(chǔ)方法，按照以 下步驟實(shí)施：
[0006] 1)建立頻率域聲波方程；
[0007] 2)對(duì)阻抗矩陣進(jìn)行實(shí)數(shù)化處理；
[0008] 3)壓縮存儲(chǔ)阻抗矩陣；
[0009] 4)基于GPU加速的LSQR算法求解聲波方程；
[0010] 5)計(jì)算GPU中相關(guān)核函數(shù)；
[0011] 6)壓縮存儲(chǔ)后計(jì)算虛擬震源、梯度、Hessian矩陣；
[0012] 7)計(jì)算近似步長；
[0013] 8)實(shí)現(xiàn)頻域波形反演。
[0014] 本發(fā)明的有益效果是，采用GPU加速計(jì)算LSQR算法，這樣做有兩點(diǎn)好處：一是 LSQR計(jì)算速率明顯提高，二是壓縮存儲(chǔ)后使得CPU和GPU之間內(nèi)存交換的數(shù)據(jù)量減少，10 效率提高。
[0015] 本發(fā)明與傳統(tǒng)的完全LU分解法求解線性方程組相比，內(nèi)存占用量小，并且可以在 普通的PC機(jī)上輕松實(shí)現(xiàn)大型稀疏矩陣線性方程。在整個(gè)過程中對(duì)阻抗矩陣壓縮存儲(chǔ)，節(jié)省 了內(nèi)存開銷、提高了對(duì)梯度、Hessian矩陣的計(jì)算效率；近似的步長計(jì)算可以不需要額外的 正演模擬；采用GPU加速LSQR算法，使得求解線性方程組時(shí)計(jì)算效率提高數(shù)千倍；機(jī)群上 對(duì)每個(gè)節(jié)點(diǎn)單獨(dú)進(jìn)行GPU加速計(jì)算，與同等硬件條件下相比計(jì)算效率提高數(shù)千倍。
【附圖說明】
[0016] 圖1是本發(fā)明波形反演的整體流程圖。
【具體實(shí)施方式】
[0017] 下面結(jié)合附圖和【具體實(shí)施方式】對(duì)本發(fā)明進(jìn)行詳細(xì)說明。
[0018] 本發(fā)明基于LSQR法頻率域波形反演的存儲(chǔ)方法，主要是要解決構(gòu)造阻抗矩陣、求 解聲波方程、梯度、Hessian矩陣、步長等問題，其中求解聲波方程是決定算法可行性與高效 性的關(guān)鍵所在。本發(fā)明是在對(duì)阻抗矩陣壓縮存儲(chǔ)的基礎(chǔ)上，采用LSQR求解方程組的解，并 采用GPU加速實(shí)現(xiàn)LSQR算法，以達(dá)到減少內(nèi)存占有量提高計(jì)算效率的目的，進(jìn)一步使得在 普通的計(jì)算機(jī)上能夠?qū)崿F(xiàn)波形反演。
[0019] 本發(fā)明基于LSQR法頻率域波形反演的存儲(chǔ)方法，參照?qǐng)D1所示，按照以下步驟實(shí) 施：
[0020] 1)建立頻率域聲波方程
[0021] 如式（1)，根據(jù)九點(diǎn)差分格式并加入PML吸收邊界構(gòu)造阻抗矩陣，
[0023]
為拉普拉斯（Laplace)算子，U(x, z, ω)為介質(zhì)中每一點(diǎn)的 波場(chǎng)值（例如：壓強(qiáng)或位移），ω為角頻率，ν(χ，ζ)為介質(zhì)速度，F(xiàn)(c〇)為震源波場(chǎng)。將波 場(chǎng)空間劃分為NxXNz的網(wǎng)格，令Μ = NxXNz，寫成矩陣形式：
[0024] SU = F, (2)
[0025] 其中，U是復(fù)數(shù)離散波場(chǎng)值，即為所求頻率域內(nèi)波長值，排列為MX 1的列向量；S是 一個(gè)MXM的復(fù)值阻抗矩陣，F(xiàn)是震源項(xiàng)，也是一個(gè)MX 1的列向量。
[0026] 2)對(duì)阻抗矩陣進(jìn)行實(shí)數(shù)化處理
[0027] 本發(fā)明采用LSQR方法求解方程組（2)的解，首先將復(fù)系數(shù)矩陣轉(zhuǎn)化為實(shí)系數(shù)分塊 矩陣⑵式可改寫為：
[0028] (SRe+iSIn) (URe+iUIn) = FRe+iFIn, (3)
[0029] 其中，Sfe、SIni分別為阻抗矩陣的實(shí)部和虛部，U fe、UIni分別為波場(chǎng)值的實(shí)部和虛部， 分別為震源的實(shí)部和虛部，i2 = -1為虛數(shù)單位。
[0030] 將式（3)式寫成分塊矩陣的形式可得：
[0032] 轉(zhuǎn)換成實(shí)數(shù)矩陣的阻抗矩陣大小為2MX2M。
[0033] 3)壓縮存儲(chǔ)阻抗矩陣
[0034] 根據(jù)九點(diǎn)差分構(gòu)造出來的阻抗矩陣某一行代表著模型中某一點(diǎn)處對(duì)整個(gè)模型的 波阻抗，因此一行中最多有9個(gè)非零元素，轉(zhuǎn)換為實(shí)數(shù)矩陣后最多為18個(gè)非零元素。首先 將矩陣S按行構(gòu)造，將每一個(gè)非零元素的點(diǎn)設(shè)五個(gè)變量m、Val、Val_p ;KX代表這個(gè)元 素在矩陣s中的列坐標(biāo)，隊(duì)代表這個(gè)元素在整個(gè)阻抗矩陣s當(dāng)前的非零元素的累計(jì)個(gè)數(shù)，Νη 代表這個(gè)元素所在行之前整個(gè)s阻抗矩陣中非零元素的累計(jì)個(gè)數(shù)，Val是該元素的值，Val_ P為該元素對(duì)模型參數(shù)的偏導(dǎo)數(shù)。依次將s中所有元素按上述方法存儲(chǔ)，可得壓縮存儲(chǔ)后的 5個(gè)向量，其中Kx、Val、Val_p向量中每一個(gè)向量最大可能大小為36MX1，隊(duì)和N n向量的大 小都為2MX1，與未壓縮的阻抗矩陣相比，這種壓縮存儲(chǔ)阻抗矩陣的方法占內(nèi)存空間縮小 了（2MX2MV(3X18X2M+2X2M) =M/28倍，一般Μ可取到數(shù)萬到數(shù)十萬。
[0035] 4)基于GPU加速的LSQR算法求解聲波方程
[0036] 在這個(gè)計(jì)算過程中主要將矩陣與向量相乘、向量歸一化、以及更新向量在GPU中 運(yùn)算，其余在CPU中計(jì)算。并且在計(jì)算過程中，CPU內(nèi)存與GPU內(nèi)存（顯存）之間的數(shù)據(jù)交 換主要為一些向量的二范數(shù)以及正交變換過程中的系數(shù)，不存在向量和矩陣的交換，這樣 就不會(huì)由于大量的數(shù)據(jù)交換導(dǎo)致10瓶頸。在顯存中存儲(chǔ)均為壓縮存儲(chǔ)的阻抗矩陣和相關(guān) 向量，數(shù)據(jù)量不大，一般的顯卡均能滿足。
[0037] 初始化過程GPU主要計(jì)算向量（震源項(xiàng)）的二范數(shù)、單位化以及阻抗矩陣與初始 解的二范數(shù)等，交換的數(shù)據(jù)是由GPU計(jì)算得到的二范數(shù)到CPU，在由CPU到GPU計(jì)算歸一化 向量（矩陣）。對(duì)角化過程中主要涉及一些矩陣與向量的乘積，因此這步主要在GPU中完 成，再將數(shù)據(jù)交換到CPU。正交變換的過程主要涉及到一些邏輯判斷、