☰
DEMATEL-ISM模型Python实现:从专家打分到层级结构全流程
2026/10/4 5:40:17 网站建设 项目流程

做影响因素分析的时候,最怕的就是手上一堆专家打分数据,算完描述性统计之后不知道怎么往下走。你当然可以给每个因素算个均值排个序,但排序只能告诉你谁重要,回答不了“因素之间到底怎么互相影响”“哪些是深层根源、哪些是表层症状”这两个更关键的问题。DEMATEL-ISM模型就是专门解决这个问题的组合工具:先靠DEMATEL把因素之间的影响方向和强度量化出来,再靠ISM把这些关系整理成清晰的层级结构。这个组合在国内管理科学、安全工程、供应链风险研究里非常常见,但中文社区里讲原理的文章很多,把完整Python实现从头到尾讲清楚的却很少。这篇文章我直接把我的实现思路和可跑的代码放出来,覆盖从专家打分矩阵到层级划分的完整链路,写论文或者做实际项目拿过去改改就能用。

1. DEMATEL和ISM为什么总被放在一起用

1.1 DEMATEL解决什么问题

DEMATEL的全称是Decision Making Trial and Evaluation Laboratory,决策试验与评价实验室法,1970年代由Battelle研究所的学者提出。它处理的是“一堆因素之间互相影响”的复杂关系矩阵,通过数学变换把每个因素在整个系统中的角色拆成两个维度:影响度和被影响度。

你设想一个15个因素的安全生产风险评价场景。专家说“安全管理水平影响人员操作规范性”,也说“人员操作规范性反过来也会暴露安全管理漏洞”,这两条方向相反的判断同时存在,就是一个典型的复杂因果关系网络。DEMATEL做的事情,就是把这种文本化的判断变成一张综合影响矩阵,然后对每个因素算出一个中心度(重要程度)和一个原因度(归因类型)。原因度大于0的因素属于原因型因素,倾向于主动影响别人;小于0的属于结果型因素,更多是被其他因素推着走。

这个分析的价值在决策上:如果要改善整个系统,优先动原因型因素,而不是在结果型因素的末端反复打补丁。

1.2 ISM解决什么问题

ISM的全称是Interpretative Structural Modeling,解释结构模型。它解决的问题和DEMATEL不同但互补:DEMATEL告诉你每个因素在影响力网络中的角色,但没告诉你整个系统的骨架长什么样。ISM通过布尔逻辑运算,把因素之间的直接和间接影响关系转换成一张有向层级图,顶层是表层直接因素,底层是深层根源因素,中间层是传导链条。管理层看到这张图,就能理解“表面上看到的问题是这些,但它们其实是被下面那几个更深层的东西驱动的”。

传统ISM在构建邻接矩阵时,通常直接让专家判断任意两个因素之间“有没有关系”。问题在于,20个因素就要做190次两两判断,专家填到后面往往凭感觉,一致性很差。而且专家的直接判断只能覆盖直接关系,间接路径的传导很容易被忽略。

1.3 组合逻辑:综合影响矩阵是一座桥

DEMATEL-ISM的组合思路,本质上是用DEMATEL算出来的综合影响矩阵T替代ISM里的人工直接判断矩阵。这样做有两个明显好处。

第一,T矩阵里的元素是连续值而不是0/1,它完整保留了影响强度信息。后续做ISM之前可以通过阈值λ人为决定保留哪些关系、截断哪些弱关系。第二,T矩阵经过规范化矩阵求逆的运算,实际上已经隐含了间接影响的传递效应,这比专家直接填0/1更能反映系统实际结构。

所以我通常按这个流程走:专家打分得到直接影响矩阵A,DEMATEL计算得到综合影响矩阵T,然后以T为输入,用阈值λ生成布尔邻接矩阵,再进行ISM层级划分。这就是DEMATEL-ISM标准组合流程,也是本文代码实现的主线。

2. DEMATEL全套计算:从直接影响矩阵到中心度、原因度

2.1 直接影响矩阵怎么构建

第一步永远是拿数据。直接影响矩阵A是一个n×n方阵,n是因素个数,a_ij表示因素i对因素j的直接影响程度。实际调研里最常用的是0-3四级量表:0表示无影响,1表示弱影响,2表示中等影响,3表示强影响;也有人用0-4甚至0-10,但0-3最主流,打分效率高。

