SOD,全称是 Standing Order for Data,是一个可以自动筛选、下载以及预处理地震数据的工具,其具有高度可定制的特点,可以满足日常科研中数据申请的几乎全部需求。

本文是 SOD 的入门教程,目的是使读者对 SOD 有初步的认识,详细教程及用法见官方网站。

Why SOD?

从 IRIS 上申请数据的工具有很多,比如 Wilber、BREQ_FAST,为什么要用 SOD 呢? SOD 有哪些优势呢?

先回顾一下申请地震波形数据的一般流程:

  1. 获取完整的地震目录,并根据实际需求筛选出自己需要的地震目录
  2. 获取台站列表,根据需求对震中距、方位角等进行筛选
  3. 根据发震时刻、震中距、震相走时等信息确定要申请的波形数据的起始、结束时间
  4. 申请数据
  5. 下载数据
  6. 数据解压以及预处理

以上提到的每个步骤都需要用户写一个以上的脚本来实现,如果需求改变,则可能需要修改多个脚本,总的来说是件很麻烦的事情。

SOD 的优势如下:

  1. 只需要一个配置文件(SOD中称为recipe)即可完成上面提到的全部步骤

  2. 灵活的地震筛选规则:

    • 可以在 IRIS/USGS 提供的地震目录的基础上根据发震时刻、深度、震级、震源所在区域等条件进行筛选
    • 可以由用户自行提供地震目录(需要满足 SOD 对地震目录的格式要求)
    • 可以构建虚拟地震以申请连续波形数据
  3. 灵活的台站筛选规则:

    • 可以选择一个或多个指定的台网、台站或通道
    • 可以筛选在某个特定点的特定距离、方位角、反方位角内的台站
    • 可以筛选某个区域内的台站
    • 可以下载指定通道的仪器响应文件(SAC PZ 和 RESP)
  4. 灵活的波形筛选规则:

    • 仅下载震中距、方位角、反方位角满足条件的数据
    • 仅当某个震相存在时才下载该台站的数据(比如仅当震中距大于115时才有PKIKP震相)
    • 仅当地震与台站的连线中点位于某个区域是才下载数据(比如想要用PcP研究某个特定区域的核幔边界结构)
    • 仅当指定震相的前后多少秒内没有其他震相干扰时才下载数据
    • 灵活的时窗截取规则,可自动计算理论到时以设置数据窗
    • 可直接生成 BREQ_FAST 文件以通过邮件方式申请数据
    • 以 miniSEED、SAC、ASCII、wav 格式将数据保存到磁盘中
  5. 数据并行下载,速度快无延迟

  6. 将数据写入磁盘前/后,可对数据进行常见的数据预处理:

    • rmean, rtrend, taper, merge
    • decimate
    • differentiate, integrate
    • filter, transfer
    • phaseCut

SOD 除了上面提到的功能/优点外,还有一个优点:入门很容易,别人写的 recipe 一看就懂,改起来也简单,但一个很明显的缺点是:自己从头写一个 recipe 很困难,因为 SOD 其实是个很复杂的东西。

为了尽可能解决这个问题,我创建了 SODrecipes 项目,其中收集整理了 SOD 的常见用法,用户可以通过阅读、修改和使用这些 recipe 片段来理解 SOD 的用法。也欢迎更多用户贡献自己的recipe。

安装 SOD

由于 SOD 是 Java 写的,因而可以跨平台直接使用。

  1. 下载:

     $ wget http://www.seis.sc.edu/downloads/sod/3.2.8/sod-3.2.8.tgz
    
  2. 解压并安装

     $ tar -xvf sod-3.2.8.tgz
     $ sudo mv sod-3.2.8 /usr/local
    
  3. 修改环境变量 PATH

     $ echo 'export PATH=${PATH}:/usr/local/sod-3.2.8/bin'>> ~/.bashrc
    
  4. 重启终端,并执行 sod -h 以确认 SOD 是否正常安装

SOD 的使用

SOD 的 bin 目录下提供了8个命令:

  • sod
  • psod
  • bgsod
  • find_events
  • find_stations
  • find_channels
  • find_responses
  • find_seismograms

