动态匹配追踪中利用连续相位求取瞬时频率
刘汉卿, 张繁昌, 代荣获, 印兴耀
中国石油大学 地球科学与技术学院,山东 青岛 266555

作者简介: 刘汉卿(1988-),女,中国石油大学(华东)在读硕士研究生,主要研究方向为地球物理反演。

摘要

动态匹配追踪算法将信号的瞬时特征作为先验信息引入到信号分解中,缩小了时频原子库的搜索范围,提高了匹配追踪效率。但在利用瞬时相位获取瞬时频率信息时会得到无物理意义的负频率,不能直接应用于匹配追踪,已有文献并没有给出瞬时频率存在负异常时的解决方案。基于此,笔者通过研究主值相位与解缠相位之间的差异,提出了利用连续相位求取瞬时频率,并推导了连续相位的求取公式;然后通过实例分析,获得了反映三种不同相位数值变化的趋势图;最后通过对比三种相位求取的瞬时频率,证明由连续相位求取的瞬时频率不存在负异常,可直接应用于匹配追踪算法。

关键词: 匹配追踪; 瞬时频率; 负异常; 瞬时相位; 连续相位
中图分类号:P631.4 文献标志码:A 文章编号:1000-8918(2015)01-0211-06
The calculation of the instantaneous frequency using the continuous phase in dynamic matching pursuit algorithm
LIU Han-Qing, ZHANG Fan-Chang, DAI Rong-Huo, YIN Xing-Yao
College of Geoscience and Technology,China University of Petroleum,Qingdao 266555,China
Abstract

The instantaneous characteristics of signals,as the priori information,have been introduced into the signal decomposition which can reduce the scanning range of the atomic library and improve the efficiency of the dynamic matching pursuit algorithm.However,when we use the instantaneous phase to obtain the instantaneous frequency,there will be negative frequency without physical meaning and the negative instantaneous frequency can't be applied to matching pursuit.As far as we know,there is no effective solution to overcome the negative instantaneous frequency in matching pursuit. In view of such a situation,the authors,through studying the difference between the principal phase and the unwrapping phase,propose a method to obtain the instantaneous frequency by introducing the continuous phase.Then the formula of calculating continuous phase can be deduced.In the numerical data tests,the authors present the trend graphs that reflect the changes of the three different phases.Finally,the comparison of the instantaneous frequencies obtained by the three phases proves that the instantaneous frequency obtained by the continuous phase avoids the negative anomalies,and hence it is a good choice for matching pursuit algorithm and can be used directly.

Keyword: matching pursuit; instantaneous frequency; negative anomaly; instantaneous phase; continuous phase

理想的信号分解方式应根据信号特点, 能够自适应选择合适的基函数进行信号分解, Mallat等[1]于1993年提出的匹配追踪正是实现信号自适应稀疏分解[2]算法的一种。匹配追踪算法虽然精度高, 但计算效率低, 限制了它在大型数据处理中的应用前景。Liu等人[3, 4]于2004年提出动态匹配追踪算法, 将信号的瞬时特征[5, 6, 7]引入到匹配追踪中, 缩小了时频原子库在搜索过程中控制参数的扫描范围, 在一定程度上提高了匹配追踪算法的效率。其基本原理是先求取瞬时振幅、瞬时相位[8, 9]、瞬时频率[10, 11, 12], 然后以这些瞬时属性为先验信息进行扫描。虽然瞬时频率是瞬时相位的导数, 但该方法求得的瞬时频率存在负值[13], 不能直接用于匹配追踪分解。

笔者基于瞬时频率存在负值的情况展开了一系列的研究, 发现瞬时频率存在负值是由于瞬时相位异常引起的。由于直接求取的瞬时相位是真实相位被周期折叠后位于区间(-π , π ]的缠绕相位, 存在明显的突变点, 求取瞬时频率时, 不可简单地直接对瞬时相位求导。解缠相位[14, 15, 16, 17, 18]虽然解决了主值相位存在不连续突变点的问题, 相位差小于π , 但它并非单调函数, 通过求导得到的瞬时频率仍存在负值。基于此, 笔者提出了利用连续相位[19]求取瞬时频率的方案, 解决了由主值相位与解缠相位求导获得的瞬时频率存在负值这一问题。

