使用Sufer软件求取高程异常值的方法

　2008年2月　　　　　　　　　　　　　　物　探　装　备

第18卷　第1期

・重磁电技术・

使用Sufer 软件求取高程异常值的方法

王福民　刘文华

(东方地球物理公司综合物化探事业部, 河北涿州072751)

3

摘　　　要

王福民, 刘文华. 使用Su fer 软件求取高程异常值的方法. 48　　Sufer , , 已经成为处理三维数据的首选软件, 。Sufer 软件操作简单和绘图能利用已知点的值求取待定点的高程异常值的方法。关键词　GPS 　高程异常值　拟合　插值

ABSTRACT

W ang Fumin and Liu Wenhu a. Using Su fer softw are to compute abnorm al elevation value. EGP ,2008,18(1) :48～51

　　The Sufer software is mainly used to draw isoline map and relevant 32D image. Its strong interpolation f unction

and the ability of drawing figure has become first 2selected software used for 32D data processing , which can rapidly transform the measured data at discrete points into continuous data curved surface. The Sufer characters of simple operation and strongly drawing capacity have been known by many technicians. The paper introduce in detail the method how to use Sufer to utilize known points to compute the abnormal elevation values at unknown points.

K ey w ords 　GPS , abnormal elevation value , fitting , interpolation

　　GPS 定位精度高、速度快、操作简单、布网灵活。对于GPS 用户来说, 最关心的问题是将获得的大地高转换成我国1956黄海高程系下的正常高(海拔高) 。高程异常值即似大地水准面到参考椭球面的距离是否准确, 将直接关系到最终的高程精度。

标转换后得到地方椭球的大地高。

由于参考椭球面可以用数学方法表示出来, 而似大地水准面则不能, 这样H 与H g 就不能建立严密的关系表达式, 这也就是GPS 测量正常高的困难所在。传统研究大地水准面起伏的方法有天文水准或天文重力水准测量法, 但其精度较低, 我国境内似大地水准面精度为0. 3～0. 6m , 西部更低一些。实际施工时, 受到项目成本和施工期限限制, 天文和重力资料不易获得, 很难将该方法应用到实际生产中。

利用已有的高程异常资料和水准资料, 与GPS 测量相结合, 采用高程拟合的方法是一种行之有效的求取高程异常的方法。

1　GPS 正常高测量原理

传统的高程测量方法得到的是测点的正常高, 也就是点位到似大地水准面的高程, 正常高用H g 表示; 似大地水准面到参考椭球面的距离, 称为高程异常, 记为N 。大地高(H ) 与正常高(H g ) 之间的关系如下:

　　　　　　H =H g +N

(1)

由于GPS 测量的坐标系统是W GS -84, 要想得到地方椭球的大地高, 首先要求得W GS -84坐标系统到北京1954坐标系统的坐标转换参数, 经坐

2　常用的高程异常求取方法

2. 1　等值线图法

等值线图法就是从国家测绘局出版的高程异常

3王福民, 男, 高级工程师,1991年毕业于华东地质学院工程测量专业, 从事重力、GPS 大地测量工作多年, 现在东方地球物理公司综

合物化探事业部海外部工作。

　第18卷　第1期王福民等:使用Sufer 软件求取高程异常值的方法

n

49　　

图上量出测点的高程异常值, 然后计算出该点的正常高。2. 2　地球模型法

地球模型法本质上是一种数字化的等值线图, 目前国际上较常用的地球模型有EGM96、OSU91A 等, 这些模型比较适合于美国本土。在一些没有高程异常资料的国家施工时, 可考虑使用这些模型。

在我国部分城市和地区, 已做过很多大地水准面精化工作, 如在新疆塔里木地区施工时, 可以使用该地区的大地水准面精化模型。2. 3　高程拟合法

　　　N (z ) =

n

z u

[d i (x , y ) ]

i =1

[d (x , y ) ]

i

(2)

u

84, 过GPS , 便可得到这些水准点的W GS -84坐标和经过坐标转换后得到的地方椭球平面坐标和大地高, 结合已有的水准资料, 我们就可以计算出这些水准点的高程异常值。通过水准点的高程异常值, 便可拟合出测区的大地水准面。当描述测区的大地水准面数学模型建立后, 根据测点的位置参数, 便可计算任意点的高程异常值。

这种方法要求GPS 联测已知水准点, 且水准点分布要均匀, 要有足够的密度。当已知水准点数量不够时, 可以通过水准测量的方法来加密, 再通过GPS 联测测定其大地高。