本文只介绍命令 sod 的用法。

命令行语法

sod 的命令行语法很简单,直接

sod -f recipe.xml

其中 recipe.xml 是 sod 的配置文件。sod 的配置文件称为 recipe,是一个 xml 格式的纯文本文件,其中定义了地震数据筛选的一系列规则。

运行示例

执行 sod --demo -r > demo.xml 会生成一个名为 demo.xml 的 recipe 示例,用编辑器查看其内容以对 recipe 的结构和格式有个初步印象。

执行以下命令则会运行该示例:

$ sod -f demo.xml
Congratulations, valid recipe.
Fiji region (-20.6, -177.7) 378 km 2003/01/04 05:15:03 UTC 6.5 mwc
...
Channel: II.AAK.00.BHE
...
Got 1 seismograms for II.RPN.00.BHN for eq on 2003/01/09 02:50:45 UTC
...

输出中的 Congratulations, valid recipe. 表明当前 recipe 可以使用,紧接着会输出事件信息、通道信息以及获取地震波形的信息。下载的地震波形数据保存在 seismograms 目录下。

SOD 数据库文件

每次执行 sod 时,sod 都会在当前目录下生成一个数据库文件夹 SodDb 以及一些 log 文件。数据库文件夹 SodDb 的数据库文件中包含了单次 sod 运行时的全部信息,比如选中了哪些事件、哪些台站以及事件台站对所对应的波形,也包含了哪些数据已经下载哪些数据尚未下载的信息。

SOD 数据库的主要作用是:若 sod 在下载数据过程中意外退出,可以重新执行 sod 命令,由于 sod 数据库中记录了哪些数据已经被下载,因而 sod 可以继续下载那些尚未被下载的数据,即实现断点续传的功能(注:实测中发现的确可以断点续传,但似乎在所有数据下载完成后程序无法正常退出)。

recipe

一个 SOD recipe 中包含了一系列规则以限定要申请的数据的范围。

xml 格式简介

SOD recipe 本质上是一个 xml 格式的文本文件。

XML 文档形成了一种“树”结构,它从“根部”开始,不断扩展到“枝叶”。一个简单的 XML 文件如下:

<?xml version="1.0"?>
<root>
    <!-- this is a comment -->
    <child>
        <name>Child One</name>
        <age>10</age>
    </child>
    <child>
        <name>Child Two</name>
        <age>8</age>
    </child>
    <printline />
</root>

其中:

  • <?xml version="1.0"?> 是 XML 声明,其指定了 XML 的版本信息
  • <root></root> 是一对标签的开始标签和关闭标签,从开始标签到结束标签中的内容称之为元素
  • 所有 XML 元素都必须有关闭标签
  • 标签可以嵌套多层,但必须正确地嵌套
  • 若某个标签中没有嵌套其他标签,则一对标签可以简写为 <tagName/> 的形式,例如 <printline/>
  • <!-- --> 用于注释
  • XML 标签是区分大小写的

sod recipe 结构

SOD recipe 是 XML 格式的文件,其根标签为 <sod>,根标签下有五个直系子标签,分别是:

  • properties: SOD程序相关的属性配置(一般不用)
  • eventArm: 事件筛选的规则
  • networkArm: 台站/通道筛选的规则
  • waveformArm: 确定波形的时间段,筛选、下载并处理波形数据
  • waveformVectorArm: 处理三分量波形数据(不常用)

通常只需要用到 eventArmnetworkArmwaveformArm。显然, eventArmnetworkArm 是互相独立的,而 waveformArm 则依赖于 eventArmnetworkArm

recipe结构

因而,所有的 recipe 都符合如下基本结构:

<?xml version="1.0"?>
<sod>
    <properties>
    <!-- sod 属性,很少使用 -->
    </properties>

    <eventArm>
    <!-- eventArm 全部内容 -->
    </eventArm>

    <networkArm>
    <!-- networkArm 全部内容 -->
    </networkArm>

    <waveformArm>
    <!-- waveformArm 全部内容 -->
    </waveformArm>

    <waveformVectorArm>
    <!-- waveformVectorArm 全部内容 -->
    </waveformVectorArm>
