低采樣率時間脈沖插值分析方法及應(yīng)用

程煜博 趙立軍

摘 要：地震數(shù)據(jù)采集器采樣時間精度是評價設(shè)備性能的重要指標，按照行業(yè)標準應(yīng)使用標準時間源輸出的時間脈沖信號采樣數(shù)據(jù)進行分析計算，但由于地震數(shù)據(jù)采集器的采樣率較低，導(dǎo)致時間同步誤差測試分辨力不高。為提高低采樣率時間脈沖分析能力，本分研究并實現(xiàn)了基于柱狀圖分析的脈沖信號高低電位識別算法、基于狀態(tài)信息的脈沖信號上升沿識別提取算法和基于頻率域補零的升采樣率插值算法，基于上述算法對地震數(shù)據(jù)采集器采集的時間脈沖信號進行波形升采樣恢復(fù)，提高時間精度測試分辨力的同時，規(guī)范了該項目測試數(shù)據(jù)處理方法和流程。

關(guān)鍵詞：地震數(shù)據(jù)采集器 時間誤差 升采樣率 測試分辨力

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

地震數(shù)據(jù)采集器是將模擬電壓信號轉(zhuǎn)換成數(shù)字量并記錄的裝置，通常與地震傳感器配合進行地震監(jiān)測[1-3]。在臺站，地震數(shù)據(jù)采集器一般通過外接GNSS天線授時，實現(xiàn)采樣點與UTC時間的同步，并對每個采樣點標記時間戳[4-5]。對于地震監(jiān)測而言，地震數(shù)據(jù)采集器采集信號的時間精度是一個非常重要的技術(shù)指標，按照行業(yè)標準應(yīng)使用標準時間源輸出的時間脈沖信號作為標準信號，由地震數(shù)據(jù)采集器對時間脈沖信號進行采樣，計算地震數(shù)據(jù)采集器采樣信號的脈沖上升沿時刻與整秒的時間間隔作為地震數(shù)據(jù)采集器的時間同步誤差，該指標的測試能力需要達到μs級。

地震數(shù)據(jù)采集器一般為連續(xù)記錄模式，為了避免數(shù)據(jù)冗余，通常利用數(shù)據(jù)抽取技術(shù)降低其采樣率，實際數(shù)據(jù)采樣率一般為100sps（sample per second，每秒采樣次數(shù)）和200sps。對于如此低的采樣率，在時間精度測試時，時鐘脈沖上升沿經(jīng)過地震數(shù)據(jù)采集器的數(shù)據(jù)抽取和低通濾波后已經(jīng)不是標準的階躍信號，也不是線性上升的直線，如果直接讀取上升沿脈沖幅度的50%作為采樣時間，分辨力不能滿足測試需求。為解決上述問題，美國SANDIA實驗室使用升采樣插值算法提高脈沖信號分辨力的方法開展了大量理論和實驗[6]，國內(nèi)謝劍波[7-8]等人也對該方法進行了進一步的研究和應(yīng)用。

圖1時間脈沖采樣信號

本文針對地震數(shù)據(jù)采集器時間精度測試分辨力不足的問題，使用低采樣率測試數(shù)據(jù)，采用升采樣率插值算法進行波形恢復(fù)，從而實現(xiàn)提高時間精度測試分辨力的目的。為實現(xiàn)上述目標，基于matlab分別實現(xiàn)了高低電位識別算法、脈沖信號上升沿識別提取算法、升采樣率插值算法，實現(xiàn)了測試數(shù)據(jù)的自動處理，規(guī)范了測試流程。

1測試流程設(shè)計

在地震數(shù)據(jù)采集器時間誤差測試中，一般采集標準時間源輸出的分脈沖信號[9]，連續(xù)記錄至少20min的數(shù)據(jù)。在數(shù)據(jù)分析處理中，需要先識別出采集數(shù)據(jù)中的上升沿數(shù)據(jù)段，再對每個數(shù)據(jù)段的數(shù)據(jù)進行升采樣處理，最后計算該段數(shù)據(jù)與標準時間的時間差。

