用 MapGIS繪制擬斷面 (切片)圖

彭亮華

(湖南省有色地質(zhì)勘查局二四七隊,湖南長沙 410129)

0 前言

在礦產(chǎn)地質(zhì)的普查與詳查過程中,需要進(jìn)行一系列的地球物理測量工作,從而會采集到大量的地球物理數(shù)據(jù)信息。如何對采集到的數(shù)據(jù)作相應(yīng)處理,并采用 MapGIS軟件繪制對應(yīng)的地球物理圖件,作者對此有一點(diǎn)經(jīng)驗。注意到王聲喜、康寶林[1]的文章后 (以下簡稱《化探圖件》),考慮物探圖件與化探圖件的編制基本相似,因此在這里只談用MapGIS繪制物探方法中測深裝置的擬斷面圖,供同行們參考。

1 數(shù)學(xué)模型

這里介紹一種逆時針旋轉(zhuǎn)θ°、等面積的變換,其數(shù)學(xué)模型如圖 1所示。ABFG多邊形是變換前的擬斷面圖區(qū)域,ACKH是逆時針旋轉(zhuǎn)θ°后等面積的擬斷面圖區(qū)域。

以 A點(diǎn)為原點(diǎn),AG為旋轉(zhuǎn)軸逆時針旋轉(zhuǎn)θ°后 ,E點(diǎn)變換到 R點(diǎn) ,∠EAR=θ°,A、B、C、R、E各點(diǎn)的坐標(biāo)如圖 1所示。不難知道,Xr=AR＊ cosθ°、Yr=AR ＊ sinθ°、AR=AE=(Xe-X0)。B點(diǎn)變換到 C點(diǎn),其坐標(biāo)變換為:

圖1 等面積變換的數(shù)學(xué)模型Fig.1 Mathematical model of equal Area tranformation

從式 (1)可知,旋轉(zhuǎn)后點(diǎn)的 X軸長只有原來的cosθ°倍。為了使旋轉(zhuǎn)后的面積與旋轉(zhuǎn)前的面積相等,那么 Y軸長需擴(kuò)大 (1/cosθ°)倍 (見下頁圖 2),從式 (2)就反映了這一點(diǎn)。從式 (2)可知,旋轉(zhuǎn)后的點(diǎn)在 Y軸上己向上平移了 (X56-X0)＊sinθ°段距離。下面把式 (1)、式 (2)改寫成通用的公式:

圖2 旋轉(zhuǎn)變換前、后面積相等Fig.2 Equal area of before and after rotation

其中 θ°為旋轉(zhuǎn)角;X′n為旋轉(zhuǎn)后第N個點(diǎn)的X坐標(biāo);Y′n為旋轉(zhuǎn)后第N個點(diǎn)的Y坐標(biāo);Xn為旋轉(zhuǎn)前第 N個點(diǎn)的 X坐標(biāo);Yn為旋轉(zhuǎn)前第 N個點(diǎn)的 Y坐標(biāo);X0為旋轉(zhuǎn)前左下角原點(diǎn)的 X坐標(biāo);Y0為旋轉(zhuǎn)前左下角原點(diǎn)的 Y坐標(biāo)。

式 (3)、式 (4)是逆時針旋轉(zhuǎn)θ°等面積變換擬斷面圖區(qū)域的數(shù)學(xué)公式。其中 X0、Y0可根據(jù)情況作取舍,如需把旋轉(zhuǎn)后的擬斷面圖放在旋轉(zhuǎn)前的擬斷面圖一塊進(jìn)行比較,則要輸入 X0、Y0,可能還要在式 (3)中加一平移常數(shù),然后才能進(jìn)行比較。

2 數(shù)據(jù)整理

2.1 GPS數(shù)據(jù)整理

GPS的參數(shù)設(shè)置及坐標(biāo)數(shù)據(jù) (另存為制表符分隔“＊.txt”的文本文件)與《化探圖件》中所描述的一致,下面作者將對 DX、DY、DZ、DA、DF參數(shù)的計算作補(bǔ)充說明。

2.1.1 坐標(biāo)系統(tǒng)參數(shù)計算

