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

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

2. 中國科學院大學，北京 100049

數據庫（集）基本信息簡介

引 言

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

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

基于數據驅動的方法的一個重要建模思路是機器學習，近年來被廣泛用于地質災害的評價模擬研究。針對地質災害的評價模擬，雖然基于規(guī)則的方法在樣本數據未知的情況下具有一定的優(yōu)勢，但是機器學習方法的優(yōu)化問題求解能力和對隨機過程的模擬能力遠遠優(yōu)于基于規(guī)則的方法。溜石坡調查資料是珍貴的現場地質調查資料，該數據利用GPS測量在蓋孜河谷采集。本文基于溜石坡調查資料進行數據驅動建模得到溜石坡的空間分布數據集。本數據集在災害評價研究領域和工程地質領域都具有潛在的可重用性。

1 數據采集和處理方法

1.1 數據來源與預處理

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

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

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

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

文中采用的數據見表1。

表1 數據來源列表

1.2 數據處理步驟

1.2.1 孕災因子遴選與特征矩陣構建

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

對于任意的連續(xù)數值變量A和B，使用經典的皮爾遜相關性檢驗進行判別，相關統計量為rA,B或者其平方形式r2，在統計機器學習模型中要避免rA,B大的連續(xù)變量對同時輸入到模型中。圖3是經過以上冗余檢測后得到的孕災因子集合。

圖3 冗余檢測后的孕災因子集合

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

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

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

特征遴選、特征整合、特征編碼、特征轉換是特征工程的完整流程。特征編碼和轉換的完整處理流程見圖4。

圖4 特征編碼和轉換的流程

1.2.2 算法選擇和溜石坡調查數據擴充

野外地質調查很難采集足量的樣本來建立監(jiān)督學習模型，小樣本驅動監(jiān)督學習模型是一個挑戰(zhàn)，本文嘗試采用半監(jiān)督（Semi-Supervised）學習的方法來解決這個問題。

半監(jiān)督學習基于少量的帶標簽訓練樣本便可得到不錯的預測效果，其原理是標簽傳播，即對輸入數據集中所有樣本構建一個相似性圖來給無標簽數據賦予類別標簽。相比Zhu等[9]提出的標簽傳播方法，Zhou等[10]的方法最小化一個包含正則屬性的損失函數，因此該方法更具魯棒性。本文選擇Zhou的方法，算法實現基于Scikit-Learn[11]，可供選擇的核方法有徑向基函數（RBF）和K近鄰（KNN），這兩種核函數都是在歐式空間計算距離，徑向基函數能夠將特征映射到高維空間，但是時間復雜度和空間復雜度都很高，K近鄰具有更高的計算效率和更低的空間復雜度。

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

圖5 溜石坡調查數據疊加GF-1號影像

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

1.2.3 機器學習模型訓練與參數選擇

為了權衡預測精度和算法的復雜度（包括時間復雜度和空間復雜度），K近鄰被選用為標簽傳播模型的核函數，K進鄰是典型的基于距離的機器學習算法，又被稱作懶惰學習。K近鄰對每一個無標簽數據尋找其在訓練集中K個最近的鄰居，將鄰居中出現次數最多的類別標簽賦給無標簽數據，迭代這個過程直到所有樣本都賦予標簽。K近鄰的計算流程如圖6。

圖6 基于KNN的半監(jiān)督學習計算過程

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

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

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

2 數據樣本描述

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

圖7中展示了溜石坡調查點的分布情況和溜石坡易發(fā)性預測結果的空間分布，紫色區(qū)域是預測為溜石坡的區(qū)域。

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

3 數據質量控制和評估

3.1 原始數據質量控制

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

3.2 溜石坡預測結果評價

實驗證明，少量樣本訓練出來的半監(jiān)督學習模型還是比較理想的，在驗證集上正樣本精度能夠達到57%以上，負樣本精度在99%以上，總體精度接近80%（圖8）。首先，負樣本的預測精度拉高了模型的總體精度；其次，基于精度的分析可以得到，如果一個格網點不是溜石坡就一定不會被預測為溜石坡，但如果是溜石坡便極有可能被預測為非溜石坡，也就是說即使溜石坡的預測精度不是很高，但是非溜石坡不會錯判為溜石坡。

圖8 模型精度評價

4 數據使用方法和建議

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

致 謝

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