做影响因素分析的时候最怕的就是手上一堆专家打分数据算完描述性统计之后不知道怎么往下走。你当然可以给每个因素算个均值排个序但排序只能告诉你谁重要回答不了“因素之间到底怎么互相影响”“哪些是深层根源、哪些是表层症状”这两个更关键的问题。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更能反映系统实际结构。所以我通常按这个流程走专家打分得到直接影响矩阵ADEMATEL计算得到综合影响矩阵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(axis1, keepdimsTrue) 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(axis1) c T.sum(axis0) 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_{k1} 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的值比如既有直接路径又有两条间接路径逻辑布尔运算后结果应该仍为1clip也能用但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] ], dtypefloat)这个矩阵的含义是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原因因素1S20.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矩阵中对角线之外所有非零元素的均值作为baselinebaseline 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的模拟矩阵把全流程走一遍验证代码没有报错再换成真实数据。代码跑通之后把阈值敏感性分析也一并做了这样出来的结果无论写论文还是做汇报都经得起追问。