潘黎黎，林龍斌，李細(xì)光

(1.廣東省珠海工程勘察院，廣東 珠海 519002；2.山西省第三地質(zhì)工程勘察院，山西 晉中 030620；3.廣西壯族自治區(qū)地震局，廣西 南寧 530022)

0 引言

近年來(lái)，數(shù)字地質(zhì)填圖技術(shù)在地質(zhì)行業(yè)，特別是區(qū)域地質(zhì)調(diào)查填圖、城市地質(zhì)調(diào)查工作中發(fā)揮著重要作用[1-5]。國(guó)內(nèi)地礦行業(yè)主要使用數(shù)字地質(zhì)調(diào)查系統(tǒng)DGSS[6-9]進(jìn)行填圖等相關(guān)工作。在城市地質(zhì)調(diào)查三維基礎(chǔ)地質(zhì)結(jié)構(gòu)調(diào)查[10-11]野外工作中利用DGSS子系統(tǒng)數(shù)字地質(zhì)填圖系統(tǒng)RGMap[12-14](移動(dòng)版)采集路線PRB(地質(zhì)調(diào)查點(diǎn)、點(diǎn)間路線、地質(zhì)界線)數(shù)據(jù)[15-17]，利用RGMap(PC版)對(duì)PRB野外數(shù)據(jù)進(jìn)行整理、成圖和管理等。移動(dòng)版RGMap除能記錄地質(zhì)調(diào)查點(diǎn)上的描述信息，也能記錄調(diào)查點(diǎn)的位置信息，如經(jīng)緯度、投影平面XY坐標(biāo)、高程等信息。

通常，為整體了解路線上地層、接觸關(guān)系及構(gòu)造特征，建立宏觀認(rèn)識(shí)，需要繪制路線信手地質(zhì)剖面圖。傳統(tǒng)地質(zhì)調(diào)查中，路線信手地質(zhì)剖面圖多是在野外與調(diào)查同步進(jìn)行繪制，而在數(shù)字地質(zhì)填圖工作中往往利用數(shù)字地質(zhì)填圖系統(tǒng)繪制。RGMap系統(tǒng)(PC版)可以利用數(shù)字化地形圖自動(dòng)生成路線信手剖面框架。但通常在城市地質(zhì)調(diào)查工作前期僅收集到大比例尺掃描版地形圖，而沒有數(shù)字化地形圖，這使得RGMap系統(tǒng)無(wú)法自動(dòng)生成路線信手剖面框架圖中的地形線，而需從掃描版地形圖上讀取路線高程數(shù)據(jù)來(lái)繪制信手剖面地形線，極大地增加了工作量。因此，快速獲取路線地形數(shù)據(jù)是繪制路線信手剖面首先需要解決的問題。隨著人工智能、大數(shù)據(jù)等的快速發(fā)展，Python編程語(yǔ)言依靠其語(yǔ)言簡(jiǎn)單、入門快、功能強(qiáng)大、跨平臺(tái)等優(yōu)勢(shì)，越來(lái)越受到大眾的喜愛。本文擬介紹利用Python腳本批量獲取地形線上地質(zhì)點(diǎn)橫縱坐標(biāo)，并結(jié)合Section軟件自動(dòng)繪制路線信手剖面地形線的原理和方法等。

1 原理分析

