0 引言

地震台网记录的数据质量是进行地震定位、震级测定及背景噪声层析成像的基础（孙业君，2015）。目前，山西测震台网仍主要依靠人工巡检监控台站数据状态、修改JOPENS系统内的台站仪器参数，不能实时关注各台站的背景噪声，并且台站维护人员在更换地震计或数据采集器后，需要由专业的系统运维人员在JOPENS系统的Web页面依次修改台站的仪器型号、电压灵敏度及各级仪器响应参数，这会导致台站故障修复不及时，从而影响地震波形的可用率。

许多研究者利用软件分别实现了对台站状态、数据质量、系统参数的监控及检查。陈智勇（2017）基于NetSeis/IP协议进行了JOPENS数据流状态实时监控；王俊等（2013）利用MATLAB平台研发CWQL软件实现准实时监控台站数据质量；安全等（2017）利用CWQL软件计算台站背景噪声，从而判断JOPENS系统内台站仪器参数的正确性。本研究利用JAVA开发山西地震台网运维智能管理系统，重点解决了台站运行监控和参数同步一体化的问题，实现了台站数据状态准实时监控及一键修改JOPENS系统台站参数的目标。

1 设计思路

山西地震台网运维智能管理系统实现实时更新JOPENS台站参数及数据状态监控的设计思路如下。

（1）通过对台站IP地址链路监控和连接JOPENS获取的实时数据流状态监控，判定台站是否有故障；若有故障，则维护人员填报完整的台站故障日志。

（2）根据台站故障日志的内容进行判定，是否有仪器类型和序列号的变化，在设备参数表中查找对应仪器的参数及仪器响应，并结合JOPENS系统中各类型仪器对应的仪器ID、Blockette、CoffBlob、ComplexErrorBlob、Response、Log_station等数据表的设置，自动更新JOPENS系统数据库中仪器响应表Response、通道表Channel、通道信息表Channel_info、台站信息表Station_info等表的各项参数值或ID值。

（3）对获取到的台站连续波形，根据McNamara等（2004，2005）提出的方法对数据进行预处理、加速度功率谱密度估算、加速度功率谱概率密度估算。

2 系统功能实现及展示

系统采用B/S结构（Browser/Server，浏览器/服务器模式），在MyEclipse开发环境中构建，以JAVA作为后台开发语言，MySQL作为后台数据库管理系统，采用经典3层架构框架进行模型建立和交互处理控制，通过com.entity建立实体，com.dao中包含实现系统主要功能所需方法类的实现，com.ping中包含多线程检测台站IP通信链路通断的方法类的实现，com.servlet通过request、response处理浏览器的请求及响应，调用相应的方法类并将结果返回jsp页面，com.util实现与数据库建立连接及查询、修改等操作，com.filter实现过滤器功能。台站噪声计算时主要通过python调用obspy包实现连续波形的解析来获取头文件和数据波形。

2.1 运行监控

运行监控包括通信链路、实时数据流的监控。

2.1.1 通信链路监控

通信链路监控的基本功能是通过JAVA调用InetAddress.isReachable（）方法判定台站IP的连通状态，同时，采用Executors.newFixedThreadPool（）创建与台站数量相同的线程池，从而实现一次性ping所有台站IP。实现的主要代码如下：

public static void main(String[] args) {

StaDaosd =new StaDaoImpl();

      List<sta>staAll = sd.getstaAll();

inttotalIps =staAll.size();

ExecutorService executor = Executors.newFixedThreadPool(totalIps);

      List<Future<PingResult>>list = new ArrayList<Future<PingResult>>();

      Callable<PingResult>callable = null;

for(int i=0; i<totalIps; i++){

callable = new PingTask(staAll.get(i).getip());

          Future<PingResult>future = executor.submit(callable);

list.add(future);

      }

for(Future<PingResult>fut : list){

try {

System.out.println(new Date()+ “: : ”+fut.get());

          } catch (Exception e) {

e.printStackTrace();

          }

      }

executor.shutdown();

}

publicPingResult call() {

InetAddressinet = null;

try {

inet = InetAddress.getByName(ipAddress);

intresultCode = inet.isReachable(5000)? 0 : -1;

return new PingResult(ipAddress, resultCode);

   } catch (Exception e) {

e.printStackTrace();

return new PingResult(ipAddress, -1);

   }

}

