邵祥奇 何兵壽* 史才旺

(①中國(guó)海洋大學(xué)海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，山東青島 266100； ②青島海洋科學(xué)與技術(shù)國(guó)家實(shí)驗(yàn)室海洋礦產(chǎn)資源評(píng)價(jià)與探測(cè)技術(shù)功能實(shí)驗(yàn)室，山東青島 266100； ③南方科技大學(xué)地球與空間科學(xué)系，廣東深圳 518000)

0 引言

近年來(lái)，全波形反演(FWI)方法在理論和實(shí)際應(yīng)用方面均取得了較大進(jìn)展[1-5]。但是，F(xiàn)WI的龐大計(jì)算量和過(guò)低的計(jì)算效率是限制該方法進(jìn)入三維地震勘探領(lǐng)域的瓶頸之一。目前，業(yè)界提高反演效率的思路主要有兩種： 一是提高收斂速度； 二是提升正演(反演需要不斷進(jìn)行多次正演)計(jì)算速度。前者通過(guò)改進(jìn)反演算法減少迭代次數(shù)，降低計(jì)算量[6-10]； 后者采用高性能的計(jì)算機(jī)硬件提高計(jì)算效率或采用降維思路降低計(jì)算量[11-26]。

本文主要討論第二種思路。FWI大部分運(yùn)算量集中在正演計(jì)算過(guò)程中(FWI的每次迭代至少需要三次正演)，因此業(yè)界普遍認(rèn)為高效的正演計(jì)算是解決FWI反演效率問(wèn)題的有效途徑，采用降維思路減少運(yùn)算量是提高正演效率的有效方法。降維技術(shù)的本質(zhì)是采用某種壓縮方式減少參與反演的炮數(shù)，主要方法包括震源編碼法[16-20]、炮采樣法[21-23]和主成分分析法[24-28]等。

Romero等[16]在地震資料的逆時(shí)偏移處理中提出了相位編碼技術(shù)，將多炮數(shù)據(jù)壓縮為幾個(gè)超級(jí)炮進(jìn)行計(jì)算，但編碼算法會(huì)在偏移結(jié)果中產(chǎn)生串?dāng)_噪聲。Krebs等[17]將相位編碼技術(shù)引入FWI領(lǐng)域，利用隨機(jī)的極性編碼函數(shù)壓制反演中的串?dāng)_噪聲，將所有炮轉(zhuǎn)換成一個(gè)超級(jí)炮集記錄并進(jìn)行反演，提升了運(yùn)算效率。此外，學(xué)者們還提出了平面波編碼[18]、頻域隨機(jī)相位編碼[19]、頻率組編碼[20]、頻域分割編碼[27]等方法，完善并豐富了編碼方案。

炮采樣法也是一種高效的加速策略。Díaz等[21]提出了一種隨機(jī)抽取炮集進(jìn)行反演的方法，每次迭代時(shí)僅使用部分炮，在保證效果的同時(shí)能夠減少運(yùn)算。Ha等[22]提出了在Laplace域循環(huán)炮采樣，不同于隨機(jī)炮采樣，該方法循環(huán)使用炮集中的所有炮，使反演中每一炮的使用次數(shù)較均衡。Shi等[23]將炮采樣技術(shù)應(yīng)用到頻率多尺度FWI中，在每個(gè)頻率段均隨機(jī)抽取炮集參與反演，該方法不僅能加速運(yùn)算，而且不會(huì)引入額外的串?dāng)_噪聲。

主成分分析法也是一種經(jīng)典的降維方法，Liu等[24]將主成分分析法應(yīng)用于FWI，減少了炮集的數(shù)量，提高了計(jì)算效率。段超然等[25]在高頻部分使用主成分分析法，使FWI在缺少低頻成分時(shí)也能高效收斂。

