航空物探数据的网格化和成图

原文：http://goo.gl/OAhmWH

比较流行的航空物探勘查数据是 USGS 的航磁数据， 可以从 http://mrdata.usgs.gov/magnetic/surveys.php 下载。 有些数据集相当大， 有数以百万计的数据点， 如果电脑没有足够内存的话， 处理起来是相当费劲的。

选一个适中的数据， 比如内华达州咯普莱特的航磁数据。 对于插值计算这个数据集不够理想， 因为其点距和线距差别巨大。

要为手上的数据集选一个最佳的网格化方法， 一般按下面几步来评估各网格化选项：

用分类张贴图查看数据的分布情况

* 用经度、纬度、剩磁分别做 XYZ 列绘制分类张贴图。
* 分组方法选择 **Equal Intervals**。
* 组数增加到 7。
* 估算到组间距大约为 23， 第一组的最小值大约是 -288， 最后一组的最小值大约是 -154。
* 把分类符号的线色设置为黑色并加 10% 透明，符号改为实心正方形。
* 符号大小全部改为一个较小的值，如 0.005，方便后面估算点线距。
* 分别设置各组符号的填充颜色，如深蓝、青色、绿色、黄色、红色、洋红、粉红等。
* 将分类定义保存为 ClassedPostClasses.cls 以备后用。
* 将图形放大到可以分得清每个测点，用 **Map | Measure** 测量工具，估算到点距（X）约 0.0001°，线距（Y）约 0.005°。
* 将符号调大获得更好的视觉观感，如 0.2 或 0.5。

网格化方法海选

运行 GridData_Comparison.bas 脚本， 比较各网格化方法， 将选择范围缩小。

* 脚本只是用默认参数对数据进行网格化，然后绘制等值线图，做直观的比对。
* 一眼看去，立刻就可以剔除最近邻点法、径向基函数法和改良的谢别德方法。
* 检查剩下的 5 种方法所得的等值线图与原始数据匹配度。
* 将前面做的分类张贴图复制 5 份，每个都覆盖到一个等值线图之上。
* 等值线图的等级设置为最小 -288，间距 23，每个等级的填充色与分类张贴图对应。
* 5 个网格化方法看起来都不错，但是三角剖分法在放大后看起来略诡异，说明该线性插值方法在此处不太适用，弃之。
* 剩下 4 个网格化方法。

高分辨率比对

用上一步选出的网格化方法手工对数据进行网格化， 对更多的网格化参数进行微调。

* 用剩下的距离倒数乘方法、克里格法、自然邻点法和最小曲率法分别网格化数据。
* 用估算到的点距作为 X 方向的间距，为 0.0001。
* 用估算到的线距作为 Y 方向的间距，因为相比点距过于巨大，去1/5， 为 0.001。
* 用这些高分网格再次绘图，发现最小曲率法好污，弃之。
* 剩下 3 个网格化方法。

确定最终选用的网格化方法

* 剩下的 3 个看起来都那么好，真是乱花渐欲迷人眼，都不知道怎么选择了。
* 只好祭出终极大招：数学方法。
* 每个方法都计算对应的残差，残差列取平方，再用统计功能求和。
* 比较各方法的残差平方和，发现克里格方法最小，说明该方法在高分网格时更加忠实于原始数据。
* 认定在目前的参数下，克里格是最佳的网格化方法。
* 大多数时候我们都在求上苍赐予我们最佳的网格化参数，但其实啊，最浪费生命的是剔除不适用的网格化方法和网格化参数。

对于物探佬，在没有那么多生命可浪费的情况下，直接用克里格方法吧，因为该方法对于大多数物探数据都是适用的。

用脚本完成网格化方法比对

前面的手工步骤实在繁琐，即使当中使用了 Surfer 附带的网格化比较脚本，为了尽可能减少鼠标的移动和点击次数，把各步骤都集成到一个脚本中。

