from http://hi.baidu.com/ustcxfzeng/blog/item/dc0ea8fd041cc64bd6887d8a.html
SAC是平常处理地震记录的软件,但是不单可以用来处理地震记录,实际上可以实现信号处理的很多功能,可惜没有汉化破解版:)。
推荐的两个
sac manual:
http://www.iris.edu/software/sac/manual.html
GaTech Prof. Peng zhigang 的一个短小的入门手册
http://geophysics.eas.gatech.edu/classes/SAC/
3.1 怎么开始
在服务器上,sac已经装好了,但是需要添加一点点内容到你的home目录下.bashrc里面
export SACAUX=/opt/sac/aux
export SACDIR=/opt/sac
export SACXWINDOWS=X11
export PATH=$PATH:/opt/sac/bin:
如果已经有export PATH,那么把/opt/sac/bin:加到后面
然后source .basrc就可以了,在server1下可以用sac100,sac2000,s1下用sac就可以了.
如果是在自己机器上,可以拷贝服务器上的/opt/sac,多数情况下不需要编译。 用windows的同志们就不要想了。
3.2 什么是sac文件
sac文件包括两部分, 一部分是头文件,一部分是记录。
头文件描述的是这个记录的相关信息,比如有几个采样点,采样间隔,起始时间,台站位置等
详细的见这里:
http://hi.baidu.com/ustcxfzeng/blog/item/b82fa51b014a8fffae513346.html
3.3 开始使用
首先我们可以找一个现成的sac文件x.sac,下面就是怎么打开一个文件,绘制图形,保存文件
sac100sac>rx.sac//读入sac文件sac>p//画出来sac>wy.sac//保存为y.sacsac>q//退出sac要是你觉得这个波形很pp,想保存下来显摆一下,该怎么办呢?方法1.摁键盘上的PrintScreen,然后打开你个画图,ctrl+v,当然我们没有这么老土啦。所以都有方法二,方法san啦方法2.
sac
sac> r x.sac
sac> p
sac> bg sgf //开启sgf作为绘图终端
sac>p //在sgf里面画个图
sac> ed sgf //退出sgf
sac> q
sgftopsf001.sgfx.ps//将sgf文件转换为ps文件
ps2pdf x.pdf x.ps // 把ps文件专程pdf
然后你就可以用pdf文件做BSOl啦。。
3.4 标上到时
$ sac
sac > r x.sac
sac > p
sac > ppk // 这就可以标到时了,将鼠标移到,图形窗口,在十字光标,输入t1,即标记为t1, 再摁q,退出
sac > w over //写入
sac > q
上面只是比较简单的标记,很多时候不是这么容易,那么下面的处理也许有帮助。
1。 波形记录很长的话,光标怎么移动造成的误差会比较大,那么可以利用xlim 来处理
xlim会调整绘图范围,比如P初动在110, 而整个记录从0 到1000,
sac > xlim 100 120 //将绘图范围限制在100s 到 120s
sac > ppk
另一个方法是ppk然后,在图形窗口移动光标,输入x,然后在希望的另一端点击左键。比如 希望窗口为 100 200,那么可以先将光标移动到100,摁x,然后移动到200,点击左键。那么相当于xlim 100 200. 这个步骤可以重复,如果希望回到上一步的窗口,可以直接摁o(range)。
很多时候原始记录中,信噪比SNR比较低,那么可以做一些滤波处理,比如对远震信号
sac > r x.sac
sac > bp c 0.05 0.2 n 4 p 2 //带通滤波,拐角频率是0.05Hz, 0.2 Hz,可以滤掉台站附近的高频噪声
sac > ppk
但是滤波参数应该根据震相,比如PKiKP的频率就比较高,再用上面的滤波就合适。三分量结合
比如在识别S波的时候,要避免P-SV波的影响,那么需要用到三分量波形比较,利用T(切向)分量的SH来标记。那么需要用到rot, rot的作用是讲两个分量旋转到另一个方向,比如大圆的径向和切向。
sac > cut 100 200 //只读入100s和200s间的波形
sac > r x_ew.sac x_ns.sac //读入ew, ns向
sac > rot to gcp //旋转到 gcp(大圆)
sac > w x_r.sac w x_t.sac //另存为
sac > cut off
sac > r x_r.sac x_t.sac x_ud.sac //读入三分量
sac > ppk
sac > q
需要注意的是,rot需要根据sac头文件中的cmpinc, cmpaz计算, 如果有报错,需要检查这两个头文件信息。
3.5 老板教导我们要去掉仪器响应,怎么办
首先要从seed文件中解压出仪器响应文件, 有两种,一个是零极点,一个evalresp文件.
rdseed -pf x.seed 出来的SAC_PZs_XX_x_BHZ_00 就是 XX_x_BHZ_00.SAC 的零极点响应文件
rdseed -Rf x.seed 出来的RESP.XX.x.00.BHZ 就是完整的evalresp文件
以下假设仪器响应文件和sac文件在同一个目录下
sac > r x.bhz.sac
sac > rtr
sac > rmean
sac > taper //这个可以抑制由滤波引起的记录两端虚假信号
sac > trans from evalresp to vel freq 0.003 0.004 5 10 // 去除仪器响应到速度记录. vel,acc,none都是可选项, 推荐用vel
sac > w over
如果是零极点型则是 trans from subtype polozeros s SAC_PZs… to vel freq 0.003 0.004 5 10
3.6 其他处理。
1. 如果需要知道某一个波形的频谱怎么办?
sac > r x.sac
sac > fft //做快速傅立叶变换
sac > plotsp //画强度(amplitude)和相位(phase)
sac > keepam //只保留amplitude信息
sac > w x_sp.sac //另存为
但是这时候读入x_sp.sac显得很难看,为什么? 因为用普通线性坐标造成比例不协调,可以改用指数坐标
sac > r x_sp.sac
sac > xlog //令x轴为log坐标
sac > ylog
sac > p
2. 据说Rayleigh的振动轨迹为逆进椭圆,怎么表示出来? 用plotpm
3. 采样率不合我口味怎么办?
decimate 2 //将采样率降为一半
interpolate delta 0.5 //插值为0.5s等间隔采样。
4. 对数据进行数学运算
add sub div mul abs sqrt log log10 int dif等等