針對脈沖幅度50%的這個參考點，時間標簽到記錄時鐘脈沖上升沿時間差的計算流程如下：

（1）分析數(shù)采記錄的參考時鐘脈沖串信號，得到記錄時間段內(nèi)的所有時鐘脈沖上升沿位置和波形高低電位位置;

（2）定位一個時鐘脈沖上升沿，取上升沿前后至少32 個樣點的數(shù)據(jù)，用插值的方法提升采樣率16倍，假設(shè)采樣率提升16倍后波形在相鄰樣點間是線性變化的，于是可以插值計算波形上升到脈沖幅度參考點電位處的時間;

（3）用當(dāng)前上升沿對應(yīng)的時鐘時間減參考點時間得到該時間標簽的時間差;

（4）重復(fù)步驟2、3 直至記錄時間段內(nèi)所有時鐘脈沖都處理完;

（5）計算時間差平均值;

（6）計算記錄時間段內(nèi)對時間差平均值的標準偏差。

2 ?插值方法的選擇

插值可以有線性插值和非線性插值等多種方法[10]。對于地震數(shù)據(jù)采集器采集到的階躍響應(yīng)分析來講，受采樣率限制，地震數(shù)據(jù)采集器等價為一個理想的低通濾波器，采集的波形會及產(chǎn)生吉布斯效應(yīng)，即輸出波形會存在吉布斯波紋，為了更好地恢復(fù)輸出波形，本文使用惠特克插值的方法實現(xiàn)采樣信號的升采樣率插值處理?；萏乜瞬逯悼梢栽跁r域或頻域執(zhí)行，是一種從信號的采樣值重建帶限信號的方法。這種插值方法的優(yōu)點是不增加任何額外的頻率內(nèi)容，就能重建帶限信號。

2.1 時域惠特克插值

在時域，惠特克插值在概念上等價于在現(xiàn)有數(shù)據(jù)點之間插入零值。在頻域中，插入零值的作用是生成現(xiàn)有頻率內(nèi)容的頻譜副本。此時，使用sinc濾波器可用來對信號進行低通濾波并去除插入零值產(chǎn)生的頻譜副本。對于長度為N的序列y[k]，可以通過下列公式進行因子D的插值。

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

2.2 頻域惠特克插值

在頻域，先將數(shù)據(jù)序列填充為0，使其長度為2的冪，然后對數(shù)據(jù)進行傅里葉變換，在傅里葉變換后的數(shù)據(jù)中間插入0，以模擬采樣率的增加。這種處理方法在原理上相當(dāng)于在樣本之間插入零，并在時域中執(zhí)行一個理想的sinc濾波器，具體算法將在后文中做詳細介紹。

3數(shù)據(jù)處理算法的實現(xiàn)

3.1 ?基于柱狀圖分析的脈沖信號高低電位識別算法及實現(xiàn)

對脈沖串信號的分析首先按設(shè)定的柱狀圖區(qū)間數(shù)M（取256）遍歷記錄數(shù)據(jù)以計算脈沖波形高低電位均值。用變量i遍歷柱狀圖區(qū)間，變量j遍歷波形數(shù)據(jù)，即Ai表示柱狀圖區(qū)間i邊界，Bi表示波形數(shù)據(jù)落在第i個區(qū)間內(nèi)的數(shù)量，Ci表示波形數(shù)據(jù)落在第i個區(qū)間的幅度累加，Sj表示第j個波形數(shù)據(jù)，以下是生成柱狀圖步驟：

a）搜索波形數(shù)據(jù)的最大、最小幅度值A(chǔ)max、Amin;

b）計算波形幅度范圍：AR=Amax – Amin;

c）計算柱狀圖區(qū)間寬度：△A=AR/M=（Amax-Amin）/M;

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

e）幅值A(chǔ)0=Amin;

f）對i循環(huán);

g）Ai+1=Ai+△A;

h）對j循環(huán);

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

j）j循環(huán)結(jié)束;

k）i尋款結(jié)束;

根據(jù)以上步驟生成的柱狀圖（Bi～Ai關(guān)系），可在下半段區(qū)間查找最大的Bi，則波形低電位均值：

