基于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]。