Perl 中调用 SAC

简单示例

下面的脚本展示了如何在 Perl 中调用 SAC。

下载地址： 0.simple-script.pl

#!/usr/bin/env perl
use strict;
use warnings;
$ENV{SAC_DISPLAY_COPYRIGHT}=0;

open(SAC, "| sac ") or die "Error opening sac\n";
print SAC "fg seismo \n";
print SAC "lh evla kstnm \n";
print SAC "q \n";
close(SAC);

Perl 中调用 SAC 本质上是使用 open(SAC, "| sac ") 语句定义了一个名为 SAC 指向 sac 的句柄，然后通过 print SAC 语句将要执行的 SAC 命令传递给 SAC。

数据转换

首先要将 SEED 格式的数据转换成 SAC 格式。

• rdseed 一次只能处理一个 SEED 数据

• rdseed 的 -pdf 选项会提取出 SAC 波形数据和 PZ 格式的仪器响应文件

下载地址： 1.rdseed.pl

警告

IRIS 自 2020 年 1 月开始不再提供 SEED 格式的数据下载支持，详情可阅读 IRIS 数据服务通讯。IRIS 也已不再维护 rdseed 软件。

#!/usr/bin/env perl
use strict;
use warnings;

@ARGV == 1 or die "Usage: perl $0 dirname\n";
my ($dir) = @ARGV;

chdir $dir;  # cd 进数据所在目录，以避免处理路径

# rdseed 一次只能处理一个 SEED 文件
foreach my $seed (glob "*.seed") {
    system "rdseed -pdf $seed";
}

chdir "..";

文件合并

SEED 文件的波形数据可能会因为多种原因而出现间断，导致同一个通道会解压出来多个 SAC 文件，因而需要将属于同一个通道的 SAC 数据合并起来。

• 此处使用了新版 merge 命令的语法，要求 SAC 版本大于 v101.6

• merge 命令还有更多选项可以控制数据合并的细节，见命令的语法介绍

• 合并后的数据，以最早的数据段的文件名命名

• 多余的数据段均被删除，以保证每个通道只有一个 SAC 文件

• 由于脚本运行速度比 SAC 运行速度快，因而应先退出 SAC 再删除多余的数据段

下载地址： 2.merge.pl

#!/usr/bin/env perl
use strict;
use warnings;
$ENV{SAC_DISPLAY_COPYRIGHT}=0;

@ARGV == 1 or die "Usage: perl $0 dirname\n";
my ($dir) = @ARGV;

chdir $dir;

# 利用 hash 的 key 的不可重复性构建集合:
#   hash 的 key 定义为 NET.STA.LOC.CHN
#   hash 的 value 是文件名与 key 匹配的 SAC 文件数目
my %sets;
foreach (glob "*.SAC") {
    # 将文件名用'.'分割，并取其中的第7到10个元素
    my ($net, $sta, $loc, $chn) = (split /\./)[6..9];
    $sets{"$net.$sta.$loc.$chn"}++;
}

# 在循环体外部调用 SAC，若没有数据需要 merge，则此次调用是做无用功
# 若将 SAC 调用放在循环体内，则可能会多次调用 SAC，而影响运行效率
open(SAC, "| sac") or die "Error in opening sac\n";
print SAC "wild echo off \n";
my @to_del;
while (my ($key, $value) = each %sets) {
    # 检查每个 key 所对应的 value，只有 value 大于1时才需要 merge
    next if $value == 1;

    print STDERR "merge $key: $value traces\n";
    my @traces = sort glob "*.$key.?.SAC";
    # 对 glob 的返回值做 sort，以保证第一个元素是最早的数据段

    # 在 SAC 中使用通配符而不是使用 @traces 以避免命令行过长的问题
    # merge 不支持通配符
    print SAC "r *.$key.?.SAC \n";  # SAC v101.6 or later only
    # 缺数据则补0，重叠则做平均
    print SAC "merge gap zero overlap average \n";
    print SAC "w $traces[0] \n";  # 以最早数据段的文件名保存

    # 将要删除的数据段文件名保存到数组中，待退出 SAC 后再删除
    push @to_del, splice(@traces, 1);
}
print SAC "q \n";
close(SAC);
unlink(@to_del);

chdir "..";

文件重命名

从 SEED 解压出的 SAC 文件名较长，因而对其重命名以简化。

• SEED 解压出的默认文件名格式为 yyyy.ddd.hh.mm.ss.ffff.NN.SSSSS.LL.CCC.Q.SAC

• 重命名后的文件名为 NET.STA.LOC.CHN.SAC

下载地址： 3.rename.pl

#!/usr/bin/env perl
use strict;
use warnings;

@ARGV == 1 or die "Usage: perl $0 dirname\n";
my ($dir) = @ARGV;

