張偉,胡蕾,張釗搏
(1.中國地質(zhì)調(diào)查局 成都地質(zhì)調(diào)查中心,四川 成都 610082; 2.電子科技大學 資源與環(huán)境學院,四川 成都 611731; 3.中科院成都信息技術(shù)股份有限公司,四川 成都 610041)
長周期大地電磁法目前是探測地球深部(>150 km)上地幔尺度電性結(jié)構(gòu)特征的主要方法,由烏克蘭Lviv Centre of Institute of Space Research研究所研制的LEMI-417儀器是最早引入到我國用于長時間、穩(wěn)定、連續(xù)觀測超低頻大地電磁場的一套成熟裝備。2007年起,中國科學院地質(zhì)與地球物理研究所[1]、成都理工大學[2]、中國地質(zhì)大學(北京)[3-4]以及中國地震局地質(zhì)研究所[5]等科研院所開始使用該儀器探測青藏高原及鄰區(qū)、龍門山構(gòu)造帶以及川滇黔相鄰區(qū)地震活動帶的深部地質(zhì)結(jié)構(gòu),用于揭示其關(guān)鍵地質(zhì)帶的深部動力學特征。與該儀器配套的數(shù)據(jù)處理軟件是由俄羅斯Ivan M. Varantsov等編寫的PRC_MTMV商業(yè)軟件,用于實現(xiàn)將大地電磁場從時間域到頻率域的轉(zhuǎn)換,通過傅里葉變換計算得到功率譜,并通過功率譜估計得到大地電磁阻抗要素。
眾所周知,計算機對于二進制格式數(shù)據(jù)的讀取效率遠遠大于文本格式數(shù)據(jù),因此在進行傅里葉變換時直接讀取二進制數(shù)據(jù)進行計算將大大提高程序的運行效率。但由于烏克蘭儀器生產(chǎn)廠家沒有公開給出其原始數(shù)據(jù)的二進制格式,因此在處理LEMI-417數(shù)據(jù)時,必須首先使用廠家提供的轉(zhuǎn)換程序?qū)x器二進制記錄數(shù)據(jù)轉(zhuǎn)換為文本格式,再使用PRC_MTMV軟件的LEMI2PRC子模塊將ASCII文本格式數(shù)據(jù)再次轉(zhuǎn)換為二進制數(shù)據(jù)后方能進行后續(xù)的傅里葉變換處理[3],這無疑明顯地增加了一步不必要的處理操作。為此,作者對照V5-2000儀器的二進制數(shù)據(jù)格式的特點[6],對LEMI-417儀器的原始二進制數(shù)據(jù)格式進行了深入研究,在不斷試錯后,成功剖析出了該儀器二進制記錄數(shù)據(jù)的數(shù)據(jù)結(jié)構(gòu)。
LEMI-417是專門針對長時間采集超低頻段天然電磁場而研制的一款觀測系統(tǒng),其主要特點是:①采用磁通門傳感器來采集磁場數(shù)據(jù),觀測范圍為±70 000 nT,可以記錄到幾十萬秒的長周期信號;②電道采集盒子通過專門的電路設計能夠保證在長時間范圍內(nèi)保持很低的漂移,提供了4組電道接口,可以同時采集4組電極的電場數(shù)據(jù);③采用高精度GPS同步授時能夠保證在長時間范圍內(nèi)具有很小的時移;④整體功耗很低,能夠保證在1塊12V45ah的鉛酸電池供電下連續(xù)工作一周以上。
該儀器以1 Hz的采樣率來連續(xù)觀測天然電磁場,采集的數(shù)據(jù)以二進制格式實時儲存到儀器內(nèi)部的CF卡中,文件名以L417_×××.B××的規(guī)則來進行命名,其中:L417為儀器的型號名稱;_×××為單點的數(shù)據(jù)存放文件序號,第一個數(shù)據(jù)文件編號為001,當GPS時間為格林威治時間的零點時,儀器系統(tǒng)自動新建一個新文件來存放數(shù)據(jù),文件序號自動在前一個文件序號的基礎(chǔ)上累加1;B××的首字母B是二進制(binary)的縮寫,××是采集儀器的出廠編號。
通常與測點數(shù)據(jù)采集相關(guān)的記錄信息(如儀器型號、采集日期、經(jīng)緯度、高程、電極距等)一般以文件頭的形式存放于二進制文件的首端。LEMI-417I型儀器共分配了32個字節(jié)的大小空間來存放信息(表1),各字節(jié)位置信息闡述如下。
表1 LEMI-417原始二進制記錄數(shù)據(jù)的頭文件格式(tag)
1~4字節(jié):4個字節(jié)長度,存放采集儀器的型號信息,對于LEMI-417I型儀器來說,第1個字節(jié)為ASCII字符‘4’(16進制為34),第2字節(jié)為ASCII字符‘1’(16進制為31),第3字節(jié)為ASCII字符‘7’(16進制為37),第4字節(jié)為ASCII字符‘I’(16進制為49)。
5字節(jié):1個字節(jié)長度,存放采集儀器的出廠編號。
6字節(jié):1個字節(jié)長度,存放開始采集數(shù)據(jù)時的當前年信息。
7字節(jié):1個字節(jié)長度,存放開始采集數(shù)據(jù)時的當前月信息。
8字節(jié):1個字節(jié)長度,存放開始采集數(shù)據(jù)時的當前日信息。
9字節(jié):1個字節(jié)長度,存放開始采集數(shù)據(jù)時的當前小時信息(24小時制)。
10字節(jié):1個字節(jié)長度,存放開始采集數(shù)據(jù)時的當前分信息。
11字節(jié):1個字節(jié)長度,存放開始采集數(shù)據(jù)時的當前秒信息。
12~15字節(jié):4個字節(jié)長度,存放當前測點GPS同步的緯度。
16字節(jié):1個字節(jié)長度,存放當前測點坐標是位于北緯還是南緯,以我國位于北半球為例,在我國測量數(shù)據(jù)時應為ASCII字符‘N’( 16進制為4E)。
17~21字節(jié):5個字節(jié)長度,存放當前測點GPS同步記錄的經(jīng)度信息。
22字節(jié):1個字節(jié)長度,存放當前測點坐標是位于東經(jīng)還是西經(jīng),以我國位于東半球為例,在我國測量數(shù)據(jù)時應為ASCII字符‘E’( 16進制為45)。
23字節(jié):1個字節(jié)長度,為空字符(對應的16進制為00)。
24~25字節(jié):2個字節(jié)長度,存放當前測點GPS同步的高程(單位:m)。
26字節(jié):1個字節(jié)長度,存放儀器的采樣頻率。
27字節(jié):1個字節(jié)長度,存放儀器內(nèi)部的VIN(外部供電電源電壓)信息。
28字節(jié):1個字節(jié)長度,存放儀器內(nèi)部的VBAT(內(nèi)部紐扣電池電壓)信息。
29字節(jié):1個字節(jié)長度,存放電道盒子上第1對電極的極距長度(單位:m)。
30字節(jié):1個字節(jié)長度,存放電道盒子上第2對電極的極距長度(單位:m)。
31字節(jié):1個字節(jié)長度,存放電道盒子上第3對電極的極距長度(單位:m)。
32字節(jié):1個字節(jié)長度,存放電道盒子上第4對電極的極距長度(單位:m)。
如圖1所示,在頭文件(Tag)后緊接著的是時間序列數(shù)據(jù),LEMI-417以30個字節(jié)的大小空間來存放一個采樣時間片采集到的電磁場數(shù)據(jù)段(segment,表2所示),由17個時間片組成一個記錄塊(record),一個數(shù)據(jù)段結(jié)束后再接下一個文件頭數(shù)據(jù),依次循環(huán)、直到采集結(jié)束。其中,單個時間片數(shù)據(jù)段中各字節(jié)位置信息闡述如下。
圖1 LEMI-417原始數(shù)據(jù)結(jié)構(gòu)分解示意Fig.1 Schematic diagram of original data structure decomposition
表2 LEMI-417原始二進制記錄數(shù)據(jù)的時間序列格式(segment)
1~3字節(jié):3個字節(jié)長度,存放磁通門傳感器X方向的測量值(單位:nT)。
4~6字節(jié):3個字節(jié)長度,存放磁通門傳感器Y方向的測量值(單位:nT)。
7~9字節(jié):3個字節(jié)長度,存放磁通門傳感器Z方向的測量值(單位:nT)。
10~13字節(jié):4個字節(jié)長度,存放第1道電極間的測量值(單位:mkV/m)。
14~17字節(jié):4個字節(jié)長度,存放第2道電極間的測量值(單位:mkV/m)。
18~21字節(jié):4個字節(jié)長度,存放第3道電極間的測量值(單位:mkV/m)。
22~25字節(jié):4個字節(jié)長度,存放第4道電極間的測量值(單位:mkV/m)。
26~27字節(jié):2個字節(jié)長度,存放磁通門傳感器溫度(單位:℃)。
28~29字節(jié):2個字節(jié)長度,存放電道采集盒子的溫度(單位:℃)。
30字節(jié):1個字節(jié)長度,存放時間片采樣間隔(單位:m)。
在C++語言中可利用struct類型變量來結(jié)構(gòu)化地讀取和存儲頭文件信息,再按表1所示各信息的字節(jié)位置和字節(jié)長度,定義其內(nèi)部的子變量數(shù)據(jù)類型,其核心代碼如下:
// 頭文件數(shù)據(jù)結(jié)構(gòu)定義
structTagData
{
unsignedcharstr_1; //4
unsignedcharstr_2; //1
unsignedcharstr_3; //7
unsignedcharstr_4; //儀器型號
unsignedcharno; //儀器編號
unsignedcharyear; //年
unsignedcharmonth; //月
unsignedcharday; //日
unsignedcharhour; //時
unsignedcharminute; //分
unsignedcharsecond; //秒
unsignedcharlatitude[4]; //緯度
charN_S; //北緯還是南緯
unsignedcharlongitude[5]; //經(jīng)度
charE_W; //東經(jīng)還是西經(jīng)
charseparator; //空字符
unsignedcharaltitude[2]; //高程
unsignedcharaverage; //采樣率
unsignedcharvin; // 外部電壓
unsignedcharvbat; // 內(nèi)部電壓
unsignedcharL1; //第一對極距長度
unsignedcharL2; //第二對極距長度
unsignedcharL3; //第三對極距長度
unsignedcharL4; //第四對極距長度
}
// … …由于篇幅有限省略部分非關(guān)鍵代碼
// 讀取儀器原始數(shù)據(jù)文件
TagDatatag;
FILE*fp=fopen(strFilePath,"rb");
chardataBuf[32] = "";
fread(dataBuf,sizeof(dataBuf),1,fp);
fclose(fp);
// 結(jié)構(gòu)體各子變量賦值
tag.str_1=dataBuf[0];
tag.str_2=dataBuf[1];
// … …由于篇幅有限省略部分代碼
tag.latitude[0]=dataBuf[11];
tag.latitude[1]=dataBuf[12];
tag.latitude[2]=dataBuf[13];
tag.latitude[3]=dataBuf[14];
// … …由于篇幅有限省略部分代碼
tag.L4=dataBuf[31];
同理,利用struct類型變量來結(jié)構(gòu)化地讀取和存儲時間片數(shù)據(jù)段信息,再按表2所示各信息的字節(jié)位置和字節(jié)長度,定義其內(nèi)部的子變量數(shù)據(jù)類型,其核心代碼如下:
// 單個時間片數(shù)據(jù)段結(jié)構(gòu)定義
structSegmentData
{
charHx[3]; //Hx
charHy[3]; //Hy
charHz[3]; //Hz
charE1[4]; //E1
charE2[4]; //E2
charE3[4]; //E3
charE4[4]; //E4
charTF[2]; //TF
charTE[2]; //TE
charNo; //
}
// … …由于篇幅有限省略部分非關(guān)鍵代碼
// 定義變量
intHx=0,Hy=0,Hz=0,E1=0,E2=0,E3=0,E4=0;shortintTf=0,Te=0;
charNo;
TagDatatag;
SegmentDatasegment;
// 讀取儀器原始數(shù)據(jù)文件
fseek(infile_origin,0L,SEEK_SET);
cur_rc= 0;
intNo= 0;
// 如果當前文件指針讀取位置不是結(jié)尾
while(cur_rc { if(i% 17 == 0) { fread(tag,sizeof(unsignedchar),32L,infile_origin); cur_rc=ftell(infile_origin); i++; continue; } fread(segment,sizeof(unsignedchar),30L,infile_origin); cur_rc=ftell(infile_origin); i++; j++; // 將二進制數(shù)值轉(zhuǎn)換為整型變量,詳見下文 convertSegmentData(segment,&Hx,&Hy,&Hz,&E1,&E2,&E3,&E4,&Tf,&Te,&No); pre_rc=cur_rc; } LEMI-417原始數(shù)據(jù)是在儀器采集時存儲的二進制數(shù)據(jù),將它賦值到計算程序語言的整型變量時必須進行二進制位操作轉(zhuǎn)換。以磁通門傳感器的x方向磁場數(shù)據(jù)(Hx)為例,由于整型變量(int)無論是在32位還是64位操作系統(tǒng)其字節(jié)長度均為32位,幸運地避開了不同操作系統(tǒng)程序代碼的特殊處理問題,位轉(zhuǎn)換的實質(zhì)就是將原始數(shù)據(jù)中各字節(jié)位置的數(shù)值完整地賦值到整型變量的各字節(jié)位置上,其核心代碼如下: inttemp1=0,temp2=0,temp3=0; temp1 =segment.Hx[2]; temp1=temp1 ? 16; temp1 =temp1 & 0x00FF0000; temp2 =segment.Hx[1]; temp2 =temp2 ? 8; temp2 =temp2 & 0x0000FF00; temp3 =segment.Hx[0]; temp3 =temp3 & 0x000000FF; Hx=temp1|temp2|temp3; if(Hx>=0) { Hx=Hx& 0x00FFFFFF; } else { Hx=Hx| 0xFF000000; } 在對原始數(shù)據(jù)格式成功剖析的基礎(chǔ)上,筆者針對該LEMI-417儀器在實際工作中存在的一些問題,補充性地開發(fā)出一套基于MFC的、可以在任何Windows系統(tǒng)運行的界面程序,輔助一線野外工作人員對個別偶然會出現(xiàn)的數(shù)據(jù)問題進行人為糾正。 實際應用中,LEMI-417型儀器在無GPS信號時會出現(xiàn)記錄紊亂的現(xiàn)象[4]。在峽谷地貌或林木茂密區(qū)開展工作時可能會出現(xiàn)個別時段GPS衛(wèi)星信號偶然解鎖的現(xiàn)象,使得儀器無法同步當前的GPS時鐘和測點坐標,在對應的文件頭信息中記錄下多組差別迥異的信息。例如,野外工作中儀器部署在深切山谷地貌區(qū)時,在部分時段會出現(xiàn)儀器GPS無法鎖定到≥3顆衛(wèi)星的情況,此時儀器記錄的測點經(jīng)緯度和高程均變?yōu)?,個別情況甚至會出現(xiàn)1顆衛(wèi)星都無法鎖定的情況,此時便會出現(xiàn)衛(wèi)星失鎖問題,無法同步到正確的GPS時鐘,使采集數(shù)據(jù)時序出現(xiàn)混亂,從而導致在后期數(shù)據(jù)導出時無法順利形成連續(xù)的時間序列數(shù)據(jù),需要進行2次返工,造成不必要的成本支出。 針對該常見問題,開發(fā)出了Lemi-417數(shù)據(jù)修改器(圖2),可以將當前測點采集的連續(xù)多天的數(shù)據(jù)加載、合并到一起,遍歷每個文件頭中的GPS時鐘、測點經(jīng)度、緯度、高程信息,判斷每段數(shù)據(jù)對應的GPS時鐘是否是以設置的采樣率進行連續(xù)觀測,會缺失的、不連續(xù)的時間序列進行三次樣條插值修復,此外判斷每段數(shù)據(jù)對應的經(jīng)度、緯度和高程信息是否出現(xiàn)較大的數(shù)據(jù)偏差,通過取眾數(shù)方法,識別出正常狀態(tài)下和故障狀態(tài)下的記錄,對故障狀態(tài)下的記錄進行自動修復。最后通過最小二乘法對正常狀態(tài)下的GPS經(jīng)度、緯度和高程進行線性回歸擬合得到高精度的測點定位信息。由此,通過計算機程序解決了以往由于GPS偶然失鎖而需要重復返工的問題。 圖2 窗口程序界面截圖Fig.2 Screenshot of the window program interface 隨著我國廣大學者對大地電磁時頻轉(zhuǎn)換算法研究的不斷深入[7-9],國產(chǎn)化數(shù)據(jù)處理軟件[10]、儀器設備[11]的研制工作也已取得長足進步。本文深入剖析了LEMI-417I型地球深部探測主流觀測系統(tǒng)的原始數(shù)據(jù)格式,較詳細地給出了對應的核心程序代碼,為廣大從事天然電磁場現(xiàn)代數(shù)字信號處理方面研究的學者提供一定的前期基礎(chǔ)支撐,對今后同類型國產(chǎn)儀器、軟件的研發(fā)工作具有一定的借鑒意義。此外,也針對LEMI-417I型儀器存在的GPS解鎖時的記錄紊亂問題,提出了一個經(jīng)濟有效的計算機輔助糾正技術(shù)方案。2.3 二進制數(shù)值到整型變量的轉(zhuǎn)換
3 應用開發(fā)實例
4 結(jié)論