G856和GSM—19T質(zhì)子磁力儀觀測數(shù)據(jù)日變改正處理新方法

王榮華+李文堯

摘要： 使用G856和GSM-19T質(zhì)子磁力儀在同一個測區(qū)進(jìn)行磁測數(shù)據(jù)采集，由于這兩種型號的質(zhì)子磁力儀觀測數(shù)據(jù)文件格式不一樣，用其配套軟件相互日變改正特別困難和繁瑣，并且用其配套軟件日變改正沒有五點(diǎn)或七點(diǎn)滑動平均。為了解決這些問題，在Microsoft Office Excel平臺編寫五點(diǎn)或七點(diǎn)滑動平均和日變改正VBA程序。該程序經(jīng)實(shí)驗(yàn)驗(yàn)證和一個測區(qū)20971個磁測數(shù)據(jù)生產(chǎn)應(yīng)用驗(yàn)證表明VBA程序的正確性。這個VBA程序適用不同型號的質(zhì)子磁力儀觀測數(shù)據(jù)之間日變改正，只要把觀測數(shù)據(jù)轉(zhuǎn)換到Microsoft Office Excel中并插入VBA程序代碼即可進(jìn)行日變改正，非常實(shí)用、快捷。

Abstract： When using G856 proton magnetometer and GSM-19T proton magnetometer to collect magnetic data in the same area， the mutual diurnal correction is particularly difficult and complicated with its supporting software， and its supporting software diurnal correction does not get five-point or seven-point moving average， cause these two types of proton magnetometer data file format is different. In order to solve these problems， the paper writes the five-point or seven-point moving average and the diurnal variation of VBA program in Microsoft Office Excel platform. The program is verified by experiment and a survey area 20971 applied magnetic data proved the correctness of VBA programs. The VBA program for proton magnetometer data can be applied in different types of correction for diurnal variation， it just need to convert the data into Microsoft Office Excel and insert the VBA program code to do diurnal correction， it's very practical and fast.

關(guān)鍵詞： G856質(zhì)子磁力儀；GSM-19T質(zhì)子磁力儀；日變改正；VBA程序

Key words： G856 proton magnetometer；GSM-19T proton magnetometer；correction for diurnal variation；the VBA program

中圖分類號：P716+.82；P413 文獻(xiàn)標(biāo)識碼：A 文章編號：1006-4311（2014）10-0001-03

0 引言

對同一測區(qū)使用G856和GSM-19T兩種型號的質(zhì)子磁力儀進(jìn)行磁測數(shù)據(jù)采集，兩種磁力儀觀測數(shù)據(jù)輸出文件格式不一樣，并且兩種磁力儀數(shù)據(jù)日變改正[1]軟件不一樣。G856質(zhì)子磁力儀數(shù)據(jù)輸出文件格式為stn[2]，日變改正軟件為MagMap2000[2]；GSM-19T質(zhì)子磁力儀數(shù)據(jù)輸出文件格式為txt[3]，日變改正軟件為GEMLinkW[3]。因此，這兩種型號的磁測觀測數(shù)據(jù)用配套軟件相互日變改正特別困難和繁瑣。此外，根據(jù)《地面高精度磁測技術(shù)規(guī)程》（DZ/T0071-93），為了提高日變改正精度[4]，先要對日變觀測原始數(shù)據(jù)作五點(diǎn)或七點(diǎn)滑動平均，壓制噪聲后，再對測點(diǎn)觀測數(shù)據(jù)進(jìn)行日變改正，但G856和GSM-19T質(zhì)子磁力儀配套日變改正軟件日變改正過程中沒有五點(diǎn)或七點(diǎn)滑動平均。針對這些問題，查閱了相關(guān)的文獻(xiàn)和軟件，管志寧[1]提出日變改正值從日變曲線上查得，在日變曲線上量得某時刻相對早基點(diǎn)時間的日變值并取反號，即為該時刻的日變改正值，這種日變改正方法速度慢、精度不高。于寶顯[5]等提出不同型號質(zhì)子磁力儀之間數(shù)據(jù)共享，經(jīng)轉(zhuǎn)換后的數(shù)據(jù)可以利用其中一種型號質(zhì)子磁力儀所配套軟件做日變改正，這種方法有些繁瑣，并且沒有做五點(diǎn)或七點(diǎn)滑動平均。中國地質(zhì)大學(xué)（武漢）劉天佑教授開發(fā)的磁法勘探軟件系統(tǒng)MAGS 2.0沒有實(shí)現(xiàn)這兩種型號質(zhì)子磁力儀觀測數(shù)據(jù)之間的日變改正處理。對此，通過Microsoft Office Excel平臺VBA[6-7]編程解決。

