补充二、希尔伯特变换及其应用

1、1、连续时间信号 实连续时间信号 的Hilbert变换 定义为： )( tx)(tx t txd tx d t x tx 1 *)( )(1)(1 ) ( tth/1)( 0 0 )sgn()( j j jjH 补充内容2： 希尔伯特（Hilbert）变换及其使用 Hilbert变换与解析信号 结论：结论：Hilbert变换器是幅频特性为1的全通滤波器。 信号x(t)通过Hilbert变换器后，其负频率成分作 +90相移，而正频率成分作90相移。其幅频、 相频特性为： Hilbert反变换 HilbertHilbert变换与解析信号变换与解析信号 )sgn()()sgn()( )()()( j

2、jXjjX jHjXjX 由时域卷积定理有： t th t txtx 1 )(, 1 )() ( 由此可得：HilbertHilbert反变换的公式反变换的公式 )( )sgn()(jXjjX d t x tx t tx ) ( 1 ) ( * 1 )( 解析信号 ) ( )()(tx jtxtz设 为 的Hilbert变换，定义 为信号 的解析信号解析信号(analytic signal)。 )( tx )(tx )(tx HilbertHilbert变换与解析信号变换与解析信号 )()()()( )()(jXjjHjXjXjjXjZ 0 0 )sgn()( j j jjH又 00 0)(2

3、 )( jX jZ ()X j ()Z j ()X j 结论：结论：由Hilbert变换构成的解析信号，只含有正频率 成分，且是原信号正频率分量的2倍。 例 HilbertHilbert变换与解析信号变换与解析信号 设 ，求其Hilbert变换及解析信号)2cos()( 0t fAtx 解：解： )()( 2 )()( 2 )( )()( 2 )( 00 00 00 A j jj A jX A jX )()( )()( 0 AjXjjXjZ又 tfj Aetz 0 2 )( 解析函数为: )2sin()( 0t fAtx可以证明，若 )2cos()( 0t fAtx则其Hilbert变换 00

4、 2 f令 )2sin() ( 0t fAtx 的Hilbert变换为：)(tx 可实现系统是因果系统，其冲激响应可实现系统是因果系统，其冲激响应 tuthth 00 tth即即: 其傅里叶变换其傅里叶变换 j 1 j 2 1 jHH )j (jjej)j ( j XRHH 又又 则则 jj)j (XR j 1 jjj 2 1 XR 1 jj 2 1 XR 1 jj 2 j RX d j 2 1 j 2 1 jjj X RXR所所以以 d j 2 1 2 j j RX 根据实部与实部相等，虚部与虚部相等，解得根据实部与实部相等，虚部与虚部相等，解得 d j 1 )j ( X R d j 1 j

5、 R X 因果系统系统函数因果系统系统函数 )j ( H的实部与虚部满足希尔的实部与虚部满足希尔 伯特变换约束关系。伯特变换约束关系。 伯特变换的约束关系。伯特变换的约束关系。 的实部与虚部满足希尔的实部与虚部满足希尔，证明，证明已知已知)()()(thFtueth ta a 因为因为 a a a a j tueFthF t 1 )( 即系统函数即系统函数 a a a a a a jjXjRjjH 2222 式中实部式中实部 22 a a a a jR 虚部虚部 22 a a jX a a a a a a a a a a d 2 1 2 1 1 22 jjjj jXH 现在求现在求 jX的希尔

6、伯特变换的希尔伯特变换 d 1 jX jXH a d 1 22 aaa C j B j A 22 令 可求出各分式系数可求出各分式系数 22 , 2 1 , 2 1 a a a a a a C j B j A 则则 a a a a a a a a a a d 2 1 2 1 1 22 jjjj jXH a a a a a a d 1 22 2 22 a a a a a a a a d 1 2222 2 22 a a a a a a a a a a lnlnarctg 1 22 22 00 22 1 22 a a a a 22 a a a a R X HilbertHilbert变换与解析信号变

7、换与解析信号 离散时间信号 0 0 )( j j eH j 设Hilbert变换器的单位抽样响应为h(n)，与连续信 号Hilbert变换器的频率响应H(j)对应，h(n)的频 率响应H(ej)为： 作H(ej)的傅里叶反变换，求得： de jdjedeeHnh njnjnjj 0 0 2 1 2 1 )( 2 1 )( 为奇数 为偶数 n n n n nh n 2 0 ) 1(1 )( )( nx)(nx离散时间信号 的Hilbert变换记为： HilbertHilbert变换与解析信号变换与解析信号 的Hilbert变换 为：)( nx)(nx m m mnx nhnxnx ) 12( )