1 相位分析

对瞬时相位进行差分运算能够获得瞬时频率信息, 从而应用于匹配追踪等其他方面的研究, 所以准确的相位是求取瞬时频率的关键。

1.1 主值相位

由于复数信号相位的周期性, 直接求取的瞬时相位位于(-π , π ], 是被周期折叠之后的相位, 该区间的相位值是主值相位, 它与真实相位之间存在着π 的整数倍差异。

1.2 解缠相位

相位谱中含有丰富的信息, 但由于周期性使得相位谱发生不连续突变, 难以提取其中所含的有用信息; 又因为瞬时频率是瞬时相位的微分, 所以主值相位使瞬时频率在不连续点处突变, 且存在负异常现象, 不能直接应用于匹配追踪, 因此, 需要对瞬时相位进行解缠, 恢复被折叠掉的周期[20, 21, 22, 23]

一般认为在解缠相位中, 相邻两个相位之间的差值不应超过半个周期, 即不超过π 。若相邻相位差超过π , 则通过积分进行解缠绕。具体求取过程如下。

ϕ˙ϕ 分别为解缠相位和主值相位, ff-1分别为缠绕算子和解缠算子, 则

ϕ˙=f-1ϕ, (1)ϕ=f(ϕ˙)=ϕ˙+2, -π< ϕ=f(ϕ˙)π,  (2)

其中:n为整数。

根据解缠相位的定义, 可知相邻相位差满足

Δϕ˙(i)=ϕ˙(i+1)-ϕ˙(i), -π< Δϕ˙π(3)

根据式(2)、式(3)得主值相位差分

Δϕi=ϕi+1-ϕi=Δϕ˙i+2πΔni, (4)

由于主值相位的差分存在突变, 故对其进行缠绕运算, 得

fΔϕ(i)]=fΔϕ˙(i)+2πΔn(i)]=Δϕ˙(i)+2πΔn(i)+2πm(i)(5)

由于fϕ (i)]和Δ ϕ˙(i)均位于区间(-π , π ], 故

2πΔn(i)+2πm(i)=0(6)

所以, 式(5)可简化为

Δϕ˙(i)=fΔϕ(i)], (7)

将上式代入式(3)得解缠相位计算公式

ϕ˙(i+1)=ϕ˙(i)+fΔϕ(i)), i2; ϕ˙(1)=ϕ(1)(8)

1.3 连续相位

虽然解缠相位解决了相位不连续问题, 但由于解缠相位只是相邻相位差小于π , 并不确定是否为递增序列, 所以由解缠相位求得的瞬时频率仍可能存在负异常, 不能直接应用于匹配追踪。而连续相位是递增序列, 其相邻相位差为正值, 故通过对连续相位求导得到的瞬时频率不存在负异常, 求取公式如下。

φ 为连续相位, χ 为连续算子, 则

ϕ=fφ=φ+2, (9)φ=χϕ=ϕ+(10)

其中:n, k均为整数。

由式(9)、式(10)得主值相位的差分为

Δϕ(i)=Δφ(i)+2πΔn(11)

由连续相位的定义可知

Δφ(i)=φ(i+1)-φ(i)(0, π], (12)

对式(11)两端取缠绕算子和连续算子后得

χ{fΔϕ(i)]}=χΔφ(i)+2πΔn+2]=  Δφ(i)+2πΔn+2+(13)

由于χ {fϕ (i)]}和 Δ φ (i)于(0, π ], 故 2π Δ n+2mπ +kπ =0, 从而得到连续相位公式

φ(i+1)=χ{fΔϕ(i)]}+φ(i)(14)

2 相位实例分析

为了对比主值相位、解缠相位以及连续相位之间的差异, 分别利用理论信号与实际信号进行分析。

