最小平方约束反演谱分析方法的应用效果分析

刘炳杨，李胜军，高建虎，刘军迎

(1.中国石油天然气股份有限公司勘探开发研究院西北分院，甘肃兰州730020;2.中国石油天然气集团公司油藏描述重点实验室，甘肃兰州730020)

解析信号的时频特性分析研究最早可追溯到1946年Gabor[1]提出的Gabor展开。现今常用的时频分析方法大体可以分为3类。第1类是线性变换方法，包括短时傅里叶变换[2]、小波变换[3]、S变换[4]和广义S变换[5-6]等，这类方法采用线性变换形式描述信号的频谱随时间的变化规律，算法较为简单且计算效率较高，但也存在一定问题，如短时傅里叶变换方法用固定的时窗取一段信号进行分析，无法调节时间和频率分辨率；小波变换方法用小波库对信号进行分解，具有多分辨率特点，但其时间尺度的概念与频率关系不直接；S变换以及广义S变换引入了滑动高斯窗函数，其时窗大小随频率而变，能够在低频端获得较高的频率分辨率而在高频端获得较高的时间分辨率，但其窗函数形式固定，一定程度上限制了它的发展。第2类常用的时频分析方法是二次型时频变换法，如Wigner分布。这种方法具有较高的时频分辨率，但美中不足是信号的时频分量间存在交叉项，干扰时频图中的有效信号。第3类时频分析方法是近年来比较热门的贪婪算法，包括匹配追踪(MP)[7]以及Hilbert-Huang变换(HHT)[8]方法，这两种方法都是通过搜索从信号中减去一个最优的基信号，不断循环直到满足停止条件，这类方法能够获取稀疏的高分辨率时频谱，但计算量较大，且都是在单个循环中寻求最优解而不是全局最优解。可以说，每一种方法都有其优点和缺点，难以满足所有需求。近年来，随着分频地震解释[9]、分频AVO[10-11]以及吸收衰减[12-13]等其它基于地震数据体频率变化的储层分析检测方法的发展，对高精度时频分析算法的需求越来越迫切。

传统反演问题的求解方法是构造误差函数的某一范数(通常是2范数)并使其为最小，在这一条件下对未知量进行迭代求解(目前常用Tikhonov回归方法)。此外，为了使求解过程更加准确、稳定，求解过程中一般还要引入其它约束条件，如L1范数、正则约束等。

最近，发展了一类新兴的、基于反演的谱分析方法。这类方法将谱分解问题描述成一个线性反问题，并通过加入其它约束条件来对频谱进行迭代求解，可以获得高时频分辨率的谱。韩利[14]介绍了一种以L1范数为约束条件的稀疏反演谱分解方法，能够获得信号的稀疏解，得到高分辨率频谱，但其稳定性方面则略有欠缺。Puryear等[15]借鉴传统反演的思想，提出了一种最小平方约束谱分析(CLSSA)方法，这种方法实际上可以说是Portniaguine等[16]所提出的反演谱分解(Inverse Spectra Decomposition)方法中的一种。我们对最小平方约束反演谱分析(Constrained Least-Square Spectra Analysis，CLSSA)[15]这一高精度时频分析方法的应用效果进行了分析研究，发现该方法能够准确地拾取小时窗范围内地震信号的频率。应用最小平方约束反演谱分析方法对国内某油田实际地震数据进行频谱分析的结果表明，其应用效果较传统分频方法有明显改善，更为准确、清晰地刻画了研究区内目标层位发育的辫状河道。

1 基本原理

假设时窗地震信号为d，对应的频谱为m，F为傅氏反变换矩阵，可得到形如(1)式的正问题：

(1)

式中：d可以是实地震信号，也可以是由原时窗信号及其Hilbert变换得到的虚地震信号共同形成的复地震信号。试验证明，用复信号进行计算效果更好一些。

为了保证解的唯一性以及稳定性，Puryear等[15]引入了模型权重Wm以及数据权重Wd。Wm初始为单位矩阵，随迭代求解过程而改变；而对于窗函数Wd，我们选用的是Hanning窗。

(2)

可以推导得到形如(3)式的以L2范数进行约束的成本函数

(3)

(3)式的最小范数解为

(4)

这里，α的取值为

(5)

式中：αF为一常系数，一般为0.010～0.001，用以保证算法的稳定性。那么很容易得到所求频谱为

(6)

为了得到更加收敛的频谱，可以进行多次迭代求解，只需更新模型权重Wm即可，其更新公式为

(7)

采用本方法对频谱进行求解，可以直接利用(4)式和(6)式一次计算求得，还可以引入(7)式进行迭代求解。迭代次数越多，获取的频谱越发的收敛，甚至可以获得脉冲状的频谱。但这样也会带来更大的计算量，需要依实际需要进行均衡选择。原则上，当选取信号的时窗长度较小时，可以适当增加迭代次数，以获取更高的频率分辨率。

