喬玉，陳凱，陽琴

(1.中國地質(zhì)大學(xué)(北京) 地球物理與信息技術(shù)學(xué)院,北京 100083；2.中國電子科技集團(tuán)公司 第五十八研究所，江蘇 無錫 214000)

0 引言

海底大地電磁方法是通過測量海底大地電磁信號，經(jīng)數(shù)據(jù)處理及反演得到海底以下地層的電性信息，進(jìn)而獲取海底深部電性結(jié)構(gòu)。20世紀(jì)末，Steven等[1]已設(shè)計出可在1000 m深水區(qū)獲得良好的電磁響應(yīng)的海底采集系統(tǒng)。經(jīng)過近些年的發(fā)展，在設(shè)備研制、資料處理方法上取得了長足進(jìn)展，但國外在設(shè)備研制方面一直對我國進(jìn)行技術(shù)封鎖。為了打破國外在OBEM設(shè)備方面的壟斷，國內(nèi)同行從1998年開始進(jìn)行海洋大地電磁相關(guān)的方法理論、儀器設(shè)備、資料處理、反演解釋等研究工作，近年來獲得了顯著的成功[2-7]。中國地質(zhì)大學(xué)(北京)在多期“863計劃”課題及國家專項的支持下，建立了基于ARM 的海底大地電磁數(shù)據(jù)采集系統(tǒng)——OBEM-Ⅲ，初步實現(xiàn)了低功耗、小型化以及智能化[8]，在油氣勘探、水合物調(diào)查、深部構(gòu)造研究等領(lǐng)域取得了一定的成果[9-10]，儀器的主要技術(shù)指標(biāo)已經(jīng)與以美國SIO為代表的國際先進(jìn)同行齊平。在國內(nèi)外同行的推進(jìn)下，海底大地電磁方法技術(shù)在海域油氣勘探[11]、地下水調(diào)查[12]等資源環(huán)境領(lǐng)域和洋脊擴(kuò)張、板塊俯沖[13]、海底火山[14]等構(gòu)造地質(zhì)研究領(lǐng)域取得了顯著的成果。

在野外測量之前[2]，為確認(rèn)儀器性能并校正測量數(shù)據(jù)中包含的儀器通道響應(yīng)，通常需要對儀器開展通道標(biāo)定。通道標(biāo)定的過程包括產(chǎn)生標(biāo)準(zhǔn)方波并采集存儲和標(biāo)定計算兩個過程，涉及硬件電路控制以及軟件數(shù)據(jù)處理：① 在硬件電路方面，通過相關(guān)命令控制集成在采集艙內(nèi)部的標(biāo)準(zhǔn)信號發(fā)生模塊產(chǎn)生3組不同頻率的方波，并送至采集通道的最前端進(jìn)行采集存儲，生成通道標(biāo)定文件(簡稱標(biāo)定文件)；② 在數(shù)據(jù)處理程序中，讀取標(biāo)定文件中的數(shù)據(jù)，將該數(shù)據(jù)與原標(biāo)準(zhǔn)方波(標(biāo)定信號輸入)進(jìn)行比較計算，得到通道的標(biāo)定結(jié)果，比較計算的過程，稱為標(biāo)定計算。在無需外接任何信號發(fā)生器的情況下，由標(biāo)定結(jié)果可判斷儀器工作狀態(tài)，并為后期資料處理提供通道響應(yīng)改正數(shù)據(jù)。

OBEM-Ⅲ是基于ARM-Linux平臺開發(fā)的數(shù)據(jù)采集系統(tǒng)，但現(xiàn)有的標(biāo)定計算需借助Windows平臺開展，每次標(biāo)定計算前需將標(biāo)定文件下載到本地計算機(jī)，再通過Matlab程序進(jìn)行計算，存在用戶操作繁瑣、作業(yè)效率低等缺點(diǎn)。為提高海上作業(yè)效率，亟需開發(fā)一個標(biāo)定計算程序，并將其封裝為可在ARM-Linux平臺上運(yùn)行的APP，可在海底電磁接收機(jī)本地端自主完成標(biāo)定計算。

1 硬件電路原理

標(biāo)定過程的硬件電路原理框如圖1所示，主要由選擇開關(guān)、采集電路、標(biāo)定信號產(chǎn)生電路、標(biāo)定計算電路組成。

