蓋孜河谷溜石坡調(diào)查和易發(fā)性分布數(shù)據(jù)集

田德宇，張耀南，韓立欽，康建芳，羅立輝，艾鳴浩，敏玉芳

1. 中國(guó)科學(xué)院西北生態(tài)環(huán)境研究院，科學(xué)大數(shù)據(jù)中心，蘭州 730000

2. 中國(guó)科學(xué)院大學(xué)，北京 100049

數(shù)據(jù)庫(kù)（集）基本信息簡(jiǎn)介

引 言

中巴公路從帕米爾高原東北緣穿過(guò)了公格爾-慕士塔格-昆蓋山的冰緣地帶（Preglacial Region）（圖1），這里屬于塔里木河的支流喀什葛爾河流域的主要組成部分。崩塌及其次生災(zāi)害溜石坡在蓋孜河谷的公路兩側(cè)廣泛發(fā)育，堆積物有砂石也有大的礫石，其緩慢前進(jìn)對(duì)公路邊坡?lián)p害較大，堆積物甚至可以蔓延到路面，在強(qiáng)大外力作用下的迅速推進(jìn)容易徹底堵塞交通[1]。特別是，溜石坡是中巴公路北部地區(qū)廣泛發(fā)育的一類天然崩積坡，呈錐形或面坡形態(tài)沿公路延伸，其形態(tài)、顆粒、沉積等特征有別于川藏公路沿線的典型溜砂坡[2]。王宗盛[1]對(duì)中巴公路（國(guó)內(nèi)段）的溜石坡災(zāi)害有過(guò)調(diào)查研究，但是少見(jiàn)針對(duì)這類災(zāi)害的評(píng)價(jià)模擬，也未形成溜石坡的分布數(shù)據(jù)集，因此本工作對(duì)于溜石坡災(zāi)害的評(píng)價(jià)模擬研究具有較為重要的意義，在工程地質(zhì)防治方面具有啟迪和應(yīng)用價(jià)值。

圖1 蓋孜河谷所在流域示意圖

基于數(shù)據(jù)驅(qū)動(dòng)的方法的一個(gè)重要建模思路是機(jī)器學(xué)習(xí)，近年來(lái)被廣泛用于地質(zhì)災(zāi)害的評(píng)價(jià)模擬研究。針對(duì)地質(zhì)災(zāi)害的評(píng)價(jià)模擬，雖然基于規(guī)則的方法在樣本數(shù)據(jù)未知的情況下具有一定的優(yōu)勢(shì)，但是機(jī)器學(xué)習(xí)方法的優(yōu)化問(wèn)題求解能力和對(duì)隨機(jī)過(guò)程的模擬能力遠(yuǎn)遠(yuǎn)優(yōu)于基于規(guī)則的方法。溜石坡調(diào)查資料是珍貴的現(xiàn)場(chǎng)地質(zhì)調(diào)查資料，該數(shù)據(jù)利用GPS測(cè)量在蓋孜河谷采集。本文基于溜石坡調(diào)查資料進(jìn)行數(shù)據(jù)驅(qū)動(dòng)建模得到溜石坡的空間分布數(shù)據(jù)集。本數(shù)據(jù)集在災(zāi)害評(píng)價(jià)研究領(lǐng)域和工程地質(zhì)領(lǐng)域都具有潛在的可重用性。

1 數(shù)據(jù)采集和處理方法

1.1 數(shù)據(jù)來(lái)源與預(yù)處理

在同一個(gè)坐標(biāo)空間對(duì)慕士塔格山周圍的ASTER GDEM產(chǎn)品、STRM DEM產(chǎn)品以及ASTER L1T單景DEM制作三維散點(diǎn)圖（圖2），發(fā)現(xiàn)ASTER L1T單景DEM（圖中紅色）有明顯異常，ASTER GDEM（圖中綠色）和STRM DEM（圖中黃色）無(wú)明顯異常，但是STRM DEM最大值比慕士塔格峰主峰海拔多出200 m以上，而ASTER GDEM的最大值最接近最高值，最終選擇ASTER GDEM參與到整個(gè)研究中。