1 觀測數(shù)據(jù)格式統(tǒng)一

為了統(tǒng)一G856和GSM-19T質(zhì)子磁力儀觀測數(shù)據(jù)格式，用MagMap2000導(dǎo)入G856質(zhì)子磁力儀stn格式文件數(shù)據(jù)，然后輸出dat格式文件數(shù)據(jù)，最后把dat文件格式數(shù)據(jù)復(fù)制到Microsoft Office Excel文件中；GSM-19T質(zhì)子磁力儀txt文件格式數(shù)據(jù)可以直接導(dǎo)入Microsoft Office Excel文件中。

2 日變觀測數(shù)據(jù)五點(diǎn)或七點(diǎn)滑動平均

根據(jù)《地面高精度磁測技術(shù)規(guī)程》（DZ/T0071-93），地磁場短周期變化的振幅與質(zhì)子磁力儀噪聲均方根值是近于同一數(shù)量級，噪聲是隨機(jī)的，地磁場脈動變化是有規(guī)律的，先要對日變觀測原始數(shù)據(jù)先作五點(diǎn)或七點(diǎn)滑動平均，壓制噪聲水平后，再對測點(diǎn)觀測數(shù)據(jù)進(jìn)行日變改正，即可提高日變改正的精度。下面介紹日變觀測數(shù)據(jù)五點(diǎn)和七點(diǎn)滑動平均原理[7]。設(shè)在一個范圍內(nèi)某點(diǎn)經(jīng)平滑后的日變數(shù)值可表示：T（x）=a0+a1x （1）endprint

（1）式中：a0和a1為待定系數(shù)。設(shè)日變觀測的數(shù)值為T（xi），它的平滑值為T（xi），則可得出：

δ=■a■+a■x■-T（x■）■=min （2）

（2）式中δ為偏差的平方和。利用微分求極值的方法將（2）式對a■和a■求導(dǎo)數(shù)，然后令其為零得：

■=■2a■+a■x■-T（x■）=0

■=■2a■+a■x■-T（x■）x■=0 （3）

若x■以剖面上的點(diǎn)距為單位，步長為1，則式（3）中的x■=0，±1，±2…±m(xù)。代入式（3）中，可解得a■和a■分別為

a■=■，a■=■ （4）

由（1）式可知x=0時，T（0）=a■，T（0）=■ （5）

由此可見，（5）式就是平滑公式，當(dāng)m=±1，三點(diǎn)平滑公式為T（0）=■[T（-1）+T（0）+T（1）] （6）

同理可得五點(diǎn)平滑公式為

T（0）=■[T（-2）+T（-1）+T（0）+T（1）+T（2）] （7）

七點(diǎn)平滑公式為

T（0）=■[T（-3）+T（-2）+T（-1）+T（0）+T（1）+T（2）+T（3）] （8）

根據(jù)五點(diǎn)和七點(diǎn)滑動平均公式（7）和（8）編寫Microsoft Office Excel VBA程序，代碼如下，利用此程序代碼把日變原始觀測數(shù)據(jù)（見圖1）五點(diǎn)和七點(diǎn)滑動平均，見圖2和圖3。

Sub 日變原始數(shù)據(jù)圓滑（）

Dim i As Integer

Sheet1.Cells（1， 1） = "日變原始數(shù)據(jù)"

Sheet1.Cells（1， 2） = "日變觀測時間"

Sheet1.Cells（1， 3） = "日變五點(diǎn)圓滑后數(shù)據(jù)"

Sheet1.Cells（1， 4） = "日變七點(diǎn)圓滑后數(shù)據(jù)"

n = 649 ′此值為原始日變數(shù)據(jù)對應(yīng)的行數(shù)′

m = 5 ′選擇5點(diǎn)或7點(diǎn)滑動平均圓滑′

If 5 = m Then

For i = 4 To n - 2