使用这个脚本之前，要确定一些参数，可以直接查看数据也可以生成分类张贴图来查看。

首先要确定的是 X、Y、Z 所在的列，把列号分别赋予脚本的 xCol, yCol, zCol 变量，选一个空白列，把列号赋予 rCol 变量，用来临时存放残差值。

然后要确定做高分辨率网格化时使用的网格间距，可以通过统计、计数或者测量的方式，估算数据 X、Y 方向的间距，把他们赋予脚本的 xSize, ySize 变量。

第一次运行脚本，应该是做粗略的比对，所有网格化方法都用默认参数参与进来，脚本的对应设置如下：

    bHighRes = False    '用默认参数做网格化
    For i = 1 To 12
        algSelected(i) = True    'True 表示该算法被选中
    Next

脚本运行后自动生成比对的等值线图。剔除掉不需要的网格化方法，把剩下的网格化方法序号记下来，比如经过比对，认为距离倒数乘方法、克里格法和自然邻点法可以做进一步的比对筛选，他们对应的序号为1、2、5，则脚本设置如下：

    bHighRes = True    'True 表示将进行高分辨率网格化
    For i = 1 To 12
        algSelected(i) = False    'False 表示该算法旁观
    Next
    algSelected(1) = True    'True 表示该网格化方法参加比对
    algSelected(2) = True
    algSelected(5) = True

再次运行脚本，生成3个叠加了分类张贴图的高分辨率等值线图，并且根据残差平方和推荐一个最佳网格化方法。

当然，所谓的最佳网格化方法，只是对应于当前所使用的网格化参数而言的，换一组参数，最佳就可能是另一种方法。

'=====================================================
'用 Surfer 13 的 12 种网格化方法对数据进行网格化，
'并绘制等值线图进行比较。
'运行此脚本之前，需要手工通过分类张贴图确定分组情况，
'并将此保存为 ClassedPostClasses.cls 文件。
'在 Scripter + Golden Software Surfer 13.x 测试通过。
'holz [AT] 21cn.com
'2016-4-12
'=====================================================

