梁娟娟,楊耀民,朱志偉,杜德文*,王春娟,葉 俊,閆仕娟,楊 剛
(1.國家海洋局 第一海洋研究所,山東 青島266061;2.國家深海基地管理中心,山東 青島266061)
多波束測深數據是獲取海底微地形地貌的主要手段[1]。它相對重力、磁力、地震等地球物理勘探手段,數據覆蓋率更高、獲取速度更快、預算成本更低,并能生成高精度的數字高程模型(DEM),提取地形特征,進行海底地形分類。
基于DEM提取地形特征是遙感信息進入地學數據庫的關鍵技術之一[2]。1964年,Hammond[3]提出通過DEM統計單元內的坡度、相對起伏度來劃分地形類型的想法,之后由Dikau等[4]、Brabyn[5]、Morgan和Lesh[6]、Dragut和 Blaschke[7]、Prima等[8]、Iwahashi和 Pike[9]進行了實現、發展、修正以及完善。Buea等[10]也曾基于DEM提取地形特征對火星地形進行過非監督分類。
前人對海底地形地貌也做過諸多研究。耿秀山[11]、林美華[12]、潘定安等[13]對近海及海灣地形進行了深入研究;陶春輝[14]用分形技術對多波束測量的南海大陸坡海底火山地形進行了研究;王英等[15]從多波束資料提取水深和坡度信息,探討地形與海底多金屬結核的分布關系;劉忠臣等[16]根據多波束調查對東海的地形進行了研究;楊剛等[17]提出了基于多波束測量數據進行海底地形分割的“三態值模型”法,可對數字海底地形進行快速準確地分割。Herzfeld和Higginson[18]利用統計學方法對大西洋中部海脊進行分類,建立了一套海底地形自動分類系統。Bishwajit[19]基于多波束原始水深數據,利用神經網絡對海底的海脊、裂谷和平原地形進行分類。這說明,基于DEM對海底進行地形統計分類是學者們認可的方法。前人工作多單純針對地形分類,少有利用地形分類技術進行礦產資源勘探靶區的圈定。
本文基于DEM的地形統計分類,將南大西洋中脊示范區進行地形分類,并將分類結果與熱液硫化物航次調查信息進行關聯,尋找見礦大概率地形區,將其作為下一步找礦詳查區域,實現找礦勘探靶區逐步縮小的目的。
選取南大西洋中脊13°~14°S區段,針對每個水深點提取地形特征;劃分統計單元,提取每個單元地形特征的統計參數;然后對這些統計單元進行非監督分類;將地形類型圖層與硫化物礦點進行關聯,篩選出熱液硫化物礦點的大概率地形單元,并將其視為后繼調查工作的重點靶區。
南大西洋中脊是慢速擴張洋脊,構造地質背景與北大西洋脊較為相似。我國在13°S,15°S和26°S附近發現并成功抓取了熱液硫化物,表明南大西洋中脊具有極大的成礦潛力。選取南大西洋13°~14°S脊段作為研究區(圖1),該區域位于洋中脊兩組轉換斷層之間,南側是洋中脊與轉換斷層的交界地帶,北端位于非轉換偏移區域,區域整體位于洋中脊裂谷東側。