</sod>

下面会介绍 eventArm、networkArm 和 waveformArm 的写法和功能。

eventArm

eventArm 的作用是定义一定的规则以筛选所需要的地震事件。

下面看一个 eventArm 的简单示例:

<?xml version="1.0"?>
<sod>
    <eventArm>
        <fdsnEvent>
            <originTimeRange>
                <startTime>2010-01-01T00:00:00.000Z</startTime>
                <endTime>2011-01-01T00:00:00.000Z</endTime>
            </originTimeRange>
            <magnitudeRange>
                <min>5.0</min>
                <max>8.5</max>
            </magnitudeRange>
            <originDepthRange>
                <unit>KILOMETER</unit>
                <min>30</min>
            </originDepthRange>
        </fdsnEvent>
        <printlineEventProcess/>
        <CSVEventPrinter>
            <filename>events.csv</filename>
        </CSVEventPrinter>
    </eventArm>
</sod>

在每个 Arm 内,SOD 会按照各个标签出现的先后顺序去执行相关操作。那么,如何去读一个recipe呢?其实很简单,由于整个XML文档是一个“树”结构,因而只要 一级一级读 即可。

以上面的例子为例, eventArm 的下一级标签有三个,SOD 会按照三个标签出现的顺序依次执行。三个标签是:

  • fdsnEvent 的作用是从 USGS 的 FDSN web service 中获取全部地震目录
  • printlineEventProcess 输出地震目录中每个事件的信息
  • CSVEventPrinter 将事件信息保存到CSV文件中

读完了这一级标签后,再读下一级标签:

  • fdsnEvent 的下一级标签有: originTimeRangemagnitudeRangeoriginDepthRange, 分别用于限制地震的发震事件、震级和发震时刻。
  • CSVEventPrinter 下一级标签 filename 则表明要将事件信息写到 filename 所指定的 CSV 文件中。

eventArm 里可以放哪些标签? fdsnEvent 标签里又可以嵌套哪些标签?这就是自己写 recipe 的难点,因为 SOD 的标签实在太多了。完整的标签列表在 Ingredient Listing 中。

networkArm

networkArm 用于筛选出需要的台站。

<?xml version="1.0"?>
<sod>
    <networkArm>
        <fdsnStation/>
        <networkCode>II</networkCode>
        <channelCode>BH*</channelCode>
        <printlineChannelProcess/>
    </networkArm>
</sod>

解释一下:

  • fdsnStation 是从 IRIS 上获取了全部的台站信息
  • networkCode 从中筛选除了 II 台网
  • channelCode 从中进一步筛选除了 BH* 通道

最终筛选出了 II 台网所有 BH* 通道。

SOD中的Arm就像一个流水线,Arm中的每个标签就像一个质检员。首先获取所有的原材料,若材料满足标签所指定的条件,则进入流水线的下一步,否则会直接被丢弃。最终留下来的材料则是满足所有条件的材料。

waveformArm

在确定了地震事件和台站之后,就可以对二者进行组合,并进一步筛选所需要的波形数据。

waveformArm流程

waveformArm 相对比较复杂,其包含了如下几个步骤(对应上图中的1-6):

  1. 根据事件台站信息对震中距、方位角进行筛选
  2. 根据地震发震时刻或者震相走时确定要申请的数据时间范围
  3. 向数据服务器发送请求,检查对应时间段内是否有数据
  4. 检查服务器的返回,确认是否存在数据
  5. 从服务器中获取数据
  6. 对数据进行预处理

下面的示例展示了一个简单的 waveformArm:

<?xml version="1.0"?>
<sod>
    <eventArm>eventArm同之前的示例,出于篇幅考虑,此处省略</eventArm>
    <networkArm>networkArm同之前的示例,出于篇幅考虑,此处省略</networkArm>
    <waveformArm>
        <distanceRange>
            <unit>DEGREE</unit>
            <min>30</min>
            <max>90</max>
        </distanceRange>
        <phaseRequest>
            <model>prem</model>
            <beginPhase>ttp</beginPhase>
            <beginOffset>
                <unit>SECOND</unit>
                <value>-120</value>
            </beginOffset>
            <endPhase>tts</endPhase>
            <endOffset>
                <unit>SECOND</unit>
                <value>360</value>
            </endOffset>
        </phaseRequest>
        <fdsnDataSelect/>

        <printlineSeismogramProcess/>
        <rMean/>
        <rTrend/>
        <taper/>
        <sacWriter/>
    </waveformArm>