震源編碼技術(shù)實(shí)現(xiàn)簡(jiǎn)單、效果明顯，是目前應(yīng)用最廣泛的加速方法，但同時(shí)也存在一些缺陷，主要表現(xiàn)在傳統(tǒng)震源編碼技術(shù)嚴(yán)重依賴于觀測(cè)系統(tǒng)。震源編碼技術(shù)應(yīng)用的前提是采用檢波點(diǎn)固定的方式觀測(cè)，因而在模擬同時(shí)激發(fā)震源時(shí)，無(wú)法對(duì)單炮獲取特定炮檢距的地震記錄。對(duì)于滾動(dòng)排列觀測(cè)系統(tǒng)，震源編碼會(huì)造成模擬與觀測(cè)資料的不匹配。鑒于此，Choi 等[29]提出了基于歸一化互相關(guān)目標(biāo)函數(shù)的波形反演，提升了震源編碼對(duì)海上拖纜資料的反演精度。夏冬明等[30]在此基礎(chǔ)上發(fā)展了基于歸一化局部互相關(guān)目標(biāo)函數(shù)的反演方法，進(jìn)一步改善了對(duì)拖纜數(shù)據(jù)的反演效果。Huang等[31]提出了基于分頻編碼的最小二乘偏移方法，能夠真正實(shí)現(xiàn)拖纜資料中多炮數(shù)據(jù)的完全分離。Zhang等[32]將分頻編碼應(yīng)用于聲波FWI，并且討論了其在計(jì)算效率、內(nèi)存消耗、抗噪性方面的優(yōu)點(diǎn)。

在前人研究的基礎(chǔ)上，本文主要解決目前編碼方法只適用于特殊觀測(cè)系統(tǒng)的缺陷，提出了一種基于震源分頻編碼的混合域彈性波FWI方法。首先為同時(shí)進(jìn)行正向延拓的每一炮賦以不同的頻率，并利用該頻率對(duì)應(yīng)的諧波代替原時(shí)間域子波進(jìn)行正演，通過(guò)相位靈敏度檢測(cè)(Phase Sensitivity Detection，PSD)提取單頻波場(chǎng)，從而實(shí)現(xiàn)各炮頻率域波場(chǎng)的完全分離；然后實(shí)現(xiàn)多炮編碼反演。這樣既能壓制炮間串?dāng)_噪聲，又能對(duì)每一炮波場(chǎng)和記錄實(shí)現(xiàn)靈活處理，可以適應(yīng)排列滾動(dòng)的觀測(cè)系統(tǒng)。

1 震源編碼技術(shù)

混合域彈性波FWI的運(yùn)算量和炮集數(shù)成正比，運(yùn)用震源編碼技術(shù)將多炮數(shù)據(jù)進(jìn)行組合、正演，可有效減少計(jì)算量。震源編碼技術(shù)借助于編碼函數(shù)，將炮集疊加、組合得到超級(jí)炮記錄，利用常規(guī)方法對(duì)該超級(jí)炮記錄進(jìn)行反演。

在頻率域，超級(jí)炮記錄可表示為

(1)

式中： 上標(biāo)“～”表示超級(jí)炮；dn是第n炮得到的觀測(cè)地震記錄；γn表示第n炮的編碼函數(shù)序列；NS為總炮數(shù)。

相對(duì)應(yīng)的超級(jí)炮震源可表示為

(2)

(3)

(4)

式中：Pn是第n炮的正傳波場(chǎng)；Bn是第n炮的殘差波場(chǎng)。

推導(dǎo)得出超級(jí)炮對(duì)應(yīng)的梯度

(5)

式中：E(m)為目標(biāo)函數(shù)，m為待反演參數(shù)，本文指速度； R(·)表示取實(shí)部；上角標(biāo)“*”表示取共軛；A是波動(dòng)方程的阻抗矩陣，滿足AP=s，s為震源子波。

2 分頻編碼技術(shù)

2.1 分頻編碼策略

