用 MapGIS绘制拟断面 (切片)图

彭亮华

(湖南省有色地质勘查局二四七队,湖南长沙 410129)

0 前言

在矿产地质的普查与详查过程中,需要进行一系列的地球物理测量工作,从而会采集到大量的地球物理数据信息。如何对采集到的数据作相应处理,并采用 MapGIS软件绘制对应的地球物理图件,作者对此有一点经验。注意到王声喜、康宝林[1]的文章后 (以下简称《化探图件》),考虑物探图件与化探图件的编制基本相似,因此在这里只谈用MapGIS绘制物探方法中测深装置的拟断面图,供同行们参考。

1 数学模型

这里介绍一种逆时针旋转θ°、等面积的变换,其数学模型如图 1所示。ABFG多边形是变换前的拟断面图区域,ACKH是逆时针旋转θ°后等面积的拟断面图区域。

以 A点为原点,AG为旋转轴逆时针旋转θ°后 ,E点变换到 R点 ,∠EAR=θ°,A、B、C、R、E各点的坐标如图 1所示。不难知道,Xr=AR＊ cosθ°、Yr=AR ＊ sinθ°、AR=AE=(Xe-X0)。B点变换到 C点,其坐标变换为:

图1 等面积变换的数学模型Fig.1 Mathematical model of equal Area tranformation

从式 (1)可知,旋转后点的 X轴长只有原来的cosθ°倍。为了使旋转后的面积与旋转前的面积相等,那么 Y轴长需扩大 (1/cosθ°)倍 (见下页图 2),从式 (2)就反映了这一点。从式 (2)可知,旋转后的点在 Y轴上己向上平移了 (X56-X0)＊sinθ°段距离。下面把式 (1)、式 (2)改写成通用的公式:

图2 旋转变换前、后面积相等Fig.2 Equal area of before and after rotation

其中 θ°为旋转角;X′n为旋转后第N个点的X坐标;Y′n为旋转后第N个点的Y坐标;Xn为旋转前第 N个点的 X坐标;Yn为旋转前第 N个点的 Y坐标;X0为旋转前左下角原点的 X坐标;Y0为旋转前左下角原点的 Y坐标。

式 (3)、式 (4)是逆时针旋转θ°等面积变换拟断面图区域的数学公式。其中 X0、Y0可根据情况作取舍,如需把旋转后的拟断面图放在旋转前的拟断面图一块进行比较,则要输入 X0、Y0,可能还要在式 (3)中加一平移常数,然后才能进行比较。

2 数据整理

2.1 GPS数据整理

GPS的参数设置及坐标数据 (另存为制表符分隔“＊.txt”的文本文件)与《化探图件》中所描述的一致,下面作者将对 DX、DY、DZ、DA、DF参数的计算作补充说明。

2.1.1 坐标系统参数计算

搜集应用区域内 GPS“B”级网三个以上网点WGS84坐标系 B、L、H值,以及我国坐标系 (BJ54或西安 80)B、L、h、x值。(注:B、L、H分别为大地坐标系中的大地纬度、大地经度及大地高,h、x分别为大地坐标系中的高程及高程异常。各参数可以通过各省级测绘局或测绘院具有“A”级、“B”级网的单位获得。)

计算不同坐标系三维直角坐标值,计算公式如下:

X=(N+H)＊cosB＊cosL

Y=(N+H)＊cosB＊ sinL

Z=[N(1-E＊ E)+H]＊ sinB

不同坐标系对应椭球的有关常数详见表 1。

表 1 不同坐标系统对应椭球的有关数据Tab.1 Relativedatatableofellipsoidcorrespondingtodifferentcoordinationsystem

在表 1中,X、Y、Z为大地坐标系中的三维直角坐标:A为大地坐标系对应椭球之长半轴;E＊E为大地坐标系对应椭球第一偏听偏信心率;F为对应椭球之扁率;N为该点的卯酉圈曲率半径,N=A/(1-E＊ E＊ sinB＊ sinB)∧0.5;H=h+x,该处 H为BJ54或西安 80坐标系中的大地高。

2.1.2 求出 DX、DY、DZ、DA、DF

利用WGS84坐标系的 X、Y、Z及 A、F值,减去我国坐标系的对应值,可以得出实现坐标系统转换的五个参数 (DX、DY、DZ、DA、DF),还应算出WGS84与北京和西安坐标系二套参数。

2.1.3 参数验证

参数计算之后必须进行验证。验证的方法是在应用区域内选择五个以上水准点进行实测,然后进行参数计算之后必须进行的验证。实测值与测绘部门提供的理论值对比,如果最大误差不大于15m,平均误差不大于 10m,则计算出的参数可以使用;否则,要重新计算或查找出现问题的原因。

例如:某 8号剖面井口的WGS84坐标是 B=25.7809°;L=113.16472°;H=222。

BJ54坐标是 X=2852880;Y=416160;H=222。

B=25.7809°=25∶46∶51.200

L=113.16472°=113∶09∶53.000