土壤類型數(shù)據(jù)依據(jù)聯(lián)合國(guó)糧農(nóng)組織制作的 FAO土壤分類體系，擁有優(yōu)于 1公里的水平分辨率[3]。在這個(gè)分類體系下，不同類型土壤數(shù)據(jù)定義可能不同：比如淺層土（Leptosol，LP）是依據(jù)土壤形成的地形條件定義，主要指侵蝕的高地。鈣積土（Calcisol，CL）、石膏土（Gypsisol，GY）、栗鈣土（Kastanozem，KS）、黑土（Phaeozem，PH）等是依據(jù)土壤形成的氣候、時(shí)間以及有機(jī)體等條件定義。該數(shù)據(jù)可在寒區(qū)旱區(qū)科學(xué)數(shù)據(jù)中心免費(fèi)獲?。╤ttp://westdc.westgis.ac.cn/data/844010ba-d359-4020-bf76-2b58806f9205），中國(guó)境內(nèi)的數(shù)據(jù)源為第二次全國(guó)土地調(diào)查的1∶100萬(wàn)土壤數(shù)據(jù)。將研究區(qū)的數(shù)據(jù)重投影為UTM Zone 43N World WGS-84坐標(biāo)系，重采樣為30 m的分辨率。根據(jù)這份數(shù)據(jù)，研究區(qū)內(nèi)除冰川外共有土壤類型12種?；谶@份數(shù)據(jù)的對(duì)比分析發(fā)現(xiàn)：其一，這份數(shù)據(jù)的冰川（代碼11930）范圍與Landsat 8 運(yùn)行陸地成像儀（Operational Land Imager，OLI）融雪末影像提取的最新冰川的邊緣差異很大；其二，冰川末梢冰緣地帶的土壤類型均為薄層土。因此對(duì)獲得研究區(qū)的這份土壤類型數(shù)據(jù)進(jìn)行修正，調(diào)整后的土壤類型包含5類，分別是淺層土、黑土、石膏土、鈣積土、栗鈣土，此外還有冰川、水體2種水體。

區(qū)域巖性地圖來(lái)自漢堡大學(xué)的Hartmann等[4]制作的全球巖性矢量數(shù)據(jù)，這份巖性數(shù)據(jù)在研究區(qū)范圍內(nèi)的巖性圖由新疆地質(zhì)礦產(chǎn)局采集自1992年，數(shù)據(jù)格式是shapefile文件，比例尺1∶150萬(wàn)，在冰川區(qū)內(nèi)無(wú)數(shù)據(jù)。依據(jù)此概略巖性圖，研究區(qū)內(nèi)以碳酸鹽沉積巖（Carbonate Sedimentary Rocks，SC）和松散沉積物（Unconsolidated Sediments，SU）為主，此外還有混合沉積巖（Mixed Sedimentary Rocks，SM）和酸性火成巖（Acid Plutonic Rocks，PA）。以1∶150萬(wàn)的區(qū)域巖性圖為地面真值概覽底圖，基于地質(zhì)學(xué)的理論知識(shí)，利用ASTER TIR和Landsat OLI傳感器的數(shù)據(jù)提取巖性指標(biāo)，來(lái)對(duì)區(qū)域巖性圖的巖性類型做出更細(xì)致的解釋，最終制備一份30 m空間分辨率的巖性地質(zhì)填圖。