選擇開關(guān)為單刀雙擲繼電器，可通過微控制單元(MCU，microcontroller unit)控制采集電路接入信號，通道標(biāo)定時切換到標(biāo)準(zhǔn)信號產(chǎn)生電路，實際信號測量時切換到外部傳感器信號輸入。

標(biāo)定信號產(chǎn)生電路由復(fù)雜可編程邏輯器件(CPLD, complex programmable logic device)、晶振、斬波電路以及衰減電路組成。標(biāo)準(zhǔn)信號的產(chǎn)生過程如下：CPLD 對晶振產(chǎn)生的信號進(jìn)行分頻，輸出峰峰值為3.3 V的單極性方波，再經(jīng)過斬波電路和衰減電路，得到峰峰值為100 mVpp 的雙極性方波，即標(biāo)定信號。斬波電路運(yùn)用比較器的原理，將CPLD輸出的單極性方波轉(zhuǎn)換為與其頻率一致的雙極性方波，通過改變CPLD對晶振信號的分頻系數(shù)，即可輸出不同頻率的標(biāo)定信號，晶振的頻率準(zhǔn)確性與斬波電路的性能，直接影響產(chǎn)生標(biāo)定信號的準(zhǔn)確性。

采集電路包括運(yùn)放(amp, amplifier)、模數(shù)轉(zhuǎn)換器(ADC，analogy-to-digital converter)、CPLD、MCU、GPS、SD 卡等，對輸入信號進(jìn)行采集、存儲。電路進(jìn)行通道之前，由MCU獲取當(dāng)前GPS時間并記錄，同時儀器本地時間以此為準(zhǔn)繼續(xù)運(yùn)行。通道標(biāo)定開始后，由標(biāo)定信號產(chǎn)生電路生成的信號通過運(yùn)放進(jìn)行濾波、放大，經(jīng)ADC進(jìn)行模數(shù)轉(zhuǎn)換、 CPLD 完成數(shù)據(jù)整合，最后由 MCU將采集到的數(shù)據(jù)寫入SD卡內(nèi)，生成對應(yīng)的標(biāo)定文件。后期用Matlab 對采集到的標(biāo)定文件進(jìn)行復(fù)原，其在時域上表現(xiàn)為方波且無明顯畸變，幅值準(zhǔn)確性達(dá)±1%；頻域上基本滿足標(biāo)準(zhǔn)方波的頻譜分布，相位精度小于±0.37152°。

標(biāo)定計算電路包括MCU 和SD 卡，MCU 讀取SD卡中的標(biāo)定文件，通過在頻域上與標(biāo)準(zhǔn)信號進(jìn)行比較計算，輸出標(biāo)定結(jié)果文件。

圖1 標(biāo)定原理框Fig.1 Block diagram of calibration

2 計算原理

標(biāo)定計算程序以二進(jìn)制原始時間序列為輸入，以十進(jìn)制標(biāo)定計算結(jié)果為輸出，程序設(shè)計包括3部分：讀取標(biāo)定文件、對標(biāo)定數(shù)據(jù)進(jìn)行頻譜分析、比較計算得到輸出結(jié)果文件。

2.1 讀取標(biāo)定文件

標(biāo)定文件包含TSH、TSM、TSL這3組數(shù)據(jù)，分別對應(yīng)不同的采樣率以及標(biāo)定信號頻率，詳細(xì)的參數(shù)見表1，其中基頻代表輸入標(biāo)定信號的頻率，保留諧波數(shù)是頻譜分析時保留的最高次諧波。在進(jìn)行數(shù)據(jù)讀取之前，需判斷輸入標(biāo)定文件的類型并設(shè)置對應(yīng)的參數(shù)。

表1 標(biāo)定文件的部分參數(shù)

每個標(biāo)定文件包括16 B的塊頭以及1MB的數(shù)據(jù)體。塊頭中包含數(shù)據(jù)塊 ID 和相對時間δt，ID 表示當(dāng)前數(shù)據(jù)塊的編號，δt代表與對鐘時間相間隔的時間。數(shù)據(jù)體中包含8CH×32768個采樣點(diǎn)的數(shù)據(jù)，每個采樣點(diǎn)數(shù)據(jù)占4 B，以小端模式存儲，組織形式為：