搜集應(yīng)用區(qū)域內(nèi) GPS“B”級網(wǎng)三個以上網(wǎng)點(diǎn)WGS84坐標(biāo)系 B、L、H值,以及我國坐標(biāo)系 (BJ54或西安 80)B、L、h、x值。(注:B、L、H分別為大地坐標(biāo)系中的大地緯度、大地經(jīng)度及大地高,h、x分別為大地坐標(biāo)系中的高程及高程異常。各參數(shù)可以通過各省級測繪局或測繪院具有“A”級、“B”級網(wǎng)的單位獲得。)

計算不同坐標(biāo)系三維直角坐標(biāo)值,計算公式如下:

X=(N+H)＊cosB＊cosL

Y=(N+H)＊cosB＊ sinL

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

不同坐標(biāo)系對應(yīng)橢球的有關(guān)常數(shù)詳見表 1。

表 1 不同坐標(biāo)系統(tǒng)對應(yīng)橢球的有關(guān)數(shù)據(jù)Tab.1 Relativedatatableofellipsoidcorrespondingtodifferentcoordinationsystem

在表 1中,X、Y、Z為大地坐標(biāo)系中的三維直角坐標(biāo):A為大地坐標(biāo)系對應(yīng)橢球之長半軸;E＊E為大地坐標(biāo)系對應(yīng)橢球第一偏聽偏信心率;F為對應(yīng)橢球之扁率;N為該點(diǎn)的卯酉圈曲率半徑,N=A/(1-E＊ E＊ sinB＊ sinB)∧0.5;H=h+x,該處 H為BJ54或西安 80坐標(biāo)系中的大地高。

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

利用WGS84坐標(biāo)系的 X、Y、Z及 A、F值,減去我國坐標(biāo)系的對應(yīng)值,可以得出實(shí)現(xiàn)坐標(biāo)系統(tǒng)轉(zhuǎn)換的五個參數(shù) (DX、DY、DZ、DA、DF),還應(yīng)算出WGS84與北京和西安坐標(biāo)系二套參數(shù)。

2.1.3 參數(shù)驗證

參數(shù)計算之后必須進(jìn)行驗證。驗證的方法是在應(yīng)用區(qū)域內(nèi)選擇五個以上水準(zhǔn)點(diǎn)進(jìn)行實(shí)測,然后進(jìn)行參數(shù)計算之后必須進(jìn)行的驗證。實(shí)測值與測繪部門提供的理論值對比,如果最大誤差不大于15m,平均誤差不大于 10m,則計算出的參數(shù)可以使用;否則,要重新計算或查找出現(xiàn)問題的原因。

例如:某 8號剖面井口的WGS84坐標(biāo)是 B=25.7809°;L=113.16472°;H=222。

BJ54坐標(biāo)是 X=2852880;Y=416160;H=222。

B=25.7809°=25∶46∶51.200

L=113.16472°=113∶09∶53.000

把B值與 L值輸入程序三參數(shù)計算中 (見圖3),右邊輸入 BJ54坐標(biāo),經(jīng)確認(rèn)后得 DX、DY、DZ(見圖 4)。選三至五個水準(zhǔn)點(diǎn)進(jìn)行實(shí)測,驗證參數(shù)。

DX=-14.065,DY= -109.798

DZ= -51.947

圖3 三參數(shù)之計算Fig.3 Calculating the three parameter

圖4 經(jīng)計算后的三參數(shù)值Fig.4 After calculated the three parameter

DA=-108(WGS84之長半軸 -BJ54之長半軸)

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

一般地 DA= -108,DF= +0.00000050以不變值輸入 GPS與Mapsource軟件中,不過你會發(fā)現(xiàn)Mapsource軟件中位置設(shè)置界面之 Delta會發(fā)生變化 (4.99999998002099e-007),見圖 5。

圖5 系統(tǒng)自動修正的 DF值Fig.5 System automatically correction of the DFvalue

2.2 剖面測深點(diǎn)的數(shù)據(jù)整理

2.2.1 測深點(diǎn)物探數(shù)據(jù)的整理

