G856和GSM—19T质子磁力仪观测数据日变改正处理新方法

王荣华+李文尧

摘要： 使用G856和GSM-19T质子磁力仪在同一个测区进行磁测数据采集，由于这两种型号的质子磁力仪观测数据文件格式不一样，用其配套软件相互日变改正特别困难和繁琐，并且用其配套软件日变改正没有五点或七点滑动平均。为了解决这些问题，在Microsoft Office Excel平台编写五点或七点滑动平均和日变改正VBA程序。该程序经实验验证和一个测区20971个磁测数据生产应用验证表明VBA程序的正确性。这个VBA程序适用不同型号的质子磁力仪观测数据之间日变改正，只要把观测数据转换到Microsoft Office Excel中并插入VBA程序代码即可进行日变改正，非常实用、快捷。

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.

关键词： G856质子磁力仪；GSM-19T质子磁力仪；日变改正；VBA程序

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

中图分类号：P716+.82；P413 文献标识码：A 文章编号：1006-4311（2014）10-0001-03

0 引言

对同一测区使用G856和GSM-19T两种型号的质子磁力仪进行磁测数据采集，两种磁力仪观测数据输出文件格式不一样，并且两种磁力仪数据日变改正[1]软件不一样。G856质子磁力仪数据输出文件格式为stn[2]，日变改正软件为MagMap2000[2]；GSM-19T质子磁力仪数据输出文件格式为txt[3]，日变改正软件为GEMLinkW[3]。因此，这两种型号的磁测观测数据用配套软件相互日变改正特别困难和繁琐。此外，根据《地面高精度磁测技术规程》（DZ/T0071-93），为了提高日变改正精度[4]，先要对日变观测原始数据作五点或七点滑动平均，压制噪声后，再对测点观测数据进行日变改正，但G856和GSM-19T质子磁力仪配套日变改正软件日变改正过程中没有五点或七点滑动平均。针对这些问题，查阅了相关的文献和软件，管志宁[1]提出日变改正值从日变曲线上查得，在日变曲线上量得某时刻相对早基点时间的日变值并取反号，即为该时刻的日变改正值，这种日变改正方法速度慢、精度不高。于宝显[5]等提出不同型号质子磁力仪之间数据共享，经转换后的数据可以利用其中一种型号质子磁力仪所配套软件做日变改正，这种方法有些繁琐，并且没有做五点或七点滑动平均。中国地质大学（武汉）刘天佑教授开发的磁法勘探软件系统MAGS 2.0没有实现这两种型号质子磁力仪观测数据之间的日变改正处理。对此，通过Microsoft Office Excel平台VBA[6-7]编程解决。

1 观测数据格式统一

为了统一G856和GSM-19T质子磁力仪观测数据格式，用MagMap2000导入G856质子磁力仪stn格式文件数据，然后输出dat格式文件数据，最后把dat文件格式数据复制到Microsoft Office Excel文件中；GSM-19T质子磁力仪txt文件格式数据可以直接导入Microsoft Office Excel文件中。

2 日变观测数据五点或七点滑动平均

根据《地面高精度磁测技术规程》（DZ/T0071-93），地磁场短周期变化的振幅与质子磁力仪噪声均方根值是近于同一数量级，噪声是随机的，地磁场脉动变化是有规律的，先要对日变观测原始数据先作五点或七点滑动平均，压制噪声水平后，再对测点观测数据进行日变改正，即可提高日变改正的精度。下面介绍日变观测数据五点和七点滑动平均原理[7]。设在一个范围内某点经平滑后的日变数值可表示：T（x）=a0+a1x （1）endprint

（1）式中：a0和a1为待定系数。设日变观测的数值为T（xi），它的平滑值为T（xi），则可得出：

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

（2）式中δ为偏差的平方和。利用微分求极值的方法将（2）式对a■和a■求导数，然后令其为零得：

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

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

若x■以剖面上的点距为单位，步长为1，则式（3）中的x■=0，±1，±2…±m。代入式（3）中，可解得a■和a■分别为

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

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

由此可见，（5）式就是平滑公式，当m=±1，三点平滑公式为T（0）=■[T（-1）+T（0）+T（1）] （6）

同理可得五点平滑公式为

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

七点平滑公式为

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

根据五点和七点滑动平均公式（7）和（8）编写Microsoft Office Excel VBA程序，代码如下，利用此程序代码把日变原始观测数据（见图1）五点和七点滑动平均，见图2和图3。

Sub 日变原始数据圆滑（）

Dim i As Integer

Sheet1.Cells（1， 1） = "日变原始数据"

Sheet1.Cells（1， 2） = "日变观测时间"

Sheet1.Cells（1， 3） = "日变五点圆滑后数据"

Sheet1.Cells（1， 4） = "日变七点圆滑后数据"

n = 649 ′此值为原始日变数据对应的行数′

m = 5 ′选择5点或7点滑动平均圆滑′

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 日变观测数据线性插值