Sheet1.Cells（i， 3）.Value = （Sheet1.Cells（i - 2， 1）.Value + Sheet1.Cells（i - 1， 1）.Value + Sheet1.Cells（i， 1）.Value + Sheet1.Cells（i + 1， 1）.Value + Sheet1.Cells（i + 2， 1）.Value）/m

Next i

End If

If 7 = m Then

For i = 5 To n - 3

Sheet1.Cells（i， 4）.Value = （Sheet1.Cells（i - 3， 1）.Value + Sheet1.Cells（i - 2， 1）.Value + Sheet1.Cells（i - 1， 1）.Value + Sheet1.Cells（i， 1）.Value + Sheet1.Cells（i + 1， 1）.Value + Sheet1.Cells（i + 2， 1）.Value + Sheet1.Cells（i + 3， 1）.Value）/m

Next i

End If

End Sub

3 日變觀測數(shù)據(jù)線性插值

根據(jù)《地面高精度磁測技術(shù)規(guī)程》（DZ/T0071-93）質(zhì)子磁力儀日變觀測時間間隔為5秒至20秒，為了求出測點(diǎn)觀測與日變觀測同一時刻的日變數(shù)據(jù)，先要求出日變觀測數(shù)據(jù)每一秒的數(shù)值，因此對日變數(shù)據(jù)進(jìn)行分段線性插值[9-10]，公式為（9）式，公式中T（n）為前一時刻的日變觀測數(shù)據(jù)，T（n+m）為后一時刻的日變觀測數(shù)據(jù)，m為日變觀測時間間隔，j為日變觀測第j秒，T（j）為分段線性插值后第j秒數(shù)值，Microsoft Office Excel VBA 程序代碼如下：

T（j）=■·j+T（n） （9）

Sub 日變分段線性插值（）

Dim i As Integer

Dim j As Integer

Dim n As Integer

Dim m As Integer

Sheet1.Cells（1， 1） = "日變原始數(shù)據(jù)"

Sheet1.Cells（1， 2） = "日變插值后數(shù)據(jù)"

n = 1173 ′此值為原始的日變觀測數(shù)據(jù)對應(yīng)的行數(shù)′

m=20′此值為日變觀測數(shù)據(jù)時間間隔′

For i = 2 To n - 1

For j = 0 To m-1

Sheet1.Cells（m* （i - 2） + j + 2， 2）.Value = （Sheet1.Cells（i + 1， 1）.Value - Sheet1.Cells（i， 1）.Value） / m * j + Sheet1.Cells（i， 1）.Value

Next j

Next i

Sheet1.Cells（m * （i - 2） + 2， 2）.Value=Sheet1.Cells（i， 1）

End Sub

4 提取測點(diǎn)與日變觀測同一時刻的日變數(shù)據(jù)

求出日變觀測每一秒數(shù)值后，要提取測點(diǎn)觀測與日變觀測同一時刻的日變數(shù)據(jù)，Microsoft Office Excel VBA 程序如下：

Sub 提取數(shù)據(jù)（）endprint

Sheet3.Cells（1， 1）.Value = "日變每一秒值"

Sheet3.Cells（1， 2）.Value = "日變時間"

Sheet3.Cells（1， 3）.Value = "觀測點(diǎn)時間"

Sheet3.Cells（1， 4）.Value = "觀測點(diǎn)觀測值"

Dim i As Integer

Dim j As Integer

For i = 2 To 108 ′此值為觀測點(diǎn)時間對應(yīng)的行數(shù)′

For j = 2 To 24142 ′此值為通過插值得到的數(shù)據(jù)對應(yīng)的行數(shù)′

If Sheet3.Cells（i， 3）.Value=Sheet3.Cells（j， 2）.Value Then

Sheet3.Cells（i， 4）.Value = Sheet3.Cells（j， 1）.Value

End If

Next j

Next i

End Sub

5 日變改正

提取測點(diǎn)與日變觀測同一時刻的日變數(shù)據(jù)后，用測點(diǎn)觀測的數(shù)據(jù)減去同一時刻的日變數(shù)值就是日變改正后的值。

6 驗(yàn)證