2 模型试算分析

为了分析CLSSA方法进行谱分解的效果，分别对单个子波、混叠余弦信号进行了试算，并同短时傅里叶变换(STFT)以及S变换方法进行了对比。

2.1 单子波试算分析

取一30Hz主频的Ricker理论子波进行试算，并与傅里叶方法进行对比(子波信号采样间隔为1ms)。以子波峰值点为中心，用Hanning窗分别取40，20ms时窗数据(图1a)进行频谱试算，用来比较不同方法在小时窗情况下的频率分辨率，结果如图1b所示。其中，用CLSSA方法进行频谱分析时均采用了一次迭代。为了便于对比，对各频谱做了归一化处理，以便叠合在一起显示。图1b 中黑线是对整个子波采用傅氏变换得到的频谱，假定为真实频谱；绿线为STFT方法的分析结果，可明显看到其频率分辨率较差，40ms时窗情况下(绿色实线)勉强可以识别出子波主频，但当时窗长度缩小为20ms时(绿色虚线)，主频完全无法识别；蓝线为仅采用原始信号的CLSSA方法分析结果，可见40ms时窗情况下(蓝色实线)所得频谱的频率分辨率相比于STFT明显得到改善，但当时窗缩小为20ms时(蓝色虚线)，依然无法正确识别子波主频；红线为对原始信号进行Hilbert变换得到虚部道，然后与原信号一起合成复数道，再采用CLSSA方法对复数道进行谱分析的结果，可以看到40ms时窗情况下(红色实线)可以得到高频率分辨率的频谱，当时窗缩小为20ms时，依旧能够准确地识别子波的主频，且频率分辨率较40ms时窗的STFT计算结果高得多。可见，相较于STFT方法，采用复数道的CLSSA方法可以用较小的分析时窗得到更为精确的分频结果，在提高时间分辨率的同时还能改善频率分辨率。

图1 子波频谱分析a 子波及计算时窗; b 不同方法频谱对比

2.2 模拟道试算分析

为了进一步对CLSSA方法进行评价分析，设计了一个用不同频率余弦信号混叠而成的模拟记录道。模拟记录道长300ms，1ms间隔采样，其中，0～100ms时窗内为20Hz余弦信号；100～200ms时窗内为50Hz的余弦信号；200～300ms时窗内为(20+50)Hz混叠的余弦信号(图2a)。分别采用STFT，S变换和CLSSA方法进行时频分析，结果分别如图2b，图2c和图2d所示(其中STFT方法和CLSSA方法采用的分析时窗为40ms)。对比可以看出，STFT方法频率分辨率最差，尤其是对于两频率混合的余弦信号，几乎无法分辨；S变换方法获得的时频谱相较于STFT方法有所改善，基本可以分辨混叠信号，但效果不太理想；而用CLSSA方法经过15次迭代后得到的时频谱效果好得多，可清晰、准确地识别出包括混合信号在内的各时窗内信号的频率，分频效果理想。

地震记录通常用地震子波与反射系数的褶积来描述，为了分析CLSSA方法对地震信号的时频分辨效果，对一个由单子波构成的信号道进行了试算分析。子波选用的是30Hz主频的Ricker子波，记录道时长200ms，1ms间隔采样，如图3a所示。分别用STFT，S变换和CLSSA方法进行时频分析(STFT方法和CLSSA方法的分析时窗长度为20ms)，结果分别如图3b，图3c和图3d所示。分析可见，STFT方法得到的分频结果时间分辨率还可以，但频率分辨率较差；S变换方法得到的分频结果的频率分辨率较STFT方法有所提高，但时间分辨率变得很差，尤其是低频段；而采用CLSSA方法所得的分频结果，无论是时间分辨率还是频率分辨率要远远优于STFT方法和S变换方法，分频效果较为理想。

图2 合成余弦信号时频分析a 合成余弦信号; b STFT时频切片; c S变换时频切片; d CLSSA时频切片

图3 单子波道时频分析a 单子波道; b STFT方法时频切片; c S变换时频切片; d CLSSA时频切片

3 实际资料应用效果分析

通过模型数据分析结果可知，CLSSA方法具有十分理想的时频分辨率，因此该方法获得了广泛的应用[17-18]。我们分别采用STFT，S变换和CLSSA方法对国内某油田实际地震资料进行分频分析，提取最大振幅属性以及主频属性。

3.1 单频剖面分析

