基于投影寻踪-博弈论-云模型的滑坡风险评价MATLAB实现

简介：这份资料围绕MATLAB实现投影寻踪博弈论-云模型的滑坡风险评价展开，面向地质灾害研究者、城市规划与环境科学从业者及灾害应急管理人员，提供从数据采集与预处理、投影寻踪博弈论建模、云模型处理到风险评估输出的一体化项目范例。资源包仅含1个docx文档，约55KB，正文以项目实例与代码详解为主线，覆盖项目背景、目标意义、挑战与解决方案、特点创新和应用领域等章节，并包含GUI设计与完整程序解析，便于按目录逐模块研读和复现。目前已有78人学习。读者可获得可直接对照的算法实现思路、多学科融合的评估流程、面向滑坡风险预测与防治建议的决策支持框架，也可将其迁移至城市规划、农林开发、环境保护及应急预案制定等场景，为后续深度学习、大数据融合等改进提供参考。

1. 从一组滑坡编录数据说起：投影寻踪、博弈论和云模型为什么要一起用

去年帮一个山区公路项目做边坡排查，手上是 24 个边坡点的编录表：坡高、坡度、岩性、结构面发育程度、年降雨量、地震动峰值加速度、人类工程活动强度，一共 8 列指标，外加一列专家给的初步等级。真正卡住人的不是算不出来，而是算出来的结论没人敢签字：换成 AHP 让专家打分，权重稍微调一调，"中风险"就跳到"高风险"；换成熵权法让数据自己说话，某一年降雨异常就把整条路的等级全顶上去了。滑坡风险评价的麻烦在于，它同时有三个毛病——指标维度太高、权重来源说不清、等级边界本身是模糊的。

投影寻踪负责第一个毛病，把 8 维指标按数据自身的聚类结构压到一维，投影方向的分量平方天然就是一组客观权重；博弈论负责第二个毛病，让主观权重和客观权重"谈判"出一个离两者偏差都最小的组合权重；云模型负责第三个毛病，用期望、熵、超熵三个数字把"低风险""高风险"这种带过渡带的定性概念量化，输出的是对每个等级的确定度而不是一刀切。MATLAB 在这条链路里承担两件事：跑遗传算法求投影方向、算云滴和确定度，再用 App Designer 把整条流程包成一个能交付给项目组点按钮的工具。

2. 投影寻踪把 8 个滑坡指标压成一维

2.1 滑坡指标的三类标准化公式与代码

投影寻踪对量纲极其敏感，坡高是几十米、降雨量是上千毫米，不标准化的话投影方向会被数值大的指标绑架。滑坡指标按危险性方向分三类：越大越危险（坡度、坡高、降雨量、地震加速度）、越大越安全（岩体完整性系数、内摩擦角）、区间型（高程、距断层距离常存在一个最不利区间）。三套公式不能混用，混用之后投影值的经济含义就乱了。

function Xn = normalize_index(X, type) % X : n×m 原始指标矩阵，n 为边坡样本数，m 为指标数 % type : 1×m 元胞，取值 'pos'（越大越危险）、'neg'（越大越安全）、'mid'（区间型） [n, m] = size(X); Xn = zeros(n, m); for j = 1:m col = X(:, j); switch type{j} case 'pos' Xn(:, j) = (col - min(col)) / (max(col) - min(col)); case 'neg' Xn(:, j) = (max(col) - col) / (max(col) - min(col)); case 'mid' a = 0.9 * mean(col); b = 1.1 * mean(col); % 最不利区间上下界 Xn(:, j) = 1 - max(a - col, col - b) ./ max(a - min(col), max(col) - b); end end Xn(Xn < 0) = 0; % 区间型可能出现负值，截断到 0 end

逻辑上，前两类是极差归一化，输出严格落在 [0,1]；区间型用"到最不利区间的相对偏离"衡量，落在区间内取 1，越远越小。a、b的取法需要结合具体指标，比如距断层距离的最不利区间通常由规范给定，而不是用均值推，我在实际项目里会把它做成normalize_index的第三个入参传进来，避免硬编码。

2.2 投影指标函数：类间散度乘类内密度

投影寻踪的核心思想是找一个方向a，让样本投影到这条直线上以后，整体尽量散开、局部尽量抱团。散开用投影值的标准差衡量，抱团用窗宽R内的点对距离累计衡量，两者相乘就是投影指标函数：