值得注意的是, 在小范围和平原地区施工时, 高程异常变化较为平缓, 其拟合的准确度在1dm 以内。对于高程异常变化剧烈的地区(如山区) , 这种方法的准确度有限, 这主要是因为在这些地区, 已知点上的高程异常值很难将似大地水准面的特征表示出来。

当测区的已知水准点数量少或其它基础资料不易获得时, 可在工区范围内选择一定数量且分布均匀的测点, 其高程异常值可在高程异常等值线图上量取, 然后再通过这些已知测点的高程异常值对待定点进行拟合。常用拟合方法有以下几种。2. 3. 1　反距离加权插值法

又称谢别德法(Shepard ) 。设平面上分布一系列离散点(x i , y i ) , 已知其位置坐标和高程异常值; N (x , y ) 为任一待定点, 根据其周围离散点的高程异常值, 通过距离加权插值法求出N 点的高程异常值。周围离散点因与N 点位置差异, 对N (z ) 点影响也不相同。我们把这种影响称为权函数, 它主要与距离相关。

式中, d i (x , y ) =x -x i ) 2+(y -y i ) 2, 表示离散点到N 点的距离; N (z ) 为要求的待定点的高程异常值; z i 为离散点的高程异常值。

由式(2) N 点距离, 权重迅速减小, 权u 。

, u 可以取2, 但待定点到已知点的距离不宜过长。2. 3. 2　平面拟合法

平面拟合法也称一次多项式拟合。它是一种线性拟合法, 适合平原地区施工, 其模型比较简单:

ζ=a 0+a 1d B +a 2d L (3) 　　　　

根据已知数据的高程异常, 组成观测值的误差

方程式, 按最小二乘法求解拟合系数a 0, a 1, a 2。2. 3. 3　曲面拟合法

曲面拟合包括二次曲面拟合、三次曲面拟合、多面函数拟合。

二次曲面拟合模型如下(二次多项式) :ζ=a 0+a 1d B +a 2d L +a 3d B 2　　

(4) +a 4d L 2+a 5d B d L

其中:

d B =B -B 0; d L =L -L 0;

B 0=

n B ; ∑L 0=

n

L ∑

;

n 为已知高程异常的点数。

通过最小二乘法可以求解出多项式的系数。

在工区范围较大、高程异常变化剧烈的地区施工时, 用低阶多项式拟合, 其精度必然很差; 而高次多项式拟合, 又可能产生解的不稳定性。因此, 最好把区域分成若干块, 对各分块采用不同的函数进行拟合, 并且要考虑相邻分块函数间的连续性。

3　使用Sufer 软件求取高程异常值

使用Sufer 软件求取高程异常值可以看作是前面介绍的等值线图法的进一步细化和延伸。我们知

50　　　　　　　　　　　　　　　　　　　　　物　探　装　备

2008年

道, 等值线图法在求取少量测点的高程异常值时, 操作简便; 但当点数较多时, 不仅工作量大且量取误差较大, 同时不可避免的伴随着粗差的产生。而使用Sufer 软件求取高程异常值的方法正好能避免上述缺点。

Sufer 软件具有强大的插值功能和绘制图件的能力, 它采用一定的格网化方法(即数学模型) 对不

Data (数据) ”, 出现“生成格网文件”对话窗, 如图1

所示

。

规则分布的原始数据点进行插值, 生成在原始数据

分布范围内规则间距的数据点分布。下面详细介绍该方法的使用过程。3. 1　准备工作

(1) 行扫描, 1∶250万全国高程异常图或1∶100万高程异常图(W GS -84坐标系统) 等, 选定的扫描范围应包含整个施工区域。

(2) 对扫描后的高程异常图进行数字化。

首先将扫描后的图像文件调入数字化软件, 然后进行坐标校正工作; 顺时针方向选择四个(或四个以上) 经纬网交点, 然后输入交点的真实坐标进行坐标校正; 选择交点时, 应将图形尽可能地放大。

坐标校正后, 即可沿着高程异常等值线数字化, 数字化时, 也应将图形尽可能地放大。采点时, 十字丝放到等值线的中心部位, 采点的密度可灵活掌握, 高程异常变化较大的地方, 应尽可能增加采点的密度。最后形成以纵坐标、横坐标、高程异常值(x , y , z ) 三维格式的数字化文档。

数字化软件可采用Ozi Explorer 或其他专业数字化软件。工区范围较大时, 建议使用Ozi Explorer 软件。该软件可直接输入大地坐标对扫描图像进行坐标校正(高程异常图的坐标为大地坐标) , 能有效减少坐标校正带来的误差。3. 2　使用Sufer 软件的插值功能, 求取高程异常值