為了驗(yàn)證VBA程序的正確性，在野外磁測過程中做了實(shí)驗(yàn)驗(yàn)證。用一臺GSM-19T質(zhì)子磁力儀做日變觀測，兩臺GSM-19T質(zhì)子磁力儀做測點(diǎn)觀測，點(diǎn)距約為2米，點(diǎn)數(shù)為64個。測點(diǎn)觀測數(shù)據(jù)用GEMLinkW軟件和VBA程序進(jìn)行日變改正，用VBA程序日變改正時沒有做五點(diǎn)和七點(diǎn)滑動平均。日變改正后的數(shù)據(jù)見表1。從表1可以看出：GEMLinkW軟件和VBA程序日變改正后的值相等（用GEMLinkW軟件日變改正后的值軟件默認(rèn)進(jìn)行四舍五入保留了兩位小數(shù)，用VBA程序日變改正后的值進(jìn)行四舍五入保留了三位小數(shù)），證明了VBA程序的正確性。

此外，對VBA程序還進(jìn)行生產(chǎn)應(yīng)用驗(yàn)證，在云南省鳳慶縣某鐵礦區(qū)開展高精度磁測工作中，同時使用G856和GSM-19T質(zhì)子磁力儀進(jìn)行磁測數(shù)據(jù)采集，用VBA程序?qū)y區(qū)的20971個磁測數(shù)據(jù)進(jìn)行5%抽查人工日變改正驗(yàn)證，結(jié)果也表明VBA程序的正確性。

7 結(jié)論

在高精度磁測野外生產(chǎn)過程中經(jīng)常會遇到用不同型號的磁力儀進(jìn)行數(shù)據(jù)采集，數(shù)據(jù)采集完后必須進(jìn)行日變改正，然而不同種磁力儀之間往往不能用配套的軟件進(jìn)行日變改正。本文在Microsoft Office Excel平臺編寫五點(diǎn)或七點(diǎn)滑動平均程序和日變改正VBA程序，并且實(shí)驗(yàn)和生產(chǎn)應(yīng)用驗(yàn)證了此程序的正確性，適合野外現(xiàn)場計算，方便、快捷。

參考文獻(xiàn)：

[1]管志寧.地磁場與磁力勘探[M].北京：地質(zhì)出版社，2005：75.

[2]美國勞雷工業(yè)有限公司.G856質(zhì)子磁力儀用戶手冊[S].美國勞雷工業(yè)有限公司，1990：20-22.

[3]加拿大GEM公司.GSM-19T 操作手冊[S].加拿大GEM公司，2011：46-49.

[4]李才明，李軍，余舟等. 提高磁測日變改正精度的方法[J].物探化探計算技術(shù)，2004，26（3）.

[5]于寶顯，杜利明，王懷坤等.不同型號質(zhì)子磁力儀之間數(shù)據(jù)共享[J].內(nèi)蒙古煤炭經(jīng)濟(jì)，2013，5.

[6]劉炳文.Visual Basic 程序設(shè)計教程[M].北京：清華大學(xué)出版社，2003：216-223.

[7]John Walkenbach. Excel 2007 VBA Programming For Dummies[M].USA：Wiley Publishing，Inc，2007：36-52.

[8]曾華霖.重力場與重力勘探 [M].北京：地質(zhì)出版社，2005：196-197.

[9]李慶揚(yáng)，王能超，易大義.數(shù)值分析[M].北京：清華大學(xué)出版社，2008：23-24.

[10]張?zhí)炝粒c梅，許明夫等.基于分段線性插值法的高精度測溫研究[J].信息技術(shù)與信息化，2012（1）.endprint

Sheet3.Cells（1， 1）.Value = "日變每一秒值"

Sheet3.Cells（1， 2）.Value = "日變時間"

Sheet3.Cells（1， 3）.Value = "觀測點(diǎn)時間"

Sheet3.Cells（1， 4）.Value = "觀測點(diǎn)觀測值"

Dim i As Integer

Dim j As Integer

For i = 2 To 108 ′此值為觀測點(diǎn)時間對應(yīng)的行數(shù)′

For j = 2 To 24142 ′此值為通過插值得到的數(shù)據(jù)對應(yīng)的行數(shù)′

If Sheet3.Cells（i， 3）.Value=Sheet3.Cells（j， 2）.Value Then

Sheet3.Cells（i， 4）.Value = Sheet3.Cells（j， 1）.Value

End If

Next j

Next i

End Sub

5 日變改正