后台运行结果见图 1。

2.1.2 实时数据流监控

实时数据流监控主要是通过连接JOPENS系统的流服务器，调用stat monitor命令，根据返回值判断每个台站的数据流状态，返回值1代表台站有数据，0代表台站无数据，进而排除通信链路正常但数据采集器无输出的情况。主要实现代码如下：

Socket sk =new Socket(“xx.xx.xx.xx”, xxx);

PrintWriter pw =new PrintWriter(new BufferedWriter(new OutputStreamWriter(sk.getOutputStream())), true);

BufferedReader in =new BufferedReader(new InputStreamReader(sk.getInputStream()));

System.out.println(“~~连接成功~~!”);

pw.write(“user ld\n”);

pw.flush();

pw.write(“pass ld\n”);

pw.flush();

pw.write(“stat monitor\n”);

pw.flush();

String reply=null;

while(!((reply=in.readLine())==null)){

        System.out.println(“接收服务器的信息: ”+reply);

        intlen=reply.length();

        if(reply.substring(0, 2).equals(“SX”) & & reply.substring(len-1, len).equals(“0”)){

              String sta =””;

              Streamstatusss=new Streamstatus();

              sta=reply.substring(3, 6);

              System.out.println(sta+”台站断记！”);

              ss.setstacode(sta);

              ss.setstaip(null);

              ss.setstart_time(new Timestamp(System.currentTimeMillis()));

              ss.setend_time(new Timestamp(System.currentTimeMillis()));

              stalist.add(ss);

        }

}

后台运行结果见图 2。

2.2 参数同步

参数同步包括在JOPENS中对台站参数和台站维护日志进行更新。

2.2.1 更新台站日志

根据系统中维护人员上报的台站日志信息，自动同步更新JOPENS系统中台站日志表Log_station，其主要代码如下：

publicbooleanaddjopensstalog(Logstalogsta) {

          boolean flag = false;

          DBconjopens.init();

          int i =DBconjopens.addUpdDel(“insert into Log_station(Net_code, Sta_code, Net_cname, Sta_cname, Operator, Start_time, End_time”+

            “Occur_type, Treat_type, Log_content, Create_time, Remark) values(‘”+”SX”+”’, ’”+logsta.getstacode()+”’, ’”+”山西”+”’+” +”’”+logsta.getstaname()+”’, ’”+logsta.getoperator()+”’, ’”+logsta.getstart_time()+”’, ’”+logsta.getend_time()+”’, ’”+logsta.getoccur_type()+”’, ’”+logsta.gettreat_type()+”’, ’”+logsta.getlog_con()+”’, ’”+logsta.getcreat_time()+”’, ’”+logsta.getremark()+”’)”);

          if(i>0){

                    flag = true;

          }

          DBconjopens.closeConn();

          return flag;

   }

2.2.2 更新台站参数及仪器响应

通过台站维护日志信息中备注Remark的信息，判断是否有仪器更换。若序列号发生变化、型号未变，则需要根据序列号查找地震计电压灵敏度，在JOPENS系统数据库中更新0级响应中Blockette058的系统响应灵敏度、1级响应中Blockette058的地震计电压灵敏度；若地震计型号发生变化，则需要将JOPENS中台站信息表Station_info中对应的Instrid字段进行更新，同时更新0级响应中Blockette058的系统响应灵敏度、1级响应中Blockette053的地震计零极点、Blockette058的电压灵敏度；若数据采集器型号发生变化，则需要将JOPENS中台站信息表Station_info中对应的datareid字段进行更新，同时更新2、3、4级响应中Blockette054、Blockette057、Blockette058对应的仪器响应。

各级仪器响应更新后，获取最新的各级响应对应的Blockette表的ID号，更新仪器响应表Response中三通道相应0、1、2、3、4级数对应的Blockette053、Blockette054、Blockette057、Blockette058的ID，更新Channel中对应的Response的ID、Channel_info中对应的Channel的ID、Station_info中对应的Channel_info的ID，并将最新的响应参数按照固定格式写成新的xml格式响应文件来替换Station_info中Response字段值。

部分功能实现代码如下：