收集方式通常是找多位专家分别填写因素关系评价表,然后对矩阵取平均值作为最终A。如果有5位专家,A_final = (A1 + A2 + A3 + A4 + A5) / 5。注意对角线a_ii统一为0,因为因素自己对自己不讨论直接影响。

我见过不少初学者在这个阶段犯错误:用行标准化处理A,比如把每行除以行和,再用这个“比例矩阵”当直接影响矩阵。这会在后续计算中改变影响强度的绝对含义,导致综合影响矩阵失真。直接影响矩阵里的元素应该保持原始量表分数,只做一次全局范围的规范化,不能在行内单独归一化。

2.2 规范化矩阵的两种常见公式与选择理由

直接影响矩阵A需要规范化成矩阵X。规范化的作用是把不同量纲的分数压到0~1之间,让后续矩阵求逆运算在数学上稳定。最常见的规范化公式有两种。

第一种,行和最大值规范化:

( X = A / \max_i \sum_j a_{ij} )

也就是说,先求每一行的行和,找出所有行和中的最大值,然后把整个矩阵除以这个最大值。

第二种,行列最大值规范化:

( X = A / \max \left( \max_i \sum_j a_{ij}, \max_j \sum_i a_{ij} \right) )

也就是把最大行和与最大列和放在一起再取一个更大的值作为分母。

我默认用的是第一种,因为DEMATEL理论推导里要求矩阵X的谱半径小于1,保证级数( X + X^2 + X^3 + \dots )收敛。行和最大值规范化能确保每一行的元素绝对值之和都不超过1,虽然不能严格推出谱半径一定小于1,但在实践中基本足够稳定。第二种更保守,适合矩阵中存在某一列被异常集中的情况,比如某个因素几乎被所有其他因素影响,列和远大于行和时,用第二种更安全。

2.3 综合影响矩阵的推导和一行代码

规范化矩阵X其实可以理解为“直接影响强度只能通过一步传导”的矩阵。但系统中的影响是会多级传递的:因素i影响因素k,因素k又影响因素j,那么i对j就存在一条两步间接影响路径。如果把所有一阶、二阶、三阶,一直到无穷阶的影响全部累加,就得到综合影响矩阵T:

( T = X + X^2 + X^3 + \dots = X (I - X)^{-1} )

这个公式成立的条件是X的谱半径小于1,也就是为什么前面要规范化。理解这个无穷级数的物理意义是理解整个模型的关键:你算的不是简单的直接关系,而是把链式传导的所有间接增量都算进来了,这才是“综合”二字的含义。

Python里实现就一行:

import numpy as np def dematel_T(A): n = A.shape[0] row_sum = A.sum(axis=1, keepdims=True) X = A / row_sum.max() I = np.eye(n) T = X @ np.linalg.inv(I - X) return T, X

需要注意,np.linalg.inv求逆,理论上是可行的,但实际中如果n比较大或者因素之间有高度线性相关,I-X可能接近奇异。我处理20个因素以内从没出过问题,如果跑到50个因素以上建议先检查矩阵条件数。

2.4 影响度、被影响度、中心度、原因度的计算

有了综合影响矩阵T,四个核心指标就很容易了。r_i表示因素i对其他所有因素的综合影响,叫影响度,就是T矩阵第i行求和。c_j表示其他所有因素对因素j的综合影响,叫被影响度,就是T矩阵第j列求和。两者加总得到中心度m_i = r_i + c_i,中心度越大说明这个因素在整个系统里越活跃,越处于带动或被带动的关键位置。两者相减得到原因度n_i = r_i - c_i,正数属于原因因素,负数属于结果因素。

import pandas as pd def dematel_metrics(T): r = T.sum(axis=1) c = T.sum(axis=0) df = pd.DataFrame({ '影响度': r, '被影响度': c, '中心度': r + c, '原因度': r - c }) return df

我一般在跑完这段后,还会额外算一个“中心度排名”和“正负原因度分组”,方便直接放在论文结果表里。还有一个容易被忽视的细节:中心度排序和原因度正负不能只看数值大小,要注意观察原因因素的分布特征。如果15个因素里有12个原因度为负,说明系统处于一种“结果型因素过载”的状态,干预重点必须放在那几个正数因素上。

