1∶50 000水系沉積物地球化學(xué)勘查采樣點(diǎn)位自動(dòng)設(shè)計(jì)的方法研究

廖國(guó)忠， 張 偉， 吳文賢， 梁生賢

(中國(guó)地質(zhì)調(diào)查局 成都地質(zhì)調(diào)查中心，成都 610081)

0 引言

地球化學(xué)勘查是一種科學(xué)的找礦方法，在找礦實(shí)踐中已經(jīng)取得了許多重大成果。在區(qū)域地球化學(xué)掃面發(fā)現(xiàn)的成礦預(yù)測(cè)區(qū)內(nèi)，基于1∶50 000水系沉積物測(cè)量圈定找礦靶區(qū),已成為地質(zhì)礦產(chǎn)普查的常用手段，且在我國(guó)大部分地區(qū)取得了較好的應(yīng)用效果，尤其在金礦勘查方面已成為首選手段[1-3]。自該方法引入國(guó)內(nèi)以來，在眾多找礦方法中，一直處于舉足輕重的地位，得到廣泛地應(yīng)用。

水系沉積物測(cè)量是以水系沉積物為采樣對(duì)象所進(jìn)行的地球化學(xué)勘查工作,采樣點(diǎn)位的設(shè)計(jì)是否合理，是否能最大限度控制工作區(qū)匯水面積，直接影響著水系沉積物地球化學(xué)勘查的效果。

當(dāng)前，水系沉積物地球化學(xué)勘查采樣點(diǎn)位的設(shè)計(jì)，主要是根據(jù)地形圖上的水系和等高線信息，手動(dòng)勾繪水系分布圖，然后再根據(jù)水系分布圖手工確定采樣點(diǎn)位。長(zhǎng)期以來都是純手工作業(yè)的方式，工作強(qiáng)度大，對(duì)工作人員工作經(jīng)驗(yàn)要求較高，即使有相當(dāng)工作經(jīng)驗(yàn)的不同工作人員所設(shè)計(jì)的采樣點(diǎn)位也有較大差別。讓工作人員從繁重的純手工設(shè)計(jì)采樣點(diǎn)位的工作中解放出來,是筆者力圖解決的問題。

前人針對(duì)這一純手工作業(yè)方式進(jìn)行了些許改進(jìn)：王海鵬等[4]提出的“MAPGIS條件下完美實(shí)現(xiàn)不規(guī)則網(wǎng)化探采樣點(diǎn)自動(dòng)編號(hào)”； 陳可忠等[5]提出的“化探采樣點(diǎn)自動(dòng)編號(hào)及航點(diǎn)自動(dòng)生成的方法”；馮偉華等[6]提出的“基于VBA實(shí)現(xiàn)設(shè)計(jì)點(diǎn)位編號(hào)和坐標(biāo)提取”。通過研究這些文獻(xiàn)發(fā)現(xiàn)，目前研究?jī)H局限于對(duì)采樣點(diǎn)點(diǎn)位進(jìn)行自動(dòng)編號(hào)，盡管此類研究對(duì)提高工作效率有了一定的改觀，但是對(duì)于采樣點(diǎn)位如何自動(dòng)合理的布置這個(gè)關(guān)鍵問題完全沒有涉足。

筆者提出了基于GIS自動(dòng)提取水系的1∶50 000 水系沉積物采樣點(diǎn)位自動(dòng)布置的工作方法，極大地提高了工作效率，具有重要的實(shí)用價(jià)值。

1 研究思路

通過分析研究水系沉積物采樣點(diǎn)定位的工作流程可知，最基本的步驟是通過已知水系和地形地貌的特征，勾繪出水系分布圖。因此水系地提取是該項(xiàng)工作最基本的步驟。

水系自動(dòng)提取是地理信息系統(tǒng)GIS和水文科學(xué)相結(jié)合的交叉研究領(lǐng)域，其中水文特征提取是一個(gè)較為常見的功能，已經(jīng)在生產(chǎn)中得到了廣泛地應(yīng)用[7-9]。因此要實(shí)現(xiàn)水系沉積物地球化學(xué)勘查中最基礎(chǔ)的水系自動(dòng)提取，引入水文特征提取這一功能便可解決。