① 地震计型号改变，修改地震计ID及响应参数

   InstrdicDaoinsd=new InstrdicDaoImpl();

   List<Instrdic>inslist=insd.getinsAll();

   Instrdic ins=new Instrdic();

   int ID=0;

   for(int k=0;k<inslist.size(); k++){

   if(m2.equals(inslist.get(k).getInstr_model())){

                  ID=inslist.get(k).getId();

   }

}

jsid.update1(stacode, ID);

jsd.update1(staid, ID);

② 获得当前地震计型号对应的零极点

s153=b53d.getb53(s[2]);

String zero_str=s153.getZero();

String pole_str=s153.getPole();

Long zid=b53d.insterzreo(zero_str);

Long pid=b53d.insterpole(pole_str);

zero z=new zero();

pole p=new pole();

List<zero>zero_list=new ArrayList<zero>();

List<pole>pole_list=new ArrayList<pole>();

String[] read_zeroLF=rd.read_LF(zero_str);

String[] read_poleLF=rd.read_LF(pole_str);

for(int x=0;x<read_zeroLF.length; x++){

      z=rd.readzero(read_zeroLF[x]);

      zero_list.add(z);

}

for(int x=0;x<read_poleLF.length; x++){

      p=rd.readpole(read_poleLF[x]);

      pole_list.add(p);

}

③ 向xml格式的response文件中写入内容

stage0.setStage_sequence(0);

stage0.setBlockette058(b058_s0);

stagelist.add(stage0);

stage1.setStage_sequence(1);

stage1.setBlockette053(b053_s1);

stage1.setBlockette058(b058_s1);

stagelist.add(stage1);

stage2.setStage_sequence(2);

stage2.setBlockette054(b054_s2);

stage2.setBlockette057(b057_s2);

stage2.setBlockette058(b058_s2);

stagelist.add(stage2);

stage3.setStage_sequence(3);

stage3.setBlockette054(b054_s3);

stage3.setBlockette057(b057_s3);

stage3.setBlockette058(b058_s3);

stagelist.add(stage3);

if(b054_s4!=null){

stage4.setStage_sequence(4);

stage4.setBlockette054(b054_s4);

stage4.setBlockette057(b057_s4);

stage4.setBlockette058(b058_s4);

stagelist.add(stage4);

}

response.setNet_code(netcode);

response.setLoc_id(locid);

response.setSta_code(stacode);

response.setChn_code(chncode);

response.setStage(stagelist);

String xmlres=rd.creatxml(response);

cid.updateRes(id_cinfo, xmlres);

2.3 噪声计算

噪声计算主要是通过对原始背景噪声波形数据按照特定长度分段，进行数据预处理，采用直接傅里叶变换法估算加速度功率谱密度PSD，进而得出加速度功率谱概率密度函数PDF的估计。计算方法和流程如下（McNamara et al，2004, 2005）。

2.3.1 数据预处理

将获取到的台站连续波形解析为1 h的时间序列段，并将时间段进行50%的重叠，每小时的波形再次分为13个段，进行75%的重叠，并将13个段中每个片段的样本数均截断为2的次幂；将数据转换为零平均值，通过平均斜率方法消除长期线性趋势；在每个被截断和去趋势的时间序列段的末端采用10%的余弦渐变抑制傅里叶变换结果的旁瓣泄露；通过原始傅里叶变换中的总功率与平滑滤波器中的总功率之比1.142 857进行量化，用于校正最终频谱中的绝对功率。

2.3.2 加速度功率谱密度

通过原始数据的有限范围傅里叶变换计算PSD，周期时间序列y(t)的有限范围傅里叶变换如下

$ Y\left(f, T_r\right)=\int\limits_0^T y(t) \mathrm{e}^{-i 2 \pi f t} \mathrm{d} t $

(1)

其中，T_r为时间序列的长度，总长为2¹⁵ = 819.2 s；f为频率。对于离散频率值f_k，傅里叶分量定义为

$ Y_k=\frac{Y\left(f_k, T_r\right)}{\Delta t} $

(2)

其中，f_k= k/NΔt，k = 1，2，⋯，N-1；采样间隔Δt为0.025 s；N为每个时间序列段的采样数，N=T_r/Δt。

因此使用式（2），总功率谱密度估计值PSD为