</sod>
  • distanceRange 筛选出满足一定震中距的数据
  • phaseRequest 定义了数据时间窗,具体时间范围是P波到时前120秒,S波到时后360秒
  • fdsnDataSelect 表示从IRIS获取数据
  • printlineSeismogramProcess 打印获取的地震波形的信息
  • rMeanrTrendtaper 数据预处理
  • sacWriter 以 SAC 格式将数据写入磁盘

高级用法、技巧与陷阱

eventArm 错误示例

上面的 eventArm 示例还有另一种等效的写法:

<?xml version="1.0"?>
<sod>
    <eventArm>
        <fdsnEvent>
            <originTimeRange>
                <startTime>2010-05-19T00:00:00.000Z</startTime>
                <endTime>2010-10-20T00:00:00.000Z</endTime>
            </originTimeRange>
        </fdsnEvent>
        <magnitudeRange>
            <min>5.0</min>
            <max>8.5</max>
        </magnitudeRange>
        <originDepthRange>
            <unit>KILOMETER</unit>
            <min>30</min>
        </originDepthRange>
        <printlineEventProcess/>
        <CSVEventPrinter>
            <filename>events.csv</filename>
        </CSVEventPrinter>
    </eventArm>
</sod>

这种写法也是正确的,但是这种写法的速度却比前一种写法慢很多。其原因在于,SOD 在执行 fdsnEvent 时,前一种写法是获取了2010年、震级在5.0到8.5、深度大于30千米的地震,而后一种写法是先获取了2010年一整年的地震目录,然后筛选出其中震级在5.0到8.5范围内的地震,再筛选出深度大于30km的地震,虽然最终效果是一样的,但是后一种写法中 fdsnEvent 请求了更多的数据,因而速度更慢。

逻辑运算

前面 eventArm 的示例已经在 fdsnEvent 中获取了2010年整年的震级在5.0到8.5、深度大于30km的地震。如果想要更复杂的筛选逻辑, fdsnEvent 就无法实现了,此时就需要用到 SOD 提供的逻辑运算功能。( fdsnEvent 本质上是调用了 IRIS 或 USGS 提供的 web service 接口,因而其功能由 web service 所限制,而无法实现复杂的筛选逻辑,所以需要在执行 fdsnEvent 之后再由 SOD 对地震目录作进一步筛选)。

比如在已有地震的基础上想要筛选出满足如下条件的地震事件:震级大于6级或者震级大于5级且深度在100到200千米之间。这一筛选条件写成一般的伪代码大概是这样:

if ((mag>5 and 100<depth<200) or mag>6.0):
    pass
else
    fail

这一筛选条件可以写成以下recipe 片段并放在 printlineEventProcess 之前即可:

<originOR>
    <originAND>
        <magnitudeRange>
            <min>5</min>
        </magnitudeRange>
        <originDepthRange>
            <unit>KILOMETER</unit>
            <min>100</min>
            <max>200</max>
        </originDepthRange>
    </originAND>
    <magnitudeRange>
        <min>6</min>
    </magnitudeRange>
</originOR>

不单单是 eventArm 中可以使用OR、AND、NOT逻辑运算,networkArm 中也有很多标签可以使用这些逻辑运算以实现灵活筛选的目的,比如想要获取两个台网的数据,其伪代码是:

if (network == IU or network == II)
    pass
else
    fail

所以 recipe 片段可以写成:

<networkOR>
    <networkCode>II</networkCode>
    <networkCode>IU</networkCode>
</networkOR>

而不能写成:

<networkAND>
    <networkCode>II</networkCode>
    <networkCode>IU</networkCode>
</networkAND>

从头开始写 recipe