chdir $dir;

foreach my $file (glob "*.SAC") {
    my ($net, $sta, $loc, $chn) = (split /\./, $file)[6..9];
    rename $file, "$net.$sta.$loc.$chn.SAC";
}

chdir "..";

添加事件信息

若 SEED 中不包含事件信息，则解压得到的 SAC 文件中也不会包含事件信息。因而需要用户手动添加事件的发震时刻、经纬度、深度和震级信息。

• 输入参数包括：目录名、发震时刻、经度、纬度、深度、震级

• 发震时刻的格式为 yyyy-mm-ddThh:mm:ss.xxx，其中 T 用于分隔日期和时间

下载地址： 4.eventinfo.pl

#!/usr/bin/env perl
use strict;
use warnings;
use POSIX qw(strftime);

@ARGV == 6 or die
    "Usage: perl $0 dirname yyyy-mm-ddThh:mm:ss.xxx evlo evla evdp mag\n";
my ($dir, $origin, $evlo, $evla, $evdp, $mag) = @ARGV;

# 对发震时刻做解析
my ($date, $time) = split "T", $origin;
my ($year, $month, $day) = split "-", $date;
my ($hour, $minute, $second) = split ":", $time;
# 秒和毫秒均为整数
my ($sec, $msec) = split /\./, $second;
$msec = int(($second - $sec) * 1000 + 0.5);

# 计算发震日期是一年中的第几天
my $jday = strftime("%j", $second, $minute, $hour, $day, $month-1, $year-1900);

chdir $dir;
open(SAC, "| sac") or die "Error in opening SAC\n";
print SAC "wild echo off \n";
print SAC "r *.SAC \n";
print SAC "synchronize \n";   # 同步所有文件的参考时刻
print SAC "ch o gmt $year $jday $hour $minute $sec $msec \n";
print SAC "ch allt (0 - &1,o&) iztype IO \n";
print SAC "ch evlo $evlo evla $evla evdp $evdp mag $mag \n";
print SAC "wh \n";
print SAC "q \n";
close(SAC);

chdir "..";

去仪器响应

使用 PZ 文件去仪器响应。若数据的时间跨度太长，在该时间跨度内可能仪器响应会发生变化，因而会存在一个通道有多个 PZ 文件的情况。目前该脚本在遇到这种情况时会自动退出。

下载地址： 5.transfer.pl

#!/usr/bin/env perl
use strict;
use warnings;
$ENV{SAC_DISPLAY_COPYRIGHT} = 0;

@ARGV == 1 or die "Usage: perl $0 dirname\n";
my ($dir) = @ARGV;

chdir $dir;

open(SAC, "| sac") or die "Error in opening sac\n";
foreach my $sacfile (glob "*.SAC") {
    my ($net, $sta, $loc, $chn) = split /\./, $sacfile;

    # 找到对应的PZ文件
    # PZ文件名为: SAC_PZs_NET_STA_CHN_LOC_BeginTime_EndTime
    my @pz = glob "SAC_PZs_${net}_${sta}_${chn}_${loc}_*_*";
    die "PZ file error for $sacfile \n" if @pz != 1;

    print SAC "r $sacfile \n";
    print SAC "rmean; rtr; taper \n";
    print SAC "trans from pol s $pz[0] to none freq 0.01 0.02 5 10\n";
    print SAC "mul 1.0e9 \n";
    print SAC "w over\n";
}
print SAC "q\n";
close(SAC);

chdir "..";

分量旋转

将成对的水平分量旋转到大圆路径。

• 检查三分量是否缺失

• 检查 delta 是否相等

• 取三分量中的最大 b 和最小 e 值作为数据窗口，此操作要求三分量的 kzdate 和 kztime 必须相同，这一点 在添加事件信息时使用 synchronize 已经实现

• 检查两个水平分量是否正交，若不正交，则无法旋转

下载地址： 6.rotate.pl

#!/usr/bin/env perl
use strict;
use warnings;
use List::Util qw(min max);

@ARGV == 1 or die "Usage: perl $0 dirname\n";
my ($dir) = @ARGV;

chdir $dir;

# 利用hash的key的不可重复性构建集合:
#   hash的key定义为 NET.STA.LOC.CH (分量信息不是BHZ而是BH)
#   hash的value是文件名与key匹配的SAC文件数目，正常情况下是3的整数倍
my %sets;
foreach (glob "*.SAC") {
    my ($net, $sta, $loc, $chn) = split /\./;
    my $chn2 = substr $chn, 0, 2;
    $sets{"$net.$sta.$loc.$chn2"}++;
}

