地球物理学进展  2014, Vol. 29 Issue (6): 2685-2689   PDF    
引用本文
李海英, 武国宁, 陈小民. 2014. 最小二乘约束下的频谱分析技术及其应用[J]. 地球物理学进展, 29(6): 2685-2689, doi: 10.6038/pg20140632  
LI Hai-ying, WU Guo-ning, CHEN Xiao-min. 2014. Least square constrained time frequency method and its applications. Progress in Geophysics, 29(6): 2685-2689, doi: 10.6038/pg20140632   
最小二乘约束下的频谱分析技术及其应用
李海英1, 武国宁2, 陈小民2    
1. 中国石油化工股份有限公司西北油田分公司, 乌鲁木齐 830011;
2. 中国石油大学(北京)理学院, 北京 102249
收稿日期: 2014-07-17 修回日期: 2014-11-10 .
基金项目：中国石油大学(北京)基础科学研究基金(JCXK-2011-12)、科研基金(KYJJ2012-06-13)、中石化项目(KY2012-S-077)和国家自然科学基金(41140033)联合资助.
作者简介：李海英, 男, 1981年生, 汉族, 工程师, 2004年毕业于中国地质大学(武汉), 主要研究领域:油气勘探, 地震资料解释与处理.(E-mail:lhy.1029@163.com)
摘要：时频分析方法是地震数据处理和解释领域的一个重要工具.本文基于最小二乘原理和正则化方法,提出了约束最小二乘时频分析方法.该方法是在有限支撑、非正交基函数空间下,通过迭代整形正则化方法求取信号傅里叶系数的一种时频分析方法.该方法是短时傅里叶变换时频分析方法的延伸和发展,具有较高的时频分析精度,能较准确地识别信号的频谱特征.理论分析和实际资料应用表明,该方法能够较精确的刻画储层的位置和边界,能够进一步提高储层预测的精度.
关键词： 约束最小二乘     正则化     时频分析     储层预测    
Least square constrained time frequency method and its applications
LI Hai-ying1, WU Guo-ning2, CHEN Xiao-min2    
1. China Petroleum & Chemical Corporation, Urumqi 830011, China;
2. College of Science, China University of Petroleum(Beijing), Beijing 102200, China
Abstract: Time-frequency analysis is an important technical tool in seismic data analysis area. In this paper, we proposed the least square constrained method (LSCM) based on the least square and shaping regularization principles. This new method intends for computing the Fourier coefficients of interest signal under finite supported and non-orthogonal base space according to the least square and regularization principles. This is an extension of the short time Fourier transform (STFT). Comparisons with STFT using synthetic signals and real seismic data validates that the proposed method has better time-frequency resolution, can better identify the reservoir's boundaries and locations, further improve the accuracy of reservoir prediction.
Key words: constrained least-square method     regularization     time-frequency analysis     reservoir prediction    
0 引 言

时频分解技术是将一维信号分解为二维时频平面.地震信号由于其非平稳性，地下构造在不同频带范围显现出不同特征.在地球物理学领域，时频分析技术经常被用在地震数据解释和处理领域.时频分解技术是将一维信号分解为二维时频平面(Partyka et al., 1999).地震信号由于其非平稳性，地下构造在不同频带范围显现出不同特征.在地球物理学领域，时频分析技术经常被用在地震数据解释和处理领域.短时傅里叶变换常被应用于河道与沙体的识别及其薄层预测等方面(Partyka et al., 1999)； Chakraborty和OKaya(1995)比较了短时傅里叶变换和小波变换方法，结果表明，小波变换在高频具有较高的频率分辨率，在低频具有较高的时间分辨率； Sinha等(2005)提出了另一种实现连续小波变换的方法，获得了信号的时频特征；Castagna等(2003)将匹配追踪算法应用于储层预测领域；Wang(2007，2010)利用Morlet小波为匹配的原子核，在匹配追踪原理下，对地震信号进行了频谱分析；袁艳华等(2013)实现了最小二乘匹配算法；武国宁等(2012a，2012b)实现了基于复数地震道的匹配追踪算法.概况起来讲，时频分析方法主要有：短时傅里叶变换、小波变换和匹配追踪算法.这些时频分析方法均受窗函数的影响(海森堡原理).如果想提高时间分辨率，不得不牺牲频率分辨率；反之亦然.1 算法原理

要提高时频分析的精度，需要重新考虑傅里叶分析的原理.傅里叶分析是在基函数空间中求信号的傅里叶系数，基函数为正弦和余弦波，在周期内基函数是正交的.所求的傅里叶系数为信号在基函数空间的最佳逼近系数 .即，所求解的傅里叶系数

其中基函数满足：

短时傅里叶变换首先利用窗函数选取部分信号，其它部分补零，然后进行傅里叶分析.短时傅里叶变换受窗函数的影响，往往会产生泄频和拖尾效应(Bracewell，1986).

