廖國忠， 張 偉， 吳文賢， 梁生賢

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

0 引言

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

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

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

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

筆者提出了基于GIS自動提取水系的1∶50 000 水系沉積物采樣點位自動布置的工作方法，極大地提高了工作效率，具有重要的實用價值。

1 研究思路

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

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

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

2 采樣點位自動設(shè)計方法的實現(xiàn)

采樣點位自動設(shè)計方法的實現(xiàn)需要以下幾個步驟:

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

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

圖1 采樣點位自動設(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ù)柵格化，可以很方便地計算柵格之間高程數(shù)據(jù)的關(guān)系，從而提取DEM中包含的地形、地貌、水文信息。

筆者采用ArcGIS 10.2軟件，首先利用工作區(qū)內(nèi)的高程點和等高線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)的水系自動提取方法中，直接通過DEM進行提取，發(fā)現(xiàn)自動提取的水系與實際水系的分布有較大出入，前人的研究也發(fā)現(xiàn)這一問題，并作出了相應(yīng)的研究[10-12]，取得了一定的成果。

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

2.1.3 水系自動提取

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

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

3) 柵格流量。填充洼地后的DEM，再計算各柵格的流向。假設(shè)每一個柵格因為大氣降水產(chǎn)生流量為“1”的水量，按柵格流向向下游流動，每向下流動一個柵格，水量將增加“1”，而每一個柵格的流量便為所有流向此柵格的柵格個數(shù)。據(jù)此計算出每一個柵格點累積流量。

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

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

表1 水系特征點的分類

2.2 水系分類以及特征點提取

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

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

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

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

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

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

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

2.4 設(shè)計點位自動編號

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

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

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

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

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

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

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

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

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

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

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

3.1 試驗區(qū)概況

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

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

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

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

3.2 采樣點位自動布置效果

3.2.1 不同方法的水系提取效果對比

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

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

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

圖3 采樣點位自動布置結(jié)果圖Fig.3 Distribution map of sampling point

3.2.2 采樣點自動布置結(jié)果

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

4 結(jié)論

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

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

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

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