水系提取成功后，繼續(xù)提取水系中的特征點(diǎn)(水系入口，出口，中間點(diǎn)等)作為備選采樣點(diǎn)。根據(jù)《地球化學(xué)普查規(guī)范》，1∶50 000水系沉積物地球化學(xué)勘查工作中，采樣點(diǎn)布設(shè)密度為4個(gè)點(diǎn)/km2～8個(gè)點(diǎn)/km2；而通過提取所獲得的備選采樣點(diǎn)的密度有可能遠(yuǎn)遠(yuǎn)大于4個(gè)點(diǎn)/km2～8個(gè)點(diǎn)/km2。因此，必須根據(jù)《地球化學(xué)普查規(guī)范》建立一個(gè)評(píng)價(jià)體系，選取最佳的備選采樣點(diǎn)作為最終的采樣設(shè)計(jì)點(diǎn)位。將篩選的采樣設(shè)計(jì)點(diǎn)位根據(jù)《地球化學(xué)普查規(guī)范》要求自動(dòng)編號(hào)。

2 采樣點(diǎn)位自動(dòng)設(shè)計(jì)方法的實(shí)現(xiàn)

采樣點(diǎn)位自動(dòng)設(shè)計(jì)方法的實(shí)現(xiàn)需要以下幾個(gè)步驟:

2.1 已知水系限制下的水系提取

以獲取的1∶50 000數(shù)字地形圖為基礎(chǔ)，開展后續(xù)各項(xiàng)研究工作。圖1展示了水系沉積物地球化學(xué)勘查采樣點(diǎn)自動(dòng)設(shè)計(jì)工作流程圖。

圖1 采樣點(diǎn)位自動(dòng)設(shè)計(jì)實(shí)現(xiàn)流程圖Fig.1 Working flow chart of automatic design of sampling location

2.1.1 建立DEM

數(shù)字高程模型(Digital Elevation Model, DEM)是水系流域提取的基礎(chǔ)，它通過將高程數(shù)據(jù)柵格化，可以很方便地計(jì)算柵格之間高程數(shù)據(jù)的關(guān)系，從而提取DEM中包含的地形、地貌、水文信息。

筆者采用ArcGIS 10.2軟件，首先利用工作區(qū)內(nèi)的高程點(diǎn)和等高線shape線文件數(shù)據(jù)創(chuàng)建TIN文件，再利用ArcToolbox工具箱內(nèi)的Conversion Tools轉(zhuǎn)換工具，將TIN數(shù)據(jù)轉(zhuǎn)換為柵格，本文所采用的柵格為10 m×10 m。

2.1.2 加入已知水系

傳統(tǒng)的水系自動(dòng)提取方法中，直接通過DEM進(jìn)行提取，發(fā)現(xiàn)自動(dòng)提取的水系與實(shí)際水系的分布有較大出入，前人的研究也發(fā)現(xiàn)這一問題，并作出了相應(yīng)的研究[10-12]，取得了一定的成果。

為了避免上述問題發(fā)生，我們提出了利用已知水系限制DEM的方法，很好地解決了這一問題。利用建立數(shù)字高程模型時(shí)的參數(shù)，對(duì)已知水系進(jìn)行柵格化，然后通過柵格計(jì)算器，將沒有水系的柵格設(shè)置為“0”，有水系的柵格設(shè)置為-100，然后將水系柵格與地形DEM進(jìn)行柵格相加，從而得到已知水系限制的DEM。

2.1.3 水系自動(dòng)提取

1) 水流方向計(jì)算。對(duì)比中心柵格與相鄰柵格之間的高程，根據(jù)水往低處流的常理，假定水流向最大坡降的柵格流動(dòng)，利用D8算法(D8算法的核心思想是：假定雨水降落在地形中某一個(gè)格子上，該格子的水流將會(huì)流向其周圍8個(gè)格子中地形最低的格子中)，計(jì)算出DEM中所有柵格的流向。ArcGIS 10.2軟件中，應(yīng)用水文分析庫下的流向確定命令,生成8方向水流流向，流向原理可參考文獻(xiàn)[11]。

2)匯檢測(cè)和洼地填補(bǔ)。流向確定后，所有的水流都向低處流動(dòng)，然而工作區(qū)內(nèi)的局部低洼點(diǎn)，周邊柵格高程均高于局部低洼點(diǎn)，僅根據(jù)流向計(jì)算，所有的水系將終止于此柵格點(diǎn)。尋找低洼點(diǎn)的過程稱為匯檢測(cè)。匯的出現(xiàn)與實(shí)際現(xiàn)象相違背，因?yàn)榈屯蔹c(diǎn)在積滿水后會(huì)向缺口的地方流動(dòng)。因此找出洼地，然后填充洼地可以解決這一問題。

3) 柵格流量。填充洼地后的DEM，再計(jì)算各柵格的流向。假設(shè)每一個(gè)柵格因?yàn)榇髿饨邓a(chǎn)生流量為“1”的水量，按柵格流向向下游流動(dòng)，每向下流動(dòng)一個(gè)柵格，水量將增加“1”，而每一個(gè)柵格的流量便為所有流向此柵格的柵格個(gè)數(shù)。據(jù)此計(jì)算出每一個(gè)柵格點(diǎn)累積流量。