为了更好的获得信号的瞬时频谱特征，人们在不断探讨新的时频分析方法.Vanìček(1969)提出了通过迭代方法，在最小二乘约束下求取信号的傅里叶系数； Xu等(2005)分析了频谱分析中存在的泄频及其相应的解决方法； Liu等(2013)，刘国昌等(2013)利用正则化方法获得了信号的瞬时时频特征.Puryear等(2012)提出了最小二乘约束下的谱分解方法.本质上说，以上方法都是在不同方法下求取信号的傅里叶系数.

本文的方法，是在前人工作的基础上，求傅里叶系数的另一种途径: 约束最小二乘方法(LSCM).


首先可将正问题表述为

其中F(n，k)=cos(2πkΔfnΔt)+isin(2πkΔfnΔt)，矩阵 F 的第k列(k=1，2，…，N)为频率为kΔf的谐波，列数为频率的取值范围；行数与窗函数具有同样长度.dw为窗函数在t0时刻与信号的乘积，即dw=H(d)*g.这里H(d)表示信号d的希尔伯特变换.注意dw与窗函数具有同样的长度，不需补零.问题进一步表述为

其中Ψk，k=1，2，…，N代表矩阵F 的列向量.该问题的解{k}作为信号在t0时刻的频谱.移动窗函数到下一个时刻t1，继续上述方法，可获得信号在t1时刻的频谱.重复以上步骤可得到信号的时频剖面.

上述问题的求解可表述为一个最小二乘问题：

因为在窗口长度范围内，基函数{Ψk}，k=1，2，…，N不一定正交，同时该方程为一欠定问题：方程的个数小于未知数的个数，所以，方程的解不唯一.需要另外附加一些限定条件来求解该方程.针对欠定方程组，Tikhonov等(1963)提出了正则化方法，上述问题可表述为

这里 D为正则化算子，ε2为参数.方程(5)的解为

其中 F*，D*分别表示F和D的共轭转置.Fomel(2007)在正则化方法的基础上，提出了整形正则化方法.上述问题的解为

其中 S为整形算子，且S-1=ε2 D*D+I.如果S=HH *，则(6)可表示为

下面我们采用上述方法来分析信号的频谱特征. 2 合成信号分析

图 1a为一道人工合成信号，频率分别为20 Hz、50 Hz的单频信号.分别采用短时傅里叶变换和约束最小二乘的方法进行分析，得到图 1b和图 1c.比较可知，约束最小二乘时频分析方法具有较高的时频分析精度，能够较准确地刻画信号的频谱特征.图 2a为Doppler信号，该信号的频率随时间变化而变化.分别采用短时傅里叶变换和约束最小二乘时频分析方法进行分析，得到图 2b和图 2c.比较可知，本文提出的方法具有较高的时频分析精度.能够较好的识别频率中心.

图 1(a)人工合成信号；(b)短时傅里叶变换时频图；(c)约束最小二乘时频图. Fig. 1(a)A synthetic signal;(b)Time-frequency map of STFT;(c)Time-frequency map of LSCM

图 2(a)人工合成信号；(b)短时傅里叶变换时频图；(c)约束最小二乘时频图. Fig. 2(a)A synthetic signal;(b)Time-frequency map of STFT;(c)Time-frequency map of LSCM

图 3为一道人工合成记录时频分析的比较.我们选择了STFT与约束最小二乘时频分析方法进行比较.(a)表示该合成地震道的生成过程，第五道为前四道的叠加；(b)为合成地震记录的STFT时频剖面；(c)为约束最小二乘时频剖面.第一个小波形为一个40 Hz雷克子波.比较时频分布图我们可以看出，约束最小二乘方法具有较高的时间和频率分辨率，STFT时频聚焦性较差.第二个小波形为一个40 Hz雷克子波与一个10 Hz雷克子波的叠加.比较得知：STFT能够识别10 Hz低频部分，对于40 Hz部分识别不是很清楚；约束最小二乘时频分析方法无论在低频还是高频，都能够较为准确地识别出信号的频谱.STFT由于受到窗函数的影响，时频聚焦性较差.第三、四个波形，分别为两个薄层的反射.前者为两个20 Hz雷克子波的叠加，后者为两个20 Hz和30 Hz雷克子波的叠加.对于薄层，STFT如果想获得较高的时间分辨率，不得不牺牲频率分辨率，反之亦然.比较两者得知：约束最小二乘方法的时间分辨率和频率分辨率都好于前者.特别是，该方法在高频能识别出薄层的位置，同时该方法具有较高的频带范围.最后一个波形为三个30 Hz雷克子波的叠加.与STFT相比较，约束最小二乘方法具有较高的频带范围，在高频能识别出薄层的位置.