L+M+H+channel_flag,

(1)

式中：L,M,H分別對應(yīng)低、中、高位數(shù)據(jù)，channel_flag為通道標(biāo)志位。單通道數(shù)據(jù)的實際值可表示為式(2)，單位：V：

(2)

2.2 對標(biāo)定數(shù)據(jù)進(jìn)行頻譜分析

從標(biāo)定文件中讀取出來的數(shù)據(jù)，均被放置在有限長的數(shù)組中，等待進(jìn)行下一步的數(shù)據(jù)處理。標(biāo)定計算是主要完成對采集信號的頻譜分析。對于有限長時間序列，通常采用離散傅里葉變換(DFT , discrete fourier transform)進(jìn)行頻譜分析。

對于N點(diǎn)序列x(n)，其DFT變換的定義為：

(3)

0≤k≤N-1(4)

從式(4)可以看出，N點(diǎn)DFT計算需要N2次復(fù)數(shù)乘法運(yùn)算、N(N-1)次復(fù)數(shù)加法運(yùn)算。因為實現(xiàn)一次復(fù)數(shù)相乘，需要四次實數(shù)相乘和兩次實數(shù)相加，所以復(fù)數(shù)乘法的計算次數(shù)直接影響程序的運(yùn)行時間，特別是當(dāng)N很大時，計算量呈指數(shù)增加。由于需要處理的數(shù)據(jù)較多，且程序的目標(biāo)運(yùn)行平臺計算速度不高，直接套用式(3)，運(yùn)算效率低，耗費(fèi)時間多，不利于實時對數(shù)據(jù)進(jìn)行處理。為了滿足現(xiàn)場作業(yè)的需求，選用快速傅里葉變換(FFT, fast fourier transform)對數(shù)據(jù)進(jìn)行處理。

2.2.1 混合基FFT

FFT計算(又稱蝶形計算)是式(3)經(jīng)過變形之后的圖形表達(dá)，其結(jié)構(gòu)可清楚的表現(xiàn)數(shù)據(jù)之間的依賴關(guān)系。每完成一組蝶形運(yùn)算，即可釋放輸入值的存儲空間給輸出值，實現(xiàn)同址計算，減少存儲空間消耗。常見的運(yùn)算對有基2、基3等，表2為長度為N時FFT 計算與DFT 計算復(fù)數(shù)乘法的運(yùn)算量比較情況，由表可以看出，隨著N的不斷增加，與直接DFT計算相比，固定基算法計算量的優(yōu)勢越來越明顯，但常見的固定基FFT運(yùn)算只能用于處理長度為Qk(Q為基底數(shù)，k為整數(shù))的數(shù)據(jù)，由表1可知，3 種文件的FFT計算長度并不能滿足這一要求。在Windows系統(tǒng)上，可通過末位補(bǔ)零法，增加數(shù)據(jù)長度，使其達(dá)到固定基FFT對輸入序列的長度要求但相同的方法在ARM-Linux 平臺上會大大增加運(yùn)算時間。

混合基FFT運(yùn)算是固定基FFT運(yùn)算推廣后的一般情況，受到計算點(diǎn)數(shù)的限制更少，可以對點(diǎn)數(shù)為非質(zhì)數(shù)的DFT 進(jìn)行快速計算，是標(biāo)定計算的最佳選擇。本程序中按照時間抽取的方法，先將X(k)逐層展開到x(n)，將x(n)進(jìn)行相應(yīng)的排序后，再從x(n)逐層計算回X(k)。

表2 DFT與FFT復(fù)數(shù)乘法運(yùn)算量比較

2.2.2 多基多進(jìn)制排序

由于對x(n)需進(jìn)行逐次抽取，必須對輸入序列倒位序，得到的次序相當(dāng)于是原始序列次序的碼位倒置，這一過程稱為多基多進(jìn)制的排序。以二進(jìn)制數(shù)為例，任意整數(shù)N都可表示為：

(5)

式中:nr,nr-1,…,n2,n1均小于2，等號左邊為二進(jìn)制表達(dá)方式，括號外的下角標(biāo)2代表每位都以2為基底；等號右邊為其對應(yīng)的十進(jìn)制數(shù)表示，經(jīng)過倒序后，其新的位序為(n1n2…nr-1nr)=n1×2r-1+n2×2r-2+…+nr-1×21+nr×20，新的位序也代表了其新的輸入順序。在多基多進(jìn)制的情況下，F(xiàn)FT計算長度N可以表示為：

