注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

notes

@-@

 
 
 

日志

 
 
 
 

关于仪器响应  

2009-11-06 18:44:26|  分类: 数据 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
zz from http://hi.baidu.com/%C4%AA%BB%F4%BD%E7%C3%E6/blog/item/d6ee5f3917fc932997ddd8c0.html
文一
来源于这里:http://www.eas.slu.edu/Earthquake_Center/CPS/TUTORIAL/RESPONSE/index.html
所需工具:rdseed, evalresp,sacgsac

rdseed主要是用来解压缩SEED格式文件的,SEED格式很好,它既包含了波形数据,有包含了一系列台站和仪器的参数。顺便说一句,IRIS和中国数字台网的都好用,但某些省局的SEED格式不那么标准,似乎没仪器响应。用下面这个语句解压缩:
rdseed -f SEED_VOLUME -R -d -o 1
嗯,就可以生成SAC格式的波形文件和RESP格式的仪器响应文件。
如果希望使用evalresp来去除仪器响应,那么rdseed就这么用:
rdseed -f PBMOout.seed -R
得到的RESP文件名称类似于这样:RESP.NETWORK.STATION.LOCATION.COMPONENT
RESP文件中包含了每个阶段的完整响应参数,包括地震动,摆,数字化仪,和所有的数字滤波。
evalresp这么用:
evalresp STATION_NAME COMPONENT_NAME YEAR DAY_OF_YEAR FMIN FMAX NFREQ -u 'vel' -f RESPFILE
STATION_NAME-台站名
COMPONENT_NAME-分量名
YEAR-年
DAY_OF_YEAR-某日是当年的第几天
FMIN-最小频率
FMAX-最大频率
NFREQ-FMIN和FMAX之间的频率点数
u-指定地震计的类型:如果为'dis'则为位移型,'vel'速度型,'acc'加速度型
f-指定仪器响应文件

比如:
evalresp PBMO BHZ 2007 292 0.001 10 500 -u 'vel' -f RESP.NM.PBMO..BH
会得到两个文件:
AMP.NM.PBMO..BHZ
PHASE.NM.PBMO..BHZ
(这个貌似是该台站仪器响应的振幅和相位,没用过)

SAC的零极点文件
rdseed这样使用就可以得到:
rdseed -f PBMOout.seed -p
会得到类似于这样的文件:
SAC_PZs_NM_PBMO_BHE__2007.173.00.00.00.0000_99999.9999.24.60.60.99999
SAC_PZs_NM_PBMO_BHN__2007.173.00.00.00.0000_99999.9999.24.60.60.99999
SAC_PZs_NM_PBMO_BHZ__2007.173.00.00.00.0000_99999.9999.24.60.60.99999
SAC_PZs_NM_PBMO_HHE__2007.173.00.00.00.0000_99999.9999.24.60.60.99999
SAC_PZs_NM_PBMO_HHN__2007.173.00.00.00.0000_99999.9999.24.60.60.99999
SAC_PZs_NM_PBMO_HHZ__2007.173.00.00.00.0000_99999.9999.24.60.60.99999
SAC_PZs_NM_PBMO_LHE__2007.173.00.00.00.0000_99999.9999.24.60.60.99999
SAC_PZs_NM_PBMO_LHN__2007.173.00.00.00.0000_99999.9999.24.60.60.99999
SAC_PZs_NM_PBMO_LHZ__2007.173.00.00.00.0000_99999.9999.24.60.60.99999

它提供地动位移的传递函数,例如某个垂直向分量的仪器响应文件时这样:
ZEROS 5
867.0800 904.7790
867.0800 -904.7790
POLES 4
-0.1480 0.1480
-0.1480 -0.1480
-314.1590 202.3190
-314.1590 -202.3190
CONSTANT 7.028933e+07
如果使用的是sac或gsac,那么该文件还可以用*来注释某些行,如
* ****
* STATION (KSTNM ): PBMO
* COMPONENT (KCMPNM): HHE
* LOCATION (KHOLE ):
* NETWORK (KNETWK): NM
* START : 2007,173,00:00:00.0000
* END : No Ending Time
* INPUT : METERS
* OUTPUT : COUNTS
* ****
ZEROS 5
867.0800 904.7790
867.0800 -904.7790
POLES 4
-0.1480 0.1480
-0.1480 -0.1480
-314.1590 202.3190
-314.1590 -202.3190
CONSTANT 7.101501e+07
如何查看仪器响应(gsac为例):
gsac << EOF
#####
# create an impulse with unit area
#####
fg impulse delta 0.05 npts 8192
w imp.sac
r imp.sac
ch KSTNM PBMO
ch KCMPNM BHZ
wh
#####
# obtain the velocity sensitivity
#####
transfer from none to eval subtype AMP.NM.PBMO..BHZ PHASE.NM.PBMO..BHZ
w velcount.sac
fft
bg plt
color list red
psp fmin 0.01 fmax 20
plotnps -F7 -W10 -EPS -K < P001.PLT > respplot.eps
convert -trim respplot.eps respplt.png


#####
# obtain the displacement sensitivity
#####
r imp.sac
transfer from none to polezero subtype SAC_PZs_NM_PBMO_BHZ__2007.173.00.00.00.0000_99999.9999.24.60.60.99999
w discount.sac
fft
color list blue
psp fmin 0.01 fmax 20
plotnps -F7 -W10 -EPS -K < P002.PLT > pzplot.eps
convert -trim pzplot.eps pzplt.png
quit
EOF