综合以上讨论我们可以看出，无论是孤立的单层，还是相连的薄层，无论在低频还是在高频，约束最小二乘方法相对于STFT来讲，都取得了较高的精度.能够较为准确的刻画信号的频谱特征.

图 3(a)人工合成单道记录；(b)短时傅里叶变换时频图；(d)约束最小二乘时频图. Fig. 3(a)A synthetic seismic signal;(b)Time-frequency map of STFT;(c)Time-frequency map of LSCM.
3 实际资料分析

地下介质由于其非完全弹性，地震波在传播过程中往往伴随有频散和衰减0-0.在油气聚集区域，往往伴随有高频能量衰减，低频能量相对加强，在地震剖面上将产生亮点.Catagna等(2003)利用匹配追踪算法对地震数据进行分频处理，很好的刻画出储层的位置和边界.Liu等(2013)利用局部谱方法对地震数据进行分频处理，较好地刻画了储层的位置和边界，同时识别出河道和沙体的展布.

为了更好的检验约束最小二乘方法的有效性，我们将此方法应用于储层预测方面. 图 4为某地地震剖面，蓝色圈定的区域为储层的位置. 该储层为含气储层. 抽取井旁道进行分析，见图 5，图 5a为井旁道地震记录，储层大约位于630~635 ms时间深度.比较STFT和约束最小二乘时频分析方法得知：STFT能够识别储层的位置，但是频率分辨率不如约束最小二乘方法.约束最小二乘方法具有较高的频带范围.

图 4 某地地震剖面 Fig. 4 The seismic profile of somewhere

图 5 单道地震记录时频分析对比图(a)单道地震记录；(b)STFT时频图；(c)约束最小二乘方法时频图. Fig. 5 Comparisons of time-frequency maps of single trace(a)Single trace;(b)Time-frequency map of STFT;(c)Time-frequency map of LSCM.

采用约束最小二乘方法对地震数据进行分频处理，得到图 6a~c的时频剖面.图 6a为25 Hz分频剖面，储层位置得到清楚成像，并且储层下方有许多阴影的存在.图 6b为45 Hz分频剖面，储层能够成像，储层下面的阴影有所减弱.图 6c为70 Hz分频剖面，储层及其阴影不再清楚.图 7a~c分别为STFT得到的25 Hz、45 Hz和70 Hz单频剖面.在较低频带范围内，储层能够清楚成像，但下方往往伴随有阴影存在(图 7a)，随着频率的增加阴影逐渐减弱(图 7b)，如果频率进一步提高，储层不再清楚成像(图 7c).对比约束最小二乘时频分析方法与STFT得到：约束最小二乘方法相对于STFT具有较好的时间和频率分辨率，能够较准确的刻画储层的位置(图中箭头所示).

图 6 LSCM单频地震剖面(a)25 Hz；(b)45 Hz；(c)70 Hz. Fig. 6 Iso-frequency slices of CLSM(a)25 Hz;(b)45 Hz;(c)70 Hz.

图 7 STFT单频地震剖面(a)25 Hz；(b)45 Hz；(c)70 Hz. Fig. 7 Iso-frequency slices of STFT(a)25 Hz;(b)45 Hz;(c)70 Hz.
4 小 结

反问题的求解受到非适定因素的影响，具有不稳定性特征.本文在最小二乘约束下，借助于整形正则化算子求取信号的频谱特征.理论和实际资料分析表明：该方法有效地解决了非适定因素的影响，获得了较高的频谱特征，提高了频谱分析的精度.本文采用了高斯函数作为整形光滑算子，不同的整形正则算子对结果的影响需要进一步研究.另外，该算法的速度相对于常规窗口傅里叶变换较慢，如何提高算法的速度也是一个需要研究的内容.

