低采樣率時間脈沖插值分析方法及應用

程煜博 趙立軍

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

關鍵詞：地震數據采集器 時間誤差 升采樣率 測試分辨力

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]。對于地震數據采集器采集到的階躍響應分析來講，受采樣率限制，地震數據采集器等價為一個理想的低通濾波器，采集的波形會及產生吉布斯效應，即輸出波形會存在吉布斯波紋，為了更好地恢復輸出波形，本文使用惠特克插值的方法實現采樣信號的升采樣率插值處理?；萏乜瞬逯悼梢栽跁r域或頻域執行，是一種從信號的采樣值重建帶限信號的方法。這種插值方法的優點是不增加任何額外的頻率內容，就能重建帶限信號。

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同步和非同步狀態下的一個小時標準時間脈沖測量數據?？梢钥吹剑贕NSS同步狀態下，地震數據采集器采集信號的時間偏差較為穩定，最大時間偏差僅為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