2016, Vol. 31
2. 禁核试北京国家数据中心, 北京 100085
2. CTBT BeiJing National Data Centre, Beijing 100085, China
全面禁止核试验条约(The comprehensive nuclear-test-ban treaty, CTBT)要求所有缔约国不能进行“任何武器及非武器核试验”.全面禁止核试验条约条约组织(The comprehensive nuclear-test-ban treaty Organization, CTBTO)临时技术秘书处(Provisional Technical Secretariat，PTS)负责监督条约的执行，其下的国际数据中心(International Data Centre, IDC)负责对来自国际监测系统(International Monitoring System, IMS)台站数据的接收、处理、分析、报告及归档.同时条约要求IDC要不断进行技术能力的提升(Comprehensive nuclear-test-ban treaty 1996, Protocol, part I, paragraph 18(b)).
不断降低在全球范围内的监测阈值，提高自动处理结果的准确性是IDC一直以来面临的挑战.IDC目前运行的波形数据(地震、水声和次声，Seismic、Hydro-acoustic、and Infrasonic, SHI)处理软件Release 3已经应用近20年(Given, 2011)，新的处理算法、软件架构不断出现，PTS也一直进行软件改进以满足监测技术发展而带来的新的需求.PTS组织召开的国际科学技术大会ISS09(2009 International Scientific Studies Conference)(CTBTO, 2009)、ISS11、ISS13及ISS15(CTBTO, 2015)为各国专家共同讨论解决CTBT面临的技术挑战提供了很好的平台，涌现出了很多好的想法，部分算法已引起PTS关注并开始在IDC进行测试评估.
本文重点关注IDC地震数据处理软件最新研究进展，首先介绍了IDC地震数据处理过程，然后从软件结构优化、信号检测与关联算法改进、区域走时模型改善等几个方面总结了目前IDC在地震数据处理软件方面开展的部分研究工作，对涉及的相关算法进行了分析，对测试、应用及重点关注问题进行了介绍.
1 IDC地震数据自动处理过程IDC地震数据处理的目标是在给定的时间期限内产生一系列描述地震事件的公报(Coyne et al., 2012).对于IMS台站获取的数据，首先通过自动处理产生一系列的标准事件公报(Standard Event Bulletin, SEL)，然后分析员对自动处理公报进行审核和更改，产生审核事件公报(Reviewed Event Bulletin，REB)公报.作为IDC的产品，REB中包含事件发生时间、深度、震级、误差和其他事件特征.最后，按照CTBT要求对审核公报的结果进行筛选处理，特征明显的天然事件被筛选出去，形成标准筛选事件报告(SSEB，Standard Screened Event Bulletin)提供给缔约国做进一步分析(Le Bras and Wuster, 2002).
自动处理包括台站处理和台网关联处理.台站处理目标是检测信号、参数计算及震相识别.目前IDC软件中，信号检测采用短时能量和长时能量比值的方法实现，台阵在信号检测前要进行聚束滤波，台阵的方位角慢度采用频率-波数(frequency-wavenumber，f-k)(Kværna and Doornbos, 1986)分析实现，三分向台站方位角采用极化分析实现，采用AIC (Akaike Information Criterion)(Kværna, 1995)方法进行到时精确估计，通过预设条件和神经网络方法进行初步的震相识别.
台网关联处理是从多台检测信号中形成事件.目前IDC软件中采用基于格点搜索的全球关联(Global Association, GA)(Bache et al., 1993; Le Bras et al., 1994)方法实现事件关联，该方法类似于广义聚束法.预先采用相互覆盖的圆在地球表面建立定位格点文件(在一些了解深度地震活动性的区域还有一些额外的深度格点)，针对每个格点，形成“首震相到达台站列表”，关联时假设每个格点为可能的事件位置，以距离格点最近台站的震相作为驱动震相来预测在格点内发生事件的时间，形成种子事件，然后依据该事件在其他台站搜索相容震相，相容震相越多则事件越可信，由驱动震相和相容震相进行事件重新定位，根据与设定的事件定义标准进行比较形成最终结果.该方法最终给出一个详尽的事件、震相列表.关联处理过程中，需要解决相邻时间段、相邻区域间形成的事件的冲突问题.
2 IDC地震数据处理软件研究进展 2.1 软件结构优化自从IDC第一版SHI波形数据处理软件发布以来，由于监测技术的发展带来的新的功能需求及基础设施的不断变化，软件系统的部分模块已经逐渐被替换.但由于缺乏现代统一的软件架构，这些变换造成了数据处理软件系统的碎片化，存在着软件代码冗余、软件架构技术过时等问题.为了提高目前SHI软件可维护和可扩展性，IDC在2011年启动了SHI软件再工程化项目(Given, 2012)，分为三个阶段，其中第一阶段(2011-2014)以提高单个模块性能为重点，第二阶段(2014-2016)重点关注建立统一的波形数据处理软件架构，为进一步软件开发做好准备.第三阶段(2016-)进行新版软件的实现.根据统一软件开发过程(Rational Unified Process, RUP)方法论(姚登峰等，2009)，再工程化的第二阶段涵盖了初始(Inception)和细化(Elaboration)阶段，第三阶段涵盖了构造(Construction)和交付(Transition)阶段.再工程化第一阶段的工作目标是基于原有处理系统，显著改进目前的SHI自动处理和交互分析软件性能.主要开展的内容包括：采用基于贝叶斯推论和机器学习技术(NET-VISA)替换SHI台网处理软件，基于开源消息中间件开发新的SHI自动处理管道系统，设计和开发新的波形数据质量控制数据模型和应用编程接口等.
再工程化第二阶段目标是研发现代的、基于模型组件的波形数据处理软件架构，同时在该阶段建立基本的IT容灾恢复系统.具体内容包括：统一波形数据处理软件架构设计、开发扩展性强的SHI交互软件、提供基本的IT容灾恢复系统以及集成包括来自国家设施捐赠、国家技术方式和其他开源可用的SHI数据源等.该阶段的工作将有助于CTBTO中期战略计划(MTS，Mid Term Strategy)中核查体系运行和维护目标的实现，特别是增强了IDC基础设施以及PTS系统和网络的可维持性.该阶段软件研发将基于统一模型语言(UML，Unified Modelling Language)和RUP方法来实施，主要涉及到RUP中的初始(Inception)和细化(Elaboration)阶段.其中的Inception阶段包括软件需求分析和初始用例模型建立, 2014年已基本完成，主要开展了系统需求文档(System Requirements Document，SRD)、系统规格文档(System Specification Document，SSD)及用例模型概览(Use Case Model Survey)的编制工作，并召开专家技术会议，邀请各签约国专家对相关文档进行了审核和建议.在Elaboration阶段需要继续开展用例模型建立及实现、软件架构文档建立以及可执行的架构原型实现.第二阶段涉及到的核心技术包括数据存储和访问、改善分析接口、先进的处理和分析技术、高可用性和故障容错以及数据和产品分发.新版的数据处理软件需求主要围绕数据获取、申请和导入，数据处理，交互分析，事件筛选，特殊事件分析，数据产品分发以及监测和控制工具研发等部分展开.数据获取、产品分发等软件已经在第一阶段更新，第二阶段将重点关注自动处理、交互分析、事件筛选及监视和控制工具等部分的设计.
2.2 数据处理算法改进 2.2.1 信号检测1)广义F检测算法
地震台阵技术由于其在信号增强和噪声抑制方面的良好性能，已广泛应用于核查地震监测中(Filson, 1975; Rost and Thomas, 2002; Douglas, 2003).在台阵信号检测方面，英国NDC Selby (Selby, 2008)针对小孔径台阵处理提出了广义F检测算法.传统的F检测方法(Melton and Bailey, 1957)主要应用于台阵信号检测，可以认为是信号功率与噪声功率的比值.构建符合F分布的检验统计量，可以将检测问题转化为针对给定先验信噪比的检测概率问题.F检测的要求是各通道间信号是相关的，噪声是不相关的.但对于小孔径台阵( < 7 km)在部分监测敏感频带噪声可能是高度相关，这使得F检测性能下降，广义F检测器就是为解决这一问题而提出，它允许信号和噪声功率谱是任意的，并考虑各子台噪声间的相关性.
台阵技术用于地震信号检测主要优势包括：通过聚束方法改善信噪比、准确的方位角估计以及统计检测理论的应用.根据阵列信号处理理论，每个台阵所设计的孔径应等于所研究的特定类型地震信号视波长λa，国际监测系统的台阵半数以上是小孔径台阵( < 7 km), 主要用于检测区域距离(Δ < 17°)的地震和爆破的高频信号.
Blandford (1974)等首先将F检测成功应用于地震监测，对于台阵信号构建如下检验统计量，公式为
|
(1) |
其中dtk是子台k在采样时刻t处的值，K是所有子台数，T是计算F值用到的采样总数，
该方法是信号检测、估计理论与最优阵列处理理论的结合.定义的广义F检验统计量为
|
(2) |
s是对信号的加权最小二乘估计，其表达式为
|
(3) |
信号检测在频域进行，首先对检测信号采用加权的最小二乘法(Tarantola, 2005)进行估计，其中信号功率谱需要考虑震源谱和路径衰减，也可进行适当简化.噪声功率谱模型由Peterson (Peterson, 1993)低噪声模型确定初值，采用相互重叠的时间窗口计算得到，通过除以长时平均对模型进行不断更新，为了减少信号影响，当F值超过检测阈值时，不更新.根据Backus等(1964)方法设置噪声相关模型，噪声相关模型由不相关的白噪声和方位角各向同性且传播速度已知的噪声联合产生，假设由方位角各向同性的平面波产生相关噪声，计算得到相关矩阵.已知传播速度Vn下，频率为wi时，台站j和k之间的噪声互相关矩阵为
|
(4) |
这里Δ=|xj-xk，xj是台站j的位置，J0是第一类0阶贝塞尔函数.
根据信号估计值构建符合F分布的检验统计量，其中由于噪声的主要方位角不是固定的，可能随时间逐步改变，因此采用长时平均(Long-term Average，LTA)方法对F值进行调整，即计算得到的F值除以长时平均进行调整，该方法可以减少误检率，计算公式为
|
(5) |
F≥Ft时，检测到信号，阈值Ft是先验信噪比SNR=2时，非中心分布函数F(n1, n2, λ)在95%处的值，即观测到信噪比为2的信号概率为95%，则检测到信号.因此，对于给定自由度和先验信噪比情况下，广义F检测可以转化为信号出现的概率.
理想情况下，地震信号和噪声都是高斯、平稳的随机过程，广义F检测将是最优的信号检测器.在实际应用中，检测器性能完全依靠信号和噪声模型与理想情况的偏离程度.
在IMS13个小孔径台阵测试表明(Selby, 2011)，广义F算法检测数量只有IDC目前的51%，但F检测中的45%与REB中初至P波一致，而IDC目前的检测中只有20%与REB一致.IDC中未关联的检测数量远高于广义F中的数量，表明使用广义F进行信号检测而形成事件将更有效，能够减少虚假事件数量，也表明广义F检测器优于传统的STA/LTA检测器.目前开展的研究工作包括针对特定台站校正、噪声和信号模型校正，振幅校正等.
2)互相关检测算法
在信号检测理论中，相关器类似匹配滤波器，它是加性高斯白噪声中已知信号检测的最佳检测器(赵树杰和赵建勋, 2005).波形互相关技术在地震信号处理中也得到广泛应用(Gibbons and Ringdal, 2006, 2012; Harris and Paik, 2006; Schaff, 2008; Selby, 2010; Harris and Dodge, 2011)，特别是在低震级事件检测、震级判定以及事件聚类等方面有显著优势.由Bobrov等人(2014)将互相关技术引入IDC地震数据处理中，证明其能够明显降低全球范围内的地震监测阈值并提高震相及其参数计算的准确性.目前，IDC主要利用互相关方法进行大事件余震序列自动和交互分析处理.在IDC建立新的数据处理管道实现互相关处理，处理过程包括信号检测、震相识别和事件关联等步骤.根据预定义的一些高质量的主事件寻找相似的信号，并计算他们的参数，包括峰峰值、均方根幅值、方位角、慢度、F统计值和f-k值等，对检测出的信号进行关联形成事件列表XSEL，所有形成的事件都符合IDC事件定义标准.
主事件主要来自IMS台阵高质量的模版信号，这些模版信号通过时间窗口和滤波参数来定义，由于核查中的地震监测感兴趣的是全球范围内的小震级事件，因此检测频带应覆盖0.8~6 Hz的整个频谱，模板窗口不宜过长，依赖频带而变化.模板中应包含信号和信号前的部分噪声，对于低频信号，选择模板长度为6.5 s，包含信号到时前1 s的噪声，对于高频信号，选择模板长度为4.5 s.
采用标准化的互相关函数计算互相关系数，公式为
|
(6) |
其中Δt为采样间隔，N为采样点数，[v(tv), w(tw)]为两个波形数据间的内积.
通常情况下，同一台站互相关系数高的两个波形来自距离较近的两个事件，但也有一些情况距离较远的两个事件互相关系数也较高，IDC采用Gibbons和Ringdal (2006)提出的基于互相关系数f-k分析方法有效区分这一特殊情况.对于有效信号，方位角、慢度估计值和理论值接近，但对于虚假信号，由于其互相关迹只是噪声，则计算的方位角和慢度有很大的不确定度.基于互相关迹的f-k分析公式为
|
(7) |
其中Fi(f)表示在频率f下互相关系数迹的傅立叶变换，J表示台阵子台数，Sn·dnorthi+Sedeasti表示第i个子台相对于参考点的理论时间延迟.
采用F统计方法实现互相关迹的信号检测，基于标准的F统计公式，采用2 s时间窗口滑动计算，到时后4 s计算得到检测信号的最大F统计值，定义Fprob>0.5作为有效检测信号的阈值.
对于具有相同源地点但不同震级的事件，在Gibbons等提出的振幅比例因子的基础上，采用相对震级方法提高检测准确性.检测事件中的每个台站计算与主事件的相对震级，若得到的每个台站相对震级结果在一定偏差范围内，则认为该检测事件可靠.同时该方法可以减少高相关性的远距离事件情况对检测性能的影响，其定义为
|
(8) |
其中x, y分别是主事件和被检测事件的矢量数据.
采用局部关联(Local Association, LA)方法实现震相关联，首先由一个主事件发现所有有效的震相，再由每个震相到时减去主事件发生位置到达该台站的经验走时得到发震时间，当有3个或3个以上台站计算的事件时间在预先定义的时间间隔内，则形成一个假设事件.当所有震相的相对震级的变化范围在预先定义的区间内，则形成一个XSEL公报.关联过程中还包括由不同主事件形成事件的冲突解决.
通过对朝鲜宣称的两次核试验以及对2008年5月中国汶川地震的余震测试表明，基于波形互相关的检测技术能够显著降低全球事件的检测阈值，其检测出的震相和生成的事件具有更好质量和一致性，基本可以覆盖REB中95%的常规地震公报，能够提供更多的真实事件，可以增加1.5到2倍的有效REB事件，特别是对于大事件余震序列的处理性能有显著提升，分析员工作量能够降低两倍以上.采用局部关联方法代替全球关联，形成的事件具有更小的误差椭圆.测试也表明，互相关方法在台阵中具有更好的应用效果，结合新的关联技术NET-VISA可以帮助其在未知区域发现新的主事件.同时，也指出互相关方法在实时计算中需要大量计算资源，可以采用并行计算方法提高效率.
在互相关处理中需要对检测和关联中的一些参数进行调整.在检测方面主要依靠主事件和台站对进行参数调整，调整的参数包括检测频带、事件窗口，滤波参数，互相关定义，互相关系数迹的STA/LTA阈值，到时，F统计阈值，相对震级范围.局部关联主要依靠主事件进行参数调整，包括关联时间窗口，形成事件的台站数，台站权重、层级聚类参数及事件最小时间间距等.
2.2.2 事件关联在事件自动关联方面，原有IDC方法，通过分步计算信号特征，然后进行关联和定位，将关联问题分成几个子问题，通过流水线方式分步解决，前一步计算结果将直接影响到后面步骤计算准确性.由美国Nimar S. Arora和Stuart Russell提出的基于贝叶斯框架的关联方法NET-VISA (Arora et al., 2013)打破了这一模式，该方法首先建立一个基于地震学和观测经验的全球尺度的概率模型，该模型综合考虑了信号到时、方位角、慢度、事件震级、正确震相、误检信号和尾波信号等事件特征，然后采用一种简单的基于贝叶斯模型的推理方法，最终得到事件经纬度、深度和时间等信息.
NET-VISA首先是构建地震事件产生、传播和检测的概率模型，该模型可分为两个部分，一个是地震事件发生概率Pθ，另一个是由事件或噪声触发的检测信号概率Pφ，X为随机的所有可能发生的地震事件，Y为在台站观测到的可能信号，则Pθ(X)描述了参数化的事件先验概率，Pφ(Y|X)描述了信号的传播和测量概率.假设观测到Y=y，我们感兴趣的是事件参数，即后验分布P(X|Y=y)，根据贝叶斯准则：
|
(9) |
在CTBT监测中，我们主要是得到使上式概率最大的x*，即
|
(10) |
NET-VISA概率模型主要由事件、真实震相、误检信号和尾波信号等部分构成，见表 1.
|
表 1 NET-VISA中地震事件产生、传播和检测的概率模型 |
事件模型主要描述事件发生频度、地理分布、深度分布和震级分布等.其中事件发生概率和时间是由随机的泊松过程产生；事件位置分布是由一个均匀的密度分布(允许爆炸和新的地震活动发生在任何地方)和基于历史事件位置密度分布混合组成.事件震级模型是由著名的震级-频度关系Gutenberg-Richter分布构成(Gutenberg and Richter, 1954).公式(11)是事件发生概率模型, 其中λe是泊松过程参数, e是一系列事件, |e|表示事件数量，T为时间周期.公式(11)为
|
(11) |
真实震相模型依靠检测概率和震相属性而建立.其中检测概率主要和台站、震相、事件震级和震中距等相关，采用逻辑回归方法获得特定台站、震相的检测分布.震相属性主要包括到时、方位角、慢度、振幅等，到时、方位角和慢度残差都采用拉普拉斯分布.而振幅类似于检测概率，依赖事件震级、震深和震中距，采用带有高斯噪声的线性回归模型建立.
误检信号模型描述了由噪声产生的误检信号属性在各台站的分布.误检信号数量由随机泊松过程产生，其时间、方位角和慢度等在各自范围内属均匀分布.
尾波信号模型，台站处理中任何震相都可能产生尾波信号，尾波信号和虚假信号模型不一样，它和触发其产生的震相有较强的相关性.尾波信号和触发信号的时间延迟符合分布，尾波的方位角、慢度和振幅对数与触发震相的之差符合拉普拉斯分布.同时尾波信号模型不是特定台站的模型，而是基于震相的通用模型.
最终形成的完整概率模型公式为
|
(12) |
其中e是假设的事件，Λ是真实震相，ξ为误检信号，η是尾波信号，A是所有检测的震相集合.
第二步是完成事件推理.根据所有台站观测的信号结果，推断出式(13)最可能的解释，形成一系列事件，式(13)定义由式(12)给出.公式为
|
(13) |
事件推理是通过一系列过程来完成，包括生成过程、改善过程和消亡过程.采用相关的概率值来计算检测和事件的评分.推理算法可以看作是梯度法的优化形式，也被称为贪婪的局部搜索算法(Cormen et al., 2009)，其结构类似马尔科夫链蒙特卡罗方法(Markov Chain Monte Carlo, MCMC)(Andrieu et al., 2003)，且在速度和准确性方面优于它.
算法初始假设震相窗口中的所有新加入震相都是虚假的，首先“生成过程”在事件窗口中形成新事件，这些事件加入到假设事件中.然后循环N次以下过程：对于每个检测在检测窗口内进行“改善检测过程”，对于每个事件在事件窗口内进行“改善事件过程”.最后“消亡过程”删除虚假事件.为了简化这些过程中形成的假设事件概率之间比较的计算量，NET-VISA定义一个事件的评分Se作为两个假设事件的概率比值，公式为
|
(14) |
当一个事件评分Se < 1时，表明当前假设事件发生概率小于被替代的假设事件.对于每个检测同样定义评分标准Sd，其定义为
|
(15) |
它是该台站的检测关联到事件i的j震相的假设与检测属于默认分类的概率比值.根据定义，任何评分小于1的检测很可能是虚假检测或尾波震相.这两个评分定义在后续推理过程中起到关键作用.
“生成过程”是在一定时间窗口内基于所有检测信号生成可能的事件.首先根据台站处理得到的每个独立震相的相关参数，然后反演得到初始备选事件列表(依据慢度确定震中距，再结合方位角，粗略估算震中和发震时刻)，然后在每个备选事件一定范围内(经纬度5°范围，时间50 s范围)利用所有检测构建最佳可能事件.只要找到事件的最佳评分大于1，则搜索过程一直持续进行.
“改善检测过程”，对于检测窗口内的每个检测，在所有发生在该检测到时之前特定时间段内的事件及其震相中，为该检测寻找最佳的事件震相对.如果找到的最佳事件震相评分小于1或最佳事件震相已经属于其他具有更高评分的检测，则这个检测将被视为误检测信号或尾波信号.
“改善事件过程”，对于每个事件，在一定范围内(经纬度2°范围，深度100 km范围，时间5 s范围以及2个震级范围)随机选择100个点，寻找具有更高评分的事件，如果这样的事件找到，则改变事件参数.
“消亡过程”，任何评分小于1的事件将被删除，它所关联的所有检测被标注为噪声.
目前NET-VISA正在IDC中进行测试评估，采用事件审核公报REB作为参考事件集，将GA台网关联结果SEL3与NET-VISA台网关联结果进行对比，结果表明NET-VISA结果一致性更好，与现有系统相比在同样误警率下漏检率减少2/3，且NET-VISA方法形成的事件具有可溯源性.目前继续开展的研究工作包括：模型训练机制的确定，先验知识的建立(Arora et al., 2015)，实现在IDC管道处理中的应用，持续研究改进基于贝叶斯的推理方法, 以及NET-VISA方法在次声监测和水声监测中的应用(Arora et al., 2015)等.
2.3 区域走时模型信号传播模型在地震核查监测中十分重要，影响到事件检测、定位、识别和当量估计.区域地震震相(定义为事件-台站距离小于2000 km范围的主要震相)对于降低全球监测阈值具有重要作用，然而目前采用的1维地球模型进行区域走时预测准确度远低于对远震走时预测(Robert et al., 1998; Flanagan et al., 2007).由美国提出的区域地震走时模型(regional seismic travel times, RSTT)(Myers et al., 2010)目前也正在IDC进行测试，以改善区域走时预测和事件定位性能.RSST参数化模型近似表示3维地壳及上地幔结构，通过若干个无缝全球密布的具有速度剖面的三角形节点构成，采用插值方法生成一个连续的速度模型，每个节点的速度剖面由一个7层的地壳速度模型和一个线性梯度的地幔速度模型构成.采用线性梯度参数化上地幔速度模型，使实时计算Pn走时更快(约1 mm)，能够有效地应用于地震日常监测分析中，对于Pn走时采用文献(Zhao and Xie, 1993; Zhao, 1993)计算公式为
|
(16) |
其中d和s是莫霍面各点之间构成的离台站和事件最近的大圆路经上第i个分段的距离和慢度，α和β分别是台站和震源下方的地壳走时，γ用于解释俯冲射线，RSTT给出了这些参数的具体定义.
采用层析成像方法调整模型速度，包括平均地壳速度、莫霍面地幔速度及梯度地幔速度.RSTT应用、扩展和调整都十分方便，不需要进行预处理.RSTT采用横向变化的莫霍面速度模型和上地幔速度梯度方法近似上地幔结构，这种近似当震中距超过15度时(Pn/Sn最大距离为15度)有效性降低.该方法还可扩展以计算Sn、Pg和Lg走时.RSTT已用于ISC (International Seismic Centre)定位程序中，目前正通过IDC已有的SSSC (Source specific station correction)框架集成到IDC数据处理中(Given, 2014).RSTT最新发布的3.0.4版本可以实现对Pn、Pg、Lg、Sn震相的走时校正，同时该校正从震中距15度延伸到20度，对软件的计算性能也进行了改善.
在欧亚大陆和北非地区测试结果表明采用RSTT模型能够提高定位精度.在走时计算方面，英国国家数据中心Nippress和Bowers (2014)对采用RSTT校正方法的处理结果和原系统进行了比较，结果表明采用RSTT校正后的Pn走时残差比采用1维的IASPEI91全球走时模型，以及目前IDC系统采用SSSC校正的结果都具有更小的离散度，效果更好.同时采用该方法可实现对目前IDC运行系统中的所有IMS台站进行校正，而原有SSSC方法只能应用于29%IMS地震台站.
RSTT为所有的IMS台站提供了准确的区域走时计算的扩展框架，通过对模型的不断调整能提高走时预测的性能.模型调整的方法包括：在初始模型中结合近震和区域震的研究成果、采用标准事件(Ground-Truth)进行模型校正、采用层析成像方法调整模型速度、通过对非校准集中事件的重定位来验证模型有效性等.
3 结束语2010年起，PTS启动了虚拟数据研发中心(virtual Data Exploitation Centre, vDEC)(Vaidya et al., 2009)项目，该虚拟中心通过丰富的数据资源、分析专家、数据处理设施以及评估机制，为全球专家共同改进IDC数据处理算法，解决IDC面临的大规模数据处理挑战提供了研发和测试环境.结合vDEC, PTS正在从信号采集、信号传播路径模型及信号分析处理等方面开展研究以提高地震监测系统性能，台站的自校准和自诊断技术、基于物理和经验的地球模型、互相关信号检测技术、基于统计分类和推理的数据处理方法等都陆续在vDEC进行测试评估，测试结果表明这些技术的应用从某些方面能够提高目前系统监测性能.本文主要分析了目前PTS在地震数据处理软件中开展的部分研究工作，包括软件架构、信号检测、信号关联及传播路径模型等，希望对我国家数据中心数据处理软件的发展，以及对核查技术相关的研究工作起到借鉴作用.致谢感谢审稿专家提出的修改意见和编辑部的大力支持！
致谢 感谢审稿专家提出的修改意见和编辑部的大力支持！| [1] | Andrieu C, de Freitas N, Doucet A, et al .2003. An introduction to MCMC for machine learning[J]. Machine Learning, 50 (1-2) : 5–43. |
| [] | Arora N S, Mialle P. 2015. Fusion of Infrasound with Hydroacoustic and Seismic Data in NET-VISA[C].//CTBT:Science and Technology conference. Vienna, Austria. |
| [3] | Arora N S, Russell S, Sudderth E .2013. NET-VISA:Network processing vertically integrated seismic analysis[J]. Bulletin of the Seismological Society of America, 103 (2A) : 709–729. DOI:10.1785/0120120107 |
| [] | Arora N S, Storchak D, Russell S. 2015. Using International Seismological Centre data to build a prior of seismicity for NET-VISA[C].//CTBT:Science and Technology Conference. Vienna, Austria. |
| [5] | Bache T C, Bratt S R, Swanger H J, et al .1993. Knowledge-based interpretation of seismic data in the Intelligent Monitoring System[J]. Bulletin of the Seismological Society of America, 83 (5) : 1507–1526. |
| [6] | Backus M, Burg J, Baldwin D, et al .1964. Wide-band extraction of mantle P waves from ambient noise[J]. Geophysics, 29 (5) : 672–692. DOI:10.1190/1.1439404 |
| [7] | Blandford R R .1974. An automatic event detector at the Tonto Forest Seismic Observatory[J]. Geophysics, 39 (5) : 633–643. DOI:10.1190/1.1440453 |
| [8] | Bobrov D, Kitov I, Zerbo L .2014. Perspectives of cross-correlation in seismic monitoring at the International Data Centre[J]. Pure and Applied Geophysics, 171 (3-5) : 439–468. DOI:10.1007/s00024-012-0626-x |
| [9] | Cormen T H, Leiserson C E, Rivest R L, et al .2009. Introduction to Algorithms[M].3rd ed.. Cambridge, Massachusetts: The MIT Press . |
| [] | Coyne J, Bobrov D, Bormann P, et al. 2012. CTBTO:Goals, networks, data analysis and data availability[A].//Bormann P ed. New Manual of Seismological Observatory Practice 2 (NMSOP-2)[M]. Potsdam:Deutsches GeoForschungsZentrum GFZ. |
| [] | CTBTO. 2009. International Scientific Studies Conference 2009[C/OL]. http://www.ctbto.org/specials/the-international-scientific-studies-project-iss/. |
| [] | CTBTO. 2015. CTBT:Science and Technology 2015 Conference:Book of Abstracts[C/OL]. https://www.ctbto.org/fileadmin/user_upload/SnT2015/SnT2015_Abstracts.pdf. |
| [] | Douglas A. 2003. Seismometer arrays-their use in earthquake and test ban seismology[A].//Lee W H K, Kanamori H, Jennings P C, et al eds. International Handbook of Earthquake & Engineering Seismology[M]. London:Academic Press, 357-367. |
| [14] | Filson J .1975. Array seismology[J]. Annual Review of Earth and Planetary Sciences, 3 : 157–181. DOI:10.1146/annurev.ea.03.050175.001105 |
| [15] | Flanagan M P, Myers S C, Koper K D .2007. Regional travel-time uncertainty and seismic location improvement using a three-dimensional a priori velocity model[J]. Bulletin of the Seismological Society of America, 97 (3) : 804–825. DOI:10.1785/0120060079 |
| [16] | Gibbons S J, Ringdal F .2006. The detection of low magnitude seismic events using array-based waveform correlation[J]. Geophysical Journal International, 165 (1) : 149–166. DOI:10.1111/gji.2006.165.issue-1 |
| [17] | Gibbons S J, Ringdal F .2012. Seismic monitoring of the North Korea nuclear test site using a multichannel correlation detector[J]. IEEE Transactions on Geoscience and Remote Sensing, 50 (5) : 1897–1909. DOI:10.1109/TGRS.2011.2170429 |
| [] | Given J. 2011. A plan to re-engineer the IDC seismic, hydroacoustic, and infrasound processing software[C].//The 36th Session of Working Group B. Vienna. |
| [] | Given J. 2012. Progress in Re-engineering the IDC Seismic, Hydro-acoustic, and infrasonic processing system[C].//The 39th Session of Working Group B. Vienna. |
| [] | Given J. 2014. Status of Integrating RSTT into IDC processing[C].//Waveform Expert Group 42nd meeting of Working Group B CTBTO PTS. |
| [21] | Gutenberg B, Richter C F .1954. Seismicity of the Earth and Associated Phenomena[M].2nd Ed.. Princeton, New Jersey: Princeton University Press . |
| [22] | Harris D B, Dodge D A .2011. An autonomous system for grouping events in a developing aftershock sequence[J]. Bulletin of the Seismological Society of America, 101 (2) : 763–774. DOI:10.1785/0120100103 |
| [] | Harris D B, Paik T. 2006. Subspace detectors:Efficient implementation[R]. Lawrence Livermore National Laboratory Technical Report UCRL-TR-223177. |
| [] | Kværna T. 1995. Automatic onset time estimation based on autoregressive processing[R]. Semiannual Technical Summary, NORSAR, Sci. Report No. 1-95/96, Kjeller, Norway, 113-133. |
| [] | Kværna T, Doornbos D J. 1986. An integrated approach to slowness analysis with arrays and three-component stations[R]. Semiannual Technical Summary, NORSAR Sci. Report 2-85/86, Kjeller, Norway. |
| [] | Le Bras R, Swanger H, Sereno T, et al. 1994. Global association[R]. Science Applications International Corp. Tech. Rept. ADA304805, San Diego, California. |
| [] | Le Bras R, Wuster J. 2002. IDC Processing of Seismic, Hydroacoustic, and Infrasonic Data, Revision 1. IDC Documentation User Guides[M]. |
| [28] | Melton B S, Bailey L F .1957. Multiple signal correlators[J]. Geophysics, 22 (3) : 565–588. DOI:10.1190/1.1438390 |
| [29] | Myers S C, Begnaud M L, Ballard S, et al .2010. A crust and upper-mantle model of Eurasia and North Africa for Pn travel-time calculation[J]. Bulletin of the Seismological Society of America, 100 (2) : 640–656. DOI:10.1785/0120090198 |
| [] | Nippress S, Bowers D. 2014. Reviewed Event Bulletin Pn time residuals:Comparing SSSC and RSTT corrections[R]. Waveform Expert Group 42nd meeting of Working Group B CTBTO PTS. |
| [] | Peterson J. 1993. Observations and modelling of seismic background noise[R]. U.S. Geological Survey, Open-File Report, 93-322. |
| [32] | Robert E E, Van Hilst R D, Buland R P .1998. Global teleseismic earthquake relocation with improved travel times and procedures for depth determination[J]. Bulletin of the Seismological Society of America, 88 (3) : 722–743. |
| [33] | Rost S, Thomas C .2002. Array seismology:Methods and applications[J]. Reviews of Geophysics, 40 (3) : 2–1. |
| [34] | Schaff D P .2008. Semiempirical statistics of correlation-detector performance[J]. Bulletin of the Seismological Society of America, 98 (3) : 1495–1507. DOI:10.1785/0120060263 |
| [35] | Selby N .2010. Relative locations of the October 2006 and May 2009 DPRK announced nuclear tests using international monitoring system seismometer arrays[J]. Bulletin of the Seismological Society of America, 100 (4) : 1779–1784. DOI:10.1785/0120100006 |
| [36] | Selby N D .2008. Application of a generalized F detector at a seismometer array[J]. Bulletin of the Seismological Society of America, 98 (5) : 2469–2481. DOI:10.1785/0120070282 |
| [37] | Selby N D .2011. Improved teleseismic signal detection at small-aperture arrays[J]. Bulletin of the Seismological Society of America, 101 (4) : 1563–1575. DOI:10.1785/0120100253 |
| [38] | Tarantola A .2005. Inverse Problem Theory. Society for Industrial and Applied Mathematics, Mathematics and its Applications[M]. Philadelphia: SIAM . |
| [] | Vaidya S, Engdahl R, Le Bras R, et al. 2009. Strategic initiative in support of CTBT data processing:vDEC (virtual Data Exploitation Centre) (abstract DM-01/A)[C].//CTBTO International Scientific Studies. Hofburg, Vienna, Austria, 10-12. |
| [40] | Yao D F, Han Y M, Qiu Y F, et al .2009. The Practice of Software Testing Based on RUP[M]. Beijing: Tsinghua University Press . |
| [41] | Zhao L S .1993. Lateral variations and azimuthal isotropy of Pn velocities beneath Basin and Range province[J]. Journal of Geophysical Research, 98 (B12) : 22109–22122. DOI:10.1029/93JB02641 |
| [42] | Zhao L S, Xie J K .1993. Lateral variations in compressional velocities beneath the Tibetan Plateau from Pn traveltime tomography[J]. Geophysical Journal International, 115 (3) : 1070–1084. DOI:10.1111/gji.1993.115.issue-3 |
| [43] | Zhao S J, Zhao J X, et al .2005. The Theory of Signal Detection and Estimation[M]. Beijing: Tsinghua University Press . |
| [44] | 姚登峰, 韩玉敏, 邱云峰, 等.2009. 基于RUP的软件测试实践[M]. 北京: 清华大学出版社 . |
| [45] | 赵树杰, 赵建勋.2005. 信号检测与估计理论[M]. 北京: 清华大学出版社 . |