根据《地面高精度磁测技术规程》（DZ/T0071-93）质子磁力仪日变观测时间间隔为5秒至20秒，为了求出测点观测与日变观测同一时刻的日变数据，先要求出日变观测数据每一秒的数值，因此对日变数据进行分段线性插值[9-10]，公式为（9）式，公式中T（n）为前一时刻的日变观测数据，T（n+m）为后一时刻的日变观测数据，m为日变观测时间间隔，j为日变观测第j秒，T（j）为分段线性插值后第j秒数值，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） = "日变原始数据"

Sheet1.Cells（1， 2） = "日变插值后数据"

n = 1173 ′此值为原始的日变观测数据对应的行数′

m=20′此值为日变观测数据时间间隔′

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 提取测点与日变观测同一时刻的日变数据

求出日变观测每一秒数值后，要提取测点观测与日变观测同一时刻的日变数据，Microsoft Office Excel VBA 程序如下：

Sub 提取数据（）endprint

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

Sheet3.Cells（1， 2）.Value = "日变时间"

Sheet3.Cells（1， 3）.Value = "观测点时间"

Sheet3.Cells（1， 4）.Value = "观测点观测值"

Dim i As Integer

Dim j As Integer

For i = 2 To 108 ′此值为观测点时间对应的行数′

For j = 2 To 24142 ′此值为通过插值得到的数据对应的行数′

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 日变改正

提取测点与日变观测同一时刻的日变数据后，用测点观测的数据减去同一时刻的日变数值就是日变改正后的值。

6 验证

为了验证VBA程序的正确性，在野外磁测过程中做了实验验证。用一台GSM-19T质子磁力仪做日变观测，两台GSM-19T质子磁力仪做测点观测，点距约为2米，点数为64个。测点观测数据用GEMLinkW软件和VBA程序进行日变改正，用VBA程序日变改正时没有做五点和七点滑动平均。日变改正后的数据见表1。从表1可以看出：GEMLinkW软件和VBA程序日变改正后的值相等（用GEMLinkW软件日变改正后的值软件默认进行四舍五入保留了两位小数，用VBA程序日变改正后的值进行四舍五入保留了三位小数），证明了VBA程序的正确性。

此外，对VBA程序还进行生产应用验证，在云南省凤庆县某铁矿区开展高精度磁测工作中，同时使用G856和GSM-19T质子磁力仪进行磁测数据采集，用VBA程序对测区的20971个磁测数据进行5%抽查人工日变改正验证，结果也表明VBA程序的正确性。

7 结论

在高精度磁测野外生产过程中经常会遇到用不同型号的磁力仪进行数据采集，数据采集完后必须进行日变改正，然而不同种磁力仪之间往往不能用配套的软件进行日变改正。本文在Microsoft Office Excel平台编写五点或七点滑动平均程序和日变改正VBA程序，并且实验和生产应用验证了此程序的正确性，适合野外现场计算，方便、快捷。

参考文献：

[1]管志宁.地磁场与磁力勘探[M].北京：地质出版社，2005：75.

[2]美国劳雷工业有限公司.G856质子磁力仪用户手册[S].美国劳雷工业有限公司，1990：20-22.

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

[4]李才明，李军，余舟等. 提高磁测日变改正精度的方法[J].物探化探计算技术，2004，26（3）.

[5]于宝显，杜利明，王怀坤等.不同型号质子磁力仪之间数据共享[J].内蒙古煤炭经济，2013，5.

[6]刘炳文.Visual Basic 程序设计教程[M].北京：清华大学出版社，2003：216-223.

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

[8]曾华霖.重力场与重力勘探 [M].北京：地质出版社，2005：196-197.

[9]李庆扬，王能超，易大义.数值分析[M].北京：清华大学出版社，2008：23-24.

[10]张天亮，姚庆梅，许明夫等.基于分段线性插值法的高精度测温研究[J].信息技术与信息化，2012（1）.endprint

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

Sheet3.Cells（1， 2）.Value = "日变时间"

Sheet3.Cells（1， 3）.Value = "观测点时间"

Sheet3.Cells（1， 4）.Value = "观测点观测值"

Dim i As Integer

Dim j As Integer

For i = 2 To 108 ′此值为观测点时间对应的行数′

For j = 2 To 24142 ′此值为通过插值得到的数据对应的行数′

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 日变改正

提取测点与日变观测同一时刻的日变数据后，用测点观测的数据减去同一时刻的日变数值就是日变改正后的值。

6 验证

为了验证VBA程序的正确性，在野外磁测过程中做了实验验证。用一台GSM-19T质子磁力仪做日变观测，两台GSM-19T质子磁力仪做测点观测，点距约为2米，点数为64个。测点观测数据用GEMLinkW软件和VBA程序进行日变改正，用VBA程序日变改正时没有做五点和七点滑动平均。日变改正后的数据见表1。从表1可以看出：GEMLinkW软件和VBA程序日变改正后的值相等（用GEMLinkW软件日变改正后的值软件默认进行四舍五入保留了两位小数，用VBA程序日变改正后的值进行四舍五入保留了三位小数），证明了VBA程序的正确性。

