基于變分模態分解的彈性參數核密度估計方法

朱鑫杰, 張宏兵, 曾繁鑫, 祝新益

(河海大學地球科學與工程學院, 南京 211100)

隨著近年來油氣資源的開采,陸上油氣資源逐漸變少,油氣開發逐漸向深層油氣藏、薄互層巖性油氣藏等領域發展[1]。由于對儲層的勘探刻畫提出了更加精細化的要求,提高地震彈性參數分辨率以及提高地震模型的精度就相應地成為了地震勘探工作中的重點和難點。

目前常見的高分辨率地震資料處理技術有:反褶積方法、反Q濾波方法、譜白化方法以及基于時頻技術的恢復方法等[2-5]。其中,時頻恢復技術通過分析信號自身特性,利用特定頻帶的信息來描述儲層的相關性質。隨著小波變換[6]以及廣義S變換[7]的使用,時頻分析方法得到了進一步發展。

時頻恢復技術中模態分解的精度對處理效果有著重要影響[8]。小波變換受制于小波函數的尺度變化及Heisenberg不確定原理,存在小波基難以選擇、時間分辨率和頻率分辨率不能兼顧等問題,而廣義S變換受容許性條件和交叉項等內容的限制,其時頻譜分辨率也無法達到最優。Huang等[9]提出的希爾伯特黃變換方法(Hilbert-Huang transform, HHT) 不受Heisenberg不確定性原理的限制,其中的經驗模態分解方法(empirical mode decomposition, EMD)被大量應用于特征提取、故障檢測、信息檢測、去噪降噪等領域[10-12],但EMD的包絡遞歸方法缺乏嚴謹的數學支撐,同時存在模態混疊和端點效應等問題。由Dragomiretskiy等[13]提出的變分模態分解方法區別于EMD,通過構造變分問題并解答的方式,以嚴格的數學推導為支撐,有效解決了模態混疊的問題并且還能有效地避免噪聲影響[14],是一種精度更高的模態分解方法。

建立概率密度模型是隨機模擬中至關重要的環節,而建立的概率密度模型是否忠于數據的真實統計規律,則直接決定了模型以及最后生成的數據的合理性。常規的參數估計方法都需要對數據的分布形式作固定性假設,如路承功等[15]進行可靠性建模時預設數據為韋伯分布;雷慶文等[16]使用極大似然估計時預設降雨量數據以皮爾遜Ⅲ型曲線分布。但實際上,經反演獲得的彈性參數數據是非平穩的時間序列數據,其統計規律并非是某一已知的分布形式,因此,將彈性參數數據限定為某一具體分布形式顯然是不合理的。核密度估計方法是一種非參數密度估計方法,該方法不依賴于密度函數形式的假定,能夠在更大程度上保留數據本身的統計特征,張衛貞等[17]和王黎等[18]的實驗表明以該方法生成的隨機數據能更加符合樣本數據的統計規律;同時針對地質模型本身具有的層狀分布特點,采用區別于總體數據建模的分層建模思想以獲得更合理的概率密度模型。

基于以上內容,現提出一種基于變分模態分解的彈性參數核密度估計方法。首先利用VMD對測井彈性參數進行分解并篩選疊加,得到測井彈性參數中的高頻成分,同時克服模態混疊問題;然后采用核密度估計方法對提取出的測井彈性參數高頻成分進行分層概率密度建模,得出不添加任何先驗條件的概率密度模型;最后以此模型進行隨機抽樣得到符合測井高頻成分概率密度分布的隨機高頻成分并疊加至井旁地震道上,以達到提高地震彈性參數分辨率的目的。

1 基本原理

1.1 VMD原理

Huang等[19]提出的EMD指出,任意的目標信號都可以被分解為一系列本征模態函數。在此基礎上提出的VMD可以將輸入信號分解為K個具有不同中心頻率的有限帶寬信號。與EMD包絡迭代的處理方式不同,VMD采用了變分求解的計算方式,以最終解出的模態的估計帶寬之和最小為約束條件,通過希爾伯特變換和高斯平滑,構建了該模態分解變分問題的模型,即

(1)

式(1)中:K為分解出的模態函數個數;uk為分解出的各模態函數;ωk為各模態函數對應的中心頻率;δ(t)為狄拉克分布函數;t為時間;?t為t時刻的偏導算子;*為卷積;f為原輸入信號。

在構建了變分問題的模型之后,通過引入拉格朗日乘數和二次懲罰項,將約束性變分問題轉化為非約束性變分問題進行求解,得到的增廣拉格朗日函數為