以正東西向測深剖面某測深點(diǎn)為例,一般以GPS觀測數(shù)據(jù)的 Y坐標(biāo)值標(biāo)注為測深圖件的 X軸坐標(biāo)值,列為第一列;以 GPS觀測數(shù)據(jù)的高程坐標(biāo)值 (也可從地形圖上取高程)為起點(diǎn),減去該測深點(diǎn)不同記錄位置的距離,并標(biāo)注為測深圖件的 Y軸坐標(biāo)值,列為第二列;以該測深點(diǎn)的視電阻率、視極化率參數(shù)值分別列為第三列和第四列。排完第一個測深點(diǎn)后,接著往下排第二個、第三個測深點(diǎn),直到排完同一剖面上的測深點(diǎn)為止。在電子表格里生成工作表后,另存為不規(guī)則網(wǎng) (前加 notgrid)逗號分隔符的“展點(diǎn).csv”文件之后,變換為“展點(diǎn).det”文件 ,“展點(diǎn).det”文件里的坐標(biāo)未旋轉(zhuǎn) 。

第一列、第二列為坐標(biāo)數(shù)據(jù),且 GPS的數(shù)據(jù)是以米為單位,因此要注意數(shù)據(jù)單位、比例尺與Map-GIS軟件規(guī)定之間的關(guān)系。

在實(shí)際工作中,對測深剖面的布設(shè),往往取決于地層界線、地質(zhì)構(gòu)造的走向。正東西向、正南北向的剖面,比較容易繪制擬斷面圖,但北東向 (一、三象限)、北西向 (二、四象限)的剖面,比較難繪制擬斷面圖 (指繪制的擬斷面圖方向與實(shí)際剖面方向一致時)。現(xiàn)介紹正東西向的測深剖面,在以西端為軸 (原點(diǎn))逆時針旋轉(zhuǎn)θ°后,其等面積擬斷面圖數(shù)據(jù)的整理。

把列 A、列 B的坐標(biāo)數(shù)據(jù),按式 (3)、式 (4)重新生成旋轉(zhuǎn)θ°角后生成的新的坐標(biāo)數(shù)據(jù);把列 C、列 D數(shù)據(jù)與新的坐標(biāo)數(shù)據(jù),配成 (新建制表符分隔的文本文件 )列 A、列 B、列 C、列 D,依次是 X、Y(軸坐標(biāo)數(shù)據(jù)),視電阻率、視極化率數(shù)據(jù),如圖 6所示。在電子表格里,把每條剖面整理成一個工作表,并另存為制表符分隔“旋轉(zhuǎn)數(shù)據(jù).txt”的文本文件。

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

在電子表格里以 GPS觀測的 Y坐標(biāo)值標(biāo)注為X軸坐標(biāo)值,列為列 A;以 GPS觀測數(shù)據(jù)的高程坐標(biāo)值 (也可從地形圖上取高程)標(biāo)注為 Y軸坐標(biāo)值,列為列 B。列 A、列 B按式 (3)、式 (4)重新生成旋轉(zhuǎn)θ°角后形成新的坐標(biāo)數(shù)據(jù)列 C、列 D,然后把列 C、列 D數(shù)據(jù)另存為制表符分隔“＊.txt”的文本文件,在投影變換中生成線,即地形剖面線“地形.WL”,見下頁圖 7。

在電子表格里,以 GPS觀測的 Y坐標(biāo)值標(biāo)注為 X軸坐標(biāo)值,列為第一列;以 GPS觀測數(shù)據(jù)的高程坐標(biāo)值 (也可從地形圖上取高程)為起點(diǎn),減去該測深點(diǎn)最近、最遠(yuǎn)記錄位置的距離,并標(biāo)注為 Y軸坐標(biāo)值,列為第二列。把第一列、第二列按式(3)、式 (4)重新生成旋轉(zhuǎn)θ°角后形成新的坐標(biāo)數(shù)據(jù)列 C、列 D,把列 C、列 D數(shù)據(jù)另存為制表符分隔“＊.txt”的文本文件。應(yīng)注意,數(shù)據(jù)要按順序排列生成一個大閉合圈,并在投影變換中生成線框。該線框即為裁剪框,注意把線框生成區(qū)文件“裁剪框.WP”。

圖6 旋轉(zhuǎn) 30°的坐標(biāo)變換Fig.6 30°rotation of the coordinatetrans for mation