把B值与 L值输入程序三参数计算中 (见图3),右边输入 BJ54坐标,经确认后得 DX、DY、DZ(见图 4)。选三至五个水准点进行实测,验证参数。

DX=-14.065,DY= -109.798

DZ= -51.947

图3 三参数之计算Fig.3 Calculating the three parameter

图4 经计算后的三参数值Fig.4 After calculated the three parameter

DA=-108(WGS84之长半轴 -BJ54之长半轴)

DF=+0.00000048079548834552(WGS84之扁率 -BJ54之扁率)

一般地 DA= -108,DF= +0.00000050以不变值输入 GPS与Mapsource软件中,不过你会发现Mapsource软件中位置设置界面之 Delta会发生变化 (4.99999998002099e-007),见图 5。

图5 系统自动修正的 DF值Fig.5 System automatically correction of the DFvalue

2.2 剖面测深点的数据整理

2.2.1 测深点物探数据的整理

以正东西向测深剖面某测深点为例,一般以GPS观测数据的 Y坐标值标注为测深图件的 X轴坐标值,列为第一列;以 GPS观测数据的高程坐标值 (也可从地形图上取高程)为起点,减去该测深点不同记录位置的距离,并标注为测深图件的 Y轴坐标值,列为第二列;以该测深点的视电阻率、视极化率参数值分别列为第三列和第四列。排完第一个测深点后,接着往下排第二个、第三个测深点,直到排完同一剖面上的测深点为止。在电子表格里生成工作表后,另存为不规则网 (前加 notgrid)逗号分隔符的“展点.csv”文件之后,变换为“展点.det”文件 ,“展点.det”文件里的坐标未旋转 。

第一列、第二列为坐标数据,且 GPS的数据是以米为单位,因此要注意数据单位、比例尺与Map-GIS软件规定之间的关系。

在实际工作中,对测深剖面的布设,往往取决于地层界线、地质构造的走向。正东西向、正南北向的剖面,比较容易绘制拟断面图,但北东向 (一、三象限)、北西向 (二、四象限)的剖面,比较难绘制拟断面图 (指绘制的拟断面图方向与实际剖面方向一致时)。现介绍正东西向的测深剖面,在以西端为轴 (原点)逆时针旋转θ°后,其等面积拟断面图数据的整理。

把列 A、列 B的坐标数据,按式 (3)、式 (4)重新生成旋转θ°角后生成的新的坐标数据;把列 C、列 D数据与新的坐标数据,配成 (新建制表符分隔的文本文件 )列 A、列 B、列 C、列 D,依次是 X、Y(轴坐标数据),视电阻率、视极化率数据,如图 6所示。在电子表格里,把每条剖面整理成一个工作表,并另存为制表符分隔“旋转数据.txt”的文本文件。

2.2.2 地形剖面和裁剪框的形成

在电子表格里以 GPS观测的 Y坐标值标注为X轴坐标值,列为列 A;以 GPS观测数据的高程坐标值 (也可从地形图上取高程)标注为 Y轴坐标值,列为列 B。列 A、列 B按式 (3)、式 (4)重新生成旋转θ°角后形成新的坐标数据列 C、列 D,然后把列 C、列 D数据另存为制表符分隔“＊.txt”的文本文件,在投影变换中生成线,即地形剖面线“地形.WL”,见下页图 7。

在电子表格里,以 GPS观测的 Y坐标值标注为 X轴坐标值,列为第一列;以 GPS观测数据的高程坐标值 (也可从地形图上取高程)为起点,减去该测深点最近、最远记录位置的距离,并标注为 Y轴坐标值,列为第二列。把第一列、第二列按式(3)、式 (4)重新生成旋转θ°角后形成新的坐标数据列 C、列 D,把列 C、列 D数据另存为制表符分隔“＊.txt”的文本文件。应注意,数据要按顺序排列生成一个大闭合圈,并在投影变换中生成线框。该线框即为裁剪框,注意把线框生成区文件“裁剪框.WP”。

图6 旋转 30°的坐标变换Fig.6 30°rotation of the coordinatetrans for mation

3 拟断面图的绘制

3.1 将 TXT文本、DET文件投影成 MapGIS中的点、线、面文件

3.1.1 生成测深数据记录点点文件

双击 MapGIS软件的主界面,[空间分析 ]→[DTM分析 ],进入MapGIS数字地面模型子系统,选择菜单项 [模型应用 (A)]→[高程点标注制图 ],打开“展点.det”文件,见下页图 8。X轴方向数据选列A,Y轴方向数据选列 B,L标注数据选列 C或列D。按所标注的内容进行选择,比如要标注视电阻率参数,那么在本例中选列 C。子图符号选 247,红色,角度选0°,符号的尺寸一般选1～3不等。标注位置则取决于个人习惯,标注字体一般选宋体、黑色、3mm(如字体大小选 0mm,可绘制实际材料图);标注格式一般取整,视情况而定。按 [确认 ]键后,即可生成“展点.WT”文件。

图7 生成旋转后的裁剪框数据Fig.7 Generation of the cutting box data after rotating

