基于Crystal Ball的泥沙來源指紋判別及不確定性評價

史忠林，陳佳村,2，嚴冬春，文安邦，龍 翼，吳 毅

(1.中國科學院/水利部 成都山地災害與環境研究所，四川 成都 610041；2.中國科學院大學，北京 100049)

流域泥沙來源的研究由來已久，方法眾多。傳統的研究方法主要有徑流小區法、水沙資料分析法和大面積人工或遙感調查法等[1]。20世紀70年代以來，基于泥沙性質的指紋識別技術興起并快速發展，目前已成為流域泥沙來源研究的熱點技術之一。指紋技術的基本原理為利用某一種或一組在不同泥沙源地間具有顯著差異的理化性質作為指示因子，建立源地物質與目標泥沙之間的對應關系，定量研究各源地對目標泥沙的相對貢獻。指紋技術的發展大體經歷了單因子識別與復合因子識別兩個階段。早期研究使用礦物性質、地質特征、磁性、顆粒顏色及植物孢粉等作為示蹤因子。20世紀90年代以來，放射性核素137Cs、210Pb、7Be等相繼用于侵蝕泥沙示蹤研究，指紋因子的種類不斷豐富[2]?？紤]到單因子在診斷多物源時固有的局限性，研究者開始嘗試采用多因子組合代替單一指紋因子，復合指紋技術應運而生[3]。復合指紋法目前在國外研究較多且較為深入[4]，在國內則相對較少且以應用性探索為主，對技術本身的研究相對欠缺[5-7]。

盡管復合指紋法歷經了20多年發展已取得重要進展，但人們在研究中逐漸認識到該技術的應用仍存在諸多不確定性，例如對潛在泥沙源地的預判、樣品采集的代表性、指紋因子在源匯轉化過程中的穩定性、物源與泥沙指紋因子間的可比性及混合模型的選擇與優化等。早前的研究多利用經統計檢驗篩選得到的指紋因子的平均值或中位數來表征不同物源與目標泥沙的差異，然而土壤本身固有的空間變異性及對侵蝕強度響應程度的不同，造成利用采集的少量樣品來代表某一物源或泥沙時存在較大的不確定性，這就極大地影響到物源貢獻判別結果的可靠性。為此，近年來研究者嘗試采用不同的分布函數(如正態分布[8]、t分布[9])來描述物源或泥沙中的指紋因子，然后通過蒙特卡洛隨機抽樣方法實現混合模型的優化求解和物源貢獻結果的不確定性評價。

已有的泥沙來源研究中，蒙特卡洛模擬運算主要通過Matlab[10]、R[11]等軟件平臺進行，要求具有計算機基礎，這對一般研究者來說具有一定挑戰性。為此，本研究介紹一種簡便快捷的分析軟件Crystal Ball (V.11.1)，并結合實例分析其在流域泥沙來源判別中的應用。

1 蒙特卡洛方法及Crystal Ball概述

蒙特卡洛方法是一種以概率統計為理論基礎的隨機抽樣模擬方法，通過估計不確定變量的概率分布，多次抽樣模擬試驗得到結果變量的概率分布，從而進行風險評價及敏感度分析。

假定函數Y=f(x1,x2,…,xn)，其中x1,x2,… ,xn為隨機變量，其概率分布已知。通過重復抽樣取出每一組隨機變量的值(x1,x2,…,xn)，然后按函數關系式確定函數值Y=f(x1,x2,…,xn)。反復獨立抽樣若干次，可得到函數Y的一批抽樣數據Y1,Y2,…,Yk，當達到合理的模擬次數時，樣本特征可以近似反映總體的特征。若誤差在可接受范圍內，這時可以得出與實際情況相近的函數Y的概率分布及其數字特征[12]。

Crystal Ball中文名為水晶球，是美國Oracle公司基于Windows平臺開發的一款商業分析和評估軟件，廣泛應用于工程投資和公共安全等領域。Crystal Ball是加載在Excel電子表格上的第三方程序，其界面友好，操作簡便，分析結果以圖表和報告的形式呈現，直觀易懂。用戶可以通過定義假設變量的概率分布，用計算機產生的隨機數進行多次模擬計算，獲得滿足約束條件下決策變量的一組解。

2 泥沙來源指紋判別分析案例