本文把分頻編碼技術(shù)引入彈性波波形反演中。其核心思想是：在正演過(guò)程中，同時(shí)對(duì)多個(gè)炮記錄進(jìn)行正向延拓，給每一炮賦以不同頻率的諧波震源，使同時(shí)正向延拓的各炮波場(chǎng)在頻譜上互不重合。同時(shí)，在炮點(diǎn)波場(chǎng)正傳和殘差反傳過(guò)程中，采用相位靈敏度檢測(cè)技術(shù)提取每炮的單頻波場(chǎng)，實(shí)現(xiàn)各炮波場(chǎng)的完全分離，避免炮間串?dāng)_。

2.2 基于相位靈敏度檢測(cè)的分頻多炮正演

分頻編碼技術(shù)的核心在于給每炮賦以不同頻帶的震源函數(shù)，避免多炮同時(shí)延拓中的炮間串?dāng)_。震源的頻帶一般由震源函數(shù)的中心頻率控制，不同炮可以選擇不同主頻的震源函數(shù)。但由于不同主頻的震源函數(shù)往往具有一定的頻帶寬度，當(dāng)炮數(shù)較多時(shí)，難以完全分離多炮頻譜，利用諧波作為震源可以避免這一缺陷。

利用諧波震源進(jìn)行正演，首先需要提取正確的波場(chǎng)信息。本文引入PSD正演方法實(shí)現(xiàn)頻率域波場(chǎng)的準(zhǔn)確提取。Nihei等[26]提出的基于PSD的正演方法用諧波作為震源，在正演過(guò)程中可以很方便地提取出不同頻率的單頻波場(chǎng)，實(shí)現(xiàn)多震源的同步模擬且互不干擾。

PSD方法是利用一個(gè)參考信號(hào)和參考信號(hào)的90°相移信號(hào)計(jì)算原諧波信號(hào)的頻譜。對(duì)于含有兩個(gè)頻率的信號(hào)，原信號(hào)εsig用振幅Esig和相位θsig表示為

εsig=Esig1cos(ω1t+θsig1)+Esig2cos(ω2t+θsig2)

(6)

式中下標(biāo)“1”和“2”分別對(duì)應(yīng)這兩個(gè)頻率。

選取兩個(gè)不同的頻率ω1和ω2， 令它們滿足以下關(guān)系

(7)

式中k為正整數(shù)。

以ω1為例，設(shè)計(jì)如下參考信號(hào)和相移信號(hào)

(8)

式中 ref1表示參考信號(hào)。

ω1對(duì)應(yīng)的單頻信號(hào)的振幅和相位[26]可以表示為

(9)

其中

(10)

Nihei等[26]證明PSD算法與傅里葉變換的功能是相同的，但是前者更適用于對(duì)諧波信號(hào)進(jìn)行處理。對(duì)只包含幾個(gè)等間隔頻率的信號(hào)，PSD算法僅需要選擇很小長(zhǎng)度的信號(hào)進(jìn)行運(yùn)算就能夠得到其頻率響應(yīng)，效率更高。

在時(shí)間域正演過(guò)程中，通常使用Ricker子波作為震源，但PSD算法需要以諧波震源作為輸入。為了得到與常規(guī)正演方法等價(jià)的結(jié)果，需要建立Ri-cker子波和諧波震源之間的數(shù)學(xué)聯(lián)系。實(shí)際上，諧波震源能夠視為Ricker子波的某一個(gè)或幾個(gè)頻率成分，所以僅需要將Ricker子波做離散傅里葉變換后，提取出某幾個(gè)頻率的響應(yīng)，再做傅里葉反變換，就能得到諧波震源。處理實(shí)際資料時(shí)，對(duì)真實(shí)的子波做類似操作即可。當(dāng)然，這樣做的前提是能夠獲取較為準(zhǔn)確的震源子波，否則使用不正確的子波將會(huì)導(dǎo)致錯(cuò)誤的反演結(jié)果。子波的提取是波形反演需要解決的一個(gè)重要問(wèn)題，但是本文的研究目標(biāo)是加快反演速度，暫不考慮如何提取子波的問(wèn)題。另外，本文方法可以與當(dāng)前的不依賴于子波的反演方法相結(jié)合[33-35]，即使利用錯(cuò)誤的震源子波進(jìn)行反演也能得到正確的反演結(jié)果。