从研究区内切出201线×201道的地震数据(道间距为20m)进行研究(数据采样间隔为2ms)。图4a为原始剖面，图中黄线标示的是本次研究的目标层位。通过频谱分析，目标层段的主频约为15Hz。因此，分别采用STFT，S变换和CLSSA方法抽取15Hz单频数据体作对比分析，结果分别如图4b，图4c和图4d所示(STFT和CLSSA方法的分析时窗为20ms)。对比可见，STFT方法和CLSSA方法提取的单频切片都具有较好的时间分辨率，相比之下S变换方法提取的单频剖面的时间分辨率显得很差；STFT方法抽取的单频剖面，其能量轴的强弱关系与原始地震记录剖面完全一致，而CLSSA方法抽取的单频剖面则更加突出了高能量异常体。

另外，图4a箭头所示为河道发育处，通过观察我们发现，沿目标层位发育的河道没有表现出明显的强反射轴，反而在河道边缘表现出弱反射特征。

图4 实际地震数据单频剖面a 原始地震剖面; b STFT方法15Hz振幅剖面; c S变换方法15Hz振幅剖面; d CLSSA方法15Hz振幅剖面

3.2 沿层属性分析

CLSSA方法的主要优势在于能用较小的时窗准确识别地震道某时刻的频率，用来识别地震数据体空间范围内的频率变化有望取得很好的效果。

首先沿图4所示研究区目的层开40ms时窗，分别采用STFT，S变换和CLSSA方法对时窗数据进行计算得到频谱，提取最大振幅值及其对应的峰值频率(主频)，结果如图5和图6所示。

图5 沿层最大振幅属性切片(40ms时窗)a 原始数据切片; b STFT方法最大振幅切片; c S变换方法最大振幅切片; d CLSSA方法最大振幅切片

图6 沿层主频属性切片(40ms时窗)a STFT方法主频切片; b S变换方法主频切片; c CLSSA方法主频切片

研究区内发育河道砂体，从图5a原始地震数据沿层切片中隐约可见。图5b为40ms时窗长度下，采用STFT方法提取的最大振幅值沿层切片，河道隐约可见，但展布特征不清晰；图5c为采用S变换方法提取的最大振幅属性，可见，对河道的识别效果较STFT方法要好得多；图5d为采用CLSSA方法提取的沿层最大振幅切片，同样可以清晰地展现研究区内的河道特征，并且同S变换提取的切片相比，CLSSA方法提取的切片中河道的低振幅特性更加突显。

当河道砂体中饱含流体时，会导致地震波主频向低频端衰减，特别是含气情况下，这种衰减现象变得尤为明显。这是利用地震数据频率变化属性进行储层及含油气性检测的理论基础。我们提取了主频属性对研究区进行分析。图6为图5中用不同方法求取的最大振幅值处所对应的峰值频率，即主频切片。从图6中可见，河道发育处地震波主频明显降低。图6a为STFT方法提取的主频切片，河道特征并不明显，很大一部分淹没在背景信息里；图6b是S变换方法提取的主频切片，相较于STFT方法的主频切片，下面的“人”字形河道以及左上部的曲流河特征都得到了较为清晰的展现，但同时也存在局部特征不连续现象；图6c为CLSSA方法提取的主频切片，河道特征更加清晰连续，尤其是红框所示区域，相较于前两种方法，雕刻效果有较大改善。

进一步缩小时窗，对比STFT方法和CLSSA方法的效果。图7为20ms时窗分析结果。STFT方法的最大振幅切片(图7a)中，仅隐约可见河道，而主频切片(图7b)出现大片零频率区域，说明在20ms时窗下，STFT方法已经无法准确地进行频率识别。而相比之下，CLSSA方法拾取的最大振幅切片(图7c)中河道特征依然很显著，同40ms时窗(图5d)相比没有大的变化；主频切片(图7d)中，对河道的刻画效果比图6c更好，分辨率有所提高，雕刻更加细微。

将分析时窗长度缩小到10ms，STFT方法以及CLSSA方法所得分频结果如图8所示。可以看到，STFT方法已经完全不能进行频率识别，主频切片(图8b)中绝大部分区域都成了零频率，而少数非零频率也拾取得完全不准确，看起来更像是奇异值点。而CLSSA方法提取的最大振幅(图8c) 和主频(图8d)依然能够准确拾取频率特征，清晰刻画了河道特征。

实际资料应用结果表明，与STFT和S变换相比，CLSSA方法能够更为准确地拾取地震记录中的频率变化，对研究区内目标层位引起局部主频下降的河道砂体进行精细雕刻。此外，后续减小时窗后STFT和CLSSA方法的应用效果对比也表明，CLSSA方法对时窗长度要求较低，能够用一个很小的时窗准确拾取信号频率，也就是说，CLSSA方法在用小时窗保证时间分辨率的同时，还能够获得很好的频率分辨率。

图7 20 ms时窗沿层属性a STFT最大振幅切片; b STFT主频切片; c CLSSA最大振幅切片; d CLSSA主频切片