L({uk},{ωk},λ)=

(2)

式(2)中:α為二次懲罰因子;λ(t)為拉格朗日乘法算子。

(3)

(4)

根據式(3)和式(4)可以更新拉格朗日乘數λ,即

(5)

式(5)中:τ為噪聲容忍度。

最后設置了一個判別精度參數ε>0,當達到以下條件時,ADMM停止更新迭代。

(6)

由VMD的計算過程可以看出,VMD具有嚴謹的數學支撐,通過傅里葉變換將頻域內的更新過程轉換到時域,不斷地更新中心頻率的估計值,最后得出所需的K個有限帶寬的IMF。分解出的IMF具有不同的中心頻率,需要對其進行篩選,挑選出中心頻率較高的項進行疊加,得到進行后續工作所需要的樣本高頻成分。

1.2 核密度估計

核密度估計是一種非參數密度估計方法。相較于極大似然估計法等參數密度估計方法,非參數密度估計方法不需要對樣本的分布形式預先做出假設,只利用樣本數據本身進行統計檢驗和判斷分析,核密度函數表示為

(7)

式(7)中:n為樣本個數;h為核密度估計的帶寬;K(u)為核函數;xi為樣本數據;x為自變量。

常用的核函數有均勻核函數、Eqanechnikov核函數和Gaussian核函數等,本文研究使用的Gaus-sian核函數為

(8)

理論上和大量的模擬已經證明,在樣本容量足夠大的情況下,核函數的選取對總體密度的估計并不是至關重要的,對于總體密度函數的任何一個估計,窗寬的選取是核函數密度估計能否成功的關鍵[20]。通常認為使得估計密度函數fn(x)的漸進積分均方誤差達到最小的帶寬為核密度估計的最優帶寬,漸進積分均方誤差表示為

(9)

由于漸近積分均方誤差中的密度函數fn(x)是未知的,所以需要通過其他方式來求得最優帶寬。Sheather等[21]提出的改進解方程插值法(solve-the equation plug-in approach, Sheather and Jones plu-gin, SJPI),通過添加一個預估計的帶寬g(h),使得滿足漸進積分均方誤差最小的最優帶寬hAMISE轉變為對于一個定點方程的解hSJPI,表示為

(10)

Botev等[22]以SJPI為基礎,設計了改進的SJ選擇器(improved Sheather-Jones, ISJ),進一步提高了該方法的性能,基于以上分析,本文研究中采用ISJ選擇器(the improve Sheather and Jones’ selector)獲得核密度估計最優帶寬h。

1.3 方法流程

基于以上原理,提出一種基于VMD的彈性參數核密度估計方法。模型的結構和計算過程如圖1所示,主要包括以下步驟。

圖1 基于VMD的彈性參數核密度估計方法Fig.1 Kernel density estimation method of elastic parameters based on VMD

(1)收集整理研究區內的測井和地震彈性參數數據。

(2)對測井中的彈性參數數據使用VMD方法進行分解,提取出該彈性參數數據中的高頻成分。

(3)將提取出的彈性參數高頻成分依照巖層劃分,依次對每一層的彈性參數高頻成分做核密度估計,確定各層的概率密度分布,構建概率密度模型。

(4) 根據構建的概率密度模型進行隨機抽樣,組成隨機模擬的彈性參數高頻成分。

(5)將模擬出的彈性參數高頻成分疊加到井旁地震彈性參數數據上得到補充了高頻成分的地震彈性參數數據。

2 方法實踐

2.1 彈性參數的VMD分解

以中國南海珠江口盆地某研究區34號井區為研究對象,該研究區位于惠州東南部一構造隆起上,油氣分布連續。通過對區內地層資料分析可知,該區底部為砂巖、粉砂巖夾少量泥巖,中部主要為碳酸鹽巖,上部以泥巖為主,夾有少量砂巖。

本次實驗采用的彈性參數數據為測井縱波波速數據和井旁地震道縱波波速數據,數據詳情如圖2所示,其中測井縱波波速數據采樣間隔為1 ms,地震縱波波速數據采樣間隔為2 ms。

圖2 測井縱波波速數據和地震縱波波速數據Fig.2 P-wave velocity data of well logging and seismic