z(i) = Σ a(j)·x(i,j), j = 1..m S(a) = sqrt( Σ (z(i) - z̄)² / (n-1) ) D(a) = ΣΣ (R - r(i,j)) · u(R - r(i,j)), r(i,j) = |z(i) - z(j)| Q(a) = S(a) · D(a), s.t. ‖a‖ = 1

u(·)是单位阶跃函数，只有距离小于R的点对才计入密度。约束‖a‖=1必须显式加，否则a无限放大会让Q无上界，优化直接跑飞。

2.3 用遗传算法在 MATLAB 里求最佳投影方向

Q(a)关于a高度非线性，梯度信息不可用，常见做法是遗传算法或粒子群。MATLAB 的ga需要 Global Optimization Toolbox，没有的话换成自己写的实数编码遗传算法也能跑，逻辑一样。

function Q = pp_objective(a, Xn, R) % a : m×1 投影方向（可正可负，靠单位化约束幅值） % Xn : n×m 标准化指标 % R : 窗宽半径，控制"局部"的尺度 [n, ~] = size(Xn); a = a(:) / norm(a); z = Xn * a; Sz = sqrt(sum((z - mean(z)).^2) / (n - 1)); % 类间散度 r = abs(z - z.'); % 成对距离矩阵 Dz = sum(sum((R - r) .* (r < R))); % 类内密度 Q = Sz * Dz; end

m = size(Xn, 2); R = 0.1 * max(pdist(Xn)); % 常用经验值：0.1 倍最大样本间距离 obj = @(a) -pp_objective(a, Xn, R); % ga 求最小，目标取负 opts = optimoptions('ga', 'PopulationSize', 80, 'MaxGenerations', 300, ... 'CrossoverFraction', 0.8, 'FunctionTolerance', 1e-8, ... 'Display', 'off'); [a_best, ~] = ga(obj, m, [], [], [], [], -ones(m,1), ones(m,1), [], opts); w_obj = a_best(:).^2 / sum(a_best(:).^2); % 分量平方归一化 = 客观权重

PopulationSize取 80 是因为 8 维决策变量下 20~30 的种群容易早熟；MaxGenerations300 配FunctionTolerance1e-8，通常在 150 代左右收敛。把a的分量平方归一化当客观权重，是因为投影方向分量的绝对值反映该指标对投影值的贡献强度，平方后消除正负号，比直接取绝对值更符合"贡献占比"的语义。

2.4 窗宽 R、符号不确定性和局部最优的坑

R是投影寻踪最需要手调的参数。R太小时密度项只统计到极少数近邻点，D(a)接近 0，Q对方向不敏感；R太大时所有点对都计入，D(a)退化成常数项，优化退化为只最大化标准差。实践区间是0.05~0.3倍最大样本距离，我会在这个区间里跑 6 次看Q的最优值和对应等级排序是否稳定。

另一个反直觉的点是符号不定：a和-a给出的Q完全相同，因为标准差和成对距离都对整体符号不敏感。这意味着遗传算法跑两次可能得到镜像方向，客观权重的平方归一化不受影响，但如果直接拿a的分量当权重就会正负颠倒。我的做法是在拿到a_best后检查与主成分第一方向的相关系数，为负就整体取反，保证投影方向符号有物理含义。

提示：ga每次运行结果有随机性，交付前用rng(2024)固定随机种子，并在报告里附上Q的收敛曲线，方便评审看到迭代过程。

3. 博弈论组合赋权：让 AHP 主观权重和投影权重谈出一个折中

3.1 单一赋权在滑坡评价里的失效场景

主观赋权（AHP、专家排序）能体现规范条文和工程经验，但一致性比例CR稍微放宽，权重就会被人为放大；客观赋权（投影寻踪、熵权）忠实于样本数据结构，但样本量小的时候（滑坡编录往往只有二三十个点）容易被一两个异常样本带偏。滑坡评价里这两种失效都会发生：专家觉得岩性最重要，数据觉得降雨量最重要，各写一份报告结论相反，项目组没法用。

博弈论组合赋权的思路不是简单加权平均，而是把两个权重向量当作博弈的两个参与者，寻找一组组合系数α₁、α₂，使得组合权重与两个单一权重的偏差之和最小。它的数学形式比"取平均"更讲道理：偏离大的那个权重会被自动降低话语权。

3.2 组合权重的一阶最优条件与线性方程组

设主观权重w₁、客观权重w₂，组合权重w = α₁w₁ᵀ + α₂w₂ᵀ。以最小化‖w - w_kᵀ‖₂为目标，对α求一阶导数并令其为零，得到矩阵方程：

[ w₁w₁ᵀ w₁w₂ᵀ ] [α₁] [ w₁w₁ᵀ ] [ w₂w₁ᵀ w₂w₂ᵀ ] [α₂] = [ w₂w₂ᵀ ]

解出α后做归一化α* = α / (α₁+α₂)，再算w* = α₁*w₁ + α₂*w₂。整个求解只是一个 2×2 线性方程组，计算量可以忽略，真正花时间的是准备两个权重向量。

3.3 MATLAB 左除求解组合系数与负值处理

w1 = w1(:); w2 = w2(:); % AHP 权重与投影寻踪客观权重，均已归一化 assert(abs(sum(w1)-1) < 1e-8 && abs(sum(w2)-1) < 1e-8, '权重未归一化'); A = [w1.'*w1, w1.'*w2; w2.'*w1, w2.'*w2]; b = [w1.'*w1; w2.'*w2]; alpha = A \ b; % 解 2×2 线性方程组 if any(alpha < 0) % 出现负系数说明该权重被反向支配 alpha = lsqnonneg(A, b); % 非负最小二乘兜底 end alpha = alpha / sum(alpha); % 归一化组合系数 w = alpha(1)*w1 + alpha(2)*w2; % 最终组合权重 w = w / sum(w);

A \ b走的是 LU 分解，两个权重向量线性相关时A接近奇异，MATLAB 会给出警告；这时改用pinv(A)*b更稳。lsqnonneg的存在是因为博弈解在极端情境下会算出负系数，物理上无法解释，非负约束是工程上的必要妥协。

3.4 组合权重的合理性检验方式

组合权重算完不能直接用，需要两类检验。第一类是一致性方向检验：组合权重与主观权重的排序是否基本一致，如果某个指标在主观里排前三、组合后掉到倒数，必须回去检查客观权重的计算，通常是标准化公式选错了。第二类是敏感性检验：把α₁人为固定成 0.3/0.5/0.7 三档，看最终风险等级排序有没有跳变。

注意：α的归一化不能省。有些资料直接拿未归一化的α去乘权重，得到的结果不满足权重和为 1，后面云模型综合确定度会整体偏大或偏小。

4. 云模型把"高风险"这种模糊说法变成三个数字

4.1 Ex、En、He 的物理含义与等级云参数计算

云模型用期望Ex、熵En、超熵He描述一个定性概念。Ex是该等级最典型的取值，En是等级边界的模糊程度，He是熵本身的离散程度，也就是"云滴有多厚"。对滑坡风险常用的五级划分（极低、低、中、高、极高），每个等级给定一个取值区间[c_min, c_max]，双边约束下：

Ex = (c_min + c_max) / 2 En = (c_max - c_min) / 6 % 3En 覆盖区间宽度，对应正态分布 99.7% 范围 He = k · En, k 一般取 0.01 ~ 0.1

首末两级是单边约束（极低等级形如[0, a]，极高等级形如[b, 1]），需要用相邻等级的En反算边界，否则Ex会落在区间端点外，导致云图形状畸形。

4.2 条件云发生器求单指标确定度矩阵

有了等级云参数，下一步是把每个边坡的每个指标值代入"X 条件云发生器"，求它属于各等级的确定度。因为云模型带有随机性，单次计算不可靠，标准做法是重复采样取均值。

function mu = cloud_membership(x, Ex, En, He, N) % x : 指标值列向量（已归一化到 [0,1]） % Ex, En, He : 某一等级的云参数 % N : 采样次数，控制随机性，经验值 200~500 if nargin < 5, N = 300; end x = x(:); En_n = En + He * randn(N, 1); % En 服从 N(En, He²) mu = zeros(numel(x), 1); for i = 1:numel(x) mu(i) = mean(exp(-(x(i) - Ex).^2 ./ (2 * En_n.^2))); % 对 N 次采样取均值 end mu(mu < 1e-6) = 1e-6; % 防止后续归一化出现除零 end

En_n = En + He*randn(...)这一步是云模型区别于普通正态隶属函数的关键：每个云滴的熵本身在波动，He越大波动越大，云图看起来越"厚"。对每个指标、每个等级调用一次这个函数，拼成n×5×m的确定度张量，再按指标维度和权重做加权求和：

U = zeros(n, 5, m); for j = 1:m for s = 1:5 U(:, s, j) = cloud_membership(Xn(:, j), Ex(s), En(s), He(s)); end end B = zeros(n, 5); for j = 1:m B = B + w(j) * U(:, :, j); % w 为博弈论组合权重 end B = B ./ sum(B, 2); % 按行归一化，得到综合确定度 [~, grade] = max(B, [], 2); % 最大确定度对应的等级即评价结果

4.3 综合确定度与等级判定的边界处理

max(B, [], 2)给出的是硬判定，但云模型的优势在于保留了软信息。当最大确定度只有 0.28，第二确定度 0.26 时，硬判成"中风险"是有风险的，我一般会额外输出"主要等级 + 次要等级 + 两者差"，差小于 0.05 就标记为"等级临界"，提示人工复核。

结果归一化那一步经常被忽略。不归一化时，各指标对不同等级的确定度加总不为 1，权重大的指标会整体抬高所有等级的确定度，最后仍然能取最大值，但数值失去可比性，不同边坡之间没法横向排序。

4.4 He 取值的边界与云滴形态判断

He是云模型里最容易被拍脑袋定的参数。取太小（比如 0.001），云图退化成一条光滑隶属曲线，云模型就白用了；取太大（比如He > En/3），云滴极度离散，同一等级内两个几乎相同的指标值会得到差异很大的确定度，评价结果不可复现。经验区间是He = (0.01~0.1)·En，滑坡评价里偏保守取 0.05·En 比较合适。

验证He是否合理有个笨办法但很管用：固定输入，把He从 0.005 调到 0.05，重复跑 30 次记录每个边坡的等级，统计等级跳变率。跳变率低于 5% 说明He可接受；超过 15% 就必须调小，或者把采样次数N从 300 提到 1000 用均值压住随机性。

提示：采样次数N和超熵He是一对相互补偿的参数，但方向相反。He大而N小，结果随机且不可复现；He大而N大，结果稳定但云图过厚，模糊性被过度放大，实际意义反而变弱。

5. App Designer 封装与结果稳定性验证

5.1 计算逻辑与界面分离的代码组织

把整条链路塞进按钮回调里是最常见的做法，也是最难维护的做法。我一般建一个+slideEval包目录，里面放normalize_index.m、pp_objective.m、solve_weights.m、cloud_membership.m、evaluate.m，evaluate.m接收原始指标矩阵和配置结构体，返回组合权重、确定度矩阵和等级向量。App Designer 的.mlapp文件只负责三件事：读数据、调evaluate、画图。这样做的好处是命令行能直接批量跑 200 组参数做敏感性分析，不用去点界面。

5.2 数据导入、运行与云图绘制的回调写法

导入用uigetfile配合readtable，把指标列名和类型映射存到app.IndicatorType属性里；运行时把配置打包成结构体传进evaluate；绘图用UIAxes画正向云图，每个等级 1000 个云滴，横轴是综合值、纵轴是确定度，重叠部分一眼就能看出等级边界的模糊带有多宽。

% 按钮回调：运行评价 function RunButtonPushed(app, ~) cfg.R = app.WindowWidthEditField.Value; % 窗宽系数 cfg.N = app.SampleCountEditField.Value; % 云滴采样次数 cfg.HeK = app.HeRatioEditField.Value; % He/En 比例 X = app.RawData; % 导入时缓存 [res] = slideEval.evaluate(X, app.IndicatorType, cfg); app.WeightTable.Data = res.w(:).'; % 权重表 app.ResultTable.Data = [res.B, res.grade]; % 确定度矩阵 + 等级 app.UIAxes_Cloud.NextPlot = 'add'; for s = 1:5 [x, y] = slideEval.forward_cloud(res.Ex(s), res.En(s), res.He(s), 1000); plot(app.UIAxes_Cloud, x, y, '.', 'MarkerSize', 3); end end

forward_cloud是正向云发生器，输出云滴坐标，和cloud_membership是同一套参数的两个方向：一个从云参数生成云滴，一个从云滴反算确定度。把两者分开写，界面上画图和数据计算互不干扰。

5.3 稳定性验证：调一个参数，等级会不会跳

交付前必做的一件事是参数敏感性扫描。以窗宽系数R为例，从 0.05 到 0.30 步长 0.05 跑 6 组，记录 24 个边坡的等级序列，统计与基准组的差异个数。

差异数是个很直观的交付指标：如果某个参数在合理区间内动一动，风险等级就大面积跳变，说明这个参数不该由人拍，要么固定成规范推荐值，要么在报告里明确写出取值依据。我通常把这张表连同He的跳变率表一起放进交付文档，评审看到"参数动过、结论没塌"，签字会快很多。

本文还有配套的精品资源，点击获取