溜石坡調(diào)查數(shù)據(jù)是重要的地面真值數(shù)據(jù)（Ground Truth），中科院成都山地災(zāi)害研究所周公旦老師提供了2017年融雪期的地質(zhì)災(zāi)害調(diào)查數(shù)據(jù)。該數(shù)據(jù)集包含多處溜石坡災(zāi)害點(diǎn)的經(jīng)緯度位置信息?；谡溆跋駡D的室內(nèi)目視解譯對(duì)這份調(diào)查數(shù)據(jù)進(jìn)行核實(shí)，并用來(lái)進(jìn)行溜石坡災(zāi)害易發(fā)性的機(jī)器學(xué)習(xí)建模。圖1包含了災(zāi)害調(diào)查點(diǎn)分布的概覽信息。利用GF-1號(hào)PMS數(shù)據(jù)融合鑲嵌的中巴公路喀什至伊斯蘭堡段正射影像圖用作調(diào)查底圖和輔助地質(zhì)調(diào)查的資料。

文中采用的數(shù)據(jù)見(jiàn)表1。

表1 數(shù)據(jù)來(lái)源列表

1.2 數(shù)據(jù)處理步驟

1.2.1 孕災(zāi)因子遴選與特征矩陣構(gòu)建

基于精度可靠的 DEM 提取的各種地形因子是地質(zhì)災(zāi)害易發(fā)性分析的重要的孕災(zāi)因子。本文基于Conrad等[5]開(kāi)源的SAGA GIS系統(tǒng)提取多種孕災(zāi)因子用于冰緣災(zāi)害的評(píng)價(jià)研究中，這些因子可以分為基本的地形因子、水文學(xué)因子、形態(tài)測(cè)量學(xué)因子（Morphometry）以及巖土水力學(xué)因子4類。這些地形因子與巖性類型、土壤類型一道作為溜石坡易發(fā)性評(píng)價(jià)的孕災(zāi)因子，全面刻畫(huà)了溜石坡災(zāi)害的孕災(zāi)環(huán)境，但是難免會(huì)有冗余。借鑒數(shù)據(jù)倉(cāng)庫(kù)的理念[6]，整合（integration）特征來(lái)消除冗余（redundancies）和不一致（inconsistencies）。離散特征的冗余檢測(cè)采用皮爾遜χ2統(tǒng)計(jì)量q的假設(shè)檢驗(yàn)，其計(jì)算如式（1）[7]。皮爾遜χ2檢驗(yàn)的原假設(shè)H0是兩個(gè)變量相互獨(dú)立，當(dāng)且僅當(dāng)q<χ21-p(df-1)時(shí)接受H0。其中，p是顯著性水平，自由度df=(c-1)×(r-1)[7]。在χ2檢驗(yàn)中，針對(duì)離散變量A和B，假設(shè)A有c個(gè)唯一值為a1,a2,…,ac，B有r個(gè)唯一值為b1,b2,…,br。A的c個(gè)值為列，B的r個(gè)值為行，構(gòu)成列聯(lián)表（contingency table）。在這個(gè)表中，每個(gè)a值會(huì)和對(duì)應(yīng)的b值相遇構(gòu)成聯(lián)合事件（Ai,Bj），占據(jù)矩陣的一個(gè)位置。相關(guān)統(tǒng)計(jì)量稱作皮爾遜χ2統(tǒng)計(jì)。式中，oij是聯(lián)合事件（Ai,Bj）的觀測(cè)頻率，eij是該聯(lián)合事件的期望頻率：

對(duì)于任意的連續(xù)數(shù)值變量A和B，使用經(jīng)典的皮爾遜相關(guān)性檢驗(yàn)進(jìn)行判別，相關(guān)統(tǒng)計(jì)量為rA,B或者其平方形式r2，在統(tǒng)計(jì)機(jī)器學(xué)習(xí)模型中要避免rA,B大的連續(xù)變量對(duì)同時(shí)輸入到模型中。圖3是經(jīng)過(guò)以上冗余檢測(cè)后得到的孕災(zāi)因子集合。

圖3 冗余檢測(cè)后的孕災(zāi)因子集合

