互相关算法在次声监测数据处理中的应用.pdf

互相关算法在次声监测数据处理中的应用 唐伟刘俊民王晓明邱宏茂盖磊王海军 禁核试北京国家数据中心，北京 100085 摘要 次声监测数据处理技术是我国履行全面禁止核试验条约必备技术之一。重点介绍逐次多通道互相关 PMCC 算法原理、 处理流程以及与之相关的数据处理系统， 并利用典型次声事件进行了实验研究。结果表明 PMCC 在处理次 声监测数据方面具有优异性能。 关键词 次声; 次声数据处理; 逐次多通道互相关 THE APPLICATION OF THE CROSS- CORRELATION ARITHMETIC IN INFRASOUND MONITORING DATA PROCESSING Tang WeiLiu JunminWang XiaomingQiu HongmaoGai LeiWang Haijun Comprehensive Nuclear-Test-Ban Treaty Beijing National Data Centre， Beijing 100085， China AbstractThe of infrasound data processing is a necessary technology in implementing the Comprehensive Nuclear- Test-Ban Treaty. The main theory of the progressive multi-channel correlation PMCC，data flow and its processing system are particularly introduced in this paper. Experimental research on ground － trouth infrasound events indicates that PMCC has a good quality in processing infrasound monitoring data. Keywordsinfrasound; infrasound data processing; PMCC 0引言 次声监测是全面禁止核试验条约 Comprehensive Nuclear-Test-Ban Treaty 简称 CTBT 规定的四种监测 技术 之 一 ［ 1］， 在 CTBT 国 际 监 测 系 统 International Monitoring System 简称 IMS 中有 60 个次声监测台 站， 我国境内的 IMS 次声台站有 2 个， 分别位于北京 和昆明， 我 国 禁 核 试 国 家 数 据 中 心 National Data Center 简称 NDC 负责接受、 转发和分析我国境内 IMS 次声监测数据。 当前国际数据中心 International Data Center 简 称 IDC 、 美国、 法国、 澳大利亚在次声监测数据处理 领域技术较为领先。IDC 拥有基于短时平均与长时 平均比 ［ 2］ Short Term Average/Long Term Average 简 称 STA/LTA 方法的 Libinfra 软件及法国根据逐次多 通道 互相关算法 ［ 3- 4］ Progressive Multi － Channal － Correlation 简称 PMCC 开发的次声数据处理软件， 最 近几年 IDC 一直在测试研究 PMCC 处理技术; 而美国 和法国也在利用 PMCC 处理技术进行相关研究并对 软件进行不断改进; 澳大利亚在次声研究方面有着较 长的历史， 在 IDC 的资助下， 目前正在开发更先进的 三维慢度空间相关处理软件 ［ 5］ 3D Slowness-Space Correlator 。 本文通过介绍 PMCC 算法、 处理流程、 数据处理 系统， 并对 PMCC 次声信号检测性能进行评估， 在此 基础上将 PMCC 引入到我国次声监测数据处理领域。 1理论 1. 1F － 统计 经带通滤波的数据可以利用 F － 统计估计信号 的空间相干性。F － 统计为相干聚束信号的单位时 间能量与单位时间内残存能力之比。其中残余能量 定义为 传感器接收信号与经适当时间延迟排列后所 得聚束信号之差的能量。如给定通道原始波形数据 与聚束信号波形数据相同， 则残余能量为零。残余能 量可以认为是传感器通道中非相干噪声信号的幅值。 F － 统计如式 1 所示 38 环境工程 2010 年 12 月第 28 卷第 6 期 FSsn， s e ≡ J － 1 [] J ∑ n0N －1 n n0 ∑ J i 1 xin mi sn， s e 2 ∑ n0 N －1 n n0 ∑ J i 1 xin mi sn， s e － 1 [] j ∑ J j 1 xjn mj sn， s e 2 1 在计算各通道 F － 统计时， 样本值的起始点为第 n0个时域样本点， 长度为 N。 mi sn， se为慢度坐标下 的时间延迟，i 为传感器序号，sn、 se分别为慢度在南 北和东西方向上的分量，J 为台站中传感器总数。如 果使用的数据中， 偏移量具有零均值， 并具有单位方 差， 此时 F － 统计可以变换为式 2～ 式 4 形式 FSsn， s e J － 1Bsn， s e J2－ Bsn， s e J － 1 Bsn， s e /J[] 2 1 －B sn， s e /J[]{} 2 2 Bsn， s e ∑ n0N －1 n n0 ∑ J i 1 xin mi sn， s e 2 3 mi sn， s e dn， isn de， is e sr 4 其 中 B sn， s e 为 总 时 间 平 均 聚 束 能 量。 mi sn， s e 慢度平面上的时间延迟为sr为数据采样 率，dn、 de传感器与台站中心间距。 当不存在噪声时，B sn， s e 趋于 J2，F － 统计值 则趋向于无穷大; 相反， 若不存在平面波信号时， 则数 据为随机噪声， B sn， s e 趋于 J， F － 统计值则趋于 1。 1. 2归一化互相关函数 标准化互相关函数 ［ 3］是对两个具有相同采样率 且按时间排列后不具有时间延迟的时序信号的估计。 当存在非随机信号时， 阵元间具有较高的互相关性。 以第 n0个样本点为起点， 长度 N 的信号 x i n， 由时域变换到频域如式 5 xi k ≡ 1 k∑ k－1 n 0 x′ in n 0 e －2π－ 槡 1n k/K 5 当 n0≤ n ≤ n0 N － 1时 x′i n xi n， 其 余情况下 x′i n 0。k 为快速傅氏变换的总长度。 由于卷积时存在时间延迟， 因此第二个信号必须具有 足够的长度以满足存在最大延迟 L 时信号卷积的需 要。最大延迟定义如式 6 L int2 SmaxΔrsr 0. 5 6 其 中 Δ r dn，i－ dn， j 2 de，i－ de， j 槡 2， d n，i、 dn， j、 de，i、 de， j分别为第 i、 j 个传感其与台站中心在南 北方向上和东西方向上的偏移量， 其单位为公里。为 Smax为第 i 和第 j 个传感器之间的距离， sr为信号的采 样率，int 为取整函数， 为了简化处理， 整个阵列使用 同一个最大的 Δr 值。第 j 个通道的信号必须有 N 2L 个数据值， 第 n0－ L 个样本点位数据段的起始点， 其频域如式 7 所示 xj k ≡ 1 k∑ K －1 n 0 x′j n n0 － Le －2π－ 槡 1n k/K 7 当 n0－ L ≤ n ≤ n0 N － 1 L 时，x′j n xj n，其他情况下 x′j n 0。 K 满足 K ≥ N 2L 且 K 2 p p 为正整数 。 式 8 为非标准化互相关函数珓χij m为频域数 据的 FFT 反变换， 其中 － L ≤ m ≤ L。 珓 χijm∑ K －1 k 0 x′i kxj ke 2π － 槡 1mk/K 8 归一化互相关函数如式 9 、 式 10 所示 χijm 珓 χijm pipj m 9 pi ∑ n0N －1 n n0 xi n 2; pj m ∑ n0N －1 n n0 xj n m 2 10 1. 3信号一致性 通过利用各传感器之间的时间延迟， 可以计算传 感器信号的一致性。算法如式 11 所示 c ∑ n－2 i 1 ∑ n－1 j i1∑ n k j1 mi， j mj， k mi， k 2 C3 槡 n 11 其中 n 表示有信号的传感器个数;C3 n表示有 3 个传 感器同时接收到信号的组合数;mi， j为第 i 个传感器 与第 j 个传感器间的时间延迟。 1. 4建立 Family 建立 Family 是在互相关检测的基础上， 将相似 检测结果进行归类， 并输出最终结果。通过该步骤可 以有效降低检测信号数量， 提高检测结果的可信性。 其过程如下 1 判 断 信 号 相 似 性， 若 相 似 信 号 数 大 于 建 立 Family 所需最小数目时， 建立 1 个 Family。 2 对于 Family 外的检测信号， 若其与 Family 中 最后一个检测结果相比时间延迟在50 s以内， 则比较 48 环境工程 2010 年 12 月第 28 卷第 6 期 该检测信号与 Family 中所有检测结果的相似性， 若 满足相似性要求则将该检测信号加入该 Family 中。 3 若相似信号数大于建立 Family 所规定的最大 信号数时， 则将大于最大数目部分的检测信号建立一 个新的 Family。 4 计算 Family 中所有检测信号的时间、 方位角、 慢度、 频率的平均值， 并作为最终检测结果。 图 1a 为 Family 聚类前信号检测点， 其中相连小 格为特征相似的检测信号， 图 1b 是经过 Family 聚类 后的信号检测结果。 1. 5PMCC 处理原理 根据处理原理以及对软件的分析， 将 PMCC 处理 流程归纳为图 2 所示。 图 1聚类前后检测信号变化情况 图 2 PMCC 处理流程 PMCC 处理过程如下 首先对原始信号进行数据 质量分析与滤波， 满足数据质量要求是进行后续处理 的前提条件。对于满足质量要求的数据， 按每 3 个监 测单元 或子台 进行组合创建子网络， 计算每个子 网络中两两传感器所记录信号的相关性、 一致性、 方 位角与慢度、 当前子网络中心与其余监测单元的间 距， 若该间距小于给定阈值， 则将对应监测单元加入 到当前子网络， 并判断是否满足一致性条件。若子网 络中监测单元数大于给定阈值， 则重新计算方位角与 慢度等波形参数， 以及 F － 统计值与信噪比， 并输出 初步检测结果。利用相似性规则对初步检测结果进 行聚类， 以聚类结果作为最终检测结果。 2实验环境建立及参数配置 建立次声监测数据处理平台用于安装、 调试和运 行 PMCC。在硬件方面， 需要具备供数据处理与存储 的服务器、 交换机、 终端工作站 或 PC 机 ; 软件运行 环境建设包括 配置系统服务器， Soalris10 操作系统 安装与调试， Clearcase、 Oracle9、 Tuxedo 等第 3 方软件 的安装配 置、 运 行 数 据 库 与 归 档 数 据 库 的 设 计 与 实现 ［ 6］。 在 PMCC 处理系统中需配置以下参数 1 检测阈值参数。该类参数直接决定 PMCC 检 测性能， 该类参数包括 子台数目阈值、 一致性阈值、 检测结果聚类阈值等。 2 数据质量参数。该类参数用于检查待处理次声 数据的质量， 若数据质量低于设定要求， 则不进行处理。 3 聚束参数。该参数用于定义聚束处理过程中 各子台组合方式。表 1 为 PMCC 部分关键检测参数。 58 环境工程 2010 年 12 月第 28 卷第 6 期 表 1 PMCC 关键检测参数 功能分类参数名称默认值参数使用说明 检测参数Nb_subnet4子网络数 windowlength50时间窗口长度， 单位为 s。 windowgap10两个连续时间窗口间隔， 单位为 s。 Threshold_consistency0. 3检测信号所允许的最大一致性阈值 时间延迟 ， 单位为 s。 Threshold_nbsensors3子网络中有效传感器数目小于该值时， 该子网络不被处理。 Q_parameter50当子网络外传感器与子网络中心距离小于 Qλ 时， 将该传感器加入当前子网络。 建立 FamilyThreshold_date 50与 Family 中最后一个信号的时间差距大于该值时， 不能将该信号加入 Family 中。 Threshold_lfam_min7Family 中包含的检测信号数少于该值时， 删除该 Family。 Threshold_lfam_max2000Family 中包含的检测信号数大于该值时， 关闭该 Family。 信号相似性判断Speed_transition2用于区分高速与低速信号， 为一速度单位为 km/s。 Time_tol50时间误差允许量， 与 Threshold_date 相一致， 单位为 s。 Freq_tol2频率误差允许量， 单位为 Hz。 Sp_tol10. 06低速信号波速误差允许量， 单位为 km/s。 Sp_tol20. 3高速信号波速误差允许量， 单位为 km/s。 Az_tol10. 1047低速信号方位角误差允许量， 单位为 rad。 Az_tol20. 4363高速信号方位角误差允许量， 单位为 rad。 3实验结果与分析 从 www. rdss. info 2006 年次声事件公报中， 选取 信噪比较高的 34 个事件作为典型次声事件， 用于 PMCC 检测性能的分析研究， 表 2 为实验中涉及的 IMS 次声检测台站列表。 表 2实验涉及的 IMS 次声台站 序号台站负责国台站地点纬度 / 经度 / 事件数 I07AU澳大利亚沃勒曼加南 19. 9东 134. 31 I08BO玻利维亚拉巴斯南 16. 2西 68. 51 I09BR巴西巴西利亚南 15. 6西 48. 01 I10CA加拿大伯耐湖北 50. 2西 96. 019 I26DE德国弗赖翁北 48. 9东 13. 72 I31KZ哈萨克斯坦阿克纠宾司克北 50. 4东 58. 06 I53US美国费尔班克斯北 64. 9西 147. 94 在次声数据处理平台中， 利用 PMCC 与 Libinfra 两种方法同时处理选定的典型次声数据， 表 3 为信号 检测结果比较。 表 3PMCC 与 Libinfra 处理结果比较 检测 方法 检测信 号数 检测典型 事件信号 检测 率 / 时间平均 误差 /s 方位角平 均误差 / PMCC15253294. 1255. 493. 52 Libinfra3711441. 18104. 497. 81 在选定的典型次声事件中， 玻利维亚的 I08BO 台站与巴西的 I09BR 台站， 同时记录到了 2006 年 4 月 18 日拉斯卡尔火山 Lascar 喷发产生的次声信号 拉 斯 卡 尔 火 山 靠 近 智 利 与 阿 根 廷 边 境，南 纬 23. 363， 西经 67. 725， 海拔5 397 m， 火山口直径 750 m， 深 300 m， 从 18 日上午 11 时 20 分开始发生了 4 次喷发 。根据 PMCC 信号检测结果， 确定事件的发 生地为 23. 38S， 67. 86W， 与火山口仅相距 14. 5 km。 表 4 为 I08BO 与 I09BR 台站参与定位的检测信号特 征参数。 表 4关联信号的特征参数 台站信号到时方位角 / 慢度 / s － 1 I08BO16 02 00186327. 3 I08BO16 03 50172324. 6 I08BO16 10 00171307. 6 I09BR17 18 40246305. 8 I09BR17 21 20245302. 3 I09BR17 24 30244306. 5 4结论 通过对 比 典 型 次 声 事 件 数 据 处 理 结 果， 结 合 PMCC 运行情况， 可以得出以下结论 PMCC 次声信号 处理方法具有较高的信号检测能力， 较小的信号检测 误差; 以 PMCC 为基础建立的次声数据处理系统， 具 有数据处理能力强， 性能稳定可靠， 可对监测数据进 行近实时快速处理的特点， 是目前性能最好的次声数 下转第 88 页 68 环境工程 2010 年 12 月第 28 卷第 6 期 现絮 状 沉 淀， 样 品 再 次 离 心， 并 将 清 液 完 全 转 入 100 mL容量瓶中， 用试剂水定容待分析。其他操作 步骤均同国标法。 2结果与讨论 2. 1精密性、 准确性 同一土壤样品国标法前处理滤膜过滤法与离心 法平行样、 加标回收结果见表 1。 表 1样品的测定和回收率 mg/kg 过滤法离心法 本底值 可溶性基 体加标 非可溶性 基体加标 本底值 可溶性基 体加标 非可溶性 基体加标 115. 1124. 4723. 2015. 4224. 5023. 94 214. 3422. 9625. 1015. 2024. 3024. 52 314. 1225. 2023. 8614. 3825. 2124. 87 413. 9224. 8024. 5014. 0223. 4623. 70 514. 8823. 3223. 4015. 2224. 7223. 92 615. 3623. 6823. 6014. 8323. 3824. 12 平均值14. 6224. 0723. 9414. 8424. 2624. 18 标准偏差 0. 530. 810. 660. 50 0. 660. 40 RSD /3. 63. 42. 83. 42. 71. 6 回收率 /94. 5 93. 294. 293. 4 由表 1 可知， 过滤法均值 14. 62 mg/kg， 离心法均 值14. 84 mg/kg， 两法相对偏差很小， 仅为 0. 75 ; 两 法可溶性基体、 非可溶性基体加标回收率基本相同; 离心法精密性好于过滤法。 2. 2不同样品对比 对不同样品用两种方法进行测定， 测定结果见 表 2。 表 2 表明 不同样品两方法相对偏差较小， 说明 离心法能适应不同土壤样品及沉积物的前处理。 表 2方法比对测定结果 mg/kg 样品编号土壤 1 号土壤 2 号土壤 3 号沉积物 过滤法13. 8621. 764. 841. 96 离心法15. 4222. 554. 322. 22 相对偏差 /5. 31. 8 5. 76. 2 2. 3样品预处理时间对比 实验了大量的不同土壤样品， 国标法一个样品消 解后， 第一次过滤时间一般在 3 h 以上， pH 调节如出 现絮状沉淀， 第二次过滤时间至少 2 h 以上， 加上样 品消解时间， 一个样品前处理过程至少 6. 5 h 以上。 采用离心法， 一次离心时间仅需 5 min， 一个样品前处 理过程只需 1. 5 h， 缩短了前处理时间。 3结论 快速、 准确的分析方法是当前环境监测技术的要 求。国标法固体废物六价铬分析样品前处理时间较 长， 一般实验室分析大批量的样品， 尤其是分析污染 纠纷、 应急监测需及时上报数据的样品具有较大的局 限性。本法固体废物样品前处理方法快速而准确， 精 密性好， 适应性广， 具有普及推广性。 参考文献 ［1 ］ GB5085. 3 － 2007 附录 T 固体废物六价铬分析的样品前处理 碱消解法［S］ . ［2 ］ GB /T1555. 4 － 1995 固体废物 六价铬的测定二苯碳酰二肼 分光光度法［S］ . ［3 ］ 欧伏平. 含泥沙水样总磷测定方法的研究［J］ . 环境工程， 2001 6 45- 46. 作者通信处易敏414000岳阳市岳阳大道环境监测中心 2010 － 05 － 07 櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅櫅 收稿 上接第 86 页 据处理系统之一。该方法可用于我国北京与昆明次 声台站实时次声数据的分析与处理， 具有很高的应用 价值。 参考文献 ［1 ］ 张利兴. 禁核试核查技术导论［M］ . 北京 国防工业出版社， 2005 156- 197. ［2 ］ Greg W Beall，David J Brown，Jerry A Carter. IDC processing of seismic， hydroacoustic， andinfrasounddata ［R ］ .Science Applications International Corporation，2001. ［3 ］ Alexis LePichon，YvesCansi.PMCCforinfrasounddata processing［R］ . The Newsletter of Subaudible Sound， 2003 1- 2. ［4 ］ Cansi Y. An automatic seismic event processing for detection and location the PMCC ［J］ . Geophy Res Lett，1995 22 1021 － 1024. ［5 ］ Kendra Biegalski，Jane Bohlin，Hallie Magyar. Database schema ［R］ . Veridian Pacific-sierra Research，2001. ［6 ］ Roger Bowman，Jerry Carter，Benjamin Kohl.Configurafion of PIDCdatabase ［ R ］ .ScienceApplicationsInternational Corporation，2001. 作者通信处唐伟100085北京市海淀区清河毛纺路 39 号禁核试 北京国家数据中心 E- mailtangwei851 126. com 2010 － 05 － 07 收稿 88 环境工程 2010 年 12 月第 28 卷第 6 期