Option Explicit
Sub Main
    Dim SurferApp As Object, PlotDoc As Object
    Dim MapFrame As Object, MapLayer As Object
    Dim wks As Object, Axis As Object
    Dim xCol As Integer, yCol As Integer, zCol As Integer
    Dim rCol As Integer
    Dim datFile As String, grdFile As String
    Dim yLength As Single
    Dim i As Integer
    Dim grdAlgorithm(12) As String
    Dim algSelected(12) As Boolean
    Dim clsFile As String
    Dim lvlFile As String
    Dim xSize As Single, ySize As Single
    Dim bHighRes As Boolean

    '==== 以下参数在运行前应核实无误 ====
    xCol = 3        '数据文件中，哪一列用做 X 列？
    yCol = 4        '数据文件中，哪一列用做 Y 列？
    zCol = 11        '数据文件中，哪一列用做 Z 列？
    rCol = 12        '残差准备暂存到数据文件哪一列？

    '高分网格间距一般分别设定为数据 X、Y 方向的平均间距即可。
    '如果 X、Y 方向的平均间距差异过大，那么大的那个方向可取
    '稍小的值，譬如 1/5
    '举例：数据的 X、Y 平均间距分别为 1 和 50，差异巨大，
    '网格化时，网格 X、Y 间距可取 1 和 10
    xSize = 0.0001    '高分网格 X 间距
    ySize = 0.001    '高分网格 Y 间距
    yLength = 2.0    '每个地图的高度

    '把 bHighRes 设为 False 并把 algSelected 全设为 True，
    '可快速比对全部的 12 种网格化方法，剔除明显不适用的。
    '把 bHighRes 设为 True 并把 algSelected 全设为 False，
    '把需要做高精度比对的方法对应的 algSelected 设为 True，
    '可进行高分辨率的网格化，并计算残差平方和，推荐最佳方法。
    bHighRes = True    'True 表示将进行高分辨率网格化
    For i = 1 To 12
        algSelected(i) = False    'False 表示该算法旁观
    Next
    algSelected(1) = True    'True 表示该网格化方法参加比对
    algSelected(2) = True
    algSelected(5) = True

    grdAlgorithm(1) = "距离倒数乘方法"
    grdAlgorithm(2) = "克里格法"
    grdAlgorithm(3) = "最小曲率法"
    grdAlgorithm(4) = "改进的谢别德法"
    grdAlgorithm(5) = "自然邻点法"
    grdAlgorithm(6) = "最近邻点法"
    grdAlgorithm(7) = "多项式回归法"
    grdAlgorithm(8) = "径向基函数法"
    grdAlgorithm(9) = "含线性插值的三角剖分法"
    grdAlgorithm(10) = "滑动平均法"
    grdAlgorithm(11) = "数据度量"
    grdAlgorithm(12) = "局部多项式法"

    Debug.Clear
    '已知当前目录下有原始数据文件
    datFile = CurDir() & "\NV_1152.csv"
    '如果没有，则弹出对话框选择原始数据文件
    If Dir(datFile) = "" Then
        datFile = GetFilePath(,"csv|dat|txt",,"选择原始数据",4)
    End If
    '如果仍然未选择任何文件，程序结束
    If datFile = "" Then End
    Debug.Print "======= " & Date & " " & Time & " ======="
    Debug.Print "这一切始于： " & datFile

    clsFile = CurDir() & "\ClassedPostClasses.cls"
    If Dir(clsFile)="" Then
        Debug.Print "创建示范性分类文件 " & clsFile
        CreateCls(clsFile)
    End If
    lvlFile = CurDir() & "\GridComparison.lvl"
    If Dir(lvlFile)="" Then
        Debug.Print "根据分类文件 " & clsFile & vbCrLf & _
            "创建示范性等级文件 " & lvlFile
        CreateLvl2(clsFile, lvlFile)
    End If

    '创建一个 Surfer 实例
    Set SurferApp = CreateObject("Surfer.Application")
    '新增一个空白图形文档
    Set PlotDoc=SurferApp.Documents.Add(srfDocPlot)

    Dim r1 As Long            '开始的行号
    Dim r2 As Long            '结束的行号
    Dim eq As String        '变换公式字符串
    Dim RMS As Double        '用来保存最小的残差平方和
    Dim ssr As Double        '保存临时的残差平方和
    Dim AlgNo As Integer    '保存最佳的网格化方法序号
    Dim bStatus As Boolean    '网格化成功与否标记
    Dim xInc As Integer        '控制地图摆放的 X 向步进值
    Dim yInc As Integer        '控制地图摆放的 Y 向步进值

    xInc = 0
    yInc = 0
    r1 = 1
    RMS = 1.70141e+038        '要找最小的，先给个巨大的
    AlgNo = 1

    For i = 1 To 12
        Debug.Print i & "、" & grdAlgorithm(i) & "： " & _
            If(algSelected(i),"开始","旁观")
        If algSelected(i) Then
            If bHighRes Then
                Debug.Print "  做高分辨率网格化……"
                grdFile = Left(datFile,InStrRev(datFile,".")-1) & "_hd" & i & ".grd"
                bStatus = SurferApp.GridData3(DataFile:= datFile, _
                    OutGrid:=grdFile, ShowReport:=False, Algorithm:=i, _
                    xCol:=xCol, yCol:=yCol, zCol:=zCol, _
                    BlankOutsideHull:=1, xSize:=xSize, ySize:=ySize)
                Debug.Print "  计算高分网格残差……"
                Set wks = SurferApp.GridResiduals2(grdFile, datFile, _
                    xCol, yCol, zCol, rCol)
                r1 = 1
                r2 = wks.Columns(rCol).Count
                eq = Chr(64 + rCol)
                eq = eq & "=" & eq & "*" & eq
                Debug.Print "  残差求平方公式： " & eq
                wks.Transform3(wksTransformWithCol, r1, r2, eq)
                ssr = wks.Columns(rCol).Statistics(wksStatsSum).Sum
                Debug.Print "  残差平方和： " & ssr
                If RMS > ssr Then
                    RMS = ssr
                    AlgNo = i
                End If
                wks.Close(srfSaveChangesNo)
                rCol = rCol + 1
            Else
                Debug.Print "  用默认参数做网格化……"
                grdFile = Left(datFile,InStrRev(datFile,".")-1) & "_m" & i & ".grd"
                bStatus = SurferApp.GridData3(DataFile:= datFile, _
                    OutGrid:=grdFile, ShowReport:=False, Algorithm:=i, _
                    xCol:=xCol, yCol:=yCol, zCol:=zCol)
            End If
            Debug.Print "  网格化数据： " & If(bStatus,"","不") & "成功。"

            '绘制等值线图
            Set MapFrame = PlotDoc.Shapes.AddContourMap(grdFile)
            Set MapLayer = MapFrame.Overlays(1)
            If bHighRes Then
                '启用等值线填充
                MapLayer.FillContours = True
                '调用等级文件对等值线做填充
                MapLayer.Levels.LoadFile(lvlFile)
                '添加分类张贴图图层
                Set MapLayer = PlotDoc.Shapes.AddClassedPostLayer(MapFrame, _
                    datFile, xCol, yCol, zCol)
                '调用分类文件对张贴图进行设置
                MapLayer.LoadClasses(clsFile)
            End If

            With MapFrame
                '地图的高度设为指定高度
                .yLength = yLength
                '让 XY 同比例
                .xMapPerPU = .yMapPerPU
                '限定地图宽度不得超出规定大小
                If .Width > 3.5 Then .xLength = 3.5
                '水平方向把地图移到指定位置
                .Left = (.xLength + 1.0) * xInc + 2.0
                '垂直方向把地图移到指定位置
                .Top = (.yLength + 1.0) * yInc + 5.0
                yInc = yInc + 1
                If yInc > 2 Then
                    yInc = 0
                    xInc = xInc + 1
                End If

                '循环设置各个坐标轴
                For Each Axis In .Axes
                    '主刻度不要
                    Axis.MajorTickType = srfTickNone
                    '次刻度也不要
                    Axis.MinorTickType = srfTickNone
                    '刻度标注还是不要
                    Axis.ShowLabels = False
                    '如果是顶坐标
                    If Axis.AxisType = srfATTop Then
                        '写个标题注明是某算法
                        Axis.Title = i & ". " & grdAlgorithm(i)
                        '标题字体改为中文字体才能显示中文
                        Axis.TitleFont.Face = "宋体"
                        '改一个标题颜色
                        Axis.TitleFont.ForeColorRGBA.Color = srfColorBlue
                    End If
                Next
            End With
        End If
    Next

    If bHighRes Then
        Debug.Print "最佳网格化方法： " & grdAlgorithm(AlgNo)
        Debug.Print "最小残差平方和： " & RMS
        MsgBox "网格化算法： " & grdAlgorithm(AlgNo) & vbCrLf & vbCrLf & _
            " 算得最小残差平方和： " & RMS & vbCrLf & vbCrLf & _
            "建议选用 " & grdAlgorithm(AlgNo) & " 对数据进行网格化。", _
            vbOkOnly + vbInformation, "网格化方法比对完成"
    End If
    '让图形缩放到适应窗口
    SurferApp.ActiveWindow.Zoom(srfZoomFitToWindow)
    '让 Surfer 主窗口可见
    SurferApp.Visible=True
    Debug.Print "======= " & Date & " " & Time & " ======="