2.1 余弦信号

以周期函数余弦信号y=cos(ω t)为例, 分别求取它的主值相位、解缠相位和连续相位, 并对比分析三种相位之间的异同点(图1)。

图1 余弦信号y=cos(ω t)的瞬时相位
a— 主值相位; b— 解缠相位; c— 连续相位

图1可看出, 余弦信号的主值相位位于(-π , π ]之间, 呈现周期性, 且存在不连续间断点。相对于此, 解缠相位和连续相位的相位差分比较小, 不存在大的间断点。由于余弦信号的特殊性, 其解缠相位和连续相位之间不存在差别。

2.2 离散信号

为说明解缠相位和连续相位之间的差异性, 以离散信号[-0.8π , -0.5π , 0.1π , 0.7π , -0.4π , -0.5π , -0.3π , 0.6π , 0.5π , 0.8π ]为例进行验证。

图2展示了离散信号的三种相位。主值相位仍位于区间(-π , π ], 但在第4和第5个采样点之间存在明显的不连续, 其差分为-1.1π 。而解缠相位避免了这一现象, 相邻相位间的差分比较小, 不超过半个周期, 变化比较平缓, 连续性明显提高。通过图2c可以看出, 连续相位也解决了相位的不连续性, 但它不同于解缠相位, 是一个递增函数。

图2 离散信号的瞬时相位
a— 主值相位; b— 解缠相位; c— 连续相位

2.3 趋势对比图

通过对比分析发现, 主值相位位于(-π , π ]之间, 相邻相位的差分可大于π ; 解缠相位并不是递增的, 只要相邻相位差分位于(-π , π ]之间即可; 而连续相位是递增的, 相位差大于零。图3表示了这三种相位的变化趋势。

图3 三种瞬时相位的变化趋势a— 主值相位; b— 解缠相位; c— 连续相位

3 瞬时频率实例分析

通过Hilbert变换[24]求取的信号瞬时属性并非对任何信号都有物理意义, 此法要求被分析信号是窄带或平稳的, 而且对噪声很敏感。而实际地震信号既非平稳又含有噪声[25], 若直接对地震信号进行Hilbert变换, 求取的瞬时属性将缺乏物理意义甚至失真。在匹配追踪过程中, 利用瞬时相位求取瞬时频率时, 也会出现无物理意义的负频率, 若继续在该瞬时频率附近扫描, 则无法求取匹配子波主频。所以如何利用瞬时相位求取准确的瞬时频率, 是亟需解决的问题之一。

为了对比不同相位获得的瞬时频率, 依次以合成地震信号和实际信号道为例, 分别利用三种相位求取瞬时频率。

3.1 合成地震信号

合成地震信号的各参数见表1, 波形如图4所示, 该地震信号由不同主频、不同相位的子波构成。

表1 合成信号参数

图4 合成地震信号

图5~图7分别给出了合成信号的主值相位、解缠相位、连续相位及其对应的频率。图5a为主值相位, 被缠绕折叠在(-π , π ]之间, 存在多个不连续突变点, 故用主值相位求导获得的瞬时频率(图5b)在突变点处, 瞬时频率呈现较大的负值。所以由主值相位求得的瞬时频率, 不能直接用来当作匹配追踪中时频原子的主频。

相对于主值相位, 解缠相位(图6)不存在明显的不连续突变点, 变化比较平缓, 相邻相位差小于π 。

图5 主值相位及其对应的瞬时频率
a— 主值相位; b— 瞬时频率

图6 解缠相位及其对应的瞬时频率
a— 解缠相位; b— 瞬时频率

图7 连续相位及其对应的瞬时频率
a— 连续相位; b— 瞬时频率

图6可看出, 在中心时间处对应的瞬时频率与表1给出的实际频率相一致, 但在0.3~0.4 s之间出现了负瞬时频率, 这是因为解缠相位并不是递增的, 相位差存在负值, 从而使得相应的瞬时频率存在负异常值。