在使用Sufer 软件进行插值计算之前, 首先要选用一定的内插法对离散数据格网化。高程异常图数字化后得到已知点的属性可以用x 、y 、z 值来表示。由于工区范围内包括无数个测点, 所以不可能对所有测点的z 值(高程异常值) 进行采集。离散数据格网化就是根据已知点上的高程异常值求出其他待定点的属性值。3. 2. 1　插值基本参数

图1　离散数据生成格网文件对话框

(1) Data Columns (数据列) 组

指定数字化后高程异常文件的x 、y 、z 数据的

列号, 根据实际情况为数据选择适当的列号。

在“Filter (滤波) ”对话窗, 设置重复数据点的处理方式。由于零距离加权格网化算法会产生错误, 因此格网化命令会对输入数据进行检查, 以确定是否出现相同位置(x , y ) 的测点。若出现位置参数相同的测点时, 用户可以有以下3种选项:

选“None ”:当用户确认使用的数据不包含重复数据点时, 可采用忽略方式;

选“Average ”:用于平均操作。将所有相同坐标测点的属性值平均, 然后由新点代替所有相同坐标的测点;

(数据排出过滤器) :选“Data Exclusion Filter ”

用一个布尔表达式来指定排出方式(如x >100或y

(2) 格网化方法(Gridding Met hod ) 组

Sufer8提供的插值方法有12种, 不同的内插方

法又有不同的“Advanced Options (高级选项) ”。本文选用克里金插值法(Kriging ) , 是一种在诸多领域应用广泛的格网化方法。

克里金插值法首先考虑待定点的空间属性在空间位置上的变异分布, 确定对一个待插点值有影响的距离范围, 然后用此范围的采样点来估计待定点的属性值。它满足无偏估计和方差最小两个条件, 又称最优内插法。在数据点多时, 其内插结果可靠性高。

(3) Grid line Geomet ry (格网线几何学)

打开Sufer 软件, 选择菜单命令“Grid (格网) /

　第18卷　第1期王福民等:使用Sufer 软件求取高程异常值的方法51　　

用来定义格网范围和格网密度。格网范围用“Maximum ”和“Minimum ”来确定, 格网密度由格网内的行、列的数量来确定。

“Spacing (间隔值) ”指定格网范围和格网密度后, 系统自动计算生成。

(4) Outp ut grid file (输出格网文件) 指定生成的格网文件保存路径和文件名, 文件名后缀为3. grd 。

(5) Cross validation (交叉验证)

用来评价某一指定插值方法的质量。果可用来评价空间变异性, 3. 2. 2　(1) 打开Sufer , “Grid (网格) /Residuals ”。

按照对话框的提示, 打开格网文件所在路径和文件名

。

　　例如, 选择待定点文件E 列为大地高程, 计算的正常高所在的列为F , 那么待定点的高程异常值N 由E 列值减F 列值计算求得(N =E -F ) , 如图2所示。

(4) 选择图2中的菜单命令“File (文件) /Save (保存) ”, 输入插值计算结果保存路径和文件名。

4, 插值计算准确、科学, 可批量处理数据; 而采用拟合的方法, 需要编写拟合软件, 其计算方法、计算模型需要验证。

(2) 可以使用国家测绘局出版的CQ G2000或者1∶100万或1∶250万的高程异常图。高程异常资料容易搜集。

(3) 数字化过程中减少了人为在图上量取的误差, 避免了粗差。同时数字化取值范围较大, 采点较密, 有利于同以后工区图的拼接。使用Sufer 软件求取的高程异常值与图上量取的结果相比较, 误差一般不大于5cm 。

(4) 高程异常变化较大的地区(山区) , 同样能够保证插值求取的精度。而采用软件拟合的方法, 由于已知高程异常的测点数量有限, 很难准确描述出高程异常值的变化规律。

参

考

文

献

图2　待定点插值计算对话框

(2) 选择要求取高程异常值的待定点文件所在

目录, 并选择该文件; 注意:此处的x 列选择横坐

标, y 列选择纵坐标, z 列选择大地高程。

(3) “Store Residuals In Column ”选择通过插值计算后, 求取的正常高所在的列。

1　王建, 白世彪, 陈晔编著. Sufer8地理信息制图. 中国地

图出版社,20042　美国G olden 公司. Sufer8软件说明书. 2002

收稿日期:2007204204