组合赋权-改进可拓云模型在磷矿山巷道围岩稳定性评价中的Matlab实现

去年做磷矿山巷道稳定性评估时甲方丢给我一摞地质报告要求把每条巷道的围岩稳定性分个等级出来。我一开始用的模糊综合评价卡在隶属度函数上——边界怎么定都像在凑数。换BP神经网络吧样本量又不够训练完自己都不敢信。后来才转向组合赋权-改进可拓云模型把主观经验和客观数据揉在一起用正态云替代刚性区间这个问题才真正落地。这套流程跑通之后我顺手把Matlab代码整理了出来就是标题里这套东西。这篇文章我想把整套思路完完整整写出来为什么磷矿山巷道评价不能直接套用煤矿的模板、组合赋权为什么比单一赋权可靠、可拓云模型的“改进”到底改在哪、以及Matlab代码每一段在干什么。适合正在做巷道稳定性评价、岩体分级、或者准备写相关论文但卡在模型实现细节上的朋友参考。1. 磷矿山巷道围岩稳定性评价的根本难点1.1 磷矿山巷道围岩跟煤矿、金属矿有什么不一样磷矿大多赋存在海相沉积岩系里围岩常见白云岩、硅质岩、碳质页岩岩层呈层状构造层理和节理面发育。跟煤矿相比磷矿山岩巷的稳定性问题更突出岩石强度不算低但层间结合力弱开挖卸荷后容易沿层理面滑移或离层跟金属矿相比磷矿山巷道断面又偏大加上部分地段有溶蚀现象围岩性质在水平方向和竖直方向变化都很快。这里有个很实际的问题煤矿巷道稳定性评价常用锚杆支护设计中的围岩分类方法直接拿来评价磷矿山岩巷经常会偏保守或偏冒进。原因是煤矿巷道要考虑瓦斯、采动压力而磷矿山岩巷的核心矛盾是岩体结构和软弱夹层。所以评价指标体系必须结合磷矿山的实际岩性特征来定不能照搬。1.2 传统评价方法的三块短板我最早试过三种主流方法各有各的别扭。第一是模糊综合评价。它的隶属度函数需要人为指定比如“单轴抗压强度属于I级的隶属度是0.8”这个0.8主观性太强换个人打分结果就可能不一样。而且模糊评价的隶属度在边界处是突变的指标值从79MPa变到81MPa评价结果可能直接从II级跳成I级工程上很难接受。第二是灰色聚类。它不依赖大样本但对指标灰类的白化函数同样需要人为设定半梯形分布和三角形分布的划分逻辑都有点生硬。用在地质参数确实很多不确定性的场景里模型本身是合适的但结果解释起来比较费劲。第三是BP神经网络。理论上它能自动学习指标与等级之间的非线性关系但磷矿山巷道样本往往只有十几个、二十几个这点数据量根本喂不饱神经网络。我试过一次训练集准确率很高验证集惨不忍睹典型的过拟合。这些方法的问题归结起来就是一句话要么主观性过强要么对样本数量要求过高要么在边界处理上不够平滑。后来我意识到需要把“主观权重客观权重”结合用云模型的随机性和模糊性去消化边界问题这就是组合赋权-改进可拓云模型的基本出发点。1.3 组合赋权改进可拓云到底解决什么问题这套模型解决的实际上是两个层面的事。权重层面AHP反映现场工程师的经验判断熵权法反映数据本身的离散程度两者组合之后权重既不会因为某个专家拍脑袋而走偏也不会因为样本数据异常而失真。评价层面可拓学提供了“物元-经典域-关联度”的框架而云模型把经典域从固定区间变成了正态云相当于把“硬边界”做成了“软边界”。从实际效果看这套组合特别适合矿山岩土类课题指标数据有限但是有工程经验、指标边界模糊但存在统计规律、需要输出直观等级但又不希望结论对边界过于敏感。2. 组合赋权AHP和熵权法的互补逻辑2.1 主观权重层次分析法到底在算什么层次分析法AHP的解题思路是既然六个指标同时评价很困难那就两两对比把复杂的多指标判断拆成若干对简单判断。具体操作是构造判断矩阵Aa_ij表示指标i相对指标j的重要程度用1-9标度打分。比如我认为岩石单轴抗压强度Rc比岩体完整性系数Kv略重要那a_12取2反过来a_21取1/2。对角线恒为1。判断矩阵构造完之后求最大特征值对应的特征向量归一化后就是权重向量。这一步Matlab直接用eig函数做。需要重点检查的是一致性比例CR。判断矩阵是两两打分凑出来的难免出现“A比B重要、B比C重要、但C比A重要”这种逻辑矛盾。CR的计算公式是CR CI / RI其中CI (λ_max - n) / (n - 1)RI是随机一致性指标6阶矩阵对应RI1.24。CR小于0.1才认为判断矩阵可以接受否则要回头调整打分。我在论文和工程报告里见过不少直接跳过一致性检验的AHP应用这是硬伤。因为CR不过关时权重向量的意义就很弱。2.2 客观权重熵权法怎么从数据里“挤”出权重熵本来是个热力学概念信息论借它来刻画数据的不确定性。熵权法的核心逻辑非常朴素指标在所有样本上取值差异越大说明它承载的区分信息越多权重就应该越高取值都差不多的指标区分不了样本权重就低。计算分三步。第一步把每个指标下的样本值转成比重p_ij第二步计算信息熵e_j第三步由差异系数1-e_j归一化得到权重。值得提醒的是熵权法不受指标量纲影响因为同一指标列的数值都会先除以该列总和量纲和尺度差异被抵消了。所以哪怕Rc是几十的单位Kv是0.x的小数直接放进去算没问题不用先做无量纲化。但熵权法有个天然的毛病它对数据质量敏感。如果某个指标的值在样本里只有极少数异常大熵值会被拉得很低权重被抬得过高这属于“数据驱动”的副作用。2.3 组合权重乘法归一化与线性加权怎么选组合方式有两种常见做法。第一种是乘法归一化w_comb w_AHP_i * w_entropy_i再归一化。这种方法的思路是“两个权重都高的指标才是真正重要的”对某个权重接近0的指标惩罚很厉害。适合AHP和熵权法结果比较一致的情况。如果两者出现明显分歧乘法会让分歧变得很极端。第二种是线性加权w_comb alpha * w_AHP (1 - alpha) * w_entropy其中alpha是主观权重占比一般取0.5-0.7。这种方法温和很多可以用alpha调节主观经验和客观数据的信任程度。我最终在项目里用的是乘法归一化因为没有出现极端分歧。实际使用中你可以两个都算一遍如果结果差异大再检查数据或判断矩阵是否出了问题。权重计算不是一次性的事需要反复核对。3. 可拓云模型的改进逻辑从刚性边界到柔性云3.1 可拓学的基本框架物元、经典域和节域可拓学这套理论是蔡文提出的核心概念是物元。一个物元用三元组R(N, C, V)表示其中N是评价对象C是特征指标V是特征值。比如“巷道Rc62MPa”就是一个物元特征。评价时把稳定性等级分成若干级每个等级对应每个指标都划分一个范围这个范围叫经典域。比如Rc在50-80MPa属于III级那III级的经典域就是区间[50,80]。全部等级的范围并集叫节域即指标可能取值的总区间。传统可拓评价的做法是计算待评样本到各等级经典域的可拓距再用关联函数算关联度。关联度哪个等级最大样本就属于哪个等级。问题就出在这个可拓距上——它是基于固定区间计算的样本值在经典域边界处可拓距会发生明显跳变。我前面提到的“79MPa是II级、81MPa变III级”的尴尬在传统可拓模型里依然存在。3.2 云模型三参数期望Ex、熵En、超熵He究竟在描述什么云模型是李德毅提出来的用三个参数描述一个具有模糊性和随机性的概念。期望Ex是概念的中心值可以理解为等级区间的代表值。熵En在云模型里不完全是“熵”的意思它表示这个概念的可接受范围En越大云越宽说明边界越模糊。超熵He则描述熵本身的随机波动程度He越大云滴越离散、随机性越强。这三个参数合起来就可以把“指标值属于某个等级”这件事从“非黑即白”变成“一个云滴落在某个范围内的概率分布”。对矿山岩土这类充满不确定性的场景这种表达比硬邦邦的区间合理得多。用生活化类比来说传统可拓区间像一个有明确围墙的院子墙内墙外泾渭分明云模型则像一片被雾气笼罩的区域中心地带浓厚、边缘地带稀薄你无法说清楚“雾的边界”到底在哪但判断一个位置是否属于这片雾是可以通过密度来估计的。3.3 改进点经典域云化和关联度计算流程所谓改进可拓云模型最关键的一步是把各等级各指标的经典域区间[a,b]转换成云参数(Ex, En, He)。转换规则通常是Ex (a b) / 2En (b - a) / 6He取En的0.1倍左右其中En除以6的依据是正态分布的3σ原则即区间长度对应6个标准差这样区间内的覆盖率接近99.7%。“改进”的意味在于传统可拓模型把[a,b]作为硬边界而这里用正态云覆盖区间边界处的隶属度是连续衰减的不再是突变。对于单边界等级比如I级“Rc≥120”和V级“Rc25”需要设定虚拟区间。我通常把I级取[120,160]V级取[0,25]再按同样的规则转成云参数。虚拟区间的大小会影响边界等级的关联度实际操作需要结合岩体力学参数经验确定不能随意拍。关联度计算采用X条件云发生器针对某个指标值x和某个等级云(Ex, En, He)生成N个满足N(En, He²)的随机数En计算每个云滴下的高斯隶属度exp(-(x-Ex)²/(2En²))然后取N次结果的平均值作为关联度。N一般取1000-2000太小结果不稳。这里有个容易踩的坑传统可拓学的关联度可以为负值样本在经典域外部时而云模型的关联度基于高斯函数始终在0到1之间。这是两种理论的差异不是代码写错了。云关联度可以理解成“隶属于该等级的概率密度强度”不能套用可拓学里负关联度的解释。4. 指标体系、分级标准与云参数生成4.1 磷矿山巷道评价指标怎么选指标选取要遵循三个原则能反映围岩稳定性本质、数据在工程现场能拿到、指标间相对独立。我最终敲定的六项指标是岩石单轴抗压强度RcMPa岩块强度的直接度量磷矿山白云岩和硅质岩的强度差异非常大岩体完整性系数Kv反映裂隙发育程度现场通过声波测试获得RQD%岩芯采取率中大于10cm的累计长度占比是国际通用的岩体质量指标单位涌水量qL/min反映地下水影响磷矿山层间裂隙水对围岩弱化明显埋深Hm反映地应力水平磷矿山大都几百米深埋深越大构造应力越明显巷道跨度Lm反映开挖尺寸效应跨度越大围岩自稳能力越弱这六个指标一个偏工程经验RQD一个偏地应力环境埋深一个偏开挖设计跨度剩下的偏岩体基本性质覆盖面比较完整。4.2 五级评价标准与云参数生成规则等级划分为五级I级稳定、II级基本稳定、III级欠稳定、IV级不稳定、V级极不稳定。分级标准参考了工程岩体分级标准结合磷矿山岩巷的特点做一些调整具体如表所示。指标I级II级III级IV级V级Rc (MPa)≥12080~12050~8025~5025Kv≥0.750.55~0.750.35~0.550.15~0.350.15RQD (%)≥9070~9050~7025~5025q (L/min)0.50.5~2.52.5~55~10≥10H (m)200200~400400~600600~800≥800L (m)33~44~66~8≥8有了这个表就能为每个指标的每个等级生成对应的云参数。以修改后的等级区间计算Ex和EnHe取En/10。例如Rc的III级区间[50,80]Ex65En(80-50)/65He0.5。4.3 He取值不当的后果He参数的调试是我在这个模型里花时间最多的地方。一开始我把He设成En/2发现云滴发散得很厉害同一个样本跑三次综合关联度的排序结果都不太稳定。后来把He降下来结果就稳定多了。这里给个经验值He建议取En的0.05~0.2倍。En本身已经代表区间边界的不确定性He再大的话随机性会盖过规律性评价结果的可复现性就差。反过来He如果取0.001那云模型就退化成了固定高斯曲线所谓“云”的随机性又消失了。我在代码里把He设成了可配置项参数调优时可以批量跑不同He值对比结果稳定性这是最笨但也最稳妥的办法。5. Matlab代码落地从数据到判定结果5.1 主程序框架设计整套代码我按模块拆成四个部分AHP权重计算、熵权法权重计算、组合赋权、可拓云关联度评价。主程序只负责数据录入和调用函数这样以后换数据、换指标、换等级标准都不用改主体逻辑。主程序的流程是输入判断矩阵、输入样本矩阵、输入待评样本、调用AHP函数得到w_AHP、调用熵权函数得到w_entropy、计算组合权重w_comb、读取云参数表、计算综合关联度、输出判定等级。%% 主程序示例 clear; clc; % 1. 判断矩阵6个指标 A [1 2 3 4 5 6; 1/2 1 2 3 4 5; 1/3 1/2 1 2 3 4; 1/4 1/3 1/2 1 2 3; 1/5 1/4 1/3 1/2 1 2; 1/6 1/5 1/4 1/3 1/2 1]; % 2. 样本矩阵5个巷道段6个指标 X [72 0.61 62 0.6 320 4.2; 62 0.52 55 0.8 380 4.5; 45 0.33 38 2.3 520 5.2; 95 0.71 78 0.3 250 3.6; 28 0.18 22 4.5 680 6.8]; % 3. 待评样本 x0 [62 0.52 55 0.8 380 4.5]; % 4. 计算权重 [CR, w_AHP] ahp_weight(A); w_entropy entropy_weight(X); w_comb w_AHP .* w_entropy; w_comb w_comb / sum(w_comb); % 5. 可拓云评价 level cloud_evaluate(x0, w_comb); disp([评价等级: , num2str(level)]);5.2 AHP权重函数AHP这块核心就是矩阵特征值分解和一致性检验。我用eig函数求最大特征值和对应特征向量特征向量归一化后就是权重。判断矩阵构造不复杂但CI、CR这些指标一定要算出来并打印。function [CR, w] ahp_weight(A) n size(A, 1); [V, D] eig(A); lambda_max max(diag(D)); [~, idx] max(diag(D)); w abs(V(:, idx)); w w / sum(w); CI (lambda_max - n) / (n - 1); RI_table [0 0 0.58 0.90 1.12 1.24 1.32 1.41 1.45]; RI RI_table(n); CR CI / RI; if nargin 1, fprintf(CR %.4f\n, CR); end endRI_table直接内置到函数里查表方便。矩阵阶数超过9时这个RI表不够用但工程评价里指标一般不会超过9个。5.3 熵权法函数熵权法实现起来很直接。唯一要留心的是p_ij里可能出现零值log(0)会得到-Inf代码里做一个防零处理用1e-10替换。这个细节不处理的话某个样本在某个指标上恰好是最小值程序就会直接崩。function w entropy_weight(X) [m, n] size(X); P X ./ sum(X, 1); P max(P, 1e-10); e -sum(P .* log(P), 1) / log(m); d 1 - e; w d / sum(d); end这里m是样本数除以log(m)是为了把信息熵归一化到0~1之间。熵值越大说明数据越均匀权重就越低。5.4 云关联度计算与等级判定云关联度实现的核心是X条件云发生器。对每个指标和每个等级生成N个正态随机数En用高斯隶属度公式计算关联度再取均值。这里N设成1000既保证稳定性又不会太慢。function k cloud_relation(x, Ex, En, He, N) if nargin 5, N 1000; end En_prime normrnd(En, He, N, 1); k mean(exp(-(x - Ex).^2 ./ (2 .* En_prime.^2))); end等级云参数的生成我写了一个独立的脚本输入等级区间数组输出Ex、En、He矩阵。这样等级标准调整时只需要改区间的几行数据不用改主程序。综合关联度方面先算待评样本每个指标对五个等级云的关联度得到6×5的矩阵然后用组合权重加权求和得到最终的1×5综合关联度向量取最大值对应的等级。function level cloud_evaluate(x0, w_comb) load(cloud_params.mat, Ex, En, He); n_ind length(x0); n_level size(Ex, 2); K zeros(n_ind, n_level); for i 1:n_ind for j 1:n_level K(i, j) cloud_relation(x0(i), Ex(i, j), En(i, j), He(i, j)); end end K_total w_comb * K; [~, level] max(K_total); end这里cloud_params.mat是我预先算好保存的云参数文件也可以用函数直接生成。实际工程中我更推荐写成函数调用因为不同矿山的分级标准可能不一样写成参数文件改起来更灵活。6. 实例演算一个磷矿山巷道段的完整判定6.1 样本数据与判断矩阵以某磷矿山5条巷道段的地质资料为样本集待评对象是2号巷道段。原始数据在前面主程序代码里已经出现这里完整列出。样本1Rc72MPaKv0.61RQD62%q0.6L/minH320mL4.2m 样本2待评Rc62MPaKv0.52RQD55%q0.8L/minH380mL4.5m 样本3Rc45MPaKv0.33RQD38%q2.3L/minH520mL5.2m 样本4Rc95MPaKv0.71RQD78%q0.3L/minH250mL3.6m 样本5Rc28MPaKv0.18RQD22%q4.5L/minH680mL6.8m判断矩阵按4.1节的相对重要性构造六阶矩阵在代码中已经给出。运行AHP函数后CR约为0.036小于0.1通过一致性检验。不同的人对“Rc比Kv重要多少”这件事看法可能不一样但判断矩阵的评分只要CR过关权重方向基本是一致的。6.2 权重计算结果AHP权重算下来Rc权重最高达到0.377Kv次之为0.246。这符合磷矿山岩巷的实际情况岩石本身强度是硬指标层理和节理裂隙是第二位的。熵权法基于5个样本的数据各指标权重相对均衡其中Rc和跨度L的权重略高因为这几个指标在样本间的差异更明显。组合赋权用乘法归一化后最终权重如表所示。指标AHP权重熵权法权重组合权重Rc0.3770.1800.399Kv0.2460.1650.239RQD0.1630.1700.163q0.0990.1520.088H0.0630.1560.058L0.0520.1770.054组合权重中Rc和Kv两指标合计超过0.63说明这个模型判定结果主要受岩石强度和岩体完整性控制。地下水涌水量和跨度的权重有所下降不是它们不重要而是样本间这些指标的离散程度没有Rc和Kv显著。6.3 综合关联度与等级判定结果按等级标准生成云参数后对2号巷道段实测值计算云关联度。以Rc62MPa为例它对III级别对应的云Ex65En5He0.5关联度最高达到0.682对II级云Ex100En6.67关联度降到0.004附近对IV级云Ex37.5En4.17关联度约0.005。边界效应的确有体现但不像可拓距那样跳变。六个指标全部计算并加权后综合关联度向量为等级I级II级III级IV级V级综合关联度0.1850.3120.4130.0730.017最大综合关联度对应III级即“欠稳定”。这个判定结果与现场实际情况吻合2号巷道段存在轻微片帮局部有掉块需要及时喷锚支护但短时间内不会整体失稳。值得关注的是II级关联度0.312也不算低这说明该巷道段处于II-III级过渡带建议支护设计时按III级偏安全考虑。这个“过渡带”信息传统可拓模型很难体现出来因为它的关联度往往只有一个等级明显占优。7. 我在实操中踩过的坑和调参经验7.1 判断矩阵一致性调不过去怎么办判断矩阵是组合赋权里最依赖人的部分。我第一次给某矿做评价时判断矩阵打分比较随意结果CR算出来0.19远远超过0.1的限值。这时候直接改权重是没有用的问题出在判断矩阵内部有环状的逻辑矛盾。我的处理方法是先打印出各指标两两比较的差值看哪个比较导致了不一致。比如指标A比B重要3倍B比C重要3倍但A只比C重要2倍这就矛盾了因为传递下来A应该比C重要9倍左右。修正这种矛盾项CR通常能很快降下来。如果调整单个元素还是不行这个判断矩阵干脆扔掉重做不要硬调。判断矩阵的意义是表达工程经验的相对关系打分的前后一致性比具体数值本身更重要。7.2 云滴数和He的调参经验云滴数N对结果稳定性的影响我已经提过补充一个实测数字N从100增加到1000同一样本综合关联度的排序结果我第一次跑就变了N到1000再跑到5000结果基本稳定。所以1000是可以接受的下限追求稳妥就用2000。He的取值则要看数据本身的离散程度。如果同一等级的样本内部差异本来就大比如某矿的III级巷道Rc从48MPa到75MPa都有那He适当取大一点更贴近实际如果样本比较集中He取太大反而会让不同等级的云大面积重叠导致关联度失去区分度。我最终的做法是He En / 10作为默认值然后用待评样本做敏感性分析即He在0.05En到0.2En之间各算一遍看等级判定结果是否稳定。如果结果对He非常敏感说明样本本身处于等级边界附近这时任何模型都不可能给出一个绝对可靠的结论需要结合现场监测数据综合判断。7.3 边界样本和多等级并列时的判定评价中一定会碰到“怎么判都觉得两个等级都合理”的样本。比如综合关联度向量是[0.35, 0.33, 0.31, 0.01, 0.00]三个等级几乎并列。这时候强行取最大关联度可能没有实际意义。我的处理方式是增加一个“可区分度”指标取最大关联度与第二大关联度的比值。比值大于1.2认为等级判定清晰小于1.1则判定为明显的过渡带报告中给出“II~III级过渡带”的结论并建议提高监测频率。工程上“模糊地把握大致风险水平”比“精确地给出一个可能错误等级”更有价值。另外单个指标云关联度低不一定是坏事。比如埋深H380m落在II级但Rc、Kv、RQD都落在III级综合结果判定III级就是合理的。真正要警惕的是某指标的关联度对所有等级都很低比如全部小于0.05这往往意味着该指标值超出了分级标准的节域范围需要回头检查数据录入或者重新标定等级区间。7.4 对代码复现的一句总结这套Matlab代码看起来不长但里面每一步都对应一个容易出错的细节AHP的特征向量要用绝对值再归一化、熵权法要防log(0)、云关联度里的normrnd要指定He作为标准差、组合权重归一化不能忘。把这些细节处理好模型跑通之后换矿山、换指标、换分级标准都是改参数的事不需要动核心代码。我个人现在的习惯是接到一个新矿山的数据后先不急着跑模型而是把该矿历年的巷道破坏记录翻出来用它来反推判断矩阵和He的取值。模型参数永远是为工程服务的而不是反过来让工程迁就模型参数。这个习惯帮我避掉了好几个“模型跑得通但现场不认账”的尴尬场景。