open(SAC, "|sac") or die "Error in opening sac\n";
# 对所有的key做循环
foreach my $key (keys %sets) {
    my ($E, $N, $Z, $R, $T, $Z0);

    # 检查Z分量是否存在
    $Z = "${key}Z.SAC";
    if (!-e $Z) {  # 若不存在，则删除该台站的所有数据
        warn "$key: Vertical component missing!\n";
        next;
    }

    # 检查水平分量是否存在
    if (-e "${key}E.SAC" and -e "${key}N.SAC") {   # E和N分量
        $E = "${key}E.SAC";
        $N = "${key}N.SAC";
    } elsif (-e "${key}1.SAC" and -e "${key}2.SAC") {  # 1和2分量
        $E = "${key}1.SAC";
        $N = "${key}2.SAC";
    } else {   # 水平分量缺失
        warn "$key: Horizontal components missing!\n";
        next;
    }

    # 检查水平分量是否正交
    my (undef, $cmpaz_E) = split m/\s+/, `saclst cmpaz f $E`;
    my (undef, $cmpaz_N) = split m/\s+/, `saclst cmpaz f $N`;
    my $cmpaz_delta = abs($cmpaz_E - $cmpaz_N);
    unless ((abs($cmpaz_delta - 90) <= 0.01) or
        (abs($cmpaz_delta - 270) <= 0.01)) {
        warn "$key: $E $N are not orthogonal!\n";
        next;
    }

    # 假定kzdate和kztime相同
    # 检查B, E, DELTA
    my (undef, $Zb, $Ze, $Zdelta) = split " ", `saclst b e delta f $Z`;
    my (undef, $Eb, $Ee, $Edelta) = split " ", `saclst b e delta f $E`;
    my (undef, $Nb, $Ne, $Ndelta) = split " ", `saclst b e delta f $N`;
    unless ( $Zdelta == $Edelta and $Zdelta == $Ndelta) {
        warn "$key: delta not equal!\n";
        next;
    }

    # 获取三分量里的最大B和最小E值作为数据窗
    my $begin = max($Zb, $Eb, $Nb) + $Zdelta;
    my $end = min($Ze, $Ee, $Ne) - $Zdelta;

    # 输出文件名为 NET.STA.LOC.[RTZ]
    my $prefix = substr $key, 0, length($key)-2;
    $R = $prefix."R";
    $T = $prefix."T";
    $Z0 = $prefix."Z";

    print SAC "cut $begin $end \n";
    print SAC "r $E $N \n";
    print SAC "rotate to gcp \n";
    print SAC "w $R $T \n";
    print SAC "r $Z \n";
    print SAC "w $Z0 \n";
}
print SAC "q\n";
close(SAC);
unlink glob "*.SAC";

chdir "..";

数据重采样

通常有两种情况下需要对数据进行重采样：

• 原始数据的采样率过高，而实际研究中不需要如此高的采样率，此时，对数据做减采样可以大大减少数据量；

• 原始数据中，不同台站的采样率不同，可能会影响到后期的数据处理，因而需要让所有数据使用统一的采样率；

下面的Perl脚本中使用 interpolate 命令将所有数据重采样到相同的采样周期。用户可以在命令行中直接指定要使用的重采样后的采样周期，若命令行中的采样周期指定为0，则以大多数数据所使用的采样周期作为重采样后的采样周期。

下载地址： 7.resample.pl

#!/usr/bin/env perl
use strict;
use warnings;

@ARGV == 2 or die "Usage: perl $0 dirname [delta | 0]\n";
#若第二个参数为0，则取数据中出现次数最多的周期为采样周期

my ($dir, $delta) = @ARGV;

chdir $dir;

if ($delta == 0) {
    # 假定所有数据已旋转到RTZ坐标系，文件名格式为NET.STA.LOC.[RTZ]
    # hash的键为出现了的采样周期，其值为出现的个数
    my %group;
    foreach (glob "*.*.*.[RTZ]") {
        my (undef, $delta0) = split /\s+/, `saclst delta f $_`;
        $group{$delta0}++;
    }
    # 将hash按value排序，找到最大value所对应的key
    my @keys = sort {$group{$b} <=> $group{$a}} keys %group;
    ($delta) = @keys;
}

# 重采样
open(SAC, "|sac") or die "Error in opening sac\n";
foreach (glob "*.*.*.[RTZ]") {
    my (undef, $delta0) = split /\s+/, `saclst delta f $_`;
    next if $delta == $delta0;  # 不需要重采样

    print SAC "r $_ \n";
    # 用interpolate实现减采样或增采样
    # 若是减采样，则需要对数据做低通滤波以防止出现混淆效应
    # 低通滤波时或许需要加上p 2以避免滤波引起的相移
    printf SAC "lp c %f\n", 0.5/$delta if $delta > $delta0;
    print SAC "interpolate delta $delta \n";

    print SAC "w over\n";
}
print SAC "q\n";
close(SAC);

chdir "..";