本文基于PSD的彈性波正演方法流程如圖1所示。

利用PSD算法實(shí)現(xiàn)分頻多炮正演實(shí)驗(yàn)可以驗(yàn)證PSD算法的有效性。以簡(jiǎn)單的均勻介質(zhì)為例，網(wǎng)格數(shù)設(shè)為150×200，網(wǎng)格間距為10m，縱、橫波速度分別為2500、1500m/s。設(shè)有4個(gè)炮點(diǎn)，炮點(diǎn)深度為20m，炮點(diǎn)等間隔分布，橫向坐標(biāo)分別為400、800、1200和1600m，原始子波是主頻8Hz的Ricker子波，時(shí)間長(zhǎng)度設(shè)為4s，四個(gè)炮點(diǎn)分別同時(shí)正演2、3、4、5Hz的單頻波場(chǎng)。

圖1 基于PSD的彈性波分頻編碼正演流程

利用PSD方法進(jìn)行正演所得到的各炮點(diǎn)頻率響應(yīng)如圖2所示，其中幾乎見(jiàn)不到多炮的串?dāng)_。

為精確驗(yàn)證PSD算法能夠有效避免炮間串?dāng)_，做了如下對(duì)照實(shí)驗(yàn)：每次只設(shè)置一個(gè)震源一個(gè)頻率，其他參數(shù)與圖2保持一致，得到不含炮間串?dāng)_的波場(chǎng)。將兩次實(shí)驗(yàn)結(jié)果中橫坐標(biāo)1km處不同深度的波場(chǎng)值繪制成圖3。對(duì)比可知，兩次實(shí)驗(yàn)得到的波場(chǎng)幾乎一致，這也說(shuō)明PSD算法具有足夠的精度，可以在分頻率的多炮正演中得到準(zhǔn)確的頻率域波場(chǎng)值。

2.3 分頻編碼的彈性波混合域FWI

從理論上講，分頻編碼技術(shù)不會(huì)引入額外噪聲[31]。本文應(yīng)用分頻編碼方法以實(shí)現(xiàn)混合域的彈性波反演。與Huang等[31]實(shí)現(xiàn)的最小二乘偏移方法不同，混合域FWI一般采用多尺度反演策略，在每一個(gè)尺度下，可供使用的頻段是受限的，所以頻率的選擇在分頻編碼混合域反演中非常重要。

圖2 分頻編碼單次正演的各炮頻率響應(yīng)(實(shí)部)

圖3 分頻編碼與單諧波震源波場(chǎng)值對(duì)比

通常而言，隨多尺度反演的遞進(jìn)，可供選擇的頻段逐漸拓寬。如果對(duì)每一尺度反演都設(shè)置相同的頻率數(shù)目，結(jié)合上述兩個(gè)限制條件，就不難得出如下結(jié)論：小尺度反演時(shí)，頻率間隔大，Δω較大而PSD計(jì)算時(shí)長(zhǎng)?。淮蟪叨确囱輹r(shí)，頻率選擇較密集，Δω較小而PSD計(jì)算時(shí)長(zhǎng)大。反之，若Δω固定不變，則在大尺度反演時(shí)，可供選擇的頻率更少，進(jìn)一步造成能同時(shí)正演的炮數(shù)變少，正演次數(shù)就會(huì)明顯增加。本文選擇用較多的頻率個(gè)數(shù)完成多炮同時(shí)正演。具體到反演過(guò)程中，首先確定反演使用的頻率個(gè)數(shù)M，把每M炮數(shù)據(jù)組合成為一個(gè)超級(jí)炮(組合方法可以有多種，本文選擇把空間上鄰近的M炮數(shù)據(jù)進(jìn)行組合)；然后每次迭代前把M個(gè)頻率的數(shù)據(jù)隨機(jī)分配給組合中的M炮，讓迭代過(guò)程中各炮使用的頻率不斷變化，這樣可以使原始數(shù)據(jù)得到充分利用。

