3.2 多线程及网络实时监控
Matlab可以通过parfeval来进行多线程运算， parfeval是Parallel Computing Toolbox提供的一种异步并行计算机制，允许在后台并行执行多个函数，而无需阻塞主线程。本软件利用parfeval进行多线程任务，保证软件后台进行自动功率谱密度分析，不影响前台数据实时监控任务。parfeval基本语法为： futureObj=parfeval （@myFunction， nargout， arg1，arg2，...），语法中@myFunction为异步执行函数，arg1，arg2，...等为myFunction的输入参数，软件中重写myFunction函数执行功率谱密度计算及时频图绘制代码。重写myFunction函数后保存，由主程序中parfeval调用执行即可实现新线程的功率谱密度数值求取计算。
对于网络实时监控，Matlab执行操作系统ping命令，利用system（command）函数，该函数通过调用系统shell来执行操作系统命令。操作后将结果返回到status变量与cmdout变量，执行代码为system（ping 10.*.*.*）返回，可以通过灵活设置系统ping命令自身参数如发送数据包大小，网络响应时间等方式实时监控VP宽频带倾斜仪网络情况。
[status cmdout]=system（˛ping-w 50 10.*.*.*˛）；%执行Ping VP倾斜仪IP地址
3.3 实时数据流可视化绘制
前兆仪器有对应网络页面，可通过HTTP协议向仪器网页站点发送请求方式获取HTML文本内容，再清洗HTML文本内容获取观测数据，最后绘制数据流曲线，目前前兆仪器监控软件均遵循该流程。但不同编程语言依各自特色有不同实现方式，各有优劣。python使用Requests命令及网络爬虫获取清洗数据， VB使用Inet网络控件获取数据[13]。本软件时域监控数据可视化借鉴沿袭以上思路，在Matlab语言中首先通过webread函数获取VP宽频带倾斜仪网页http：//10.*.*.*/cgi-bin/realtime.cgi字节数据，再通过正则表达式函数regexp（）筛选字节并转换为实数值。最后利用Matlab动态数组、plot（）绘图函数结合其循环语句来实现数据流动态实时重绘。以下仅通过给出简单单台监控代码示例阐述上述流程，多台监控拓展编程即可，另实际程序中还涉及到网络中断台站地图信息显示、数据流长度控制（防监控时间过长，数组溢出）、阈值报警等方面，在此不再赘述。
NSstream=[]；EWstream=[]；%定义初始空数组用于存储VP数据
ConnectIp=[˛http：//10.*.*.*/cgi-bin/realtime. cgi˛]；%定义监控VP仪器DNS地址字符串
RealTime=inf；
while RealTime>0 %while循环持续获取数据并绘制曲线
try
[status cmdout]=system（˛ping -w 50 10.*.*.*˛）；%ping命令设置50ms网络响应时间
webdata=webread（ConnectIp）；
catch err %获取网络中断error信息
fprintf（˛Error：%s\n˛，err.message）；
end
if status==0；
%通过正则表达式函数regexp（）筛选字节数据分别赋予NStmpStream、EWtmpStream数组%
index=strfind（webdata，˛tr˛）；
str1=webdata（index（3）：index（4））；
str1key=regexp（str1，˛size="6">（\S+）mV</font>˛，˛tokens˛）；
tmp=str1key{1，1}；
NSvalue=str2num（tmp{1，1}）；
str2=webdata（index（5）：index（6））；
str2key=regexp（str2，˛size="6"><b>（\S+）mV</b>˛，˛tokens˛）；
tmp=str2key{1，1}；
EWvalue=str2num（tmp{1，1}）；
NSstream（end+1）=NSvalue；%数组尾部加入数据
EWstream（end+1）=EWvalue；%数组尾部加入数据
%通过正则表达式函数regexp（）筛选字节数据分别赋予NSstream、EWstream数组%
subplot（2， 1， 1）， plot（NSstream）， title（˛北南向数据˛）；
subplot（2，1，2），plot（EWstream），title（˛东西向数据˛）；
else status==1；%网络中断则设置数据尾部加入NaN空值
NSstream（end+1）=NaN；
EWstream（end+1）=NaN；
end
end
3.4 PSD功率谱密度计算去趋势项及固体潮
如前所述，VP宽频带倾斜仪数据PSD功率谱密度计算需进行去趋势项预处理，去掉一般长趋势可以使用最小二乘拟合方法，首先对数据最小二乘法拟合得到趋势曲线BaseLine，然后再对原始数据减去拟合趋势曲线BaseLine。在Matlab中最小二乘拟合去信号趋势项有相应的库函数，该库函数为detrend（），直接调用即可。对于非线性趋势可通过poly fit（）与 polyval（）进行拟合得到趋势曲线BaseLine值，再用原信号减掉趋势曲线。
LengthOfVPseconddata=length（VPseconddata）；
t=[1：1：LengthOfVPseconddata]；
p=polyfit（t， VPseconddata， 5）；% 5为拟合阶数
xtrend=polyval（p，t）；%拟合获取数据长周期趋势项
VPseconddataOfDetrend=VPseconddata-xtrend；%原始数据去趋势
除了长趋势项，VP数据还包括固体潮信息，可采用小波分解去固体潮[14]。小波分解可将原始信号S分解成趋势部分A和细节部分D，要去掉VP数据中固体潮信息，可用原始数据扣除小波分解后的日波、半日波、1/4日波等固体潮所对应的细节部分。亦可以直接通过扣除小波趋势部分A及细节部分D来达到同时去除趋势项及固体潮，笔者对比小波去趋势项与最小二乘拟合去趋势项，小波去趋势更彻底，但也会使谱密度处理结果较之最小二乘拟合方式丢失部分低频信息，软件设置了预处理选择控件，可以根据功率谱密度计算需要自行选择。在Matlab中提供了小波分解和重构函数，代码如下。
[C，L]=wavedec（VPseconddata，16，˛db4˛）；% VPseconddata数据进行db4小波基分解
A16=wrcoef（˛a˛，C，L，˛db4˛，16）；%重建小波16层趋势
D16=wrcoef（˛d˛，C，L，˛db4˛，16）；%重建小波16层细节