圖1 研究區地理位置Fig.1 Location of the study area
多波束測深數據:主要來源于“大洋一號”科學考察船22航次(2012-12-2011-12)調查多波束全覆蓋聲納測深數據,文中所用數據是由隨船科學家對現場實測多波束數據插值所得的GRD數據,精度約為0.001°×0.001°。文中基于GRD數據生成同精度的DEM圖層和網格化的水深點圖層。
站位調查數據:我國在本研究區31個地質站位進行了取樣,其中有7站取到硫化物樣品,24站未取到硫化物樣品,調查工作主要由“大洋一號”科學考察船21航次(2009-07-2010-05),22航次(2010-12-2011-12)和26航次(2012-04-2012-12)完成。國外報道中未見有本研究區取樣記錄。
通過地形特征的定性分析,選取水深、坡度、坡度變率、總曲率等微觀地形特征以及地形起伏度、地表切割度、高程變異系數等宏觀地形特征,在ArcGIS中基于DEM提取各地形特征并生成對應的特征圖層。將網格化的水深點圖層與各特征圖層疊加并提取各特征值。
統計單元劃分是礦產資源定量評價工作的重要步驟之一,按照不同礦種、不同調查數據精度,確定統計單元的劃分方案,地質體單元和網格單元是最常用的劃分方法[20]。針對海底多金屬硫化物勘探實際工作,無法選擇地質體單元。結合多波束地形勘測精度,選擇擁有適當數目網格的網格組作為統計單元。
2.2.1 統計單元劃分的依據
本文劃分統計單元的目的是為了將點參數轉化為統計單元參數,優點是增加了變量,降低了以點參數進行分類造成的隨機性;缺點是降低了數據精度。統計單元大小取決于所獲取數據資料的精度,關系到潛在礦產資源的空間位置和數量的估計精度。統計單元越小,分類精度越高,但統計參數的代表性越弱;統計單元越大,分類精度越低,但統計參數的代表性越強[21]。所以,劃分統計單元時要考慮2點:1)統計單元有足夠多的數據,保證統計參數具有代表性;2)統計單元的規格要盡量小,保證地形分類的精度。
2.2.2 統計單元的劃分
1)研究區網格整體劃分
考慮到研究區的實際情況以及所獲取數據的精度(0.001°×0.001°),設定該區統計單元大小為近似0.01°×0.01°的正方形區域,將其劃分為8 475個統計單元,每個統計單元包含100個網格單元。
2)選取有點參數的統計單元
劃分統計單元時是將研究區作為規則的矩形區,但是由于研究區的不規則性,在區域外圍有很多統計單元內沒有水深數據點,因此我們要選取包含水深數據點的統計單元。
3)選取有效的統計單元
在篩選出來包含水深數據點的統計單元中,為了削弱邊界統計單元數據點太少對統計結果的影響,剔除點參數不足(少于10個點)的統計單元,最后運用ArcGIS的空間分析選取了符合統計條件的4 267個統計單元。
每個統計單元包含約100個網格化水深數據點(不考慮邊界效應),對每個統計單元中的7個地形特征進行統計參數提取,主要包括平均值、方差、最大值、最小值。最后得到每個統計單元中水深、坡度、坡度變率、總曲率、地表切割深度、地表起伏度、高程變異系數七個地形特征的平均值、方差、最大值、最小值,共計28個變量。
因涉及的變量較多,數據處理麻煩且會造成較大的誤差,因此要對變量進行篩選,去除相關性較強以及對分析結果貢獻不大的變量[22]。為了消除量綱的影響,首先將28個變量進行標準化。傳統選取變量的方法是因子分析和相關分析相結合,但由此方法篩選出的變量分類結果并不理想,因此作者采用不同因子組合的方法選取最佳的變量組合。
經過反復驗證變量組合并與研究區實際海底地形相對比,最終得到水深與坡度的統計變量組合分類效果最佳,因此,最終參與地形分類的變量為水深均值、方差、最大值、最小值,坡度均值、方差、最大值、最小值,共8個變量。
由于相同地形類型實體具有最大的相似性、最小的差異性,不同地形類型實體具有最小的相似性、最大的差異性[23],因此從相同地形類型所提取的同種變量參數,其值也更為接近,而從不同類型的地形實體提取的地形因子則存在較大的差異。這種相似性和差異性即是地形形態類型劃分的依據[23]。
根據這種劃分依據,我們將篩選出的8個變量在統計軟件SPSS中實現非監督分類,運用的方法是K-均值聚類,設置迭代次數為10次,聚類數為5類,得到的分類結果如表1。從表1可知,地形5只包含一個統計單元,因此忽略。