3 模型算例

3.1 分頻編碼對(duì)FWI的加速效果

為證明本文方法的加速效果，用部分Overthrust模型的正演記錄進(jìn)行測(cè)試。圖4為部分Overthrust縱、橫波速度模型，該模型網(wǎng)格數(shù)為500×160，網(wǎng)格間隔為25m×25m。正演記錄的炮點(diǎn)均勻分布于地表，第一炮位于0m處，炮間距為125m，共得到100炮合成地震記錄，全排列接收，道間距為25m。反演所用的初始模型(圖5)是將圖4進(jìn)行高斯平滑得到的。

圖4 部分Overthrust速度模型

圖5 反演使用的初始速度模型

采用多尺度反演策略，每個(gè)頻率組包含10個(gè)頻率(表1)，反演分5個(gè)頻段，每個(gè)頻段迭代15次。分別采用常規(guī)全部炮反演、傳統(tǒng)震源編碼(極性編碼)反演和分頻編碼反演，并將三者反演結(jié)果(圖6)進(jìn)行對(duì)比。

表1 分頻編碼反演的頻率組信息

圖7為三種方法的誤差對(duì)比， 反演誤差定義為[19]

(11)

式中：mFWI和mTRUE分別為反演結(jié)果和真實(shí)模型的彈性參數(shù)； ‖·‖2表示L2范數(shù)。

從圖6、圖7可以看出，分頻編碼方法和極性編碼方法的反演精度與常規(guī)方法相近。但就計(jì)算量而言，由于每次迭代需要4次正演(梯度計(jì)算需要兩次，步長(zhǎng)估計(jì)需要兩次)，因此常規(guī)全部炮反演需要計(jì)算30000(100×5×15×4)次正演，而本文方法與極性編碼一樣，將10炮數(shù)據(jù)組合為1炮進(jìn)行反演，因此正演次數(shù)僅為常規(guī)方法的十分之一。

當(dāng)然，計(jì)算效率的提升與組合的炮數(shù)有關(guān)。通常反演中能夠進(jìn)行組合的炮數(shù)受頻帶寬度的影響，不能把任意多的炮集進(jìn)行組合。此算例中，每個(gè)超級(jí)炮都由固定的10炮數(shù)據(jù)組成。而在一般的多尺度反演中，大尺度反演時(shí)給定的頻帶寬度往往較小，能夠選擇的頻率點(diǎn)(亦即能夠進(jìn)行組合的炮數(shù))也比較少；隨著尺度變小，頻帶可以更寬，從而可以用更多的炮進(jìn)行組合。但是由于混合域算法中需要存儲(chǔ)每個(gè)頻率的波場(chǎng)，更多的頻率點(diǎn)意味著存儲(chǔ)量的增加，同時(shí)也意味著需更多的計(jì)算完成PSD的頻率分離，因此最終選擇多少個(gè)頻率需視情況而定。

圖6 全排列接收時(shí)不同方法縱波速度(左)、橫波速度(右)反演結(jié)果

圖7 誤差分析對(duì)比

3.2 滾動(dòng)排列觀測(cè)系統(tǒng)反演測(cè)試

應(yīng)用傳統(tǒng)震源編碼方法反演時(shí)，每炮的接收排列必須相同且覆蓋整個(gè)工區(qū)，但是排列滾動(dòng)的數(shù)據(jù)只在炮點(diǎn)附近一定范圍內(nèi)接收，遠(yuǎn)炮點(diǎn)并沒(méi)有數(shù)據(jù)，所以面對(duì)滾動(dòng)排列觀測(cè)系統(tǒng)傳統(tǒng)編碼方法反演效果很差。