信手剖面地形線常用展開法和投影法[18]繪制，RGMap系統(tǒng)默認(rèn)的繪制方法為展開法，因此，本文也僅討論展開法繪制的地形線。利用展開法繪制路線信手剖面地形線的一般步驟(圖1)：① 確定信手剖面地形線方向和起止點(diǎn)；② 確定地形線比例尺；③ 確定各點(diǎn)間的水平距離和各點(diǎn)的高程。為完整的反映路線上地質(zhì)概況，將路線首、尾地質(zhì)點(diǎn)作為信手剖面地形線的起止點(diǎn)，方向一般由起點(diǎn)指向終點(diǎn)，若路線整體方向變化較大，可在剖面上標(biāo)注多個(gè)方向。本文以比例尺為1∶1萬(wàn)的三維基礎(chǔ)地質(zhì)結(jié)構(gòu)調(diào)查為例，設(shè)定信手剖面地形線比例尺為1∶1萬(wàn)，繪制時(shí)先按設(shè)定比例尺將第一個(gè)點(diǎn)繪制在方格紙上，計(jì)算第二點(diǎn)與第一點(diǎn)的水平距離和高程，確定第2點(diǎn)的位置，以此類推，在方格紙上按順序先后繪出所有地質(zhì)點(diǎn)，最后將各地質(zhì)點(diǎn)用平滑曲線依次連接起來(lái)，即為路線信手剖面地形線。因此，繪制地形線最關(guān)鍵的是確定相鄰地質(zhì)點(diǎn)間的水平距離和各點(diǎn)高程。移動(dòng)版RGMap系統(tǒng)記錄了地質(zhì)調(diào)查點(diǎn)的平面投影x、y坐標(biāo)和高程值，在PC版RGMap系統(tǒng)中可以通過(guò)GPOINT.WT文件查看地質(zhì)調(diào)查點(diǎn)的位置信息，利用Python腳本提取、處理地質(zhì)點(diǎn)位置信息依次獲取相鄰地質(zhì)點(diǎn)間的水平距離和高程，結(jié)合Section軟件的批量點(diǎn)數(shù)據(jù)導(dǎo)入功能，即可自動(dòng)繪制信手剖面地形線。

2 數(shù)據(jù)準(zhǔn)備

PC版RGMap系統(tǒng)中GPOINT.WT文件記錄了野外地質(zhì)調(diào)查點(diǎn)的所有信息，并按屬性字段分別存儲(chǔ)，主要有ID(系統(tǒng)自增點(diǎn)記錄)、ROUTECODE(自定義路線號(hào))、GEOPOINT(自定義點(diǎn)號(hào))、LONGITUDE(經(jīng)度)、LATITUDE(緯度)、ALTITUDE(高程)、XX(投影平面X坐標(biāo))、YY(投影平面Y坐標(biāo))等(圖2a)等屬性字段。為快速提取地質(zhì)點(diǎn)的ALTITUDE、XX、YY屬性內(nèi)容，計(jì)算各點(diǎn)在地形線上的橫縱坐標(biāo)，繪制信手剖面地形線，先利用RGMap系統(tǒng)中“屬性數(shù)據(jù)導(dǎo)出Excel”的功能將這些屬性內(nèi)容導(dǎo)出為Excel文件。在RGMap系統(tǒng)中野外手圖界面將GPOINT.WT文件設(shè)置為編輯狀態(tài)，點(diǎn)擊菜單欄“工具”—“圖層屬性導(dǎo)出Excel”(圖2b)，在彈出的窗口“處理類型”默認(rèn)選擇點(diǎn)文件(圖2c)，“是否導(dǎo)出PRB描述信息”的窗口選擇“否”，系統(tǒng)自動(dòng)調(diào)用MS-Office等工具打開導(dǎo)出的屬性數(shù)據(jù)(圖2d)，將其保存為Excel文件。

3 數(shù)據(jù)提取

3.1 導(dǎo)入第三方模塊

利用Python腳本從Excel文件中提取地質(zhì)調(diào)查點(diǎn)的高程和平面投影XX、YY坐標(biāo)并計(jì)算相鄰地質(zhì)點(diǎn)之間的水平距離，首先導(dǎo)入Python腳本處理數(shù)據(jù)所需的pandas、openpyxl、math、os、tkinter、subprocess等模塊。其中pandas模塊中的read_excel和ExcelWriter用于向Excel文件中讀取和寫入數(shù)據(jù)，openpyxl模塊用于處理“.xlsx”格式的Excel文件，math模塊中的sqrt用于計(jì)算兩點(diǎn)間的水平距離，os模塊中的path用于獲取Excel文件路徑相關(guān)參數(shù)，tkinter模塊中的filedialog和Tk用于從磁盤選擇要處理的Excel文件并返回文件的絕對(duì)路徑。代碼如下：

from pandas import read_excel，ExcelWriter

import openpyxl

from math import sqrt

from os import path

from tkinter import filedialog，Tk

from subprocess import Popen