8、 12(2 )(*)() ( 0)( 0)( )()()( j j jjj ejX ejX eHeXeX 求得 后，即可构成 的解析信号 )( nx ) ( )()(nx jnxnz )(nx 使用DFT求信号的Hilbert变换步骤： Step 1. 对 做 DFT, 得：( )x n( ),0,1X kkN Step 2. 令 1, 2 0 1 2 , 2 , 1)(2 0)( )( N N k N kkX kkX kZ Step 3. 对 做逆 DFT, 得( )Z k( )z n Step 4. 由( ) IDFT( ( )( )x nj Z kX k 得)()() ( nxnzjnx

9、 实连续信号的包络、瞬时相位、瞬时频率 ) ( )()(tx jtxtz )()(| )(|)( 22 txtxtzte )( ) ( arctan)( tx tx t dt td t )( )( X（t）的包络 ,瞬时振幅 X（t）的瞬时相位 X（t）的瞬时频率 实离散信号的包络、瞬时相位、瞬时频率 ) ( )()(nx jnxnz )()(| )(|)( 22 nxnxnzne )( ) ( arctan)( nx nx n X（n）的包络 X（n）的瞬时相位 18 help hilbert HILBERT Discrete-time analytic signal via Hilbert

10、 transform. X = HILBERT(Xr) computes the so- called discrete-time analytic signal X = Xr + i*Xi such that Xi is the Hilbert transform of real vector Xr. If the input Xr is complex, then only the real part is used: Xr=real(Xr). If Xr is a matrix, then HILBERT operates along the columns of Xr. 19 n=0:

11、1:50; dt=0.001; a=0.1; x=exp(-a.*n).*sin(2*pi*0.4375.*n) subplot(2,2,1); plot(x); y=hilbert(x); rx=real(y); Ix=image(y); 20 z=sqrt(rx.*rx+ix.*ix);%求信号x的包络 %z=sqrt(abs(x).2+abs(y).2); subplot(222); plot(z); thet=atan(ix./rx);%求信号x的瞬时相位 subplot(223); plot(thet); real(z) 实部函数 imag(z) 虚部函数 abs(z)求复数z的模

12、angle(z)求复数z的辐角，其范围是（ ， conj(z) 求复数z的共轭复数 0204060 -1 -0.5 0 0.5 1 0204060 0 0.5 1 1.5 0204060 -2 -1 0 1 2 单道地震信号数值模拟单道地震信号数值模拟 1.实际地震记录信号瞬时属性提取 2.合成记录瞬时属性提取 3.分析瞬时属性的分辨率及地质意义 1.实际地震记录信号瞬时属性提取实际地震记录信号瞬时属性提取 三条同相轴，最上一条和中间一条间距小，在瞬时振幅曲线上未分辨出来， 中间一条和最下一条间距相对较远，在瞬时振幅曲线上可以分辨出来，但 也不明显。 600ms 461 700ms 单道地震信

13、号数值模拟单道地震信号数值模拟 2. 合成记录瞬时属性提取合成记录瞬时属性提取 合成地震记录参数：30Hz的Ricker子波，三个反射层位且相距较远 ， 瞬时振幅可以很清晰的分辨出三个层位 -30-20-100102030 -1 0 1 the Ricker wavelet 050100150200250300 -0.5 0 0.5 the reflection coefficient 050100150200250300 -0.5 0 0.5 the synthetic seismograms 050100150200250300 -0.5 0 0.5 time/ms amplitude p

14、art of the original signal 050100150200250300 0 0.2 0.4 0.6 0.8 time/ms amplitude the instantaneous amplitude of the synthetic seismograms 单道地震信号数值模拟单道地震信号数值模拟 单道地震信号数值模拟单道地震信号数值模拟 2. 合成记录瞬时属性提取合成记录瞬时属性提取 合成地震记录参数：100Hz的Ricker子波，四个反射层位，其中第二个和第三个相距较近，合 成地震信号有明显的四个层位，但是Hilbert变换只能清晰的辨别三个层位，这说明Hilbert变

15、换提 取瞬时振幅的分辨率是有限的，特别是对于薄层、薄互层。 -30-20-100102030 -1 0 1 time amplitue the riker wavelet of 100Hz 050100150200250300 -0.5 0 0.5 time amplitude the reflection coefficient 050100150200250300 -0.5 0 0.5 time amp the synthetic seismogram 050100150200250300 -0.5 0 0.5 time amplitude the synthetic seismogram

16、s 50100150200250300 0.2 0.4 0.6 0.8 time amplitude the instantaneous amplitude of the synthetic seismograms 单道地震信号数值模拟单道地震信号数值模拟 3. 瞬时属性的分辨率及地质意义瞬时属性的分辨率及地质意义 通过单道信号的瞬时属性的分析，可知利用瞬时属性可以反映同相轴的 局部或细微变换，但其分辨率也是有限的，而且不同瞬时属性反映的信息也 不同。 属性类别 物理意义 主要地质意义 瞬时振幅 地震反射波强度 的量度 反映岩性、不整合、断层、流体及储集层 的孔隙率变化，识别河流和湖泊相砂体、

17、 局部断层、薄层、缝（洞）等。 瞬时相位 同相轴连续性的 量度 反映岩性、地层层序变化，识别河流和湖 泊相砂体、不整合断层、缝（洞）等。 瞬时频率 相位的时间变化 率 反映岩性、地层层序变化，识别冲积扇和 三角洲砂体、异常衰减和薄层调谐变化、 缝（洞）等。 1. 瞬时地震属性提取方法综述瞬时地震属性提取方法综述 “？”表示尚无文献有相关方面研究 由上框图可知，不论小波分解还是经验模态分解，最后的分量信号依然要 利用Hilbert变换处理提取瞬时属性。 瞬时地震属性提取算法研究瞬时地震属性提取算法研究 2.频率域频率域Hilbert变换算法实现变换算法实现 )()()(1)( )()()()(

18、)()( Hsjhs hjsss jsS 0 0 )( j j h 00 02 )( H )( )()(ts jtstS 瞬时地震属性提取算法研究瞬时地震属性提取算法研究 2.频率域频率域Hilbert变换算法实现变换算法实现 频率域希氏变换流程图频率域希氏变换流程图 瞬时地震属性提取算法研究瞬时地震属性提取算法研究 瞬时地震属性提取算法研究瞬时地震属性提取算法研究 3. Hilbert变换提取瞬时地震属性的算法验证变换提取瞬时地震属性的算法验证 1、单道信号验证、单道信号验证 02004006008001000120014001600 -2 0 2 x 10 4 the original s

19、ignal 02004006008001000120014001600 0 1 2 x 10 4 the intantaneous amplitude calculated by the matlab 02004006008001000120014001600 0 1 2 x 10 4 the instantaneous amplitude calculated by the C program 瞬时地震属性提取算法研究瞬时地震属性提取算法研究 3. Hilbert变换提取瞬时地震属性的算法验证变换提取瞬时地震属性的算法验证 2、瞬时振幅剖面验证、瞬时振幅剖面验证 瞬时地震属性应用瞬时地震属性应用 1. 圈定火成岩侵入煤层范围 2. 天然气水合物BSR空间展布 瞬时地震属性应用瞬时地震属性应用 1. 圈定火成岩侵入煤层范围圈定火成岩侵入煤层范围 克拉美丽气田发现层系属石炭系火山 岩，气藏地质条件复杂，火山岩岩性与分 布的识别是制约火山岩勘探的一个“瓶颈” 难题。此次数值试验取用克拉美丽气田资 料，尝试利用瞬时属性进行火层岩侵入煤 田范围。 瞬时地震属性应用瞬时地震属性应用 1 7 0 4 线叠后剖面线叠后剖面 瞬时地震属性应用 1704线瞬时振幅剖面