State_low =Ci/Bi（formaximumBiinthelowerpartofhistogram）

同理，在上半段區(qū)間查找最大的Bi，則波形高電位均值：

State_high =Ci/Bi（formaximumBiinthelowerpartofhistogram）

則時鐘脈沖波形幅度A=State_high – State_low。

實驗中，我們使用銣原子頻率標準設(shè)備輸出占空比為1：1的時鐘分脈沖，處理從整點開始1200s的實驗數(shù)據(jù)。利用柱狀圖分析方法，可識別出脈沖信號的低電位均值為-8.2950×10-4，高電位均值為1.1293×108。

圖2 時鐘脈沖信號原始波形與直方圖

3.2 ?基于狀態(tài)信息的脈沖信號上升沿識別提取算法及實現(xiàn)

下面討論提取脈沖上升沿的方法。該算法的核心思想是將時鐘脈沖信號分為高電平、低電平和轉(zhuǎn)換態(tài)3個狀態(tài)，其中考慮到吉布斯效應(yīng)的存在，高電平和低電平的定義區(qū)間為高電位均值和低電位均值的±3%，其余為轉(zhuǎn)換態(tài)。遍歷全部數(shù)據(jù)，將數(shù)據(jù)按照高電平、低電平和轉(zhuǎn)換態(tài)進行分段。實現(xiàn)該算法的matlab語句如下：

%定義狀態(tài)和變量，State_low為狀態(tài)1，State_high為狀態(tài)2

y = S; N_samples = length（y）; ? ? ? ? ? ? ? ? ? ? ? ? ?%S為脈沖波形

N_states=2; ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?%波形狀態(tài)數(shù)

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; ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?%濾波參數(shù)

%建立波形樣點的狀態(tài)序列

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;

處理后，共將原始數(shù)據(jù)分為N_sub段，每段的起始位置和終止位置分別存在數(shù)組sub_start和sub_end中，該段的狀態(tài)信息存在數(shù)組sub_type中，其中2表示高電平、1表示低電平、0表示轉(zhuǎn)換態(tài)。波形數(shù)據(jù)分析完成后，識別時鐘脈沖上升沿的工作變成了搜索sub_class[]數(shù)組，查找到分類為transition的子片段后，再判斷其是否前一個子片段分類為狀態(tài)1而后一個子片段分類為狀態(tài)2，滿足該條件的子片段就是上升沿。

（a）狀態(tài)分段 （b）上升沿提取

圖3 狀態(tài)分段與上升沿識別提取結(jié)果

一般來講，由于時鐘信號信噪比較高，經(jīng)上述處理后即可將全部脈沖分段。但當(dāng)信號信噪比較低時，經(jīng)上述處理生成的子片段可能還會包含一些連續(xù)的不確定段，它們可能是轉(zhuǎn)換態(tài)也可能是暫態(tài)，需要進一步分析合并這些不確定段。該算法相對容易，在此不過多描述。

3.3 ?基于頻率域補零的升采樣率插值算法及實現(xiàn)

先討論采樣數(shù)據(jù)參考點處時間的計算方法。假設(shè)選定的上升沿參考點為脈沖幅度的x%處（一般選擇50%），根據(jù)柱狀圖分析結(jié)果，脈沖幅度為：

A=State_high-State_low

則參考點電位為：

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

于是可根據(jù)最接近參考點的2個相鄰采樣點的采樣數(shù)據(jù)和采樣時間，按下式計算參考點處的時間：

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īng)y_（x%-）和y_（x%+）兩點的時間。

最后，對應(yīng)所選擇上升沿的數(shù)采時標t_tag的時間同步誤差為：

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

由于數(shù)采的采樣數(shù)據(jù)是經(jīng)過抽取濾波后的離散數(shù)據(jù)，雖然按采樣定理可以恢復(fù)Nyquist頻率以內(nèi)的信號，但信號在相鄰采樣點間的實際變化不是線性的。因此，若直接對原始采樣數(shù)據(jù)按上述公式計算時間，則會存在較大的誤差。為提高測試精度，需采用升采樣插值算法先將數(shù)據(jù)采樣率提高。