(6)

式中nr,nr-1,…,n2,n1分別小于xr,xr-1,…,x2,x1，等號左邊為多進(jìn)制表達(dá)方式，括號外的下角標(biāo)xr,xr-1,…,x2,x1分別代表每一位的基底；等號右邊為其對應(yīng)的十進(jìn)制數(shù)表示。同理，其倒敘后的新位序N’表示為:

(7)

倒序之后的序列按照對應(yīng)的十進(jìn)制大小從小到大進(jìn)行排序。經(jīng)過排序后的新序列逐層抽取計算，便可得到最終的FFT計算結(jié)果。由于旋轉(zhuǎn)因子WNnk具有周期性，計算量可進(jìn)一步縮小。圖2為N=12=3×2×2時的逐層抽取的過程，其中x’(n)為經(jīng)過倒序后新序列，X11(k)為x’(n)經(jīng)過第一次基3抽取的結(jié)果，X1(k)序列為X11(k)經(jīng)過基二抽取的結(jié)果，X(k)為X1(k)經(jīng)過基二抽取得到的最終結(jié)果。圖3為整個混合基FFT的流程框，其中index2, index3, index5分別為2,3,5的冪，即混合基計算的長度LEN=2index2×3index3×5index5，基2FFT、基3FFT與基5FFT計算采用式(3)進(jìn)行。

2.3 比較計算

對標(biāo)定文件完成FFT計算之后，需與輸入的標(biāo)定信號進(jìn)行比較分析得到通道輸入響應(yīng)。在程序中定義幅值分別為1250以及1075000的理想方波來代替輸入的標(biāo)定信號，其長度與FFT計算長度相同，也進(jìn)行同樣的頻譜分析。進(jìn)行完FFT計算后的理想方波與標(biāo)定文件在相同頻點(diǎn)處進(jìn)行幅值相除、相位相減，再存儲到SD卡中，便生成最后的結(jié)果文件。

圖2 12點(diǎn)混合基FFT流程Fig.2 12-point flow chart for Mixed-radix

圖3 混合基FFT程序框Fig.3 Block diagram for Mixed-radix

3 室內(nèi)測試

編寫好的程序通過交叉編譯后上傳至SD卡，輸入相關(guān)命令，即可在OBEM-Ⅲ進(jìn)行本地標(biāo)定計算。圖4為本地標(biāo)定計算的實際過程，從計算開始到結(jié)束只需11 s，計算結(jié)果存儲到SD卡的新建文件中,即圖4中的PM21L.asc文件。

圖5為標(biāo)定計算結(jié)果的幅值相位圖，由圖可以看出，電道與磁道在通道增益以及相位方面的一致性較好，標(biāo)定計算的結(jié)果有效可靠。

圖4 標(biāo)定計算過程Fig.4 Progress of calibration

4 海上實測結(jié)果

2020年7～8月，在南海西南次海盆海底大地電磁探測科學(xué)任務(wù)中，投入了23臺海底電磁接收機(jī)，工作水深2 508～4 443 m，完成了46站位海底大地電磁采集任務(wù)。海上任務(wù)開展前，均完成了通道標(biāo)定工作，圖6給出了現(xiàn)場經(jīng)預(yù)處理得到的測深曲線(S33站位)，有效觀測頻段為30～30 000 s，獲得了高質(zhì)量海底大地電磁測深數(shù)據(jù)。

5 結(jié)論

針對海底電磁接收機(jī)標(biāo)定計算程序依賴于 PC 端的 Matlab 平臺的不足，開展了基于 ARM-Linux 平臺的計算程序開發(fā)，設(shè)計最優(yōu)計算參數(shù)，獲得了通道頻率響應(yīng)，得到了標(biāo)定文件。經(jīng)室內(nèi)測試驗證了實測結(jié)果的正確性，計算耗時約為11 s。該計算程序應(yīng)用于海上數(shù)據(jù)采集任務(wù)中，獲得了高質(zhì)量海底大地電磁測深數(shù)據(jù)，有效提升了海上現(xiàn)場作業(yè)效率。