End Sub

Sub CreateCls(clsFile As String)
    Open clsFile For Output As #1
    Print #1, "ClassData 2 1"
    Dim sSym As String
    Dim sFill As String
    Dim sLine As String
    Dim s As String

    sSym = Chr(34) & "GSI Default Symbols" & Chr(34)
    sFill = " " & Chr(34) & "Navy Blue" & Chr(34)
    sLine = " " & Chr(34) & "R0 G0 B0 A26" & Chr(34) & " 1 10 0.019685"
    s = "-288.23 -265.94 " & sSym & sFill & sLine
    Print #1, s
    sFill = " " & Chr(34) & "Cyan" & Chr(34)
    s = "-265.94 -243.65 " & sSym & sFill & sLine
    Print #1, s
    sFill = " " & Chr(34) & "Green" & Chr(34)
    s = "-243.65 -221.36 " & sSym & sFill & sLine
    Print #1, s
    sFill = " " & Chr(34) & "Yellow" & Chr(34)
    s = "-221.36 -199.07 " & sSym & sFill & sLine
    Print #1, s
    sFill = " " & Chr(34) & "Red" & Chr(34)
    s = "-199.07 -176.78 " & sSym & sFill & sLine
    Print #1, s
    sFill = " " & Chr(34) & "Magenta" & Chr(34)
    s = "-176.78 -154.49 " & sSym & sFill & sLine
    Print #1, s
    sFill = " " & Chr(34) & "Pink" & Chr(34)
    s = "-154.49 -132.1 " & sSym & sFill & sLine
    Print #1, s
    Close #1
