方位角检测程序的设计及应用
2021-07-24陈家樑张宝剑
陈家樑，张宝剑
(福建省地震局,福建福州，350003)
0 引言
福建地震台网中心从1971年开始建设，从最初3个测震台站的模拟观测及人工震中定位，发展到2020年底88个测震台的数字观测台网，得益于国家的发展壮大。福建地震台网包含了测震台、强震台、烈度台及GNSS台网，台站密度走在了全国前列。其中测震台网中部署有宽频带地震计的台站占据84%，为组建虚拟台网我中心接入了邻省台站及16个台湾台，福建本省测震台站的平均台间距约37千米，震级的观测范围为0-6级。台站数量的增加提高了地震定位的精度，而地震数据的可靠性决定了地位精度的准确性，因此有必要对地震数据质量进行多方面的监测。方位角应用于包括SKS剪切波分裂，接收函数研究等多项科学研究，如果方位角存在偏差将会影响相关研究结果的可靠性。以往都是在相关研究出现异常时才会发现台站有故障，这样的故障维护就会出现延迟，且会造成该时间段内的数据不可用。因此需要有专门软件对台站方位角进行实时监测，发现问题及时处理，这样不仅可以提高故障修复的时效性，且可以让仪器维护人员的工作更具有针对性。
图1 福建台网及台湾交换台站分布
1 程序设计
方位角的检测需要用到震中距30度以上的远震，所以该程序为准实时系统。根据具体的操作步骤，程序分为地震事件选取模块、数据预处理模块、台站属性读取模块、最小能量计算模块、相关系数计算模块及虚拟台站构建模块。为简化代码数量，减少程序运行中所占用内存，本程序使用Python中的Obspy这个开源的地震数据处理框架。
1.1 地震事件选取模块
本程序中计算方位角偏差的方法为利用P波质点的运动特性法，需要震中距在30度以上的远震，因此在准备地震事件时需要计算震中距是否达到30度。福建本省面积较小，对于远震而言省内不同台站对于同一地震事件的震中不会有太大差别，因此我们选取位于福建较为中心的永安燕西作为所有台站的参照台，只要计算该台与地震事件的震中距即可，这样可以减少运算量又不影响数据的可靠性。计算震中距我们引入Obspy.geodetics中的gps2dist_azimuth及kilometers2degrees两个方法，前者用于计算震中距，其返回值是一个数组，震中距为第一个元素且单位为千米，所以需要使用kilometers2degrees将单位从千米转化为度。程序将把符合震中距要求的地震事件三要素及震中距输出至指定路径下的TXT文件中。
图2 程序架构
1.2 数据预处理
因宽频带地震计中有60秒与120秒甚至更高的频带，为保证带宽一致所以需要进行带通滤波。台站的采样率为100HZ，为减少计算量程序采取降采样处理。利用Obspy.Stream类循环读取地震离线波形数据，全部读取后利用该类中的resample（sampling_rate）方法即可进行降采样处理。Obspy.Stream.filter可以对数据进行带通滤波处理，用该方法离线波形数据进行带通滤波处理的具体参数设置如下所示 filter(‘bandpass’,freqmin=0.02,freqmax=0.2),参数表示使用带通滤波处理数据且设置了低频与高频的具体数值。
1.3 台站属性读取
台站属性指的是台站名、经纬度及高程，这些数据均存储于MySQL数据库中。对于省级地震台网中心，有固定的数据库用于配置各台站经纬度、高程及灵敏度等参数，所以我们用xml文件配置好数据库的地址等，程序就可以通过xml文件进入数据库进行台站属性的读取。Xml文件的内容形式如下所示：
该模块中我们创建一个专门读取xml文件的类，在Python中引用xml.dom.minidom类包，就可以按照该类读取xml文件的步骤方法，将xml文件中的各配置内容读取到指定的数组中并返回该数组。
获取上述数据库的配置信息后连接MySQL数据库，宽频带地震计的通道命名有显著要求，即开头字母为B，其三通道的命名为BHN、BHE、BHZ。因此在选取宽频带台站的相关属性时对数据库的操作语句可写作“select Sta_code from Channel_info where Chn_code like ‘%BHE%’”,其 中Sta_code为台站名，Channel_info为表名称，Chn_code为通道存储字段名。通过以上操作可以获取全部宽频带地震计台站名，而后运用循环方法获取所有的台站属性后封装到dictionary中返回。
1.4 最小能量计算
通过以上以上步骤计算每个台站对应每个地震事件的P波到时，选取到时前后五秒内经过旋转的能量，再以1度为间隔进行180度的旋转计算每一度对应的能量，再进行对比选出最小能量对应的旋转度数，认为此时对应的角度为与正北方向的角度偏差。再将此时的偏差角与经过经纬度经过计算出的理论偏差值进行对比，即可得出此台对应该地震事件时间点上的角度偏差。将每个地震事件的偏差值输出至String类型的数组，计算完所有事件的偏差值后，再利用“SNR”法计算所有事件中每一度的能量值后选出最小值对应的偏差角，即认为是最终的地震计与正北方向的偏差值。
图3 时间序列图
为更清晰的显示每个事件对应的方位角偏差，在此步骤中将所得到的数据进行绘图显示。因利用了GMT绘图软件所以需要在主机上安装GMT5.0以上版本，否则会出现绘图命令错误或不兼容的情况。时间序列图的纵轴单位为的形式，所以需要用到GMT中的转义字符“@～\152\244\050\260\051@～”。将GMT绘图命令以String形式写入bat文件后，执行该文件即可进行绘图操作，绘制后的图片将按xml文件中配置的路径存储。
1.5 虚拟台阵
为更加精准的分析异常方位角偏差台站，该程序还设计了异常台站与正常台站所构建的虚拟台阵对应通道的相关系数分析模块。此模块以故障台站为中心，计算所有台站与此台的台间距，并选取台间距最近的15个正常台站构成虚拟台阵。线性叠加台阵内所有台站相同通道的数据后，与故障台相同通道进行相关系数计算，如果故障台不出现仪器方向摆放错误或极性反转等重大问题，则其相关系数会接近1相反则会接近0或负数。其相关系数的计算引入scipy.stats类中的pearsonr方法，可以迅速计算出两个相同通道的相关系数。为了便于观察，此模块同样具备输出波形形状及相关系数值的功能。
图4 上图为故障台，下图为正常台与虚拟台阵的相关系数波形图
2 程序在福建台网中的应用
通过该程序对福建地震台网所属宽频带地震计进行计算，地震事件为2015年1月到2019年8月之间的远震，得到74个宽频带地震计台站中有72个台的方位角偏差值在-7°～7°之间，这样的偏差值在该方法的理论中可以忽略不计，因此福建台网宽频带地震计台站的方位角总体合格率在97%如图5所示。其中两个存在较大偏差值的台站均为BBVS-60地震计,此类型地震计没有明显指北标志，所以导致维护人员在更换过程中无法准确的进行指北操作，经过实地校正后已使方位角的偏差值减小到要求值1°以内。
图5 方位角偏差台站分布
3 结论
通过对程序模块化开发，可以理顺每个模块的具体功能，使程序整体研发更清晰，更易调试出现问题的部分代码进行小范围的修复，而不至于出现牵一发动全身式的更改。因程序为处理地震数据而开发，所以使用到Obspy库，因其几乎支持地震学界内通常使用的所有波形格式的读写，可简化大量的数据格式之间转换的代码，同时他还集成了大量的地震学及数学所用的专有库，对于开发者只需要对其进行引入就可使用。程序应用在福建台网中心后，经过实地校核发现该程序的计算结果可靠。该程序的应用弥补了福建台网中心在方位角检测方面的不足，且提高了方位角偏差的检测效率及修复故障的时效性。