與EMD輸入信號數據即可得出IMF分量的處理方式不同,使用VMD對信號進行分解時需要預先確定幾個輸入參數,其中最為關鍵的參數為懲罰因子α和輸出的模態個數K。懲罰因子α主要與得到的模態分量的帶寬以及ADMM方法的收斂速度有關。和Dragomiretskiy等[13]在原文中設置的參數值一樣,本文將輸入數據的采樣頻率作為懲罰因子α輸入。另一個重要的參數,模態個數K的選擇對VMD的分解效果非常重要,如果設定的K過小,會使得信號分解得不夠充分,導致模態混疊,即欠分解,如果設定的K過大,則會分解出一些虛假的分量,即過分解。

對于一些已知信號組成成分的情況,K的選取較為簡單,可以根據組成信號的成分個數直接設定。但由于測井彈性參數數據中沒有明確預知的模態信息,所以需要通過其他方法選取合適的模態個數K。采用對比各模態分量中心頻率的方法來選取合適的K,即通過依次代入不同K并統計該K下各模態的中心頻率,分析各模態中心頻率的分布情況以確定最合適的K值。統計結果如表1所示。

表1 不同K下各模態的中心頻率Table 1 Center frequencies of each mode underdifferent K values

從表1可以看出當K達到6時,VMD分解出的模態中已經出現了中心頻率相近的分量IMF2和IMF3,K達到7時過分解情況更加明顯,故選取K=5作為最合適的分解模態個數。

在確定了輸入參數后,對測井縱波波速數據進行VMD分解,所得結果如圖3所示,從圖3可以看出,分解出的5個IMF分量在頻域上圍繞各自中心頻率分布,不同IMF分量間頻帶只有少量重疊,說明VMD分解效果良好,在合理分離出各模態的情況下,有效地避免了模態混疊現象的發生。

圖3 VMD分解出的各IMF及其頻譜Fig.3 The IMF of VMD decomposition and its spectrum

通常的地震資料優勢頻率介于 10～70 Hz[23],實驗里VMD分解出的IMF分量中,IMF1和IMF2兩個分量的頻率較低,均小于70 Hz,將這兩個分量合成作為測井縱波波速數據中的低頻成分,并將該低頻成分與井旁地震道測得的地震縱波波速數據進行對比,如圖4所示。

圖4 VMD分解的低頻成分與地震縱波波速數據對比Fig.4 Comparison of low frequency components of VMD decomposition with seismic P-wave velocity data

可以看出兩組數據除在峰值處存在一定差異,其他位置上的波速值和波動形態幾乎一致,進而認為使用VMD對測井縱波波速數據進行分解后,得到的低頻成分基本可以和井旁地震道所測得的地震縱波波速數據相匹配。其中,兩者在比對時峰值處存在的差異判斷為是由于地震彈性參數反演過程中的偏差,以及地震數據中原本含有的部分高頻成分和噪聲所導致。

2.2 彈性參數高頻成分的核密度估計建模

對VMD分解出的測井縱波波速低頻成分和井旁地震道測得的縱波波速數據進行對比后,認為VMD分解出的另外3個模態分量IMF3、IMF4、IMF5即為測井縱波波速數據中含有而地震縱波波速數據中缺失的高頻成分,將三個分量疊加而得到的縱波波速高頻成分如圖5所示。

圖5 VMD分解出的縱波波速高頻成分Fig.5 High frequency component of P-wave velocity decomposed by VMD

隨機模擬思想中的一個重要內容為獲取樣本后求解樣本的概率密度分布函數,依照概率密度分布函數進行建模,使得生成的預測結果在一定約束下盡可能地貼合實際情況或得出符合統計規律的預測結果。對地層的彈性參數而言,內在地質因素如巖石的礦物成分、結構、孔隙度等,外在環境因素如壓力場、溫度場等[24]都會對其產生影響,這也使得一定范圍內,地質條件相近的地層具有相似的彈性參數分布規律。故對提取出的高頻成分采取合理的分層劃分,將有助于概率密度建模過程的實現。根據已有的地震數據、測井數據以及鉆孔的實際情況,研究地點的地層根據巖性可分為3層,對應到文中的測井縱波波速數據上則是時深1 574～1 604 ms為泥巖層,時深1 604～1 729 ms為灰巖層,時深1 729～1 772 ms為砂巖層。從圖2可以明顯看出,在第2層灰巖層中,縱波速度有兩個明顯的峰值,故根據峰值位置,將灰巖層分為3個小層,對應的時深分別為1 604～1 645 ms、1 645～1 698 ms、1 698～1 729 ms,由此,根據地層巖性以及同層內縱波波速的明顯差異,將彈性參數數據分為了5個小層,再將VMD分解出的高頻成分作為樣本數據,依照相同形式分層后,選擇Gaussian核作為密度核,使用ISJ選擇器自適應地計算出最優帶寬,依次對每一層內的樣本數據進行核密度估計建模,得到每一層縱波波速高頻成分的概率密度曲線如圖6所示。