此外，对VBA程序还进行生产应用验证，在云南省凤庆县某铁矿区开展高精度磁测工作中，同时使用G856和GSM-19T质子磁力仪进行磁测数据采集，用VBA程序对测区的20971个磁测数据进行5%抽查人工日变改正验证，结果也表明VBA程序的正确性。

7 结论

在高精度磁测野外生产过程中经常会遇到用不同型号的磁力仪进行数据采集，数据采集完后必须进行日变改正，然而不同种磁力仪之间往往不能用配套的软件进行日变改正。本文在Microsoft Office Excel平台编写五点或七点滑动平均程序和日变改正VBA程序，并且实验和生产应用验证了此程序的正确性，适合野外现场计算，方便、快捷。

参考文献：

[1]管志宁.地磁场与磁力勘探[M].北京：地质出版社，2005：75.

[2]美国劳雷工业有限公司.G856质子磁力仪用户手册[S].美国劳雷工业有限公司，1990：20-22.

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

[4]李才明，李军，余舟等. 提高磁测日变改正精度的方法[J].物探化探计算技术，2004，26（3）.

[5]于宝显，杜利明，王怀坤等.不同型号质子磁力仪之间数据共享[J].内蒙古煤炭经济，2013，5.

[6]刘炳文.Visual Basic 程序设计教程[M].北京：清华大学出版社，2003：216-223.

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

[8]曾华霖.重力场与重力勘探 [M].北京：地质出版社，2005：196-197.

[9]李庆扬，王能超，易大义.数值分析[M].北京：清华大学出版社，2008：23-24.

[10]张天亮，姚庆梅，许明夫等.基于分段线性插值法的高精度测温研究[J].信息技术与信息化，2012（1）.endprint

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

Sheet3.Cells（1， 2）.Value = "日变时间"

Sheet3.Cells（1， 3）.Value = "观测点时间"

Sheet3.Cells（1， 4）.Value = "观测点观测值"

Dim i As Integer

Dim j As Integer

For i = 2 To 108 ′此值为观测点时间对应的行数′

For j = 2 To 24142 ′此值为通过插值得到的数据对应的行数′

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 日变改正

提取测点与日变观测同一时刻的日变数据后，用测点观测的数据减去同一时刻的日变数值就是日变改正后的值。

6 验证

为了验证VBA程序的正确性，在野外磁测过程中做了实验验证。用一台GSM-19T质子磁力仪做日变观测，两台GSM-19T质子磁力仪做测点观测，点距约为2米，点数为64个。测点观测数据用GEMLinkW软件和VBA程序进行日变改正，用VBA程序日变改正时没有做五点和七点滑动平均。日变改正后的数据见表1。从表1可以看出：GEMLinkW软件和VBA程序日变改正后的值相等（用GEMLinkW软件日变改正后的值软件默认进行四舍五入保留了两位小数，用VBA程序日变改正后的值进行四舍五入保留了三位小数），证明了VBA程序的正确性。

此外，对VBA程序还进行生产应用验证，在云南省凤庆县某铁矿区开展高精度磁测工作中，同时使用G856和GSM-19T质子磁力仪进行磁测数据采集，用VBA程序对测区的20971个磁测数据进行5%抽查人工日变改正验证，结果也表明VBA程序的正确性。

7 结论

在高精度磁测野外生产过程中经常会遇到用不同型号的磁力仪进行数据采集，数据采集完后必须进行日变改正，然而不同种磁力仪之间往往不能用配套的软件进行日变改正。本文在Microsoft Office Excel平台编写五点或七点滑动平均程序和日变改正VBA程序，并且实验和生产应用验证了此程序的正确性，适合野外现场计算，方便、快捷。

参考文献：

[1]管志宁.地磁场与磁力勘探[M].北京：地质出版社，2005：75.

[2]美国劳雷工业有限公司.G856质子磁力仪用户手册[S].美国劳雷工业有限公司，1990：20-22.

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

[4]李才明，李军，余舟等. 提高磁测日变改正精度的方法[J].物探化探计算技术，2004，26（3）.

[5]于宝显，杜利明，王怀坤等.不同型号质子磁力仪之间数据共享[J].内蒙古煤炭经济，2013，5.

[6]刘炳文.Visual Basic 程序设计教程[M].北京：清华大学出版社，2003：216-223.

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

[8]曾华霖.重力场与重力勘探 [M].北京：地质出版社，2005：196-197.

[9]李庆扬，王能超，易大义.数值分析[M].北京：清华大学出版社，2008：23-24.

[10]张天亮，姚庆梅，许明夫等.基于分段线性插值法的高精度测温研究[J].信息技术与信息化，2012（1）.endprint