4) 設(shè)置水系流量閾值。實(shí)際上只有部分柵格會(huì)形成水系，而有的僅形成季節(jié)性沖溝。判定是否形成水系或者沖溝，最關(guān)鍵的問題是設(shè)置閾值。根據(jù)區(qū)內(nèi)已知的水系資料，進(jìn)行多次的試驗(yàn)可以給出與實(shí)際情況相符的閾值。為了區(qū)分不同流量的水系，還可以通過分類設(shè)置閾值來達(dá)到水系分類的目的。

5)水系提取。以設(shè)定的閾值對(duì)累積流量進(jìn)行提取,大于閾值的柵格被定義為水系。

表1 水系特征點(diǎn)的分類

2.2 水系分類以及特征點(diǎn)提取

根據(jù)《地球化學(xué)普查規(guī)范》的要求[13]，每個(gè)采樣點(diǎn)最好能控制0.25 km2面積的匯水域，由于本方法柵格化時(shí)采用的柵格為10 m×10 m，因此0.25 km2面積的匯水流域?yàn)? 500個(gè)柵格產(chǎn)生的流量。

為了使水系的出口有意義，采用了水系分類提取的辦法。流量大于2 500為第一類；流量1 250至2 500為第二類；流量625至1 250為第三類；流量312至625為第四類(表1)。

顯然此時(shí)第二類，流量為1 250至2 500，水系出口流量正好為2 500，此點(diǎn)為最符合《地球化學(xué)普查規(guī)范》的要求的水系沉積物采樣點(diǎn)位。但是由于流量為2 500的點(diǎn)分布不均勻，對(duì)于山脊以及已知的大江大河內(nèi)，此類點(diǎn)很難存在。因此為了符合《地球化學(xué)普查規(guī)范》要求中點(diǎn)位分布均勻的特點(diǎn)，必須在其他流量非2 500的點(diǎn)位上設(shè)計(jì)點(diǎn)位。

筆者提出的解決方法是：將各水系的出口點(diǎn)和長(zhǎng)度大于500 m的第一類水系的中點(diǎn)提取出來，然后將這些提取出來的點(diǎn)按其所處的水系類別按表的規(guī)則為這些點(diǎn)賦予順序號(hào)，以作為水系深積物設(shè)計(jì)點(diǎn)位篩選的關(guān)鍵字。

2.3 水系特征點(diǎn)內(nèi)篩選設(shè)計(jì)點(diǎn)位

根據(jù)《地球化學(xué)普查規(guī)范》的要求[13]，累積流量為2 500的采樣點(diǎn)位，控制的匯水面積正好是0.25 km2。但是此類點(diǎn)較少，僅篩選此類點(diǎn)很難滿足采樣密度要求。因此必須選取其他特征點(diǎn)作為設(shè)計(jì)采樣點(diǎn)。

在同一個(gè)采樣小格內(nèi)，如果存在不同類型的特征點(diǎn)，則采取如下的處理方式進(jìn)行篩選：①以順序號(hào)為第一關(guān)鍵字升序排序，累積流量作為第二關(guān)鍵字降序排序的方式進(jìn)行排序；②選取同一小格內(nèi)的前兩個(gè)特征點(diǎn)作為最終的采樣設(shè)計(jì)點(diǎn)位。

2.4 設(shè)計(jì)點(diǎn)位自動(dòng)編號(hào)

在1∶50 000水系沉積物地球化學(xué)勘查中，樣品編號(hào)分為三步：

1)將圖幅范圍內(nèi)不完整的方里格網(wǎng)坐標(biāo)向圖幅外延伸，取最鄰近的整公里方里網(wǎng)為界，然后以1 km為間距縱向和橫向上劃分網(wǎng)格，并按從左到右，從上到下的順序按數(shù)字編號(hào)。

2)將①劃分好的網(wǎng)格以縱橫方向500 m的間距將網(wǎng)格劃分為四個(gè)小網(wǎng)格，左上角小網(wǎng)格編號(hào)為A，右上角小網(wǎng)格編號(hào)為B，左下角小網(wǎng)格編號(hào)為C，右下角小網(wǎng)格編號(hào)為D。

3)500 m×500 m的網(wǎng)格內(nèi)樣品再按從左到右，從上到下的順序按從“1”開始的自然數(shù)編號(hào)，確定小格內(nèi)樣品順序編號(hào)。