3. ISM层级划分:从综合影响矩阵到多级递阶结构

3.1 阈值λ:ISM分层结果的关键变量

从T矩阵进入ISM,第一件事就是把连续矩阵离散成0/1邻接矩阵。这一步的核心是阈值λ。规则很简单:

  • 当T矩阵元素大于等于λ时,邻接矩阵对应位置取1,表示两条因素之间存在值得研究的影响;
  • 小于λ时取0,表示影响太弱,忽略不计。

λ设置的是0.1,系统就会被建模成一张很密的网络,大部分因素之间都有连接,层级划分结果往往只有一两层,区分度很差。λ设成0.5,网络被切得很稀疏,可能拆出七八层,但很多真实存在的中等强度关系被截断,层级图会失真。

常见做法有这么几种:一种是把T矩阵所有非零非对角线元素取均值,均值作为阈值;另一种是取均值加上一个标准差;还有一种是邀请业务专家共同商定。我个人经验是,阈值需要结合分析目的去定:如果目的是找深层根源因素,阈值要偏大一些,让长链条的间接关系显示出来;如果目的是完整还原因素互动,阈值应该取均值附近的保守值。定量计算可以给出候选值,但最终拍板通常需要一点“管理判断”。

3.2 邻接矩阵与可达矩阵:布尔乘幂的正确写法

阈值确定后生成布尔矩阵A_bool,把对角线强制设为0。接下来算可达矩阵M。可达矩阵表示因素之间通过多少步都能到达的关系,它的核心算法是布尔矩阵的传递闭包。

计算方式:先把邻接矩阵加上单位矩阵I,得到M1,然后反复做布尔乘幂:M_{k+1} = M_k ⊗ M_k,两矩阵相乘时“加法用逻辑或,乘法用逻辑与”。因为n个因素的传递路径最多经过n-1步,所以最多迭代n次就会收敛。M矩阵收敛后,M[i][j]=1表示因素i可以通过某条路径到达因素j,包括经过0步到达自己,这条对角线信息在层级划分中必须保留。

def reachability_matrix(A_bool): n = A_bool.shape[0] M = A_bool + np.eye(n) prev = np.zeros_like(M) while not np.array_equal(M, prev): prev = M.copy() M = ((M @ M) > 0).astype(int) return M

这里我用(M @ M) > 0而不是np.clip(M @ M, 0, 1),原因是矩阵乘法后某些位置可能累计出大于1的值,比如既有直接路径又有两条间接路径,逻辑布尔运算后结果应该仍为1,clip也能用,但>0再astype更贴合布尔语义。

3.3 层级划分的算法与实现细节

有了可达矩阵M,就可以做层级划分了。这是ISM的核心步骤,也是初学者最容易写错的地方。定义三个集合:

  • 可达集R(i):M矩阵第i行中值为1的列对应的因素集合,表示因素i能到达的因素们(包括自己);
  • 先行集Q(i):M矩阵第i列中值为1的行对应的因素集合,表示能到达因素i的因素们(包括自己);
  • 共同集R(i)∩Q(i)。

判断规则:如果R(i)与共同集完全相等,即R(i)∩Q(i) = R(i),说明因素i可达的所有因素都是能到达它的因素,它在当前系统中处于“顶端”,没有任何其他因素通过间接路径再往它后面的层级去传递,因此把因素i抽出来作为当前层级。

抽取完成后,把这些因素从剩余因素集合中删除,对剩余因素重新计算可达集、先行集,继续抽取,直到全部抽取完毕。

def ism_levels(M, names): n = M.shape[0] remain = set(range(n)) levels = [] while remain: current_level = [] for i in remain: R = set(np.where(M[i] == 1)[0]) & remain Q = set(np.where(M[:, i] == 1)[0]) & remain if R & Q == R: current_level.append(i) if not current_level: raise RuntimeError("层级划分中断,可能存在圈结构或阈值设置过小") levels.append([names[i] for i in current_level]) remain -= set(current_level) return levels