表1 最終聚類中心Table 1 The final clustering center based on cluster analysis
1)地形1 水深最大值、最小值與均值都比較高,水深的方差較小;坡度的最大值、最小值、方差、均值都比較低。說明這一類地形總體水深值較大,水深變化起伏小;坡度值較小,起伏也較小。結合實際海底地形,可推知此類地形屬于“裂谷”。
2)地形2 坡度最大值、均值和方差都較高,坡度最小值較小;水深最大值、最小值為4類地形中最低的,水深的均值較小,方差的相對較大。可看出此類地形地勢相對較高,坡度起伏較大,坡度高值相對較多。結合實際海底地形,可推知此類地形屬于“裂谷壁”。
3)地形3 坡度最大值、最小值、均值都比較高,尤其是坡度最小值相比其他3類地形高出許多,但坡度的方差不是很大;水深方差相比其他地形高出很多,最大值、最小值和均值都相對較小。結合海底地形,可推知此類地形屬于“內角高地斜坡”。
4)地形4 8個變量的值都比較平均,沒有很高或很低的值,水深最大值、最小值以及均值較地形1小,地形2、3大,水深方差較小;坡度各項變量值都較小,水深變化不大,地形起伏度較小,整體趨勢較平坦。對比實際海底地形,可推知此類地形屬于“高地”。
通過研究區實際海底地形與分類結果對比(圖2),分類結果具有一定的可信度,但地形4分類效果不好,結果具有一定的誤差,因此單獨將地形4取出重新選取變量進行分類。地形4整體地勢較高,但粗糙度卻有所差異,因此,在水深-坡度的基礎上再引入粗糙度進行二次分類。綜合第一次分類結果,最終分類結果如圖3a。這樣將地形4分成2類:地形4和地形5。二次分類后的地形4共1 052個統計單元,為“洋中脊高粗糙度的高地地形”;地形5共1 093個統計單元,為“洋中脊低粗糙度的高地地形”。
綜上所述,我們將大西洋中脊研究區地形整體分為5類:地形1為“裂谷”,地形2為“裂谷壁”,地形3為“內角高地斜坡”,地形4為“高粗糙度的高地”,地形5為“低粗糙度的次高地”。

圖3 調查數據與二次地形分類結果關聯及遠景區Fig.3 Relationship between the survey data and secondary landform classification results and proposed prospecting target area
將本區所做的航次調查結果與最終所得的地形類型圖層疊加,如圖3a。
為尋找大概率地形,統計了本區各類地形的航次調查的結果(表2)。表2中站位數=地質取樣站位數;見礦數=地質取樣取得硫化物樣品的站位數;未見礦數=地質取樣未取得硫化物樣品的站位;見礦概率=(見礦數/站位數)×100%;未見礦概率=(未見礦數/站位數)×100%;100網格見礦率系數=(見礦數/單元數)×100。