图8 10ms时窗沿层属性a STFT最大振幅切片; b STFT主频切片; c CLSSA最大振幅切片; d CLSSA主频切片

4 结束语

最小平方约束谱分析(CLSSA)方法采用反演的思路对信号频谱进行迭代求解，能够获得信号的更为精准的频谱。模型数据试算结果表明，CLSSA方法相较于传统时频分析方法，能够同时获得较高的时间分辨率和频率分辨率。实际资料的应用效果也表明，CLSSA方法用很小的时窗就能够十分准确地进行主频识别，能够对研究区内引起频率变化的地质异常体(河道)进行准确识别和细致刻画。鉴于该方法对分析时窗长度要求较低，建议应用时选取较小的时窗进行分析，这样在保证高频率分辨率的同时，还能够获得更高的时间分辨率。

参 考 文 献

[1] Gabor D.Theory of communication,part 1:the analysis of information[J].Electrical Engineers-Part III:Radio and Communication Engineering,1946,93(26):429-441

[2] Nawab S H,Quatieri T F.Short-time fourier transform[J].Advanced Topics in Signal Processing,1988,6(2):289-337

[3] Daubechies I.The wavelet transform,time-frequency localization and signal analysis[J].Information Theory,IEEE Transactions on,1990,36(5):961-1005

[4] Stockwell R G,Mansinha L,Lowe R P.Localization of the complex spectrum:the S transform[J].Signal Processing,IEEE Transactions on,1996,44(4):998-1001

[5] Adams M D,Kossentini F,Ward R K.Generalized S transform[J].Signal Processing,IEEE Transactions on,2002,50(11):2831-2842

[6] 王长江,杨培杰,罗红梅,等.基于广义 S 变换的时变分频技术[J].石油物探,2013,52(5):489-494

Wang C J,Yang P J,Luo H M,et al.Time variable frequency division based on generalized S transform[J].Geophysical Prospecting for Petroleum,2013,52(5):489-494

[7] Wang Y H.Seismic time-frequency spectral decomposition by matching pursuit[J].Geophysics,2006,72(1):V13-V20

[8] Huang N E,Shen S S.Hilbert-Huang transform and its applications[M].Singapore:World Scientific Publishing Company,2005:1-311

[9] 陈学华,贺振华,黄德济,等.时频域油气储层低频阴影检测[J].地球物理学报,2009,52(1):215-221

Chen X H,He Z H,Huang D J,et a1.Low frequency shadow detection of gas reservoirs in time-frequency domain[J].Chinese Journal of Geophysics,2009,52(1):215-221

[10] 宁媛丽.频散介质基于反演谱分解的AVO方法研究[D].长春：吉林大学,2012

Ning Y L.Study on AVO dispersive medium using inverse spectral decomposition[D].Changchun:Jilin University,2012

[11] 何兵红,吴国忱,郭念民.基于叠前地震道集的 FVO 分析方法[J].石油地球物理勘探,2013，48(1):94-102

He B H,Wu G C，Guo N M.FVO analysis using pre-stack seismic data[J].Oil Geophysical Prospecting,2013,48(1):94-102

[12] 王小杰,印兴耀,吴国忱.基于 S 变换的吸收衰减技术在含气储层预测中的应用研究[J].石油物探,2012,51(1):37-42

Wang X J,Yin X Y,Wu G C.The application of an S transform based absorption and attenuation technique for prediction of gas-bearing reservoir[J].Geophysical Prospecting for Petroleum,2012,51(1):37-42

[13] 刁瑞,李振春,韩文功,等.基于广义 S 变换的吸收衰减分析技术在油气识别中的应用[J].石油物探,2011,50(3):260-265

Diao R,Li Z C,Han W G,et al.Application of absorption and attenuation analysis technique based on generalized S transform for hydrocarbon identification[J].Geophysical Prospecting for Petroleum,2011,50(3):260-265

[14] 韩利.高分辨率全谱分解方法研究[D].长春:吉林大学,2013

Han L.Research on the methods of high-resolution full spectrum decomposition[D].Changchun:Jilin University,2013

[15] Puryear Charles I,Oleg N Portniaguine,Carlos M Cobos,et al.Constrained least-squares spectral analysis:application to seismic data[J].Geophysics,2012,77(5):V143-V167

[16] Portniaguine O,Castagna J.Inverse spectral decomposition[J].Expanded Abstracts of 74thAnnual Internat SEG Mtg,2004,1786-1789

[17] Qi J,Castagna J.Application of a PCA fault-attribute and spectral decomposition in Barnett Shale fault detection[J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013,1421-1425

[18] Oyem A,Castagna J.Layer thickness estimation from the frequency spectrum of seismic reflection data[J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013,1451-1455