本發(fā)明涉及信號處理及模式識別領域，具體地，涉及一種適用于多光柵傳感探測系統(tǒng)的陵區(qū)入侵振動信號特征提取與分類識別算法。

背景技術：

聲波信號的散射衰減是描述地球介質物理化學狀態(tài)的重要參數(shù)，它表征了尾波幅度隨時間的變化趨勢。迄今為止，有關地球物理問題的彈性波散射研究主要集中在大尺度傳播的地震波研究方面，而在井孔聲場的研究中卻鮮有報道。在聲波測井常規(guī)數(shù)據(jù)處理中，散射波和尾波往往被視為噪聲而加以壓制?，F(xiàn)代聲波測井儀器已廣泛應用于油氣勘探，采集了大量包括單極、偶極和正交偶極的陣列聲波測井數(shù)據(jù)。如果能從聲波測井數(shù)據(jù)正確分析和提取彈性波的散射衰減，將會提供很有價值的地層信息，并為聲波測井的處理和解釋開辟新的應用前景。

目前主要有兩種計算散射衰減的模型：單次散射波模型和尾波能量通量模型。單次散射波模型對尾波衰減的分析沒有區(qū)別固有吸收衰減和散射衰減，而是把兩種類型的衰減歸為一項，用“Coda Q”來描述，該模型假定尾波衰減與直達波在介質中傳播的衰減相等，忽略了散射過程中直達波的能量損失對尾波衰減的影響，難于真實反映地層的散射衰減信息(Aki K,Chouet B.Origin of coda waves:Source,attenuation and scattering effects,Journal of Geophysical Research,1975,80:3322–3342)。尾波能量通量模型假定散射衰減和固有吸收衰減對尾波衰變速率的作用不一樣，隨著持續(xù)時間的增加，增加散射衰減將使更多的直達波能量轉換為尾波能量，而增加固有吸收衰減，尾波和直達波均將失去更多的彈性能，該方法在聲波測井數(shù)千赫茲頻率下固有吸收衰減品質因子的靈敏度較低，不適用于井孔聲波測井的散射衰減計算(Frankel A,Wennerberg L.Energy-flux model of seismic coda:Separation of scattering and intrinsic attenuation,Bulletin of the Seismological Society America,1987,77:1223-1251)。

地層巖石散射衰減既與頻率有關，也與尾波的持續(xù)時間有關，它與巖石的固有衰減共同作用，在反演的過程中需要綜合考量各種因素的影響。因此，本專利提出一種基于偶極聲波測井數(shù)據(jù)的多參數(shù)聯(lián)合反演方法計算巖石的散射衰減。

技術實現(xiàn)要素：

針對現(xiàn)有技術不適用于計算散射衰減的現(xiàn)狀，本發(fā)明提供了一種多參數(shù)聯(lián)合反演的方法，利用偶極聲波測井數(shù)據(jù)的尾波，綜合考慮固有衰減、散射衰減與頻率和時間的關系，為計算聲波測井頻率下散射衰減提供一種實用可行的方法。

本發(fā)明為解決上述技術問題所采用的技術方案是：

一種多參數(shù)聯(lián)合反演計算巖石散射衰減的方法，包括如下步驟：

步驟1：在深度區(qū)間內進行陣列聲波測井，得到深度區(qū)間內偶極陣列波形數(shù)據(jù)；

步驟2：獲得處理深度位置處濾波后的偶極子全波列數(shù)據(jù)v(t)，設定窗函數(shù)g(t)，利用短時傅里葉變換(STFT)計算出實測波形不同時刻的功率譜：

其中，t為時間，為角頻率，f為頻率，A(t,f)函數(shù)給出了以t時為中心的時窗g(t)內信號v(t)的傅里葉振幅譜；

步驟3：構建目標函數(shù)ΔE，它是理論的尾波振幅譜A_C(t,f)和實測的尾波振幅譜A(t,f)的函數(shù)，定義為：

ΔE(t,f,Q_I,Q_S)＝|A_C(t,f)-A(t,f)|² (2)

其中Q_I為固有衰減品質因子，Q_S為散射衰減品質因子；理論的尾波振幅譜A_C(t,f)可由三維非均勻模型的尾波幅度表達式計算得到：

其中，t_d為直達波走時，I_D是直達波振幅A_D在時間窗(t₁,t₂)的積分，表示為：