采樣設(shè)計(jì)點(diǎn)最終編號(hào)為大格編號(hào)+小格編號(hào)+小格內(nèi)樣品順序編號(hào)。自動(dòng)編號(hào)的計(jì)算相對(duì)簡(jiǎn)單，只需要簡(jiǎn)單的數(shù)據(jù)換算和排序便可以完成。

假定圖幅的四個(gè)邊界坐標(biāo)分別為left、right、top、bottom，待編號(hào)的樣品設(shè)計(jì)點(diǎn)坐標(biāo)為x、y，那么該點(diǎn)大格編號(hào)和小格編號(hào)通過運(yùn)算可以得到。首先，通過excel函數(shù)求取四個(gè)邊界的整公里坐標(biāo)。

大格編號(hào)：(INT(top / 1000) - INT(y / 1000)) * (INT(right / 1000) - INT( left / 1000)) + INT(x/1000) - INT(left/1000)

小格編號(hào)：首先計(jì)算10 * INT( MOD( x , 1000) / 500)+ INT( MOD( y ,1000)/500)，根據(jù)計(jì)算結(jié)果選定小格編號(hào)，如果結(jié)果為“0”，則小格編號(hào)為C；如果結(jié)果為“1”，則小格編號(hào)為A；如果結(jié)果為“10”，則小格編號(hào)為D；如果結(jié)果為“11”，則小格編號(hào)為B。

小格內(nèi)樣品順序編號(hào)主要是根據(jù)小格內(nèi)樣品從左到右，從上到下的順序進(jìn)行編號(hào)。首先，將計(jì)算得到的大格編號(hào)和小格編號(hào)按字符串連接起來，然后將其編號(hào)與其對(duì)應(yīng)的樣品坐標(biāo)x、y組成一行，然后按excel的排序功能，將所有的設(shè)計(jì)點(diǎn)位進(jìn)行排序，以大格編號(hào)和小格編號(hào)連接起來的字符串列為主要關(guān)鍵字，x為第一次要關(guān)鍵字，y為第二次要關(guān)鍵字。排序好之后的數(shù)據(jù)，通過判斷當(dāng)前行所對(duì)應(yīng)的大格編號(hào)和小格編號(hào)連接起來的字符串與上一行的字符串是否相同，如果不同，則小格內(nèi)樣品順序編號(hào)設(shè)置為“1”，如果相同，小格內(nèi)樣品順序編號(hào)則為上一行小格內(nèi)樣品順序編號(hào)加上“1”。

最終將大格編號(hào)，小格編號(hào)和小格內(nèi)樣品順序編號(hào)按字符串連接起來便得到了樣品編號(hào)。

3 試驗(yàn)區(qū)的應(yīng)用及效果分析

3.1 試驗(yàn)區(qū)概況

此次研究選取了貴州省黔西南1∶50 000馬嶺幅作為試驗(yàn)區(qū)。

馬嶺幅位處南盤江-右江前陸盆地大地構(gòu)造(相)之興義泥凼-貞豐者相臺(tái)緣-斜坡亞相。由師宗-彌勒斷裂帶和紫云-六盤水?dāng)嗔褞г谫F州省內(nèi)圍限區(qū)域，屬滇黔桂“金三角”的組成部分，是由揚(yáng)子被動(dòng)邊緣碳酸鹽臺(tái)地演化而成的一個(gè)中晚三疊紀(jì)周緣前陸盆地。測(cè)區(qū)地層為上古生界至中生界，其中以三疊系中上統(tǒng)的陸源碎屑復(fù)理石最引人矚目。測(cè)區(qū)構(gòu)造變形組合樣式具顯著“條”、“塊”鑲嵌特色，即北東向，北西向強(qiáng)變形帶(“條”—緊閉形褶皺及逆沖斷裂)與弱變形塊(寬緩短軸型褶皺區(qū))鑲嵌。

馬嶺幅內(nèi)目前尚未發(fā)現(xiàn)金礦或金礦點(diǎn)，但是圖幅外分布著許多金礦或者金礦點(diǎn)，如圖幅北側(cè)有泥堡金礦，樓下金礦，東側(cè)有戈塘金礦，南側(cè)有雄武金礦等。在區(qū)域地球化學(xué)金異常圖中，馬嶺幅內(nèi)不僅有幅值較小的金異常，而且區(qū)內(nèi)北東向，北西向和近東西向斷裂發(fā)育，區(qū)內(nèi)出露二疊系龍?zhí)督M地層。因此，在馬嶺幅開展1∶50 000水系沉積物地球化學(xué)勘查，對(duì)尋找以陸源硅質(zhì)碎屑巖為主要容礦巖石的微細(xì)浸染型金礦床有一定的潛力。