提取測點(diǎn)與日變觀測同一時刻的日變數(shù)據(jù)后，用測點(diǎn)觀測的數(shù)據(jù)減去同一時刻的日變數(shù)值就是日變改正后的值。

6 驗(yàn)證

為了驗(yàn)證VBA程序的正確性，在野外磁測過程中做了實(shí)驗(yàn)驗(yàn)證。用一臺GSM-19T質(zhì)子磁力儀做日變觀測，兩臺GSM-19T質(zhì)子磁力儀做測點(diǎn)觀測，點(diǎn)距約為2米，點(diǎn)數(shù)為64個。測點(diǎn)觀測數(shù)據(jù)用GEMLinkW軟件和VBA程序進(jìn)行日變改正，用VBA程序日變改正時沒有做五點(diǎn)和七點(diǎn)滑動平均。日變改正后的數(shù)據(jù)見表1。從表1可以看出：GEMLinkW軟件和VBA程序日變改正后的值相等（用GEMLinkW軟件日變改正后的值軟件默認(rèn)進(jìn)行四舍五入保留了兩位小數(shù)，用VBA程序日變改正后的值進(jìn)行四舍五入保留了三位小數(shù)），證明了VBA程序的正確性。

此外，對VBA程序還進(jìn)行生產(chǎn)應(yīng)用驗(yàn)證，在云南省鳳慶縣某鐵礦區(qū)開展高精度磁測工作中，同時使用G856和GSM-19T質(zhì)子磁力儀進(jìn)行磁測數(shù)據(jù)采集，用VBA程序?qū)y區(qū)的20971個磁測數(shù)據(jù)進(jìn)行5%抽查人工日變改正驗(yàn)證，結(jié)果也表明VBA程序的正確性。

7 結(jié)論

在高精度磁測野外生產(chǎn)過程中經(jīng)常會遇到用不同型號的磁力儀進(jìn)行數(shù)據(jù)采集，數(shù)據(jù)采集完后必須進(jìn)行日變改正，然而不同種磁力儀之間往往不能用配套的軟件進(jìn)行日變改正。本文在Microsoft Office Excel平臺編寫五點(diǎn)或七點(diǎn)滑動平均程序和日變改正VBA程序，并且實(shí)驗(yàn)和生產(chǎn)應(yīng)用驗(yàn)證了此程序的正確性，適合野外現(xiàn)場計算，方便、快捷。

參考文獻(xiàn)：

[1]管志寧.地磁場與磁力勘探[M].北京：地質(zhì)出版社，2005：75.

[2]美國勞雷工業(yè)有限公司.G856質(zhì)子磁力儀用戶手冊[S].美國勞雷工業(yè)有限公司，1990：20-22.

[3]加拿大GEM公司.GSM-19T 操作手冊[S].加拿大GEM公司，2011：46-49.

[4]李才明，李軍，余舟等. 提高磁測日變改正精度的方法[J].物探化探計算技術(shù)，2004，26（3）.

[5]于寶顯，杜利明，王懷坤等.不同型號質(zhì)子磁力儀之間數(shù)據(jù)共享[J].內(nèi)蒙古煤炭經(jīng)濟(jì)，2013，5.

[6]劉炳文.Visual Basic 程序設(shè)計教程[M].北京：清華大學(xué)出版社，2003：216-223.

[7]John Walkenbach. Excel 2007 VBA Programming For Dummies[M].USA：Wiley Publishing，Inc，2007：36-52.

[8]曾華霖.重力場與重力勘探 [M].北京：地質(zhì)出版社，2005：196-197.

[9]李慶揚(yáng)，王能超，易大義.數(shù)值分析[M].北京：清華大學(xué)出版社，2008：23-24.

[10]張?zhí)炝?，姚慶梅，許明夫等.基于分段線性插值法的高精度測溫研究[J].信息技術(shù)與信息化，2012（1）.endprint

Sheet3.Cells（1， 1）.Value = "日變每一秒值"

Sheet3.Cells（1， 2）.Value = "日變時間"

Sheet3.Cells（1， 3）.Value = "觀測點(diǎn)時間"

Sheet3.Cells（1， 4）.Value = "觀測點(diǎn)觀測值"

Dim i As Integer

Dim j As Integer

For i = 2 To 108 ′此值為觀測點(diǎn)時間對應(yīng)的行數(shù)′

