关于sac的一些笔记

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)。

1. 很多时候原始记录中，信噪比SNR比较低，那么可以做一些滤波处理，比如对远震信号
   sac > r x.sac
   sac > bp c 0.05 0.2 n 4 p 2 //带通滤波，拐角频率是0.05Hz, 0.2 Hz，可以滤掉台站附近的高频噪声
   sac > ppk
   但是滤波参数应该根据震相，比如PKiKP的频率就比较高，再用上面的滤波就合适。

2. 三分量结合
   比如在识别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等等