圖2 自動(dòng)提取水系結(jié)果對(duì)比圖Fig.2 A contrast diagram of water system extraction

3.2 采樣點(diǎn)位自動(dòng)布置效果

3.2.1 不同方法的水系提取效果對(duì)比

筆者在對(duì)水系自動(dòng)提取的處理上，相較于傳統(tǒng)的自然水系提取，采用了已知水系限制后再提取水系的方案，取得的結(jié)果與實(shí)際水系重合度更高，更好地保證了后期采樣點(diǎn)位置的精度。

沒有已知水系限制下的常規(guī)傳統(tǒng)方法提取的水系，在地勢(shì)起伏較大的地段,與已知水系分布較為吻合；而在寬闊山谷中，河道位置偏差較大，尤其在已知河道轉(zhuǎn)彎處，吻合度很低 ，提取的水系不能用作采樣點(diǎn)位布置的參考。

應(yīng)用本文提出的已知水系限制條件下自動(dòng)提取的水系，不僅在地勢(shì)起伏較大的地段吻合，而且在寬闊山谷中，水系彎道段都能有很好地吻合，大大提高了自動(dòng)提取水系的精度，使得后期點(diǎn)位設(shè)計(jì)的位置更精準(zhǔn)。

圖3 采樣點(diǎn)位自動(dòng)布置結(jié)果圖Fig.3 Distribution map of sampling point

3.2.2 采樣點(diǎn)自動(dòng)布置結(jié)果

圖3為使用本文提出的方法自動(dòng)布置的采樣點(diǎn)位(因?yàn)橐粋€(gè)圖幅較大，為便于展示，僅提取了一個(gè)流域)，采樣點(diǎn)位布置均勻，采樣點(diǎn)均位于水系或者是季節(jié)性水系上，很好地控制了采樣點(diǎn)上游的匯水域，采樣點(diǎn)編號(hào)正確，根據(jù)水系分布密度不同，采樣小格的采樣點(diǎn)最多為2個(gè)，空白小格數(shù)較少，可以滿足水系沉積物測(cè)量采樣點(diǎn)位設(shè)計(jì)的基本要求。但是，圖3中尚有不少設(shè)計(jì)點(diǎn)位的合理性不夠，如3級(jí)水系中設(shè)計(jì)點(diǎn)位太多等，需要今后開展進(jìn)一步的研究工作，完善本方法技術(shù)。

4 結(jié)論

筆者結(jié)合GIS和水文科學(xué)相關(guān)知識(shí)，為水系沉積物地球化學(xué)勘查采樣點(diǎn)位布置提供了一種自動(dòng)、快捷的新方法。以1∶50 000電子版地形圖為基礎(chǔ)數(shù)據(jù)，利用ArcGIS 10.2軟件平臺(tái)，通過已知水系的限制，自動(dòng)提取出與實(shí)際情況相吻合的水系。根據(jù)流量將水系進(jìn)行分類，并提取水系的出水口和流量大于2 500 m且長(zhǎng)度大于500 m水系的中點(diǎn)。篩選出符合《地球化學(xué)普查規(guī)范》中采樣點(diǎn)要求的特征點(diǎn)作為最終采樣點(diǎn)。通過采樣點(diǎn)坐標(biāo)換算，自動(dòng)進(jìn)行編號(hào)，實(shí)現(xiàn)了水系沉積物勘查采樣點(diǎn)位自動(dòng)布置。

新方法是多學(xué)科多領(lǐng)域交叉融合的結(jié)果，極大地提高了工作效率，大大減輕了從業(yè)人員勞動(dòng)量。同時(shí)，克服了不同經(jīng)驗(yàn)工作人員布置點(diǎn)位不同的缺點(diǎn)。

由于水系自動(dòng)提取方法自身的缺陷，在喀斯特地區(qū)和平頂山上，尤其是巖溶漏斗和“斷頭河”較多的地區(qū)，該方法誤差較大。同時(shí)，由于選擇條件設(shè)置不夠全面，篩選出來的設(shè)計(jì)采樣點(diǎn)位還有不少不符合《地球化學(xué)普查規(guī)范》的要求，需要人工干預(yù)調(diào)整。

這里僅提出了一種思路，且是初步嘗試，具有了一個(gè)很好的開端，以后還需要根據(jù)野外實(shí)際情況，提取河道更多的特征點(diǎn)。同時(shí)，按照《地球化學(xué)普查規(guī)范》中的布點(diǎn)原則，設(shè)置更多的選擇條件，最終達(dá)到篩選出完全符合實(shí)際要求的水系沉積物測(cè)量設(shè)計(jì)采樣點(diǎn)。