For j = 2 To 24142 ′此值為通過插值得到的數(shù)據(jù)對應(yīng)的行數(shù)′

If Sheet3.Cells（i， 3）.Value=Sheet3.Cells（j， 2）.Value Then

Sheet3.Cells（i， 4）.Value = Sheet3.Cells（j， 1）.Value

End If

Next j

Next i

End Sub

5 日變改正

提取測點(diǎn)與日變觀測同一時刻的日變數(shù)據(jù)后，用測點(diǎn)觀測的數(shù)據(jù)減去同一時刻的日變數(shù)值就是日變改正后的值。

6 驗(yàn)證

為了驗(yàn)證VBA程序的正確性，在野外磁測過程中做了實(shí)驗(yàn)驗(yàn)證。用一臺GSM-19T質(zhì)子磁力儀做日變觀測，兩臺GSM-19T質(zhì)子磁力儀做測點(diǎn)觀測，點(diǎn)距約為2米，點(diǎn)數(shù)為64個。測點(diǎn)觀測數(shù)據(jù)用GEMLinkW軟件和VBA程序進(jìn)行日變改正，用VBA程序日變改正時沒有做五點(diǎn)和七點(diǎn)滑動平均。日變改正后的數(shù)據(jù)見表1。從表1可以看出：GEMLinkW軟件和VBA程序日變改正后的值相等（用GEMLinkW軟件日變改正后的值軟件默認(rèn)進(jìn)行四舍五入保留了兩位小數(shù)，用VBA程序日變改正后的值進(jìn)行四舍五入保留了三位小數(shù)），證明了VBA程序的正確性。

此外，對VBA程序還進(jìn)行生產(chǎn)應(yīng)用驗(yàn)證，在云南省鳳慶縣某鐵礦區(qū)開展高精度磁測工作中，同時使用G856和GSM-19T質(zhì)子磁力儀進(jìn)行磁測數(shù)據(jù)采集，用VBA程序?qū)y區(qū)的20971個磁測數(shù)據(jù)進(jìn)行5%抽查人工日變改正驗(yàn)證，結(jié)果也表明VBA程序的正確性。

7 結(jié)論

在高精度磁測野外生產(chǎn)過程中經(jīng)常會遇到用不同型號的磁力儀進(jìn)行數(shù)據(jù)采集，數(shù)據(jù)采集完后必須進(jìn)行日變改正，然而不同種磁力儀之間往往不能用配套的軟件進(jìn)行日變改正。本文在Microsoft Office Excel平臺編寫五點(diǎn)或七點(diǎn)滑動平均程序和日變改正VBA程序，并且實(shí)驗(yàn)和生產(chǎn)應(yīng)用驗(yàn)證了此程序的正確性，適合野外現(xiàn)場計算，方便、快捷。

參考文獻(xiàn)：

[1]管志寧.地磁場與磁力勘探[M].北京：地質(zhì)出版社，2005：75.

[2]美國勞雷工業(yè)有限公司.G856質(zhì)子磁力儀用戶手冊[S].美國勞雷工業(yè)有限公司，1990：20-22.

[3]加拿大GEM公司.GSM-19T 操作手冊[S].加拿大GEM公司，2011：46-49.

[4]李才明，李軍，余舟等. 提高磁測日變改正精度的方法[J].物探化探計算技術(shù)，2004，26（3）.

[5]于寶顯，杜利明，王懷坤等.不同型號質(zhì)子磁力儀之間數(shù)據(jù)共享[J].內(nèi)蒙古煤炭經(jīng)濟(jì)，2013，5.

[6]劉炳文.Visual Basic 程序設(shè)計教程[M].北京：清華大學(xué)出版社，2003：216-223.

[7]John Walkenbach. Excel 2007 VBA Programming For Dummies[M].USA：Wiley Publishing，Inc，2007：36-52.

[8]曾華霖.重力場與重力勘探 [M].北京：地質(zhì)出版社，2005：196-197.

[9]李慶揚(yáng)，王能超，易大義.數(shù)值分析[M].北京：清華大學(xué)出版社，2008：23-24.

[10]張?zhí)炝?，姚慶梅，許明夫等.基于分段線性插值法的高精度測溫研究[J].信息技術(shù)與信息化，2012（1）.endprint