图7a显示, 连续相位存在递增趋势, 而且由它求得的瞬时频率并未出现负值; 在图7b中, 中心时间处的瞬时频率与表1的实际频率一致, 说明了由连续相位求取的瞬时频率可直接应用于匹配追踪。虽然在0.3~0.4 s处瞬时频率异常, 但在进行匹配追踪时, 是在包络能量最大值处进行扫描, 故通过一次迭代匹配出能量极值处的匹配子波, 并将其从原信号中减去后, 旁瓣消失, 进行下一次迭代时, 在该位置处的频率异常现象就会消失, 故不会追踪出假高频信息。

3.2 实际地震信号

进一步以实际地震数据检验连续相位用于匹配追踪瞬时频率求取的有效性。图8为从某工区中提取的地震道, 分别求取该道信号的主值相位、解缠相位、连续相位及其对应的瞬时频率, 如图9~图11

图8 某工区的地震道

对比三种相位及其对应的瞬时频率可知, 该实际地震信号的主值相位是被周期折叠的, 在不连续突变点处的瞬时频率呈现较大的负值; 而解缠相位虽然恢复了其折叠掉的周期, 变得平缓无突变点, 但是由其求得的瞬时频率仍存在负异常, 不能直接应用于匹配追踪算法; 而用本文的连续相位求得的瞬时频率很好地消除了负值现象, 可直接应用于匹配追踪分解。但从图中可看出, 利用连续相位直接求取的瞬时频率存在一定的问题, 且容易受到噪声影响产生波动, 如图11b中某些点处存在较大的频率值, 虽然这不会影响匹配追踪分解效果, 但在一定程度上限制了连续相位在其他方面的应用。针对这一问题, 期望后期通过研究瞬时频率的求取算法将其进行压制, 从而求得连续相位对应的稳定瞬时频率, 更好地为匹配追踪及其他应用奠定基础。

图9 某道实际地震信号的主值相位及其对应的瞬时频率
a— 主值相位; b— 瞬时频率

图10 某道实际地震信号的解缠相位及其对应的瞬时频率
a— 解缠相位; b— 瞬时频率

图11 某道实际地震信号的连续相位及其对应的瞬时频率
a— 连续相位; b— 瞬时频率

4 结论

针对动态匹配追踪中瞬时属性的求取问题, 先分析了主值相位、解缠相位之间的差异, 并在此基础上得到了连续相位的求取公式; 然后分别求取理论信号和实际信号的三种相位, 分析相位数值的变化, 得到了相应的变化趋势图; 最后, 对于合成地震信号和实际地震道, 分别利用三种相位求取瞬时频率, 证明了由主值相位和解缠相位求得的瞬时频率都存在负值, 不能直接用于匹配追踪等信号处理, 而由于连续相位的递增性, 通过它求得的瞬时频率不存在无物理意义的负频率, 可直接当作匹配追踪中时频原子的主频, 从而为后续的匹配追踪分解等研究奠定基础。

The authors have declared that no competing interests exist.