結(jié)合tkinter模塊filedialog函數(shù)返回的源Excel文件(Input_file)的絕對(duì)路徑，利用os模塊的path將結(jié)果文件(Output_file)保存于Input_file目錄下，便于查看和管理。用pandas模塊read_excel函數(shù)打開并讀取Input_file的數(shù)據(jù)并保存到內(nèi)存中。代碼如下：

root = Tk()

root.withdraw()

Input_file = filedialog.askopenfilename()

圖2 (a)地質(zhì)點(diǎn)屬性字段 (b)處理文件類型 (c)圖層屬性導(dǎo)出Excel (d)Excel屬性數(shù)據(jù)

(c)Properties of map layer to excel file (d)Excel property data

if input_file：

Output_file = path.dirname(input_file) +'/地形線_'+ path.basename(input_file)

Worksheet_data = read_excel(input_file，'Sheet1'，index_col=None)

3.2 橫坐標(biāo)的提取

定義一個(gè)list列表保存地質(zhì)調(diào)查點(diǎn)在信手剖面上的橫坐標(biāo)，并且設(shè)定第一個(gè)地質(zhì)調(diào)查點(diǎn)的橫坐標(biāo)為0(即路線信手剖面的起點(diǎn))。依次計(jì)算兩相鄰兩地質(zhì)調(diào)查點(diǎn)之間的水平距離，假定S1為第2點(diǎn)與起點(diǎn)之間的水平距離，S2為第3點(diǎn)與第2點(diǎn)之間的水平距離，以此類推，Sn為第n+1點(diǎn)與第n點(diǎn)之間的水平距離，那么第n+1點(diǎn)的橫坐標(biāo)應(yīng)等于S1，S2，S3，...，Sn之和，依此原理計(jì)算各地質(zhì)點(diǎn)的在信手剖面上的橫坐標(biāo)值。式(1)為平面投影直角坐標(biāo)系下兩點(diǎn)之間的水平距離計(jì)算公式，式(2)為地質(zhì)點(diǎn)在信手剖面上的橫坐標(biāo)值計(jì)算公式。

(1)

(2)

式(1)中，Sn為第n+1地質(zhì)點(diǎn)與第n地質(zhì)點(diǎn)之間的水平距離，Xn和Yn分別為第n個(gè)地質(zhì)點(diǎn)的平面投影坐標(biāo)XX和YY；式(2)中Tn為第n地質(zhì)點(diǎn)在信手剖面上的橫坐標(biāo)，默認(rèn)T1為0。Python中可利用math模塊內(nèi)置的sqrt計(jì)算平面直角坐標(biāo)系下兩點(diǎn)間的水平距離，并進(jìn)行循環(huán)累加求和計(jì)算各地質(zhì)點(diǎn)的橫坐標(biāo)。野外RGMap系統(tǒng)采集的XX、YY坐標(biāo)及高程ALTITUDE單位均為“m”，而路線信手剖面圖的比例尺設(shè)定為1∶1萬(wàn)，RGMap系統(tǒng)中最小單位為“mm”，所以需要將以m為單位的Tn數(shù)據(jù)轉(zhuǎn)換成1∶1萬(wàn)比例尺下以mm為單位的數(shù)據(jù)，即將Tn值除以10即得以mm為單位的橫坐標(biāo)x。代碼如下：

#計(jì)算地質(zhì)點(diǎn)的橫坐標(biāo)

List_x = [0，]

sum = 0

for i in range(1，len(Worksheet_data['XX']))：

sum+=sqrt((Worksheet_data['XX'][i]-Worksheet_data['XX'][i-1])**2+(Worksheet_data['YY'][i]-Worksheet_data['YY'][i-1])**2)/10

List_x.append(sum)

Worksheet_data['y'] = List_x

3.3 縱坐標(biāo)的提取

通常信手剖面圖橫縱比例尺保持一致，路線信手剖面上地質(zhì)點(diǎn)的縱坐標(biāo)為地質(zhì)點(diǎn)的高程，但需對(duì)RGMap系統(tǒng)采集的高程數(shù)據(jù)ALTITUDE進(jìn)行變換，與XX、YY坐標(biāo)變換方法一樣，將ALTITUDE數(shù)據(jù)除以10即為1∶1萬(wàn)比例尺下以mm為單位的數(shù)據(jù)。一幅完整的路線信手剖面圖需將比例尺及圖例擺放于圖下方，為便于制圖，將地形線整體向上平移50 mm，所以需在轉(zhuǎn)換后的高程數(shù)據(jù)(縱坐標(biāo))上加50mm得到最終縱坐標(biāo)值y。代碼如下：