SOD 具有高度可订制性的优点,同时也导致写 recipe 变得很困难,这一节说一说如何去读 SOD 官方文档。

如之前所说,一个 SOD recipe 包含了 eventArm、networkArm 和 waveformArm,文件结构如下所示:

<?xml version="1.0"?>
<sod>
    <eventArm>
    <!-- eventArm 全部内容 -->
    </eventArm>

    <networkArm>
    <!-- networkArm 全部内容 -->
    </networkArm>

    <waveformArm>
    <!-- waveformArm 全部内容 -->
    </waveformArm>
</sod>

写 recipe 的关键在于搞清楚每个 Arm 中可以放哪些标签,各个标签的先后顺序又怎样?这就需要查询 SOD 的 Ingredient Listing

以 eventArm 为例,解释一下如何读 SOD 的文档。

首先,进入 eventArm 文档页面,可以看到:

  • eventArm 的功能说明
  • Example 中给出了 eventArm 的简单示例
  • This consists of 表明 eventArm 中可以包含两部分内容,分别是是 eventSource 和 origin,其中前者必须至少出现一次(At least once),后者可以出现任意多次(Any number of times)
  • Places this can be found 表明 eventArm 只能直接出现在 sod 标签下

需要注意的是, eventSourceorigin 两边没有尖括号 < >,所以这两个 都不是标签 ,可以将这两个理解为 标签集

进入 eventSource 文档页面,可知 eventSource 中是用来获取地震事件的标签集,其包含了如下几个标签:

  • backwardsEventFinder
  • CSVEventSource
  • delayedEventSource
  • eventFinder
  • fdsnEvent
  • nowFakeEventSource
  • periodicFakeEventSource

这里我们看到了前面示例用到的 fdsnEvent 。再进入 fdsnEvent 文档页面,可知 fdsnEvent 会从USGS网站中获取地震目录。 fdsnEvent 中包含了很多标签以对事件进行筛选,比如示例中用到的 originTimeRangemagnitudeRangeoriginDepthRange ,以及示例中不曾用到的 boxAreapointDistance 等标签。

再到标签 originTimeRangemagnitudeRangeoriginDepthRange 的文档页面看一看,就很容易把 eventArm 的 eventSource 部分 recipe 片段写出来了,如下:

<?xml version="1.0"?>
<sod>
    <eventArm>
        <fdsnEvent>
            <originTimeRange>
                <startTime>2010-01-01T00:00:00.000Z</startTime>
                <endTime>2011-01-01T00:00:00.000Z</endTime>
            </originTimeRange>
            <magnitudeRange>
                <min>5.0</min>
                <max>8.5</max>
            </magnitudeRange>
            <originDepthRange>
                <unit>KILOMETER</unit>
                <min>30</min>
            </originDepthRange>
        </fdsnEvent>
    </eventArm>
</sod>

写完了 eventSource 部分,再看看 origin 部分。进入 origin 的文档页面,可以看到其包含的一堆标签,比如:

  • CSVEventPrinter 把事件以CSV格式保存到磁盘中
  • eventArea 筛选出某个规则区域的事件
  • eventPolygonFile 筛选出某个多边形内的事件
  • printlineEventProcess 终端打印出筛选出来的事件信息

再把文档一级一级看下去,即可写出本文最开始所展示的示例片段,如下:

<?xml version="1.0"?>
<sod>
    <eventArm>
        <fdsnEvent>
            <originTimeRange>
                <startTime>2010-01-01T00:00:00.000Z</startTime>
                <endTime>2011-01-01T00:00:00.000Z</endTime>
            </originTimeRange>
            <magnitudeRange>
                <min>5.0</min>
                <max>8.5</max>
            </magnitudeRange>
            <originDepthRange>
                <unit>KILOMETER</unit>
                <min>30</min>
            </originDepthRange>
        </fdsnEvent>
        <printlineEventProcess/>
        <CSVEventPrinter>
            <filename>events.csv</filename>
        </CSVEventPrinter>
    </eventArm>
</sod>

总结一下,SOD中标签很多,但由于整个 recipe 是树形结构,所以只要一级一级不断看文档,就可以定制出自己所需的任意功能的 recipe 了。