致 谢 感谢中国石油大学(北京)基础学科研究基金项目：NO.JCXK-2011-12；中国石油大学(北京)科研基金项目：NO.KYJJ2012-06-13的联合资助.
参考文献
[1] Bracewell R. 1986. The Fourier Transform and Its Applications[M]. McGraw Hill Publishing Company.
[2] Castagna J P, Sun S, Siegfried R W. 2003. Instantaneous spectral analysis: Detection of low-frequency shadows associated with hydrocarbons[J]. The leading edge, 22(2): 120-127.
[3] Chakraborty A, OKaya D. 1995. Frequency-time decomposition of seismic data using wavelet-based methods[J]. Geophysics, 60(6): 1906-1916.
[4] Chen X H, He Z H, Huang D J, et al. 2009. Low frequency shadow detection of gas reservoirs in time-frequency domain[J]. Chinese J. Geophys. (in Chinese), 52(1): 215-221.
[5] Ebrom D. 2004. The low-frequency gas shadow on seismic sections[J]. The Leading Edge, 23(8): 772.
[6] Fomel S. 2007. Shaping regularization in geophysical estimation problems[J]. Geophysics, 72(2): R29-R36.
[7] Gassmann F. 1951. Vber die Elastizität poröser Medien[J]. Vierteljahrsschrift der Naturforschenden Gesellschaft in Zürich, 96(3): 1-23.
[8] Liu G C, Fomel S, Chen X H. 2011. Time-frequency analysis of seismic data using local attributes[J]. Geophysics, 76(6): P23-P34.
[9] Liu G C, Chen X H, Song J W. 2013. Estimation of primaries and multiples by sparse inversion for OBS data with integration of streamer data[J]. Chinese J. Geophys. (in Chinese), 56(12): 4288-4296, doi: 10.6038/cjg20131231.
[10] Mavko G, Nur A. 1979. Wave propagation in partially saturated rocks[J]. Geophysics, 44(3): 161-178.
[11] Partyka G, Gridley J, Lopez J. 1999. Interpretational applications of spectral decomposition in reservoir characterization[J]. The Leading Edge, 18(3): 353-360.
[12] Puryear C I, Portniaguine O N, Cobos C M, et al. 2012. Constrained least-squares spectral analysis: Application to seismic data[J]. Geophysics, 77(5): V143-V167.
[13] Sinha S, Routh P S, Anno P D, et al. 2005. Spectral decomposition of seismic data with continuous-wavelet transform[J]. Geophysics, 70(6): P19-P25.
[14] Tikhonov A N, Arsenin V Y. 1977. Solitions of Ill-Posed Problems[M]. V. H. Winston and Sons.
[15] Vanìček P. 1969. Approximate spectral analysis by least-squares fit[J]. Astrophysics and Space Science, 4(4): 387-391.
[16] Wang Y H. 2007. Seismic time-frequency spectral decomposition by matching pursuit[J]. Geophysics, 72(1): V13-V20.
[17] Wang Y H. 2010. Multichannel matching pursuit for seismic trace decomposition[J]. Geophysics, 75(4): V61-V66.
[18] Wu G N, Cao S Y, Liu J J. 2012a. Reassigned time-frequency decomposition and its application of hydrocarbon exploration[J]. Progress in Geophys. (in Chinese), 27(2): 0596-0602, doi: 10.6038/j.issn.1004-2903.2012.02.023.
[19] Wu G N, Cao S Y, Sun N. 2012b. Matching pursuit method based on complex seismic traces and its application of hydrocarbon exploration[J]. Chinese J. Geophys. (in Chinese), 55(6): 2027-2034, doi: 10.6038/j.issn.0001-5733.2012.06.024.
[20] Xu S, Zhang Y, Pham D, et al. 2005. Anti-leakage Fourier transform for seismic data regularization[J]. Geophysics, 70(4): V87-V95.
[21] Yuan Y H, Wang Y B, Li Y K, et al. 2013. A non-dyadic curvelet transform based least-squares matching algorithm and its application[J]. Chinese J. Geophys. (in Chinese), 56(4): 1340-1349, doi: 10.6038/cjg20130428.
[22] Zhang F C, Li C H. 2012. Orthogonal time-frequency atom based fast matching pursuit for seismic signal[J]. Chinese J. Geophys. (in Chinese), 55(1): 277-283, doi: 10.6038/j.issn.0001-5733.2012.01.027.
[23] 陈学华, 贺振华, 黄德济,等. 2009. 时频域油气储层低频阴影检测[J]. 地球物理学报, 52(1): 215-221.
[24] 刘国昌, 陈小宏, 宋家文. 2013. 基于稀疏反演的OBS数据多次波压制方法[J]. 地球物理学报, 56(12): 4288-4296, doi: 10.6038/cjg20131231.
[25] 武国宁, 曹思远, 刘建军. 2012a. 谱图重排的谱分解理论及其在储层探测中的应用[J]. 地球物理学进展, 27(2): 596-602, doi: 10.6038/j.issn.1004-2903.2012.02.023.
[26] 武国宁, 曹思远, 孙娜. 2012b. 基于复数地震记录的匹配追踪算法及其在储层预测中的应用[J]. 地球物理学报, 55(6): 2027-2034, doi: 10.6038/j.issn.0001-5733.2012.06.024.
[27] 袁艳华, 王一博, 刘伊克,等. 2013. 基于非二次幂Curvelet变换的最小二乘匹配算法及其应用[J]. 地球物理学报, 56(4): 1340-1349, doi: 10.6038/cjg20130428.
[28] 张繁昌, 李传辉. 2012. 基于正交时频原子的地震信号快速匹配追踪[J]. 地球物理学报, 55(1): 277-283, doi: 10.6038/j.issn.0001-5733.2012.01.027.