#計(jì)算地質(zhì)點(diǎn)的眾坐標(biāo)

List_y = []

for i in range(len(Worksheet_data['ALTITUDE']))：

item = Worksheet_data['ALTITUDE'][i]/10 + 50

List_y.append(item)

Worksheet_data['x'] = List_y

3.4 點(diǎn)編號(hào)及數(shù)據(jù)保存

路線信手剖面上要求標(biāo)注地質(zhì)點(diǎn)號(hào)，所以需從屬性字段GEOPOINT(點(diǎn)號(hào))提取地質(zhì)點(diǎn)號(hào)GEOPOINT(圖2d)，將計(jì)算的橫縱坐標(biāo)x、y與地質(zhì)點(diǎn)號(hào)ID一并寫入Output_file，利用subprocess模塊內(nèi)置Popen自動(dòng)打開Output_file，代碼見圖3a。Output_file包含三列數(shù)據(jù)，即x、y和ID(圖3b)，分別代表地質(zhì)點(diǎn)的橫、縱坐標(biāo)和地質(zhì)點(diǎn)號(hào)。

圖3 (a)代碼段 (b)地質(zhì)點(diǎn)的橫縱坐標(biāo)和編號(hào)

4 地形線繪制

Section V4.6是在MapGIS6.7基礎(chǔ)上，利用Microsoft VC++6.0語(yǔ)言編寫的地質(zhì)圖件制作軟件?；贛apGIS 6.7輸入編輯子系統(tǒng)強(qiáng)大的點(diǎn)、線、區(qū)編輯能力，添加專業(yè)的地質(zhì)圖件制作工具，大大提高了地質(zhì)圖件的制作效率。

4.1 繪制地形線

首先創(chuàng)建路線信手剖面圖點(diǎn)、線、區(qū)文件。在PC版RGMap系統(tǒng)野外手圖界面，選擇菜單欄上“地質(zhì)填圖數(shù)據(jù)操作”—“信手剖面生成與瀏覽”—“顯示”(圖4a)，即可自動(dòng)在相應(yīng)野外路線的“素描圖”文件夾里生成信手剖面圖點(diǎn)、線、區(qū)和對(duì)應(yīng)的工程文件。打開Section軟件，找到“素描圖”文件夾，加載自動(dòng)生成的信手剖面工程文件，即可對(duì)信手剖面圖點(diǎn)、線、區(qū)文件進(jìn)行編輯。打開前述Output_file文件，點(diǎn)擊菜單欄上“1輔助工具”—“表格數(shù)據(jù)投影”—“全部數(shù)據(jù)投影(EXCEL)”(圖4b)。彈出“數(shù)據(jù)投影”窗口，去掉“線閉合”前面的勾(圖4c)，分別設(shè)置“文字圖元參數(shù)”、“子圖圖元參數(shù)”(圖4d、4e)和“線圖元參數(shù)”，為使信手剖面地形線圓滑，可以將線型設(shè)置為曲線。文字、子圖、線圖元等參數(shù)設(shè)置完成后，點(diǎn)擊“確認(rèn)”，即可自動(dòng)生成地形線，并在地形線上標(biāo)注了地質(zhì)點(diǎn)號(hào)(圖5a)。

4.2 調(diào)整點(diǎn)標(biāo)注

自動(dòng)標(biāo)注的文字和地質(zhì)點(diǎn)圖元擺放并不美觀，為便于整飾圖面，在設(shè)置文字和子圖參數(shù)時(shí)，需將文字和子圖的圖層號(hào)設(shè)置為不同的數(shù)字(圖4d、4e)，地形線生成后，點(diǎn)擊菜單欄上“A圖層”—“改層開關(guān)”—“改點(diǎn)P”將地質(zhì)點(diǎn)圖層關(guān)閉(圖4f)，移動(dòng)地質(zhì)點(diǎn)編號(hào)位置，使圖面更美觀(圖5b)，至此地形線繪制完畢，補(bǔ)充完善信手剖面圖地層、巖體界線和花紋、構(gòu)造、圖例、比例尺、方位等要素即可繪制完成信手剖面圖。