步驟4：根據(jù)步驟3計算每個接收器滿足最小二乘求解目標函數(shù)ΔE的Q_I和Q_S的全局極小值，將每個深度點的N個接收器計算得到的散射衰減品質因子Q_S的均值作為該深度點的散射衰減品質因子，取倒數(shù)得到該深度點的散射衰減值Atten：

步驟5：重復步驟二、步驟三和步驟四，獲得深度區(qū)間內的散射衰減曲線。

本發(fā)明的有益效果是：

本發(fā)明采用聲波測井的方法來計算巖石的散射衰減，可產(chǎn)生兩方面的效益：其一是為聲波測井技術指出了一個重要的應用方向，促進該技術的應用和推廣；其二是從聲波測井數(shù)據(jù)正確分析和提取彈性波的散射衰減，將會提供很有價值的地層信息，與大尺度地震波散射衰減相比，聲波測井計算的散射衰減可提供地下巖石小尺度的非均勻性信息。

附圖說明

圖1為算法流程圖。

圖2為實測單點偶極聲波數(shù)據(jù)。

圖3為實測數(shù)據(jù)的振幅譜。

圖4為理論計算的最優(yōu)振幅譜。

具體實施方式

下面結合附圖及實施例對本發(fā)明作進一步說明。以下的實施例將有助于本領域的技術人員進一步理解本發(fā)明，但不以任何形式限制本發(fā)明。

如圖1所示，本發(fā)明提供一種多參數(shù)聯(lián)合反演計算巖石散射衰減的方法，包括如下步驟：

步驟1：在深度區(qū)間內進行陣列聲波測井，得到深度區(qū)間內偶極陣列波形數(shù)據(jù)；

步驟2：獲得處理深度位置處濾波后的偶極子全波列數(shù)據(jù)v(t)，設定窗函數(shù)g(t)，利用短時傅里葉變換(STFT)即公式(1)計算出實測波形不同時刻的功率譜：

其中，t為時間，為角頻率，f為頻率，A(t,f)函數(shù)給出了以t時為中心的時窗g(t)內信號v(t)的傅里葉振幅譜；

步驟3：構建目標函數(shù)ΔE，它是理論的尾波振幅譜A_C(t,f)和實測的尾波振幅譜A(t,f)的函數(shù)，定義為：

ΔE(t,f,Q_I,Q_S)＝|A_C(t,f)-A(t,f)|² (2)

其中Q_I為固有衰減品質因子，Q_S為散射衰減品質因子；理論的尾波振幅譜A_C(t,f)可由三維非均勻模型的尾波幅度表達式計算得到：

其中，t_d為直達波走時，I_D是直達波振幅A_D在時間窗(t₁,t₂)的積分，表示為：

步驟4：根據(jù)步驟3計算每個接收器滿足最小二乘求解目標函數(shù)ΔE的Q_I和Q_S的全局極小值，將每個深度點的N個接收器計算得到的散射衰減品質因子Q_S的均值作為該深度點的散射衰減品質因子，取倒數(shù)得到該深度點的散射衰減值Atten：

步驟5：重復步驟二、步驟三和步驟四，獲得深度區(qū)間內的散射衰減曲線。

圖2給出了實測的單點偶極聲波測井數(shù)據(jù)，為了進行不同深度點的對比，首先將波形進行歸一化處理，從圖中可以看出波形序列中1ms-5ms的區(qū)間是直達波，5ms之后為地層非均勻性引起的散射波(尾波)。根據(jù)數(shù)據(jù)處理流程，將圖2中的數(shù)據(jù)進行短時傅里葉變換，得到如圖3所示的實測數(shù)據(jù)振幅譜的變密度顯示，其中橫坐標為時間，縱坐標為頻率。

首先選取時窗1ms～3ms的數(shù)據(jù)對直達波進行估計，通過改變Q_I和Q_S使得目標函數(shù)ΔE最小，即可得到在與如圖4所示實測尾波振幅譜A(t,f)最為接近的理論時頻譜，此時它對應的固有吸收衰減品質因子Q_I＝65，散射衰減品質因子Q_S＝262時，從圖中可以看到它與實測的時頻圖極為相似。

說明書中未闡述的部分均為現(xiàn)有技術或公知常識。本實施例僅用于說明該發(fā)明，而不用于限制本發(fā)明的范圍，本領域技術人員對于本發(fā)明所做的等價置換等修改均認為是落入該發(fā)明權利要求書所保護范圍內。