統(tǒng)計(jì)機(jī)器學(xué)習(xí)模型往往需要輸入樣本服從一定的分布，至少要服從一定的輸入規(guī)范。本文選擇的半監(jiān)督學(xué)習(xí)方法會(huì)涉及到歐式距離的計(jì)算，因此需要對(duì)離散特征編碼和對(duì)連續(xù)特征標(biāo)準(zhǔn)化（Standardization）處理。

特征矩陣中包含4個(gè)離散變量和9個(gè)連續(xù)變量。離散變量中包括巖性類屬、地表分類[8]、土壤類型3個(gè)類別變量和摩擦敏感失穩(wěn)指數(shù)1個(gè)二值變量。利用OneHot編碼器來(lái)對(duì)4個(gè)離散變量編碼，經(jīng)過(guò)OneHot編碼后4個(gè)離散變量轉(zhuǎn)換為了31維的二進(jìn)制變量，這個(gè)31維的二進(jìn)制特征和連續(xù)特征組成包含40個(gè)特征的特征矩陣。

新的特征矩陣中有比率標(biāo)度也有區(qū)間標(biāo)度，有的單位是弧度也有的無(wú)單位，而且取值范圍參差不齊，在建模之前要進(jìn)行標(biāo)準(zhǔn)化，這里選擇最小最大標(biāo)準(zhǔn)化方法來(lái)處理。對(duì)于任意特征X，最小最大標(biāo)準(zhǔn)化方法首先利用樣本的最大值和最小值來(lái)將其約束為?，之后繼續(xù)將?標(biāo)準(zhǔn)化到用戶定義的區(qū)間[a, b]為?。式中min和max分別為求序列的最小值和最大值的函數(shù)。

特征遴選、特征整合、特征編碼、特征轉(zhuǎn)換是特征工程的完整流程。特征編碼和轉(zhuǎn)換的完整處理流程見(jiàn)圖4。

圖4 特征編碼和轉(zhuǎn)換的流程

1.2.2 算法選擇和溜石坡調(diào)查數(shù)據(jù)擴(kuò)充

野外地質(zhì)調(diào)查很難采集足量的樣本來(lái)建立監(jiān)督學(xué)習(xí)模型，小樣本驅(qū)動(dòng)監(jiān)督學(xué)習(xí)模型是一個(gè)挑戰(zhàn)，本文嘗試采用半監(jiān)督（Semi-Supervised）學(xué)習(xí)的方法來(lái)解決這個(gè)問(wèn)題。

半監(jiān)督學(xué)習(xí)基于少量的帶標(biāo)簽訓(xùn)練樣本便可得到不錯(cuò)的預(yù)測(cè)效果，其原理是標(biāo)簽傳播，即對(duì)輸入數(shù)據(jù)集中所有樣本構(gòu)建一個(gè)相似性圖來(lái)給無(wú)標(biāo)簽數(shù)據(jù)賦予類別標(biāo)簽。相比Zhu等[9]提出的標(biāo)簽傳播方法，Zhou等[10]的方法最小化一個(gè)包含正則屬性的損失函數(shù)，因此該方法更具魯棒性。本文選擇Zhou的方法，算法實(shí)現(xiàn)基于Scikit-Learn[11]，可供選擇的核方法有徑向基函數(shù)（RBF）和K近鄰（KNN），這兩種核函數(shù)都是在歐式空間計(jì)算距離，徑向基函數(shù)能夠?qū)⑻卣饔成涞礁呔S空間，但是時(shí)間復(fù)雜度和空間復(fù)雜度都很高，K近鄰具有更高的計(jì)算效率和更低的空間復(fù)雜度。

