低采样率时间脉冲插值分析方法及应用

程煜博 赵立军

摘 要：地震数据采集器采样时间精度是评价设备性能的重要指标，按照行业标准应使用标准时间源输出的时间脉冲信号采样数据进行分析计算，但由于地震数据采集器的采样率较低，导致时间同步误差测试分辨力不高。为提高低采样率时间脉冲分析能力，本分研究并实现了基于柱状图分析的脉冲信号高低电位识别算法、基于状态信息的脉冲信号上升沿识别提取算法和基于频率域补零的升采样率插值算法，基于上述算法对地震数据采集器采集的时间脉冲信号进行波形升采样恢复，提高时间精度测试分辨力的同时，规范了该项目测试数据处理方法和流程。

关键词：地震数据采集器 时间误差 升采样率 测试分辨力

Analysis Method and Application of Low Sampling Rate Time Pulse Interpolation

CHENG Yubo

（Chongqing University-University of Cincinnati Joint Co-Op Institute， Chongqing， 400044China）

Abstract： The sampling time accuracy of seismic data collector is an important index to evaluate the performance of equipment. According to the industry standard， the sampling data of time pulse signal output by standard time source should be used for analysis and calculation. However， due to the low sampling rate of seismic data collector， which caused the poor resolution of time synchronization error test. In order to improve the ability of low sampling rate time pulse analysis， this paper studies and implements the high and low potential identification algorithm of pulse signal based on histogram analysis， the rising edge identification and extraction algorithm of pulse signal based on state information and the rising sampling rate interpolation algorithm based on zero filling in frequency domain. Based on the above algorithm， the waveform rising sampling recovery of time pulse signal collected by seismic data collector is carried out. While improving the resolution of time accuracy test， the test data processing method and process of the project are standardized.

Key Words： Seismic data collector;Time error;Rising sample rate;Test resolution

地震数据采集器是将模拟电压信号转换成数字量并记录的装置，通常与地震传感器配合进行地震监测[1-3]。在台站，地震数据采集器一般通过外接GNSS天线授时，实现采样点与UTC时间的同步，并对每个采样点标记时间戳[4-5]。对于地震监测而言，地震数据采集器采集信号的时间精度是一个非常重要的技术指标，按照行业标准应使用标准时间源输出的时间脉冲信号作为标准信号，由地震数据采集器对时间脉冲信号进行采样，计算地震数据采集器采样信号的脉冲上升沿时刻与整秒的时间间隔作为地震数据采集器的时间同步误差，该指标的测试能力需要达到μs级。

地震数据采集器一般为连续记录模式，为了避免数据冗余，通常利用数据抽取技术降低其采样率，实际数据采样率一般为100sps（sample per second，每秒采样次数）和200sps。对于如此低的采样率，在时间精度测试时，时钟脉冲上升沿经过地震数据采集器的数据抽取和低通滤波后已经不是标准的阶跃信号，也不是线性上升的直线，如果直接读取上升沿脉冲幅度的50%作为采样时间，分辨力不能满足测试需求。为解决上述问题，美国SANDIA实验室使用升采样插值算法提高脉冲信号分辨力的方法开展了大量理论和实验[6]，国内谢剑波[7-8]等人也對该方法进行了进一步的研究和应用。

图1时间脉冲采样信号

本文针对地震数据采集器时间精度测试分辨力不足的问题，使用低采样率测试数据，采用升采样率插值算法进行波形恢复，从而实现提高时间精度测试分辨力的目的。为实现上述目标，基于matlab分别实现了高低电位识别算法、脉冲信号上升沿识别提取算法、升采样率插值算法，实现了测试数据的自动处理，规范了测试流程。

1测试流程设计

在地震数据采集器时间误差测试中，一般采集标准时间源输出的分脉冲信号[9]，连续记录至少20min的数据。在数据分析处理中，需要先识别出采集数据中的上升沿数据段，再对每个数据段的数据进行升采样处理，最后计算该段数据与标准时间的时间差。

针对脉冲幅度50%的这个参考点，时间标签到记录时钟脉冲上升沿时间差的计算流程如下：

（1）分析数采记录的参考时钟脉冲串信号，得到记录时间段内的所有时钟脉冲上升沿位置和波形高低电位位置;

（2）定位一个时钟脉冲上升沿，取上升沿前后至少32 个样点的数据，用插值的方法提升采样率16倍，假设采样率提升16倍后波形在相邻样点间是线性变化的，于是可以插值计算波形上升到脉冲幅度参考点电位处的时间;

（3）用当前上升沿对应的时钟时间减参考点时间得到该时间标签的时间差;

（4）重复步骤2、3 直至记录时间段内所有时钟脉冲都处理完;

（5）计算时间差平均值;

（6）计算记录时间段内对时间差平均值的标准偏差。

2  插值方法的选择