表2 各類地形調查結果Table 2 Survey results for all landform types in the study area
分析表2可知,地形5“低粗糙度的次高地”見礦概率為28.6%,在5類地形中最高,其次是地形3“內角高地斜坡”和地形4“高粗糙度的高地”。
各類地形的面積差異很大,為了標度單位面積見礦率,引入100網格見礦率系數。統計單元的精度為0.01°×0.01°,每個統計單元的面積約為1.21km2,100網格見礦率系數即面積為121km2范圍內硫化物的產出概率。結果顯示地形3“內角高地斜坡”100網格見礦率系數為2.00,為5類地形中最高;地形4“高粗糙度的高地”與地形5“低粗糙度的次高地”分別為0.19和0.18。表明研究區相同面積區域內“內角高地斜坡”發育熱液硫化物的概率最大,為硫化物產出的大概率地形,如圖3b黑線圈定區域。考慮到地質調查站位并不是隨機分布于各地形類型內,100網格見礦率系數在指示各地形類型的成礦概率方面也只能作為參考指標。
本文綜合考慮見礦概率和100網格見礦率系數,姑且將“內角高地斜坡”作為首選后繼勘探靶區。隨著調查認識的深入,會有更完善的指標,圈定的后繼重點勘探靶區也會更合理。
本區共7個熱液硫化物站位,其中3個位于本文方法圈定出的遠景區內,4個熱液硫化物站位處于漏判狀態。漏判率=(漏判熱液硫化物站位/熱液硫化物總站位)×100%=57%。因此漏判率超過一半。
研究區共4 267個統計單元,總面積為5 163.07km2,其中圈定出的硫化物發育的大概率地形“內角高地斜坡”有150個統計單元,總面積為181.5km2。靶區縮小倍數:靶區縮小率=(勘探區域面積/總區域面積)×100%=3.5%。
地形分類與已知礦點關聯的目的是為了縮小勘探靶區,以圈定熱液硫化物的后繼勘探區。縮小靶區的過程不可避免會造成礦點的漏判。雖然漏判率超過50%,但是我們卻將靶區縮小到3.5%,即我們用超過一半的漏判率換回了3.5%的靶區縮小率。
將大西洋中脊劃分為5類地形:“裂谷”、“裂谷壁”、“內角高地斜坡”、“高粗糙度的高地”和“低粗糙度的次高地”。各地形見礦概率:“裂谷”為0,“裂谷壁”為0,“內角高地斜坡”為27.3%,“高粗糙度的高地”為20%,“低粗糙度的次高地”為28.6%;各地形的100網格見礦率系數:“裂谷”為0,“裂谷壁”為0,“內角高地斜坡”為2.0,“高粗糙度的高地”為0.19,“低粗糙度的高地”為0.18。其中,“內角高地斜坡”100網格見礦率系數最大,將其認定為發育硫化物的大概率地形,建議作為下一步的重點勘探區域。雖然該方法礦點漏判率高達57%,卻獲得了將靶區縮小到3.5%的效果。
研究區調查程度低,地質調查站位少,已知硫化物礦點亦少,本方法基于這些有限數據信息建立起來的認識,為熱液硫化物后繼勘探靶區的選擇提供一種思路,并不能排除其他地形類型也可能是熱液硫化物礦點產出的有利區域。
(References):
[1]YANG S J.Linear structural interpretation of Mid-Atlantic ridge and associate with hydrothermal mineralization[J].Acta Mineralogica Sinica,2011,(Suppl.1):706-707.閆仕娟.大西洋中脊區線性構造解譯及熱液成礦關聯[J].礦物學報,2011,(增刊1):706-707.
[2]Lü G N,QIAN Y D,CHEN Z M.Automated extraction of the characteristics of topography from grid digital elevation data[J].Journal of Geographical Sciences,1998,53(6):562.閭國年,錢亞東,陳鐘明.基于柵格數字高程模型提取特征地貌技術研究[J].地理學報,1998,53(6):562.
[3]HAMMOND E H.Analysis of properties in landform geography:an application to broad-scale landform mapping[J].Annals of Association of American Geographers,1964,54(1):11-19.
[4]DIKAU R,BRABB E E,MARK R M.Landform classification of New Mexico by computer[C].U.S.Geological Survey,Open-file Report,1991,26:91-634.
[5]BRABYN L.GIS analysis of macro landform[C]∥The 10th Colloquium of the Spatial Information Research Centre.New Zealand,1998:35-48.
[6]MORGAN J M,LESH A M.Developing landform maps using ESRI'S Model Builder[EB/OL].[2014-05-20].http:∥gis.esri.com/library/userconf/froc05/papers/pap2206.pdf,2005.
[7]DRAGUT L,BLASCHKE T.Automated classification of landform elements using object-based image analysis[J].Geomorphology,2006,81(3-4):330-344.
[8]PRIMA O D A,ECHIGO A,YOKOYAMA R,et a1.Supervised landform classification of Northeast Honshu from DEM derived the matic maps[J].Geomorphology,2006,78(3-4):373-386.
[9]IWAHASHI J,PIKE R J.Automated classifications of topography from DEMs by an unsupervised nested-means algorithm and a threepart geometric signature[J].Geomorphology,2007,86(3-4):409-440.
[10]BUEA B D,STEPINSKIB T F.Automated classification of landforms on Mars[J].Computers and Geosciences,2006,30(5):604-614.
[11]GENG X S.The geomorphic classification in the East China Sea[J].Transactions of Oceanology and Limnology,1980,2(3):32-37.耿秀山.東中國海的地貌分類[J].海洋湖沼通報,1980,2(3):32-37.
[12]LIN M H.The submarine geomorphological zones and geomorphological types in the Huanghai Sea[J].Marine Sciences,1989,13(6):7-15.林美華.黃海海底地貌分區及地貌類型[J].海洋科學,1989,13(6):7-15.
[13]PAN D A,WANG S M,SHEN H T.Approach to the formation cause of the central deep channel and the Beiniu sands in Meizhou bay[J].Acta Geographica Sinica,1994,49(1):55-63.潘定安,汪思明,沈煥庭.湄洲灣中央深槽及白牛淺灘的成因探討[J].地理學報,1994,49(1):55-63.
[14]TAO C H.The application of the fractal geometry in the study of submarine volcano[J].Donghai Marine Science,2000,18(2):9-14.陶春輝.分形幾何在海底火山地形中的應用[J].東海海洋,2000,18(2):9-14.
[15]WANG Y,LI J B,HAN X Q,et al.The influence of terrain slope on the distribution of poly metallic nodules[J].Acta Oceanologica Sinica,2001,23(1):60-64.王英,李家彪,韓喜球,等.地形坡度對海底多金屬結核的分布控制作用[J].海洋學報,2001,23(1):60-64.
[16]LIU Z C,CHE Y L,ZHOU X H,et al.Study on zoned characteristics and formation cause of the East China Sea submarine topography[J].Advances in Marine Science,2003,21(2):160-173.劉忠臣,陳義蘭,周興華,等.東海海底地形分區特征和成因研究[J].海洋科學進展,2003,21(2):160-173.
[17]YANG G,DU D W,Lü H L.Segmentation algorithm for digital sea-floor terrain[J].Advances in Marine Science,2004,22(2):204-209.楊剛,杜德文,呂海龍.數字海底地形分割算法[J].海洋科學進展,2004,22(2):204-209.
[18]HERZFELD U C,HIGGINSON C A.Automated Geostatistical seafloor classification-principles,parameters,feature vectors,and discrimination criteria[J].Computer and Geosciences,1996,22(1):35-52.
[19]BISHWAJIT C.Application of artificial netural network to segmentation and classification of topographic profiles of ridge-flank seafloor[J].Current Science,2003,85(3):306-312.
[20]ZHANG F Y,ZHANG W Y,ZHU K C,et al.Principle of resource evaluation and technique of mining area delineate for ocean polymetallic crust[M].Beijing:China Ocean Press,2001:1-150.張福元,章偉艷,朱克超,等.大洋多金屬結殼資源評價原理和礦區圈定方法[M].北京:海洋出版社,2001:1-150.
[21]LIU Y G,DU D W,LI Z S,et al.Estimation of polymetallic nodule distribution and resource quantity in the CC zone and its adjacent areas of the Pacific Ocean[J].Advances in Marine Science,2009,27(3):342-349.劉永剛,杜德文,李鐘山,等.太平洋 CC區及周邊多金屬結核分布及資源量預測[J].海洋科學進展,2009,27(3):342-349.
[22]WANG C J,DU D W,MENG X W,et al.Mineral resources prediction model and its application about gas hydrate in the Gulf of Mexico[J].Acta Oceanologica Sinica,2009,31(3):67-71.王春娟,杜德文,孟憲偉,等.墨西哥灣天然氣水合物礦產資源預測模型及其應用[J].海洋學報,2009,31(3):67-71.
[23]LIU A L,TANG G A.DEM based auto-classification of Chinese landform[J].Geo-information Science,2006,8(4):8-16.劉愛利,湯國安.中國地貌基本形態 DEM 的自動劃分研究[J].地球信息科學,2006,8(4):8-16.