溜石坡災(zāi)害調(diào)查點(diǎn)數(shù)據(jù)于2017年夏天采集自中巴公路蓋孜河谷兩側(cè)，數(shù)據(jù)格式是矢量點(diǎn)文件。將災(zāi)害調(diào)查點(diǎn)數(shù)據(jù)疊加到 GF-1號(hào)正射影像中可以清楚地看到溜石坡的形狀和分布特點(diǎn)，發(fā)育在公路兩側(cè)的溜石坡基本呈扇形延伸到河谷以及公路，有的甚至可以蔓延到對(duì)岸。圖5a是昆蓋山東坡的一個(gè)溜石坡點(diǎn)，圖5b是蓋孜河谷公格爾山西坡上的一個(gè)溜石坡災(zāi)害點(diǎn)，公路等基礎(chǔ)設(shè)施以及一些村落直接建在溜石坡上，潛在危害巨大。

圖5 溜石坡調(diào)查數(shù)據(jù)疊加GF-1號(hào)影像

本研究采取基于格網(wǎng)點(diǎn)的機(jī)器學(xué)習(xí)模型構(gòu)建策略，也就是說(shuō)每個(gè)格網(wǎng)點(diǎn)是一個(gè)樣本。選用的標(biāo)簽傳播式半監(jiān)督學(xué)習(xí)算法對(duì)樣本量的要求不高，只要將未知標(biāo)簽和已知標(biāo)簽的比例控制在60∶1左右便可，即有類別標(biāo)簽的格網(wǎng)點(diǎn)占比為1/60左右。通過(guò)溜石坡調(diào)查點(diǎn)和GF-1號(hào)正射影像圖疊加可以明顯地看到溜石坡的輪廓（圖5），基于溜石坡調(diào)查資料勾繪得到的矢量圖形作為模型訓(xùn)練的標(biāo)簽數(shù)據(jù)，特征矩陣作為模型訓(xùn)練的初始特征集合進(jìn)行模型訓(xùn)練。

1.2.3 機(jī)器學(xué)習(xí)模型訓(xùn)練與參數(shù)選擇

為了權(quán)衡預(yù)測(cè)精度和算法的復(fù)雜度（包括時(shí)間復(fù)雜度和空間復(fù)雜度），K近鄰被選用為標(biāo)簽傳播模型的核函數(shù)，K進(jìn)鄰是典型的基于距離的機(jī)器學(xué)習(xí)算法，又被稱作懶惰學(xué)習(xí)。K近鄰對(duì)每一個(gè)無(wú)標(biāo)簽數(shù)據(jù)尋找其在訓(xùn)練集中K個(gè)最近的鄰居，將鄰居中出現(xiàn)次數(shù)最多的類別標(biāo)簽賦給無(wú)標(biāo)簽數(shù)據(jù)，迭代這個(gè)過(guò)程直到所有樣本都賦予標(biāo)簽。K近鄰的計(jì)算流程如圖6。

圖6 基于KNN的半監(jiān)督學(xué)習(xí)計(jì)算過(guò)程

顯然，核函數(shù)K近鄰的一個(gè)重要的超參數(shù)是用于鄰域距離計(jì)算的鄰域點(diǎn)個(gè)數(shù)K，K值直接影響模型的0精度，因此得到最優(yōu)的K值是至關(guān)重要的參數(shù)選擇問(wèn)題。調(diào)整K值來(lái)查看驗(yàn)證集上的精度。當(dāng)K=18時(shí)，驗(yàn)證集中正樣本、負(fù)樣本以及全局精度達(dá)到最大。另一個(gè)重要的參數(shù)是迭代次數(shù)，不過(guò)迭代次數(shù)達(dá)到30就已經(jīng)收斂且是穩(wěn)定的，所以調(diào)整迭代次數(shù)到收斂即可。

1.2.4 溜石坡易發(fā)性空間區(qū)域預(yù)測(cè)

利用建立的模型遍歷公路兩側(cè)2公里的范圍得到中巴公路在蓋孜河谷沿線的潛在溜石坡易發(fā)區(qū)。雖然建模特征中沒(méi)有包含任何光學(xué)遙感數(shù)據(jù)的光譜特征和紋理特征，但是模型在蓋孜河谷公路兩側(cè)2公里的緩沖區(qū)內(nèi)的預(yù)測(cè)結(jié)果與GF-1號(hào)正射影像圖中目視解譯出來(lái)的災(zāi)害區(qū)高度匹配，即使目視看來(lái)有不少疑似錯(cuò)判的地方。