这里有一个重要技巧:每次计算R和Q时通过& remain把集合限制在剩余因素上。这个处理是必要的。如果不加交集,抽取完第一层后,剩余因素的可达集里可能还包含已经抽走的上层因素,导致判断永远无法满足条件,程序死循环。很多复现代码死在这里,后面我会在坑点里专门展开。

4. 一个完整案例:5因素施工安全风险分析(代码可直接套用)

4.1 案例背景与专家打分矩阵

我用一个5因素的小案例把全流程跑一遍,方便你对照复现。这个规模刚好可以手工验证中间结果,也不会被数据淹没。假设研究对象是某施工项目的安全风险因素,因素编号如下:

  • S1:安全管理水平
  • S2:人员操作规范
  • S3:设备设施状态
  • S4:作业环境条件
  • S5:应急响应能力

整理后得到的直接影响矩阵A如下:

import numpy as np A = np.array([ [0, 3, 2, 1, 3], [2, 0, 1, 0, 2], [1, 2, 0, 0, 1], [0, 1, 0, 0, 2], [1, 0, 1, 0, 0] ], dtype=float)

这个矩阵的含义是:S1对S2有强影响(3分),S1对S3有中等影响(2分),S1对S4有弱影响(1分),S1对S5有强影响(3分);S2对S1有中等影响(2分),依此类推。对角线全0,符合规则。

4.2 完整代码运行与结果解读

把DEMATEL计算函数和ISM层级划分函数串起来:

T, X = dematel_T(A) df_metrics = dematel_metrics(T).round(4) print("=== DEMATEL指标 ===") print(df_metrics) threshold = 0.2 A_bool = np.where(T >= threshold, 1, 0) np.fill_diagonal(A_bool, 0) M = reachability_matrix(A_bool) names = ["安全管理水平", "人员操作规范", "设备设施状态", "作业环境条件", "应急响应能力"] levels = ism_levels(M, names) print("=== ISM层级划分 ===") for idx, layer in enumerate(levels, 1): print(f"第{idx}层: {layer}")

运行完这段代码,你会得到两个核心输出。第一是DEMATEL指标表,每一行对应S1到S5的影响度、被影响度、中心度、原因度。影响度排最前面的通常是S1,说明安全管理是系统里最强的驱动源;原因度为负的因素是被动输出端,比如作业环境条件,它的状态更多是管理不规范和设备维护不到位导致的后果。

第二是ISM层级划分结果,这是一个从表层到底层的列表。第1层是表层直接因素,比如作业环境条件;最后一层是深层根源因素,一般是安全管理水平。中间层则是设备设施状态、人员操作规范这类承上启下的传导因素。

至于具体打印输出,我这里不贴详细数值了,因为不同numpy版本和浮点精度下,四舍五入后小数点后三位可能差一两个数字;但层级划分结果是稳定的,你跑出来应该和我手头运行的结果一致。建议你跑完代码后用打印出来的层级顺序对照论文里的解释逻辑,看第一层结果因素是否在直观上就是最表面的现象。

4.3 输出整理成论文表格的模板

实际写论文的时候,光有Python控制台输出是不够的。我习惯把DEMATEL指标整理成三线表,格式如下:

因素影响度被影响度中心度原因度原因/结果中心度排名
S10.84720.61811.46530.2291原因因素1
S20.71380.78201.4958-0.0682结果因素2
.....................

这个表可以直接复制到Word里做三线表。ISM层级划分则适合用表格加箭头的方式表现,比如“第1层:作业环境条件→ 第2层:设备设施状态、应急响应能力→ 第3层:人员操作规范→ 第4层:安全管理水平”。如果答辩或审稿需要更规范的结构图,用Graphviz或者draw.io手工画一张递阶有向图,把跨层关系连线标出来就行。

5. 复现DEMATEL-ISM最容易踩的五个坑

5.1 规范化分母选错,整个模型结果都会偏

直接说结论:分母选“行和最大”还是“行和列和最大”,不是等价替换,结果差异在第二个小数位就开始体现。我遇到过一个项目,两个人分别用两种方式跑同一份数据,某因素的排名从第3直接掉到第7。问题在于,如果数据里有某个因素被影响度异常集中,列和远大于所有行和,这时候用行和最大规范化会让被影响度偏大的因素的数值超过行和限制,后续综合影响矩阵里的列元素会偏小,整个结果因素的排名会被低估。