3 擬斷面圖的繪制

3.1 將 TXT文本、DET文件投影成 MapGIS中的點(diǎn)、線、面文件

3.1.1 生成測深數(shù)據(jù)記錄點(diǎn)點(diǎn)文件

雙擊 MapGIS軟件的主界面,[空間分析 ]→[DTM分析 ],進(jìn)入MapGIS數(shù)字地面模型子系統(tǒng),選擇菜單項 [模型應(yīng)用 (A)]→[高程點(diǎn)標(biāo)注制圖 ],打開“展點(diǎn).det”文件,見下頁圖 8。X軸方向數(shù)據(jù)選列A,Y軸方向數(shù)據(jù)選列 B,L標(biāo)注數(shù)據(jù)選列 C或列D。按所標(biāo)注的內(nèi)容進(jìn)行選擇,比如要標(biāo)注視電阻率參數(shù),那么在本例中選列 C。子圖符號選 247,紅色,角度選0°,符號的尺寸一般選1～3不等。標(biāo)注位置則取決于個人習(xí)慣,標(biāo)注字體一般選宋體、黑色、3mm(如字體大小選 0mm,可繪制實(shí)際材料圖);標(biāo)注格式一般取整,視情況而定。按 [確認(rèn) ]鍵后,即可生成“展點(diǎn).WT”文件。

圖7 生成旋轉(zhuǎn)后的裁剪框數(shù)據(jù)Fig.7 Generation of the cutting box data after rotating

圖8 高程點(diǎn)標(biāo)注選項Fig.8 Marked elevation point option

圖9 離散數(shù)據(jù)網(wǎng)格化選項Fig.9 Discrete data grid options

圖10 測深區(qū)域數(shù)據(jù)裁剪、無效化Fig.10 The cutting sounding data in the regional is invalid

3.1.2 等值線圖的繪制

等值線圖的繪制基本與《化探圖件》[1]相一致,因步驟與地形有關(guān),在此重復(fù)一下。

(1)第一步:先雙擊 MapGIS軟件的主界面,[空間分析 ]→[DTM分析 ],進(jìn)入 MapGIS數(shù)字地面模型子系統(tǒng),選擇菜單項 [Grd模型 ]→[離散數(shù)據(jù)網(wǎng)格化 ],打開“旋轉(zhuǎn)數(shù)據(jù).txt”文件,見圖 9。在數(shù)據(jù)列 D中 X、Y為坐標(biāo)軸數(shù)據(jù)選項,Z為所選繪圖參數(shù)的高程數(shù)據(jù)。按 [確定 ]生成“＊.GRD”文件,如圖 9所示,生成“4M.GRD”文件。

(2)第二步:在 MapGIS數(shù)字地面模型子系統(tǒng)中,單擊“打開高程數(shù)據(jù)文件”快捷圖標(biāo),選擇剛生成的“＊.GRD”文件 (如“4M.GRD”)打開后,選擇菜單項[Grd模型 ]→[規(guī)則網(wǎng)裁剪,無效化 ],見圖10。在裁剪區(qū)域選取方式中,點(diǎn)擊 [選取工作區(qū)中區(qū)域 ],并點(diǎn)擊隨后圖標(biāo),選上述“裁剪框.WP”區(qū)文件打開。操作方式選裁剪,裁剪方式選內(nèi)裁,在提示編輯裁剪范圍前打“√”,按 [確認(rèn) ],保存為“上圖數(shù)據(jù)＊.GRD”文件 (如“上圖數(shù)據(jù)Ms.GRD”文件)。

(3)第三步:在 MapGIS數(shù)字地面模型子系統(tǒng)中 ,選擇菜單項 [Grd模型 ]→[平面等值線圖繪制 ]→打開己裁剪好的“上圖數(shù)據(jù)Ms.GRD”文件,出現(xiàn)如圖 11所示 [設(shè)置等值線參數(shù) ]對話框。此時應(yīng)根據(jù)統(tǒng)計數(shù)據(jù),來確定背景值、異常范圍,以及異常級別,并設(shè)定各種參數(shù)[2]。也可以點(diǎn)擊[裝設(shè)置 ],調(diào)用己有的 VEL設(shè)置文件 (設(shè)置好后請注意保存設(shè)置)。在按[確認(rèn) ]鍵后,自動生成異常等值線圖。在認(rèn)為比較滿意時,請保存點(diǎn)、線、面文件(如“視極化率.WT、視極化率.WL、視極化率.WP”文件)。注意:[制圖幅面 ]要選取原始數(shù)據(jù)范圍。