圖6 每一層的縱波波速高頻成分概率密度曲線Fig.6 The probability density curve of the high frequency component of P-wave velocity for each layer

由概率密度曲線可以看出,不同層的縱波波速高頻成分瞬時幅值樣本分布有顯著的差異。其中,第1層泥巖的瞬時幅值相對較小,主要在-200～200 m/s,分布較為集中;第2層灰巖的第一個高速峰中,瞬時幅值的變化較為平穩,概率密度曲線在形態上近似于一個標準差較大的平闊形正態分布;第3層灰巖中低速層的瞬時幅值分布整體上向負值傾斜;第四層灰巖的第二個高速峰中,瞬時幅值出現了明顯的雙峰形態,并且整體上關于0值對稱;第5層砂巖層的高頻成分瞬時幅值則出現了3個集中分布的區域。從概率密度曲線上來看,核密度估計很好地反映出了每一層中彈性參數高頻成分的具體分布情況,體現出了核密度估計方法對樣本數據進行分析時能盡可能地保留數據原本統計特征的方法優越性。

2.3 地震縱波波速數據的高頻成分補充

得到了每一層中縱波波速高頻成分的概率密度模型之后,根據概率密度模型可以進行隨機抽樣,通過將抽樣結果按時間順序進行組合,就可以生成新的縱波波速高頻成分。根據上文所得到的結果,對VMD分解出的縱波波速高頻成分按照核密度估計的結果分層進行隨機抽樣,組合生成的高頻成分結果如圖7所示,從圖7可以看出,在每一個層內,隨機值和樣本值的波動幅度保持一致,這說明使用核密度估計方法建模后進行隨機抽樣,所得的模擬數據符合樣本數據的統計規律。通過將隨機抽樣生成的縱波波速高頻成分疊加到地震縱波波速數據上,可以得到補充了高頻成分的地震縱波波速數據如圖8所示。

圖7 隨機抽樣生成的縱波波速高頻成分Fig.7 High frequency component of P-wave velocity generated by random sampling

圖8 疊加了高頻成分的地震縱波波速數據對比圖Fig.8 Comparison of seismic P-wave velocity data with superimposed high frequency components

與原始的地震縱波波速數據對比可以看出,疊加隨機生成的高頻成分后,地震縱波波速數據明顯增添了許多細節內容,同時,通過與測井縱波波速數據進行對比,可以看出,補充高頻成分后,地震縱波波速數據出現的細節內容和測井縱波波速數據中對應的部分具有相似的形態。結合去趨勢處理后繪制的頻譜圖(圖9)可以直觀地看出,在常規地震資料優勢頻率10～70 Hz以外的區域,通過高頻成分補充,新生成的地震縱波波速數據擁有了更加豐富的高頻信息,同時也充分保留了原始地震縱波波速數據中的低頻內容。

圖9 疊加高頻成分前后的地震縱波波速數據頻譜對比Fig.9 Spectral comparison of seismic P-wave velocity data before and after superimposing high-frequencycomponents

3 結論

針對常規方法中彈性參數高頻成分提取精度不足、概率密度建模先驗條件過度約束以及彈性參數的概率密度建模分層設計等問題,提出了基于VMD的彈性參數核密度估計方法,利用VMD良好的模態分解性能提取出測井彈性參數數據中的高頻成分,并在此基礎上運用核密度估計方法建立高頻分量的概率密度模型,最后完成隨機抽樣以及對地震彈性參數高頻成分的補充。珠江口盆地34號井區縱波波速數據的實驗結果如下。

(1)使用VMD能有效地分離出測井縱波波速數據中70 Hz以上的高頻成分。各IMF圍繞其中心頻率分布,不同IMF間只有極少量的重疊,即VMD有效地避免了模態混疊現象的發生。

(2)使用核密度估計進行縱波波速樣本的概率密度建模,并根據地層特點進行分層建模,可以使得到的概率密度模型充分地保留樣本數據原本的統計特征。

(3)使用基于VMD的彈性參數核密度估計方法對井旁地震道縱波波速數據進行高頻成分疊加后,頻譜分析顯示在保留了原始低頻信息的基礎上,地震縱波波速數據70 Hz以上的高頻信息得到了明顯的補充。該方法為地震高分辨率隨機建模提供了新的思路。