仍以圖4所示的部分Overthrust模型為例，采用中間放炮的方式進(jìn)行觀測(cè)，炮點(diǎn)均勻分布于地表，第一炮位于0m處，炮間距為125m，每炮300道接收，道間距為25m，在模型兩邊排列固定不動(dòng)，炮點(diǎn)滾動(dòng)，共得到100炮兩分量記錄。首先采用傳統(tǒng)震源編碼反演方法(極性編碼)對(duì)合成記錄進(jìn)行反演，反演結(jié)果如圖8a所示。顯然，常規(guī)震源編碼方法反演得到的結(jié)果與真實(shí)模型相差很大。由圖8a可以看出，橫波的反演結(jié)果優(yōu)于縱波，這與本文的實(shí)驗(yàn)條件有關(guān)。傳統(tǒng)震源編碼方法應(yīng)用于排列滾動(dòng)的地震數(shù)據(jù)時(shí)，誤差主要源于排列范圍之外(遠(yuǎn)炮檢距)的數(shù)據(jù)缺失。具體到本文的實(shí)驗(yàn)，采用爆炸震源激發(fā)，橫波分量均來(lái)自于波的類型轉(zhuǎn)換，因此占比較少。遠(yuǎn)炮檢距存在很強(qiáng)的縱波分量，尤其是直達(dá)縱波，這一部分?jǐn)?shù)據(jù)的缺失會(huì)對(duì)縱波反演產(chǎn)生極大的不利影響。結(jié)合圖6給出的反演結(jié)果，可以看出，傳統(tǒng)編碼方法在全排列接收時(shí)可以給出正確的反演結(jié)果，但是在滾動(dòng)排列情況下會(huì)受到嚴(yán)重的干擾。

利用本文的分頻編碼方法對(duì)排列滾動(dòng)的數(shù)據(jù)進(jìn)行反演，結(jié)果如圖8b所示。由圖可見(jiàn)本文方法的反演結(jié)果結(jié)構(gòu)清晰，與真實(shí)模型基本一致，遠(yuǎn)遠(yuǎn)優(yōu)于圖8a的反演結(jié)果。計(jì)算效率方面，本文基于GPU集群環(huán)境利用MPI+CUDA實(shí)現(xiàn)多炮并行反演，采用16張Tesla K80運(yùn)算。反演共進(jìn)行75次迭代，分頻編碼方法完成反演用時(shí)70min，極性編碼方法用時(shí)66min，分頻編碼方法耗時(shí)略高。分頻編碼方法首先需要對(duì)記錄做分頻處理；其次，本文在實(shí)施相位靈敏度檢測(cè)時(shí)需要適當(dāng)延長(zhǎng)正演時(shí)間，因此產(chǎn)生了一些額外計(jì)算量。但是這些額外的計(jì)算量相比于波場(chǎng)正演和反傳而言影響較小，所以整體上兩種編碼方法計(jì)算效率相近。

圖8 流動(dòng)排列接收時(shí)不同方法縱波(左)、橫波(右)速度反演結(jié)果

4 結(jié)束語(yǔ)

分頻編碼方法在不明顯降低反演精度的基礎(chǔ)上，能夠?qū)⑦\(yùn)算效率提升數(shù)倍，有效減少正演次數(shù)，減小全波形反演的計(jì)算量。本文的分頻編碼方法首先利用諧波震源實(shí)現(xiàn)分頻的多炮同步正演，再利用基于相位靈敏度檢測(cè)的波場(chǎng)分離方法，實(shí)現(xiàn)不同炮波場(chǎng)的完全分離，從而能夠在梯度計(jì)算中避免引入多炮串?dāng)_噪聲。另外，本文使用的這種分離方法僅依靠頻率信息，與震源、接收點(diǎn)的排列布置無(wú)關(guān)，能夠適用于滾動(dòng)排列觀測(cè)系統(tǒng)，彌補(bǔ)了傳統(tǒng)震源編碼方法在觀測(cè)系統(tǒng)適應(yīng)性上存在的不足。