图8 高程点标注选项Fig.8 Marked elevation point option

图9 离散数据网格化选项Fig.9 Discrete data grid options

图10 测深区域数据裁剪、无效化Fig.10 The cutting sounding data in the regional is invalid

3.1.2 等值线图的绘制

等值线图的绘制基本与《化探图件》[1]相一致,因步骤与地形有关,在此重复一下。

(1)第一步:先双击 MapGIS软件的主界面,[空间分析 ]→[DTM分析 ],进入 MapGIS数字地面模型子系统,选择菜单项 [Grd模型 ]→[离散数据网格化 ],打开“旋转数据.txt”文件,见图 9。在数据列 D中 X、Y为坐标轴数据选项,Z为所选绘图参数的高程数据。按 [确定 ]生成“＊.GRD”文件,如图 9所示,生成“4M.GRD”文件。

(2)第二步:在 MapGIS数字地面模型子系统中,单击“打开高程数据文件”快捷图标,选择刚生成的“＊.GRD”文件 (如“4M.GRD”)打开后,选择菜单项[Grd模型 ]→[规则网裁剪,无效化 ],见图10。在裁剪区域选取方式中,点击 [选取工作区中区域 ],并点击随后图标,选上述“裁剪框.WP”区文件打开。操作方式选裁剪,裁剪方式选内裁,在提示编辑裁剪范围前打“√”,按 [确认 ],保存为“上图数据＊.GRD”文件 (如“上图数据Ms.GRD”文件)。

(3)第三步:在 MapGIS数字地面模型子系统中 ,选择菜单项 [Grd模型 ]→[平面等值线图绘制 ]→打开己裁剪好的“上图数据Ms.GRD”文件,出现如图 11所示 [设置等值线参数 ]对话框。此时应根据统计数据,来确定背景值、异常范围,以及异常级别,并设定各种参数[2]。也可以点击[装设置 ],调用己有的 VEL设置文件 (设置好后请注意保存设置)。在按[确认 ]键后,自动生成异常等值线图。在认为比较满意时,请保存点、线、面文件(如“视极化率.WT、视极化率.WL、视极化率.WP”文件)。注意:[制图幅面 ]要选取原始数据范围。

图11 设置等值线参数Fig.11 Parameter setup for contour lines

图12 7号剖面视极化率拟断面图Fig.12 Apparent polarization pseudo section map of profile No.7

3.2 绘制拟断面图

在上面的叙述中己经生成了“地形.WL、展点.WT”文件,平面等值线的“视极化率.WT、视极化率.WL、视极化率.WP”文件。双击 MapGIS软件的主界面,[图形处理 ]→[输入编辑 ]→[新建工程 ]→[确定 ]→[确定 ]→选 [不生成可编辑项 ]→[确定 ]→在工程管理窗口内右击→单击添加项目→将上述所生成的点、线、面文件打开,并逐一进行编辑。最后,保存全部点、线、面文件,并保存工程,见图 12(θ°=0)。

在绘制完一条测深剖面后,接着重复上述步骤,绘制另外一条测深剖面的拟断面图。如此反复,直到所有测深剖面的拟断面图都绘制完毕。

双击 MapGIS软件的主界面,[图形处理 ]→[输入编辑 ]→[新建工程 ]→[确定 ]→[确定 ]→选[不生成可编辑项 ]→[确定 ]→在工程管理窗口内右击→单击[添加项目 ]→把绘制旋转θ°后的拟断面图全部添加,逐一进行编辑,生成新的工程文件并保存,见下页图 13。

图13未绘制出地形线,但从整体来看立体感比较强。

图13 相邻剖面旋转θ°后的拟断面图Fig.13 Pseudo section map of adjacent cross-section with rotation angleθ°

4 结论

(1)首先先定好切片图的左下角坐标 (X0,Y0)为原点,然后利用旋转角θ°重新计算裁剪区、绘图区域的坐标,最后重新生成目标绘图区域的 GRD文件,利用VEL文件生成旋转的切片图。

(2)组合排列多条相邻剖面的拟断面图,视觉上有较好的立体感。

[1] 王声喜 ,康宝林.Mapsource与 Excel、MapGIS相结合在化探工作中的应用[J].物探化探计算技术,2009,31(2):169.

[2] 傅良魁.电法勘探教程 [M].北京:地质出版社,1983.

[3] 李金铭.地电场与电法勘探[M].北京:地质出版社,2005.

[4] 武汉中地信息工程有限公司.Mapgis地理信息系统实用教程[M].武汉:武汉中地信息工程有限公司,2002.

[5] 中华人民共和国地质矿产行业标准.DZ/T0070-93时间域激发极化法技术规定 [S].北京:中国标准出版社,1993.

[6] 中华人民共和国地质矿产行业标准.DZ/T0072-93电阻率测深法技术规程[S].北京:中国标准出版社,1993.

[7] 中华人民共和国地质矿产行业标准.DZ/T0069地球物理勘查图图式图例及用色标准[S].北京:中国标准出版社,1993.