参考文献
[1] Mallat S, Zhang Z. Matching pursuit with time-frequency dictionaries[J]. IEEE transactions on signal processing, 1993, 41(12): 3397-3415. [本文引用:1]
[2] 张春梅, 尹忠科, 肖明霞. 基于冗余字典的信号超完备表示与稀疏分解[J]. 科学通报, 2006, 51(6): 628-633. [本文引用:1]
[3] Liu J, Wu Y, Han D, et al. Time-frequency decomposition based on Ricker wavelet[C]//Expand ed Abstracts of 74th SEG Annual International Meeting, 2004: 1937-1940. [本文引用:1]
[4] Liu J, Marfurt K J. Matching pursuit decomposition using Morlet wavelets[C]//Expand ed Abstracts of 75thSEG Annual International Meeting, 2005: 786-789. [本文引用:1]
[5] 杨培杰, 印兴耀, 张广智. 希尔伯特—黄变换地震信号时频分析与属性提取[J]. 地球物理学进展, 2007, 22(5): 1585-1590. [本文引用:1]
[6] 钟韬, 肖宏跃, 雷宛. 地震瞬时属性分析技术在岩溶勘查中的应用[J]. 物探化探计算技术, 2008, 30(1): 30-33. [本文引用:1]
[7] 张繁昌, 李传辉, 印兴耀. 基于动态匹配子波库的地震数据快速匹配追踪[J]. 石油地球物理勘探, 2010, 45(5): 667-673. [本文引用:1]
[8] 孙成禹, 尚新民, 石翠翠, . 影响地震数据相位特征的因素分析[J]. 石油物探, 2011, 50(5): 444-454. [本文引用:1]
[9] 廖立坚, 杨新安, 黄小平. 高分辨率探地雷达复数道分析方法[J]. 物探与化探, 2008, 32(3): 301-303. [本文引用:1]
[10] 肖兵, 何继善. 探地雷达信号瞬态谱分析[J]. 物探与化探, 1999, 23(6): 459-461, 466. [本文引用:1]
[11] 胡明顺, 潘冬明, 徐红利, . 几种时频分析方法对比及在煤田地震勘探中的应用[J]. 物探与化探, 2009, 33(6): 691-695. [本文引用:1]
[12] 李传辉, 张繁昌, 印兴耀. 三种线性时频分析方法的影响因素及精度分析[J]. 石油天然气学报, 2010, 32(4): 239-242. [本文引用:1]
[13] 张尔华, 王晓凯, 高静怀, . 局部瞬时频率提取及其在大庆油田中的应用[J]. 地球物理学进展, 2011, 26(3): 1064-1069. [本文引用:1]
[14] 熊涛, 陈亦伦, 杨健, . 基于变相位技术的相位解缠方法[J]. 中国科学: 信息科学, 2010, 20(3): 445-457. [本文引用:1]
[15] 靳国旺, 徐国华, 佘懋勋. 基于瞬时频率估计的 In-SAR 相位解缠方法[J]. 测绘科学与技术学报, 2009, 26(1): 33-35. [本文引用:1]
[16] 程璞, 许才军, 王华. InSAR 相位解缠算法研究[J]. 大地测量与地球动力学, 2007, 27(3): 50-55. [本文引用:1]
[17] Chen C W, Zebker H A. Two-dimensional phase unwrapping with use of statistical models for cost functions in nonlinear optimization[J]. Journal of the Optical Society of America A, 2001, 18: 338-351. [本文引用:1]
[18] Goldstein R M, Zebker H A, Werner C L. Satellite radar interferometry: Two-dimensional phase unwrapping[J]. Radio Science, 1988, 23(4): 713-720. [本文引用:1]
[19] 陈静, 刘波. 连续相位调制解调技术研究[J]. 载人航天, 2014, 20(3): 278-282. [本文引用:1]
[20] Shatilo A P. Seismic phase unwrapping: methods, results, problems[J]. Geophysical Prospecting, 1992, 40: 211-225. [本文引用:1]
[21] Mcgowan R, Kuc R. A direct relation between a signal time series and its unwrapped phase[J]. IEEE Transactions on Acoustics, Speech and Signal Processing, 1982, 30: 719-726. [本文引用:1]
[22] Poggiagliolmi E, Berkhout A J, Boone M M. Phase unwrapping: possibilities and limitations[J]. Geophysical Prospecting, 1982, 30: 281-291. [本文引用:1]
[23] 杨峰涛, 罗江龙, 刘志强, . 相位展开的6种算法比较[J]. 激光技术, 2008, 32(3): 323-326. [本文引用:1]
[24] 张凤旭, 孟令顺, 张凤琴, . 利用 Hilbert 变换计算重力归一化总梯度[J]. 地球物理学报, 2005, 48(3): 704-709. [本文引用:1]
[25] 张繁昌, 李传辉. 非平稳地震信号匹配追踪时频分析[J]. 物探与化探, 2011, 35(4): 546-552. [本文引用:1]