End Sub

Sub CreateLvl(lvlFile As String)
    Open lvlFile For Output As #1
    Print #1, "LVL3"
    Print #1, Chr(39) & "Level Flags LColor LStyle LWidth " & _
        "FVersion FFGColor FBGColor FPattern OffsetX " & _
        "OffsetY ScaleX ScaleY Angle Coverage"
    Dim s1 As String
    Dim s2 As String
    s1 = " 0 " & Chr(34) & "Black" & Chr(34) & " " & Chr(34) & _
        "Solid" & Chr(34) & " 0 1 " & Chr(34)
    s2 = Chr(34) & " " & Chr(34) & "Black" & Chr(34) & " " & _
        Chr(34) & "Solid" & Chr(34) & " 0 0 1 1 0 0"
    Print #1, "-288.23" & s1 & "Navy Blue" & s2
    Print #1, "-265.94" & s1 & "Cyan" & s2
    Print #1, "-243.65" & s1 & "Green" & s2
    Print #1, "-221.36" & s1 & "Yellow" & s2
    Print #1, "-199.07" & s1 & "Red" & s2
    Print #1, "-176.78" & s1 & "Magenta" & s2
    Print #1, "-154.49" & s1 & "Pink" & s2
    Close #1
End Sub

Sub CreateLvl2(clsFile As String, lvlFile As String)
    If Dir(clsFile) = "" Then
        Debug.Print "找不到文件 " & clsFile
        Exit Sub
    End If

    Open lvlFile For Output As #1
    Print #1, "LVL3"
    Print #1, Chr(39) & "Level Flags LColor LStyle LWidth " & _
        "FVersion FFGColor FBGColor FPattern OffsetX " & _
        "OffsetY ScaleX ScaleY Angle Coverage"

    Open clsFile For Input As #2
    Dim s As String
    Line Input #2, s
    Dim s1 As String
    Dim s2 As String
    s1 = " 0 " & Chr(34) & "Black" & Chr(34) & " " & Chr(34) & _
        "Solid" & Chr(34) & " 0 1 " & Chr(34)
    s2 = Chr(34) & " " & Chr(34) & "Black" & Chr(34) & " " & _
        Chr(34) & "Solid" & Chr(34) & " 0 0 1 1 0 0"
    While Not EOF(2)
        Line Input #2, s
        Print #1, Split(s)(0) & s1 & Split(s,Chr(34))(3) & s2
    Wend
    Close #2
    Close #1
End Sub