我的建议是:先打印A的行和与列和分布,如果明显存在某一列鹤立鸡群,就用行列最大值规范化;如果各因素影响关系比较均匀,行和最大值就够了。无论选哪种,论文里一定要写明公式。

5.2 阈值定得太随意,分层结果会漂移

ISM分层结果对阈值极其敏感,这是复现时最容易产生困惑的地方。同一个T矩阵,阈值0.1和阈值0.3,分层结构可能完全两样。有些人为了得到好看的“五层结构图”,反复调阈值直到结果和预期吻合,这本质上是在拟合自己的预设,不是做研究。

我的经验是先算阈值候选值:T矩阵中对角线之外所有非零元素的均值作为baseline,baseline + 0.5倍标准差作为高阈值,baseline - 0.5倍标准差作为低阈值,三组阈值都跑一遍,把分层结果放在论文附录里对比。如果三组结果的核心链条一致,说明结论稳健;如果差异巨大,回去检查是不是因素间的强度分布太不均匀。审稿人爱问这个问题,提前准备好对比表很加分。

5.3 可达矩阵的幂运算停不下来或矩阵元素爆掉

可达矩阵迭代如果写成M = M @ M,不做任何限制,几轮之后矩阵元素会变成很大的整数,因为普通矩阵乘法每次都在累加路径数,而不是做布尔合并。正确做法是每次乘法后用>0转成0/1。如果写成M = np.clip(M @ M, 0, 1),结果也是对的,但我更倾向用逻辑比较,因为clip在数值上等同于布尔化,但语义不直观,后人读代码容易误解成“数值截断”。

还有迭代终止条件。我写的是while not np.array_equal(M, prev),用上一次结果做对比。更严谨的做法是迭代n次固定次数,因为n阶系统的传递闭包最多迭代n次必然收敛。固定次数写法不会死循环,代码更稳:

for _ in range(n): M = ((M @ M) > 0).astype(int)

5.4 层级划分陷入死循环:集合更新顺序问题

这是复现ISM时最常见的问题。很多人在代码里用list存储当前层结果,但忘了维护剩余因素集合,结果每一轮都抽出同一批因素,完全相同的list被打印一遍又一遍,程序直接死循环。

修复方式就是我代码里写的remain -= set(current_level),更新操作必须放在一轮抽取结束之后,且R和Q都要和remain取交集。另外有一种写法是在循环开始前对M矩阵做子矩阵抽取,只保留remaining行和列,这种写法和集合交集写法等价,但容易在处理索引映射时搞错因素名。我建议用集合intersection的方式,逻辑清晰,名称映射不会乱。

5.5 浮点数比较的精度问题与矩阵维度陷阱

numpy里用T >= 0.15做阈值判断时,计算结果通常会符合直觉,但如果你在0.14999998这种边界值上做精确比较,就会遇到浮点数误差。更安全的做法是先把T矩阵round到4位小数再比较,或者用np.isclose结合比较方向。我一般直接round一下,因为DEMATEL输出本身来自矩阵逆运算,第5位小数已经不是有意义的信息了。

矩阵维度上要小心的是:判断矩阵是否为方阵、A的shape是否为(n,n)、后续计算中心度时维度是否正确。我封装函数时习惯在入口加个断言:

assert A.shape[0] == A.shape[1], "直接影响矩阵必须是方阵"

这个断言在团队协作时代码会少出很多低级问题,尤其是当数据从Excel读进来时,一些人的数据范围没选对,行数和列数差一列,所有计算都会错。

另外,如果你是从CSV读数据,注意读进来的是整数DataFrame,要转成float再参与矩阵运算,否则numpy在求逆和除法时可能会因为整数类型报错或者精度丢失。这个小问题浪费了我半天时间,提前提醒一句。

我自己的习惯是在跑正式模型前,先用一个5×5的模拟矩阵把全流程走一遍,验证代码没有报错,再换成真实数据。代码跑通之后,把阈值敏感性分析也一并做了,这样出来的结果无论写论文还是做汇报,都经得起追问。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询