升采樣率插值處理可以在時域或頻域進行。由于我們的數(shù)據(jù)是已有的采樣數(shù)據(jù)，更適合用頻域方法。從實用及效率角度出發(fā)，我們假設(shè)采樣率提高16倍后信號在新的相鄰樣點間近似是線性變化的，因此做16倍升采樣處理。

定位上升沿后，截取上升沿前后至少32個樣點數(shù)據(jù)（2的整數(shù)冪長度），使之包含完整的吉布斯震蕩數(shù)據(jù)直至狀態(tài)平穩(wěn)。將截取得到的數(shù)據(jù)作狀態(tài)反相復(fù)制并與原序列拼接，得到對稱的新序列：

4 ?實驗結(jié)果

應(yīng)用本文的算法，對比插值前后的時間脈沖信號波形，可以明顯看出插值后的數(shù)據(jù)可以較好地恢復(fù)出吉布斯震蕩，時間誤差測量分辨力得到明顯提升。

圖5 升采樣前后波形比較（紅色為原采樣點，綠色為升采樣后的波形）

分別計算某地震數(shù)據(jù)采集器采集GNSS同步和非同步狀態(tài)下的一個小時標準時間脈沖測量數(shù)據(jù)。可以看到，在GNSS同步狀態(tài)下，地震數(shù)據(jù)采集器采集信號的時間偏差較為穩(wěn)定，最大時間偏差僅為0.3μs;在非GNSS同步狀態(tài)下，地震數(shù)據(jù)采集器一般使用內(nèi)部晶振提供時鐘信息，與標準時間相比時間偏差呈線性漂移，與設(shè)備原理一致。

（a）GNSS同步狀態(tài) （b）GNSS非同步狀態(tài)

圖6 數(shù)采時間偏差測試結(jié)果

5結(jié)語

本文基于matlab分別實現(xiàn)了高低電位識別算法、脈沖信號上升沿識別提取算法和升采樣率插值算法，基于上述算法對地震數(shù)據(jù)采集器采集的時間脈沖信號進行波形升采樣恢復(fù)，從而提高時間誤差測試分辨力，同時規(guī)范了該項目測試數(shù)據(jù)處理方法和流程。

參考文獻

[1]蔡經(jīng)緯.基于正弦波的數(shù)據(jù)采集器相位特性測試技術(shù)研究[D].北京：中國地震局地震預(yù)測研究所，2020.

[2]李彩華，滕云田，李小軍，等.地震數(shù)據(jù)采集器非線性指標測試方法及其影響分析[J].地震學(xué)報，2019，41（1）：92-99.

[3]王增波，黃少卿，尚民強，等.深海拖纜地震數(shù)據(jù)采集實時質(zhì)量控制[J].石油地球物理勘探，2020，55（S1）：9-14+4.

[4]劉益成，姜耕，易碧金.地震數(shù)據(jù)采集站系統(tǒng)延時測試方法[J].石油天然氣學(xué)報，2010，32（3）：55-58+407.

[5]劉益成，易碧金，鞏慶鋼.地震數(shù)據(jù)采集系統(tǒng)計時精度測試方法研究[J].地球物理學(xué)進展，2009，24（2）：759-762.

[6]劉益成，易碧金，鞏慶鋼，等.地震數(shù)據(jù)采集系統(tǒng)相位特性與測試方法[J].石油物探，2009，48（2）：168-170+174+17.

[7]謝劍波，袁松湧，葉世山，等.用時鐘方波測試地震數(shù)采時間精度及系統(tǒng)特性[J].大地測量與地球動力學(xué)，2018，38（10）：1073-1079.

[8]許衛(wèi)衛(wèi).地震數(shù)據(jù)采集器自噪聲檢測研究[J].地震學(xué)報，2017，39（5）：806-813+819.

[9]DB/T 22-2020， 地震觀測儀器進網(wǎng)技術(shù)要求 ?地震儀[S].

[10]孟慧寧. 兩類非線性插值型曲線細分方法[D].杭州：杭州電子科技大學(xué)，2018.

1008501186292