插值可以有线性插值和非线性插值等多种方法[10]。对于地震数据采集器采集到的阶跃响应分析来讲，受采样率限制，地震数据采集器等价为一个理想的低通滤波器，采集的波形会及产生吉布斯效应，即输出波形会存在吉布斯波纹，为了更好地恢复输出波形，本文使用惠特克插值的方法实现采样信号的升采样率插值处理。惠特克插值可以在时域或频域执行，是一种从信号的采样值重建带限信号的方法。这种插值方法的优点是不增加任何额外的频率内容，就能重建带限信号。

2.1 时域惠特克插值

在时域，惠特克插值在概念上等价于在现有数据点之间插入零值。在频域中，插入零值的作用是生成现有频率内容的频谱副本。此时，使用sinc滤波器可用来对信号进行低通滤波并去除插入零值产生的频谱副本。对于长度为N的序列y[k]，可以通过下列公式进行因子D的插值。

y ̃[n]=∑_（k=0）^（N-1）▒〖y[k]（sin⁡[π（n/D-k）]）/（π（n/D-k））〗

2.2 频域惠特克插值

在频域，先将数据序列填充为0，使其长度为2的幂，然后对数据进行傅里叶变换，在傅里叶变换后的数据中间插入0，以模拟采样率的增加。这种处理方法在原理上相当于在样本之间插入零，并在时域中执行一个理想的sinc滤波器，具体算法将在后文中做详细介绍。

3数据处理算法的实现

3.1  基于柱状图分析的脉冲信号高低电位识别算法及实现

对脉冲串信号的分析首先按设定的柱状图区间数M（取256）遍历记录数据以计算脉冲波形高低电位均值。用变量i遍历柱状图区间，变量j遍历波形数据，即Ai表示柱状图区间i边界，Bi表示波形数据落在第i个区间内的数量，Ci表示波形数据落在第i个区间的幅度累加，Sj表示第j个波形数据，以下是生成柱状图步骤：

a）搜索波形数据的最大、最小幅度值Amax、Amin;

b）计算波形幅度范围：AR=Amax – Amin;

c）计算柱状图区间宽度：△A=AR/M=（Amax-Amin）/M;

d）对所有i，置Bi=0，Ci=0;

e）幅值A0=Amin;

f）对i循环;

g）Ai+1=Ai+△A;

h）对j循环;

i）if（Ai+1≥Sj>Ai） thenBi=Bi+1，Ci=Ci+Sj，endif;

j）j循环结束;

k）i寻款结束;

根据以上步骤生成的柱状图（Bi～Ai关系），可在下半段区间查找最大的Bi，则波形低电位均值：

State_low =Ci/Bi（formaximumBiinthelowerpartofhistogram）

同理，在上半段区间查找最大的Bi，则波形高电位均值：

State_high =Ci/Bi（formaximumBiinthelowerpartofhistogram）

则时钟脉冲波形幅度A=State_high – State_low。

实验中，我们使用铷原子频率标准设备输出占空比为1：1的時钟分脉冲，处理从整点开始1200s的实验数据。利用柱状图分析方法，可识别出脉冲信号的低电位均值为-8.2950×10-4，高电位均值为1.1293×108。

图2 时钟脉冲信号原始波形与直方图

3.2  基于状态信息的脉冲信号上升沿识别提取算法及实现

下面讨论提取脉冲上升沿的方法。该算法的核心思想是将时钟脉冲信号分为高电平、低电平和转换态3个状态，其中考虑到吉布斯效应的存在，高电平和低电平的定义区间为高电位均值和低电位均值的±3%，其余为转换态。遍历全部数据，将数据按照高电平、低电平和转换态进行分段。实现该算法的matlab语句如下：

%定义状态和变量，State_low为状态1，State_high为状态2

y = S; N_samples = length（y）;                          %S为脉冲波形

N_states=2;                                            %波形状态数

assigned_state = [];

state_upper=[State_low+Am*0.03，State_high+Am*0.03];

state_lower=[State_low-Am*0.03，State_high-Am*0.03];

dmin = strate*10;                                      %滤波参数

%建立波形样点的状态序列

for i=1：N_samples

assigned_state（i）=0;

for j=1：N_states

if（（state_lower（j）<=y（i））&&（y（i）<=state_upper（j）））

assigned_state（i） = j;

end

end

end

%将波形序列分割成子片段

i = 1;

k = 1;

while（i<=N_samples）

current_assignment = assigned_state（i）;

sub_start（k） = i;

sub_type（k） = current_assignment;

while（（assigned_state（i）==current_assignment）&&（i<N_samples））

i = i + 1;

end

sub_end（k） = i-1;

if（（sub_end（k）-sub_start（k））<（dmin-1））

sub_type（k）=0;

end

k=k+1;

i=i+1;

end

N_sub = k-1;

处理后，共将原始数据分为N_sub段，每段的起始位置和终止位置分别存在数组sub_start和sub_end中，该段的状态信息存在数组sub_type中，其中2表示高电平、1表示低电平、0表示转换态。波形数据分析完成后，识别时钟脉冲上升沿的工作变成了搜索sub_class[]数组，查找到分类为transition的子片段后，再判断其是否前一个子片段分类为状态1而后一个子片段分类为状态2，满足该条件的子片段就是上升沿。