圖11 設(shè)置等值線參數(shù)Fig.11 Parameter setup for contour lines

圖12 7號剖面視極化率擬斷面圖Fig.12 Apparent polarization pseudo section map of profile No.7

3.2 繪制擬斷面圖

在上面的敘述中己經(jīng)生成了“地形.WL、展點(diǎn).WT”文件,平面等值線的“視極化率.WT、視極化率.WL、視極化率.WP”文件。雙擊 MapGIS軟件的主界面,[圖形處理 ]→[輸入編輯 ]→[新建工程 ]→[確定 ]→[確定 ]→選 [不生成可編輯項 ]→[確定 ]→在工程管理窗口內(nèi)右擊→單擊添加項目→將上述所生成的點(diǎn)、線、面文件打開,并逐一進(jìn)行編輯。最后,保存全部點(diǎn)、線、面文件,并保存工程,見圖 12(θ°=0)。

在繪制完一條測深剖面后,接著重復(fù)上述步驟,繪制另外一條測深剖面的擬斷面圖。如此反復(fù),直到所有測深剖面的擬斷面圖都繪制完畢。

雙擊 MapGIS軟件的主界面,[圖形處理 ]→[輸入編輯 ]→[新建工程 ]→[確定 ]→[確定 ]→選[不生成可編輯項 ]→[確定 ]→在工程管理窗口內(nèi)右擊→單擊[添加項目 ]→把繪制旋轉(zhuǎn)θ°后的擬斷面圖全部添加,逐一進(jìn)行編輯,生成新的工程文件并保存,見下頁圖 13。

圖13未繪制出地形線,但從整體來看立體感比較強(qiáng)。

圖13 相鄰剖面旋轉(zhuǎn)θ°后的擬斷面圖Fig.13 Pseudo section map of adjacent cross-section with rotation angleθ°

4 結(jié)論

(1)首先先定好切片圖的左下角坐標(biāo) (X0,Y0)為原點(diǎn),然后利用旋轉(zhuǎn)角θ°重新計算裁剪區(qū)、繪圖區(qū)域的坐標(biāo),最后重新生成目標(biāo)繪圖區(qū)域的 GRD文件,利用VEL文件生成旋轉(zhuǎn)的切片圖。

(2)組合排列多條相鄰剖面的擬斷面圖,視覺上有較好的立體感。

[1] 王聲喜 ,康寶林.Mapsource與 Excel、MapGIS相結(jié)合在化探工作中的應(yīng)用[J].物探化探計算技術(shù),2009,31(2):169.

[2] 傅良魁.電法勘探教程 [M].北京:地質(zhì)出版社,1983.

[3] 李金銘.地電場與電法勘探[M].北京:地質(zhì)出版社,2005.

[4] 武漢中地信息工程有限公司.Mapgis地理信息系統(tǒng)實(shí)用教程[M].武漢:武漢中地信息工程有限公司,2002.

[5] 中華人民共和國地質(zhì)礦產(chǎn)行業(yè)標(biāo)準(zhǔn).DZ/T0070-93時間域激發(fā)極化法技術(shù)規(guī)定 [S].北京:中國標(biāo)準(zhǔn)出版社,1993.

[6] 中華人民共和國地質(zhì)礦產(chǎn)行業(yè)標(biāo)準(zhǔn).DZ/T0072-93電阻率測深法技術(shù)規(guī)程[S].北京:中國標(biāo)準(zhǔn)出版社,1993.

[7] 中華人民共和國地質(zhì)礦產(chǎn)行業(yè)標(biāo)準(zhǔn).DZ/T0069地球物理勘查圖圖式圖例及用色標(biāo)準(zhǔn)[S].北京:中國標(biāo)準(zhǔn)出版社,1993.