2 數(shù)據(jù)樣本描述

本數(shù)據(jù)集提供蓋孜河谷的溜石坡調(diào)查GPS點(diǎn)矢量文件和基于調(diào)查點(diǎn)模擬得到的中巴公路兩側(cè)2公里范圍內(nèi)的溜石坡易發(fā)性分布柵格文件。

圖7中展示了溜石坡調(diào)查點(diǎn)的分布情況和溜石坡易發(fā)性預(yù)測(cè)結(jié)果的空間分布，紫色區(qū)域是預(yù)測(cè)為溜石坡的區(qū)域。

圖7 溜石坡易發(fā)性分布圖

3 數(shù)據(jù)質(zhì)量控制和評(píng)估

3.1 原始數(shù)據(jù)質(zhì)量控制

參與巖性填圖的Landsat OLI數(shù)據(jù)采取了嚴(yán)格的大氣校正，而ASTER TIR數(shù)據(jù)的巖性指數(shù)反演直接基于星上輻射值（radiance at sensor）進(jìn)行，這是經(jīng)過(guò)前人研究論證的[12]。DEM數(shù)據(jù)的質(zhì)量保證是經(jīng)過(guò)比較不同來(lái)源的DEM產(chǎn)品來(lái)選擇最接近研究區(qū)真實(shí)高程的DEM產(chǎn)品。土壤類型數(shù)據(jù)經(jīng)過(guò)了基于Landsat OLI提取的融雪末冰川范圍的修正來(lái)保證該數(shù)據(jù)在研究區(qū)的真實(shí)性。對(duì)研究區(qū)1∶150萬(wàn)巖性地圖的遙感填圖擴(kuò)充了巖性的類屬，遙感巖性填圖的準(zhǔn)確性通過(guò)設(shè)置合理的波段指數(shù)下限閾值來(lái)保證，即使用波段指數(shù)提出者所建議的閾值或更高的閾值進(jìn)行巖性制圖。

3.2 溜石坡預(yù)測(cè)結(jié)果評(píng)價(jià)

實(shí)驗(yàn)證明，少量樣本訓(xùn)練出來(lái)的半監(jiān)督學(xué)習(xí)模型還是比較理想的，在驗(yàn)證集上正樣本精度能夠達(dá)到57%以上，負(fù)樣本精度在99%以上，總體精度接近80%（圖8）。首先，負(fù)樣本的預(yù)測(cè)精度拉高了模型的總體精度；其次，基于精度的分析可以得到，如果一個(gè)格網(wǎng)點(diǎn)不是溜石坡就一定不會(huì)被預(yù)測(cè)為溜石坡，但如果是溜石坡便極有可能被預(yù)測(cè)為非溜石坡，也就是說(shuō)即使溜石坡的預(yù)測(cè)精度不是很高，但是非溜石坡不會(huì)錯(cuò)判為溜石坡。

圖8 模型精度評(píng)價(jià)

4 數(shù)據(jù)使用方法和建議

溜石坡調(diào)查點(diǎn)數(shù)據(jù)為矢量SHP格式，屬性表中記錄了經(jīng)度和緯度信息。溜石坡易發(fā)性分布數(shù)據(jù)保存為柵格TIF格式，二者疊加顯示最佳。ArcGIS、QGIS、ENVI、ERDAS等常用的GIS與遙感軟件可支持該數(shù)據(jù)的讀取和操作。

致 謝

感謝國(guó)家特殊環(huán)境、特殊功能觀測(cè)研究臺(tái)站共享服務(wù)平臺(tái)給與的項(xiàng)目支持。感謝中科院成都山地所周公旦研究員提供2017年的溜石坡調(diào)查資料。