采集云南省元謀縣涼山鄉某小流域農耕地、林地、溝岸等三種潛在泥沙源地的土壤樣品和流域出口處泥沙樣品，測定樣品理化性質，運用復合指紋技術分析各物源對流域出口處泥沙的相對貢獻率。測試指標經范圍檢驗、Kruskal-Wallis H檢驗和多元判別分析后，確定S、Rb、U、Y四種元素構成最佳指紋組合，統計結果見表1。

2.1 混合模型

利用多元混合模型定量計算各源地泥沙貢獻百分比，模型結構[3]為

式中：Res為殘差平方和；Cssi為流域出口處泥沙中指紋因子i的濃度；Csi為泥沙源地s中指紋因子i的濃度；Ps為泥沙源地s的泥沙貢獻百分比；m為泥沙源地數量，m=3；n為指紋因子的數量，n=4。

模型應用需滿足兩個基本前提，即各源地泥沙貢獻率Ps非負且總和為1。

2.2 變量定義

Crystal Ball軟件中需要用戶定義的變量有三個，即假設變量、決策變量和預測變量。假設變量是Crystal Ball中基本的輸入量，可以是單一的數值，也可以是服從某一分布的不確定變量。本例中假設變量為目標泥沙中指紋因子的濃度Cssi和各源地指紋因子的濃度Csi。Crystal Ball 11.1分布庫中提供了正態分布、三角形分布、均勻分布等21種預先設定的分布類型和一種自定義分布，用戶可根據需要自行選擇，也可根據實測值進行擬合。但用實測值擬合時，要求實測值的數量至少為15個。本例假設目標泥沙和源地土壤中各指紋因子濃度均服從正態分布，用均值和標準差依次定義各指紋因子的分布函數，并將其下限設置為0[13]。

決策變量為每個潛在物源的泥沙貢獻百分比(Ps)，依次定義P1、P2、P3，將其下限和上限分別設為0和1，變量類型設置為連續型。

在完成假設變量和決策變量定義后，需要定義預測變量。預測變量單元格為包含一個或多個假設變量和決策變量單元格的公式。本例中預測變量為各指紋因子濃度的殘差平方和，根據混合模型的公式編輯單元格。在運行模擬時，預測單元格會自動關聯模型中的各個單元格并計算生成最終結果。

2.3 模型運行

點擊運行Crystal Ball軟件中的OptQuest優化求解器。由于混合模型的求解要求參與運算的所有指紋因子濃度殘差平方和最小，因此將“目標”設置為最小化預測單元格的最終值。同時，在“約束”選項中將決策變量之和指定為1。

優化求解時，需要對模擬參數進行設置，主要包括試驗次數和抽樣方法等。軟件中隨機數抽樣方法有蒙特卡洛抽樣和拉丁超立方抽樣兩種，其中前者所產生的隨機數之間相互獨立，而后者把假設變量的概率分布分成幾個等概率區間，區間數量可根據需要設定，模擬試驗時將為每個區間產生不同系列的隨機數[14]。本例中將模擬次數設為1 000次，抽樣方法采用蒙特卡洛隨機抽樣。

2.4 優化求解結果與不確定性分析

將混合模型模擬運算1 000次后，統計各物源泥沙貢獻率的頻率分布，見圖1。由圖1可知，三種物源的泥沙貢獻頻率呈尖峰分布，集中趨勢顯著，表明模型計算結果的不確定性較小。

指紋因子濃度的殘差平方和(預測變量)與各物源對目標泥沙的相對貢獻百分比(決策變量)統計結果列于表2。三種源地中，溝岸對流域出口處泥沙的貢獻率最大，達到78.3%±17.1%；其次為農耕地，為17.1%±11.6%；林地的貢獻率最小，僅為4.6%±16.1%。據調查，泥沙來源判別結果符合該流域實際情況。

3 結 論

運用Crystal Ball軟件的優化求解功能，結合蒙特卡洛模擬，成功實現了某流域泥沙來源及其相對貢獻率的定量判別，模型運算結果與流域實際情況吻合。Crystal Ball通過對指紋因子的分布擬合，極大地降低了指紋因子單一取值帶來的模型求解結果的不確定性。同時，模型輸出結果給出了不同物源對流域泥沙貢獻率的平均值和標準差，相比傳統的Excel規劃求解只給出唯一解而言，提高了判別結果的科學性和合理性。此外，Crystal Ball操作簡便，不需要復雜的計算機語言，具有較強的實用性，未來將在流域泥沙來源研究中得到更廣泛的應用[15]。