$ P_k=\frac{2 \Delta t}{N}\left|Y_k\right|^2 $

(3)

从式（3）可以看出，总功率P_k只是幅度谱的平方，归一化因子为2Δt/N。在计算了所有13个段PSD估计之后，对q= 13个独立时间段的功率进行平均，其中，每个时间段的长度为T_r。PSD的最终平滑估计如下

$ P_k=\frac{1}{q}\left(P_{k, 1}+P_{k, 2}+\ldots+P_{k, q}\right) $

(4)

其中，P_{k, q}为第q个时间段的频率f_k的原始估计。

PSD的估计值将通过参数1.142 857进行校正，即P_k = P_k×1.142857，然后将PSD估计值除以仪器传递函数来消除地震计的仪器响应

$ \mathrm{PSD}_k=\frac{P_k}{\left|H\left(f_k\right)\right|^2} $

(5)

由于数字地震仪观测的物理量为地面运动速度，因此，由式（5）计算得到的PSD_k为速度功率谱密度PSD_v，其单位为（m/s）²/Hz；加速度功率谱密度PSD_a的单位为（m/s²）²/Hz，PSD_a与PSD_v间的关系如下

$ \operatorname{PSD}_{a k}\left(f_k\right)=4 \pi^2 f_k^2 \operatorname{PSD}_{v k}\left(f_k\right) $

(6)

最后将结果由（m/s²）²/Hz转为为分贝值dB

$ \mathrm{PSD}_{a k}=10 * \lg \left(\mathrm{PSD}_{a k}\right) $

(7)

2.3.3 加速度功率谱概率密度

应用概率密度函数方法将所有时间序列段的PSD值，合成为PSD的概率密度函数分布图。首先，为了对PSD进行足够的采样，以1/8倍频程间隔获取全倍频程平均值，因此，功率需要在短周期（高频）角T_s和长周期（低频）角T₁之间平均，T₁=2 * T_s，中心频率T_c = sqrt(T_s * T₁)，是倍频程内的几何平均周期，T_s以1/8倍频程递增，即以T_s = T_s * 2^0.125计算下一个时序段的平均功率。重新计算T₁和T_c，在下一个周期范围T_s至T₁中平均功率，然后继续进行处理，直到给定原始数据的时间序列窗口长度T_r/10，对于每小时的PSD估计重复此过程，最终为每个台站生成数千个平滑的PSD估计。其次，使用概率密度函数绘制每个周期的功率分布，对于给定的中心周期T_c，PDF可以估计为

$ P\left(T_{\mathrm{c}}\right)=N_{P T_{\mathrm{c}}} / N_{T_{\mathrm{c}}} $

(8)

其中，N_PTc为给定中心周期T_c落入1 dB功率仓P（–200— –80 dB）的频谱估计数量。N_Tc为具有中心周期T_c的所有功率的频谱估计总数。

概率密度函数PDF的计算使用python语言编程，调用obspy开源库的PPSD类实现，部分实现代码如下：

import obspy

from obspy.signal import PPSD

st = obspy.read(“SX.ANZ.00.BHZ.D.2019.296”) //读取波形文件

inv = obspy.read_inventory(“RESP.SX.ANZ.00.BHZ”) //读取响应文件

ppsd = PPSD(st[0].stats, inv) //创建PPSD实例

ppsd.add(st) //加入波形数据，进行计算

ppsd.plot(“SX.ANZ.00.BHZ.D.2019.296.png”) //绘图，并存入图形文件

执行结果如图 3所示。

3 总结

山西地震台网运维智能管理系统实现了对测震台站的通信链路和实时数据流的监控，可以综合判定台站数据中断的原因，并且维护人员填报维护日志后，如果有仪器更换，系统自动将对应台站的仪器和响应参数等相关信息更新至JOPENS数据库，避免了人为原因造成的参数更新不及时。系统自动计算噪声加速度功率谱密度PSD和加速度功率谱概率密度PDF估算，可对台站噪声实时监控，保障了山西地震台网台站数据质量。

正在建设的国家地震烈度速报与预警工程对台站参数的时效性和准确性提出了更高的要求，今后将重点针对该系统的台站参数模块，补充完善地震计、加速度计及数据采集器的仪器模板，优化参数更新流程。