（a）状态分段 （b）上升沿提取

图3 状态分段与上升沿识别提取结果

一般来讲，由于时钟信号信噪比较高，经上述处理后即可将全部脉冲分段。但当信号信噪比较低时，经上述处理生成的子片段可能还会包含一些连续的不确定段，它们可能是转换态也可能是暂态，需要进一步分析合并这些不确定段。该算法相对容易，在此不过多描述。

3.3  基于频率域补零的升采样率插值算法及实现

先讨论采样数据参考点处时间的计算方法。假设选定的上升沿参考点为脉冲幅度的x%处（一般选择50%），根据柱状图分析结果，脉冲幅度为：

A=State_high-State_low

则参考点电位为：

y_（x%）=State_low+A x/100

于是可根据最接近参考点的2个相邻采样点的采样数据和采样时间，按下式计算参考点处的时间：

t_（x%）=t_（x%-）+（（t_（x%+）-t_（x%-））/（y_（x%+）-y_（x%-） ））（y_（x%）-y_（x%-） ）

式中：y_（x%-）和y_（x%+）是最接近上升沿脉冲幅度的x%参考点的2个相邻样点的电位，且满足y_（x%-）≤〖y_（x%）≤y〗_（x%+）;t_（x%-）和t_（x%+）则是对应y_（x%-）和y_（x%+）两点的时间。

最后，对应所选择上升沿的数采时标t_tag的时间同步误差为：

t_（sync_err）=t_tag-t_（x%）

由于数采的采样数据是经过抽取滤波后的离散数据，虽然按采样定理可以恢复Nyquist频率以内的信号，但信号在相邻采样点间的实际变化不是线性的。因此，若直接对原始采样数据按上述公式计算时间，则会存在较大的误差。为提高测试精度，需采用升采样插值算法先将数据采样率提高。

升采样率插值处理可以在时域或频域进行。由于我们的数据是已有的采样数据，更适合用频域方法。从实用及效率角度出发，我们假设采样率提高16倍后信号在新的相邻样点间近似是线性变化的，因此做16倍升采样处理。

定位上升沿后，截取上升沿前后至少32个样点数据（2的整数幂长度），使之包含完整的吉布斯震荡数据直至状态平稳。将截取得到的数据作状态反相复制并与原序列拼接，得到对称的新序列：

4  实验结果

应用本文的算法，对比插值前后的时间脉冲信号波形，可以明显看出插值后的数据可以较好地恢复出吉布斯震荡，时间误差测量分辨力得到明显提升。

图5 升采样前后波形比较（红色为原采样点，绿色为升采样后的波形）

分别计算某地震数据采集器采集GNSS同步和非同步状态下的一個小时标准时间脉冲测量数据。可以看到，在GNSS同步状态下，地震数据采集器采集信号的时间偏差较为稳定，最大时间偏差仅为0.3μs;在非GNSS同步状态下，地震数据采集器一般使用内部晶振提供时钟信息，与标准时间相比时间偏差呈线性漂移，与设备原理一致。

（a）GNSS同步状态 （b）GNSS非同步状态

图6 数采时间偏差测试结果

5结语

本文基于matlab分别实现了高低电位识别算法、脉冲信号上升沿识别提取算法和升采样率插值算法，基于上述算法对地震数据采集器采集的时间脉冲信号进行波形升采样恢复，从而提高时间误差测试分辨力，同时规范了该项目测试数据处理方法和流程。

参考文献

[1]蔡经纬.基于正弦波的数据采集器相位特性测试技术研究[D].北京：中国地震局地震预测研究所，2020.

[2]李彩华，滕云田，李小军，等.地震数据采集器非线性指标测试方法及其影响分析[J].地震学报，2019，41（1）：92-99.

[3]王增波，黄少卿，尚民强，等.深海拖缆地震数据采集实时质量控制[J].石油地球物理勘探，2020，55（S1）：9-14+4.

[4]刘益成，姜耕，易碧金.地震数据采集站系统延时测试方法[J].石油天然气学报，2010，32（3）：55-58+407.

[5]刘益成，易碧金，巩庆钢.地震数据采集系统计时精度测试方法研究[J].地球物理学进展，2009，24（2）：759-762.

[6]刘益成，易碧金，巩庆钢，等.地震数据采集系统相位特性与测试方法[J].石油物探，2009，48（2）：168-170+174+17.

[7]谢剑波，袁松湧，叶世山，等.用时钟方波测试地震数采时间精度及系统特性[J].大地测量与地球动力学，2018，38（10）：1073-1079.

[8]许卫卫.地震数据采集器自噪声检测研究[J].地震学报，2017，39（5）：806-813+819.

[9]DB/T 22-2020， 地震观测仪器进网技术要求  地震仪[S].

[10]孟慧宁. 两类非线性插值型曲线细分方法[D].杭州：杭州电子科技大学，2018.

1008501186292