会得到这样的图像:
关于仪器响应 - xpku - notes
evalresp -u ‘vel’ 输出结果,带FIR滤波器
关于仪器响应 - xpku - notes
零极点去除后的仪器响应,位移谱,不带FIR滤波器


如何去除仪器响应:
去除仪器响应要很小心,因为真实的波形数据中包含了噪声,如果去除不当,就会造成信噪比下降,安全的做法就是在反褶积的过程中使用带通滤波器。
例如:
r imp.sac
transfer from none to eval subtype AMP.NM.PBMO..BHZ PHASE.NM.PBMO..BHZ
w velcount.sac
传递函数把高频和低频端都给去掉了,像图1中显示的那样。
如果有噪声进入了波形记录,而这噪声又和仪器无关,那么冒然进行去除仪器响应会使得噪声增强,一般采取下面的过程:
#####
# define the frequency limits for deconvolution
#####
DELTA=`saclhdr -DELTA velcount.sac`
FHH=`echo $DELTA | awk '{print 0.50/$1}' `
FHL=`echo $DELTA | awk '{print 0.25/$1}' `


#####
# now try a deconvolution with the FREQLIMITS
#####
gsac << EOF
r velcount.sac
transfer from eval subtype AMP.NM.PBMO..BHZ PHASE.NM.PBMO..BHZ to none freqlimits 0.005 0.01 ${FHL} ${FHH}
w deconvelfl.sac
quit
EOF
==================================================================================
文二
来自于这里:http://xb9he.bokee.com/6808495.html,某位地震学的大师兄,汉字的哦,好好向人家学学。

用sac去仪器响应

如果观测记录是u(t),仪器响应应该包含三部分I(t)*A*B,其中I(t)为零点极点,A为放大倍数也即是从seed文件中活得的零极点文件中constant一项,B为仪器的sensitivity.

用语法:

transfer from polezero s pz.file to none/vel/acc freq f1 f2 f3 f4

是去掉I(t)*A,其中pz.file可以从seed里面提取,文件中包含了零点极点和放大倍数A(constant).
相当于u(t)/(I(t)*A). 这之后在除上sensitivity B值,就彻底转换为了位移,速度或者加速度了。

transfer语法中:

s为subtype的简写;
pz.file为零点极点和放大倍数文件
none:位移
vel:速度
acc:加速度
none/vel/acc:每次只选择一个
freq f1 f2 f3 f4
-------------
/                \
-----------                  --------------
f1 f2           f3 f4
表示带通滤波,一般可以不在transfer这里用,而是等你去完仪器响应以后再单独进行滤波,想怎么滤波就怎么滤波。

如果地震仪为速度记录,可直接to vel,当然也可以 to none之后再微分得到vel(sac中提供微分函数dif);或者to acc之后再积分(int)成vel,效果都是一样的。

但是如果后面还有进一步去掉sensitivity, sensitivity就要考虑是速度的还是加速度的,如果是速度的
sensitvity直接先to vel,然后再除以sensitivity的值,就得到m/s;如果是加速度的
sensitvity直接先to acc,然后再除以sensitivity的值,就得到m/s^2;

但其实sensitivity其实也只是幅度影响,至于是速度,位移,加速度,都无关紧要的,但是要注意其单位,别搞错了:
-------------------------------------------------------------------------------------------

下面说说sensitivity:

用u(t)/(I(t)*A)得到的单位为count,

速度仪器的sensitivity的单位是count/kine=count/(cm/s), or: count/(m/s)

u(t)/(I(t)*A*B)最后得到的就是速度m/s, or cm/s

**********************************************

而加速度仪器的sensitivity单位是volt/g=volt/(980.65cm/s^2)

*******************************************************************

在特别强调一点,上面的用sac的transfer去掉仪器响应,其零点极点文件一定是用于to disp的,一般

从seed解压出来的都是默认的为to disp的零点极点,如果其是to vel的零点极点文件,即是用的:

transfer from polezero s PZ.file to none

其得到的也是速度地震图。切记切记!,其中to disp和to vel的零点极点文件很相识,只是零点个数不一样,to disp的比to vel的多一个。

下面这两种方式都可以得到有原始的速度记录去掉仪器响应得到速度地震图:

I have compared the transferred waveform by following two ways, they are equal.

The first way is:

I did not change the zero number, and use the command:

sac> transfer from polezero s PZ.file1 to vel

I can get the velocity seismogram directly from the original velocity record.

The second way is:

I change the zero number (eg. 3-->2), and use the command:

sac> transfer from polezero s PZ.file2 to none

I also can get the velocity seismogram from the original velocity record.

Both of the two seismograms are the same.

文件PZ.file1中零点个数比PZ.file2多一个。PZ.file1就是默认的to disp的,当然如果你想得到位移记录,可以使用:

sac> transfer from polezero s PZ.file1 to none

有点容易混淆哦。

----------------------------------------------------------------------------------------------

最后再强调一个问题是零点极点的单位,SAC默认的是弧度(rad/sec)   ;

如果不是从seed直接解压出来的文件,你得注意有可能是Hz单位,要先转换成*2*pi,转换一下再用

SAC去仪器响应。直接从seed文件解压的不会有问题。

/////////////////////////////////////////////////////////////////////////////////////////////////

如果你关注一下从seed解压出来的response文件,其中的A0 normalization factor *

gain=constant,这个constant就是零点极点文件中的constant. 增益gain,一般会有几个,乘积起来

就可以了。在response文件中的最后一项就是sensitivity。

  评论这张
 
阅读(3717)| 评论(0)
推荐 转载

历史上的今天

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2018