5 結(jié)語(yǔ)

路線信手剖面圖作為整體了解路線巖性、構(gòu)造特征的重要圖件，由于缺失數(shù)字化地形圖，無(wú)法在數(shù)字地質(zhì)填圖系統(tǒng)中自動(dòng)生成信手剖面圖框架。經(jīng)研究，利用路線上地質(zhì)點(diǎn)的平面投影坐標(biāo)計(jì)算相鄰地質(zhì)點(diǎn)間的水平距離，結(jié)合地質(zhì)調(diào)查點(diǎn)的高程，可繪制出大致的路線信手剖面地形線。

圖4 (a)自動(dòng)生成信手剖面 (b)表格數(shù)據(jù)投影 (c)數(shù)據(jù)投影

圖5 (a)自動(dòng)生成的地形線和標(biāo)注 (b)調(diào)整后的地形線和標(biāo)注Fig.5 (a)Automatically created topographic line and annotation (b)Adjusted topographic line and annotation

1)利用Python腳本批量處理RGMap系統(tǒng)中地質(zhì)點(diǎn)屬性數(shù)據(jù)，提取地質(zhì)點(diǎn)的平面投影x、y坐標(biāo)和高程數(shù)據(jù)，采用平面坐標(biāo)系中兩點(diǎn)間的距離公式計(jì)算相鄰地質(zhì)點(diǎn)間的水平距離，并進(jìn)行累加求和，即可計(jì)算出各地質(zhì)點(diǎn)在設(shè)定剖面比例尺下的橫坐標(biāo)，而地質(zhì)點(diǎn)的縱坐標(biāo)為地質(zhì)點(diǎn)高程屬性數(shù)據(jù)在設(shè)定比例尺下的值。

2)利用Section軟件的表格數(shù)據(jù)投影功能，將Python腳本計(jì)算所得的各地質(zhì)點(diǎn)橫、縱坐標(biāo)和地質(zhì)點(diǎn)號(hào)文件導(dǎo)入Section軟件，同時(shí)設(shè)置地質(zhì)點(diǎn)標(biāo)注、地質(zhì)點(diǎn)子圖、線型等參數(shù)，即可自動(dòng)繪制出路線信手剖面地形線。

3)利用Python提取、處理數(shù)據(jù)和Section批量導(dǎo)入數(shù)據(jù)繪制地形線能一定程度上提高工作效率，減少由于人為因素帶來(lái)的誤差，提高了數(shù)據(jù)的準(zhǔn)確度。

本文主要利用野外采集的數(shù)據(jù)進(jìn)行分析計(jì)算，獲得地質(zhì)點(diǎn)的橫縱坐標(biāo)。因此，通過(guò)計(jì)算獲取的橫縱坐標(biāo)數(shù)據(jù)精度主要依賴于野外采集的x、y坐標(biāo)及高程數(shù)據(jù)的精度。目前常用的野外數(shù)字填圖設(shè)備的GPS定位系統(tǒng)獲取的x、y坐標(biāo)的精度較高，基本可以滿足繪圖的需要，而通常獲取到的高程數(shù)據(jù)精度較低，這可能會(huì)影響到地形線的垂直精度，本文存在的不足之處在于直接使用了GPS采集的高程數(shù)據(jù)。因此，如需更準(zhǔn)確的高程數(shù)據(jù)，可以在RGMap系統(tǒng)中根據(jù)電子掃描地形圖對(duì)野外采集的地質(zhì)點(diǎn)高程數(shù)據(jù)進(jìn)行修正后，再利用python腳本進(jìn)行處理。Python語(yǔ)言依據(jù)其獨(dú)特的優(yōu)勢(shì)，不僅能批量處理地質(zhì)點(diǎn)數(shù)據(jù)，還可以應(yīng)用于城市地質(zhì)調(diào)查工作中的鉆孔資料整理分析、數(shù)據(jù)庫(kù)建設(shè)數(shù)據(jù)準(zhǔn)備等各階段，能極大的提高批量數(shù)據(jù)分析處理的工作效率，相信在以后的地質(zhì)工作中應(yīng)用將會(huì)越來(lái)越廣泛。