# 02.3 多准则决策 ## 核心问题 > 当多个目标相互冲突时，如何做出"最优"决策？ > 专家的判断经验如何转化为可计算的权重？ --- ## 概念讲解 ### 什么是多准则决策分析 (MCDA) 多准则决策分析是一种在多个、通常是冲突的准则下评估和选择替代方案的方法论。 ``` ┌─────────────────────────────────────────────────────────────┐ │ 多准则决策问题的结构 │ ├─────────────────────────────────────────────────────────────┤ │ │ │ ┌─────────────────────────────────────────────────────┐ │ │ │ 决策目标 │ │ │ │ "选择最适合生态修复的区域" │ │ │ └─────────────────────────────────────────────────────┘ │ │ │ │ │ ▼ │ │ ┌─────────────────────────────────────────────────────┐ │ │ │ 评估准则 │ │ │ │ ┌─────────────┐ ┌─────────────┐ ┌─────────────┐ │ │ │ │ │ 生态重要性 │ │ 实施可行性 │ │ 成本效益 │ │ │ │ │ │ (Weight) │ │ (Weight) │ │ (Weight) │ │ │ │ │ │ 0.4 │ │ 0.3 │ │ 0.3 │ │ │ │ │ └─────────────┘ └─────────────┘ └─────────────┘ │ │ │ └─────────────────────────────────────────────────────┘ │ │ │ │ │ ▼ │ │ ┌─────────────────────────────────────────────────────┐ │ │ │ 替代方案 │ │ │ │ ┌──────┐ ┌──────┐ ┌──────┐ ┌──────┐ │ │ │ │ │ 区域A │ │ 区域B │ │ 区域C │ │ 区域D │ ... │ │ │ │ └──────┘ └──────┘ └──────┘ └──────┘ │ │ │ └─────────────────────────────────────────────────────┘ │ │ │ │ │ ▼ │ │ ┌─────────────────────────────────────────────────────┐ │ │ │ 决策结果 │ │ │ │ 综合得分 + 排名 + 稳健性分析 │ │ │ └─────────────────────────────────────────────────────┘ │ │ │ └─────────────────────────────────────────────────────────────┘ ``` ### MCDA的核心组成部分 | 组成部分 | 描述 | 空间应用示例 | |---------|------|-------------| | **准则 (Criteria)** | 评估标准，反映决策目标 | 生境质量、连通性、建设成本 | | **权重 (Weights)** | 准则的相对重要性 | 生态重要性0.5，经济成本0.3，社会因素0.2 | | **得分 (Scores)** | 各方案在各准则下的表现 | 每个栅格的生境适宜性指数 | | **标准化 (Normalization)** | 将不同单位转换为可比尺度 | 0-1标准化、排名转换 | | **集结规则 (Aggregation)** | 合并多准则得分的方法 | 加权求和、加权乘积、TOPSIS | ### 常用的MCDA方法 ``` ┌─────────────────────────────────────────────────────────────┐ │ MCDA方法分类 │ ├─────────────────────────────────────────────────────────────┤ │ │ │ 1. 加权线性组合 (WLC) / 简单加权法 │ │ ┌─────────────────────────────────────────────────┐ │ │ │ Score = Σ(weight_i × score_i) │ │ │ │ │ │ │ │ 优点：简单、直观、易于理解 │ │ │ │ 缺点：允许补偿（一个准则的差可被另一个优弥补） │ │ │ └─────────────────────────────────────────────────┘ │ │ │ │ 2. 层次分析法 (AHP) │ │ ┌─────────────────────────────────────────────────┐ │ │ │ 通过成对比较确定权重 │ │ │ │ │ │ │ │ 优点：结构化、一致性检验 │ │ │ │ 缺点：比较次数多(n(n-1)/2)、可能存在不一致 │ │ │ └─────────────────────────────────────────────────┘ │ │ │ │ 3. TOPSIS (逼近理想解排序法) │ │ ┌─────────────────────────────────────────────────┐ │ │ │ 选择距正理想解最近、负理想解最远的方案 │ │ │ │ │ │ │ │ 优点：考虑方案与理想解的相对距离 │ │ │ │ 缺点：对权重敏感 │ │ │ └─────────────────────────────────────────────────┘ │ │ │ │ 4. OWA (有序加权平均) │ │ ┌─────────────────────────────────────────────────┐ │ │ │ 允许控制"风险态度"(ORness) │ │ │ │ │ │ │ │ 优点：灵活的风险偏好建模 │ │ │ │ 缺点：需要确定orness参数 │ │ │ └─────────────────────────────────────────────────┘ │ │ │ └─────────────────────────────────────────────────────────────┘ ``` ### 标准化方法 不同准则有不同的量纲，需要标准化： ``` ┌─────────────────────────────────────────────────────────────┐ │ 标准化方法 │ ├─────────────────────────────────────────────────────────────┤ │ │ │ 1. 最大-最小标准化 (Min-Max) │ │ ┌─────────────────────────────────────────────────┐ │ │ │ x_norm = (x - min) / (max - min) │ │ │ │ 适用：有明确最大最小值，线性关系 │ │ │ └─────────────────────────────────────────────────┘ │ │ │ │ 2. Z-score标准化 │ │ ┌─────────────────────────────────────────────────┐ │ │ │ x_norm = (x - mean) / std │ │ │ │ 适用：正态分布数据，异常值敏感 │ │ │ └─────────────────────────────────────────────────┘ │ │ │ │ 3. 分位数转换 (Quantile) │ │ ┌─────────────────────────────────────────────────┐ │ │ │ x_norm = percentile_rank(x) │ │ │ │ 适用：分布未知，需要稳健性 │ │ │ └─────────────────────────────────────────────────┘ │ │ │ │ 4. 目标导向标准化 │ │ ┌─────────────────────────────────────────────────┐ │ │ │ 效益型: x_norm = x / target │ │ │ │ 成本型: x_norm = target / x │ │ │ │ 适用：有明确目标值 │ │ │ └─────────────────────────────────────────────────┘ │ │ │ └─────────────────────────────────────────────────────────────┘ ``` --- ## 设计原理 ### AHP层次分析法 AHP通过成对比较确定权重，具有结构化和一致性检验的优点： ```python import numpy as np from typing import List, Dict, Tuple, Optional class AHPAnalyzer: """ 层次分析法 (Analytic Hierarchy Process) 实现 核心思想：通过成对比较构建判断矩阵，计算权重 """ # Saaty标度：1-9及其倒数 SAATY_SCALE = { 1: "同等重要", 3: "稍微重要", 5: "明显重要", 7: "强烈重要", 9: "极端重要", 2: "介于1和3之间", 4: "介于3和5之间", 6: "介于5和7之间", 8: "介于7和9之间" } def __init__(self, criteria: List[str]): """ Args: criteria: 准则列表 """ self.criteria = criteria self.n = len(criteria) self.comparison_matrix = None self.weights = None self.consistency_ratio = None def set_comparison_matrix(self, matrix: np.ndarray): """ 直接设置比较矩阵 Args: matrix: n×n的判断矩阵，matrix[i,j]表示i相对于j的重要性 """ if matrix.shape != (self.n, self.n): raise ValueError(f"矩阵大小应为{self.n}×{self.n}") # 确保对角线为1 np.fill_diagonal(matrix, 1) self.comparison_matrix = matrix def build_from_pairs(self, pairs: Dict[Tuple[str, str], float]): """ 从成对比较构建矩阵 Args: pairs: {(criterion_i, criterion_j): value} 表示i相对于j的重要性 如果j相对于i，则value应为1/value """ self.comparison_matrix = np.eye(self.n) criterion_to_idx = {c: i for i, c in enumerate(self.criteria)} for (ci, cj), value in pairs.items(): i, j = criterion_to_idx[ci], criterion_to_idx[cj] self.comparison_matrix[i, j] = value self.comparison_matrix[j, i] = 1.0 / value def compute_weights(self, method: str = 'eigenvector') -> np.ndarray: """ 计算权重 Args: method: 'eigenvector'(特征向量法) 或 'geometric'(几何平均法) Returns: 权重向量 """ if self.comparison_matrix is None: raise ValueError("请先设置比较矩阵") if method == 'eigenvector': # 特征向量法：求最大特征值对应的特征向量 eigenvalues, eigenvectors = np.linalg.eig(self.comparison_matrix) max_idx = np.argmax(eigenvalues.real) weights = eigenvectors[:, max_idx].real # 归一化 weights = weights / weights.sum() elif method == 'geometric': # 几何平均法 weights = np.exp(np.log(self.comparison_matrix).mean(axis=1)) weights = weights / weights.sum() else: raise ValueError(f"未知方法: {method}") self.weights = weights return weights def compute_consistency(self) -> Dict[str, float]: """ 计算一致性指标 Returns: 包含一致性相关指标的字典 """ if self.comparison_matrix is None or self.weights is None: raise ValueError("请先设置比较矩阵并计算权重") # 计算最大特征值 weighted_sum = self.comparison_matrix @ self.weights lambda_max = (weighted_sum / self.weights).mean() # 一致性指标 CI n = self.n ci = (lambda_max - n) / (n - 1) if n > 1 else 0 # 随机一致性指标 RI (Saaty给出的标准值) ri_table = {1: 0, 2: 0, 3: 0.58, 4: 0.90, 5: 1.12, 6: 1.24, 7: 1.32, 8: 1.41, 9: 1.45, 10: 1.49} ri = ri_table.get(n, 1.49) # 一致性比率 CR cr = ci / ri if ri > 0 else 0 self.consistency_ratio = cr return { 'lambda_max': lambda_max, 'CI': ci, 'RI': ri, 'CR': cr, 'consistent': cr < 0.1 } def get_weights_dict(self) -> Dict[str, float]: """返回准则-权重字典""" if self.weights is None: self.compute_weights() return {c: w for c, w in zip(self.criteria, self.weights)} # 使用示例 def example_ahp(): """AHP使用示例：生态源地选址""" criteria = ['生态价值', '连通性', '实施成本', '社会接受度'] ahp = AHPAnalyzer(criteria) # 设置成对比较（示例数据） pairs = { ('生态价值', '连通性'): 2, # 生态价值比连通性稍微重要 ('生态价值', '实施成本'): 5, # 生态价值比成本明显重要 ('生态价值', '社会接受度'): 3, # 生态价值比社会接受度稍微重要 ('连通性', '实施成本'): 3, # 连通性比成本稍微重要 ('连通性', '社会接受度'): 1, # 连通性与社会接受度同等重要 ('实施成本', '社会接受度'): 1/2, # 社会接受度比成本稍微重要 } ahp.build_from_pairs(pairs) weights = ahp.compute_weights() consistency = ahp.compute_consistency() print("准则权重:") for c, w in ahp.get_weights_dict().items(): print(f" {c}: {w:.3f}") print(f"\n一致性比率: {consistency['CR']:.3f}") print(f"一致性{'合格' if consistency['consistent'] else '不合格'}") return ahp if __name__ == "__main__": example_ahp() ``` ### TOPSIS方法 TOPSIS通过计算与理想解的距离进行排序： ```python import numpy as np from typing import List, Dict, Callable, Optional class TOPSISAnalyzer: """ TOPSIS (逼近理想解排序法) 实现 核心思想：选择距离正理想解最近、负理想解最远的方案 """ def __init__(self, criteria: List[str], directions: Optional[List[str]] = None): """ Args: criteria: 准则列表 directions: 每个准则的方向，'benefit'(效益型)或'cost'(成本型) """ self.criteria = criteria self.directions = directions or ['benefit'] * len(criteria) self.weights = None self.normalized_matrix = None self.ideal_positive = None self.ideal_negative = None self.scores = None def set_weights(self, weights: np.ndarray): """设置权重""" if len(weights) != len(self.criteria): raise ValueError("权重数量与准则数量不匹配") self.weights = np.array(weights) / np.sum(weights) def normalize(self, decision_matrix: np.ndarray) -> np.ndarray: """ 向量标准化 Args: decision_matrix: m×n矩阵，m个方案，n个准则 Returns: 标准化后的矩阵 """ # 向量标准化：每个列向量除以其范数 norm = np.sqrt((decision_matrix ** 2).sum(axis=0)) normalized = decision_matrix / norm # 处理除零 normalized = np.nan_to_num(normalized) self.normalized_matrix = normalized return normalized def compute_ideal_solutions(self, weighted_matrix: np.ndarray): """ 计算正理想解和负理想解 Args: weighted_matrix: 加权标准化决策矩阵 """ n = weighted_matrix.shape[1] ideal_positive = np.zeros(n) ideal_negative = np.zeros(n) for j in range(n): if self.directions[j] == 'benefit': ideal_positive[j] = weighted_matrix[:, j].max() ideal_negative[j] = weighted_matrix[:, j].min() else: # cost ideal_positive[j] = weighted_matrix[:, j].min() ideal_negative[j] = weighted_matrix[:, j].max() self.ideal_positive = ideal_positive self.ideal_negative = ideal_negative return ideal_positive, ideal_negative def compute_scores(self, decision_matrix: np.ndarray) -> np.ndarray: """ 计算TOPSIS得分 Args: decision_matrix: m×n决策矩阵 Returns: 得分向量 (0-1之间，越大越好) """ # 标准化 normalized = self.normalize(decision_matrix) # 加权 if self.weights is None: self.weights = np.ones(len(self.criteria)) / len(self.criteria) weighted = normalized * self.weights # 计算理想解 self.compute_ideal_solutions(weighted) # 计算距离 m = weighted.shape[0] d_positive = np.zeros(m) d_negative = np.zeros(m) for i in range(m): d_positive[i] = np.sqrt( ((weighted[i] - self.ideal_positive) ** 2).sum() ) d_negative[i] = np.sqrt( ((weighted[i] - self.ideal_negative) ** 2).sum() ) # 计算相对贴近度 scores = d_negative / (d_positive + d_negative) scores = np.nan_to_num(scores) self.scores = scores return scores def rank(self, decision_matrix: np.ndarray) -> List[int]: """ 返回方案排名 Returns: 排名索引列表（从最优到最差） """ scores = self.compute_scores(decision_matrix) return np.argsort(-scores).tolist() # 使用示例 def example_topsis(): """TOPSIS使用示例：生态修复区域选择""" # 准则 criteria = ['生境适宜性', '连通性指数', '实施成本', '社会效益'] directions = ['benefit', 'benefit', 'cost', 'benefit'] # 权重 weights = [0.35, 0.25, 0.20, 0.20] # 决策矩阵：5个候选区域在4个准则下的得分 # 注意：成本已经是负向的（数值越小越好） decision_matrix = np.array([ [0.85, 0.72, 500, 0.65], # 区域A [0.78, 0.85, 450, 0.70], # 区域B [0.92, 0.68, 600, 0.55], # 区域C [0.70, 0.90, 400, 0.80], # 区域D [0.88, 0.75, 550, 0.60], # 区域E ]) # 对成本型准则进行预处理（转为效益型） # 成本越高，得分越低 cost_col = decision_matrix[:, 2].copy() decision_matrix[:, 2] = cost_col.max() - cost_col topsis = TOPSISAnalyzer(criteria, directions) topsis.set_weights(weights) scores = topsis.compute_scores(decision_matrix) ranking = topsis.rank(decision_matrix) print("区域排名 (TOPSIS):") region_names = ['A', 'B', 'C', 'D', 'E'] for rank, idx in enumerate(ranking, 1): print(f" 第{rank}名: 区域{region_names[idx]} (得分: {scores[idx]:.3f})") return topsis if __name__ == "__main__": example_topsis() ``` ### 敏感性分析 敏感性分析是MCDA的重要组成部分，评估权重变化对决策结果的影响： ```python import numpy as np from typing import List, Dict, Tuple, Callable import matplotlib.pyplot as plt class SensitivityAnalyzer: """ 权重敏感性分析 评估权重变化对决策结果的影响程度 """ def __init__(self, decision_matrix: np.ndarray, weights: np.ndarray, criteria: List[str], directions: List[str] = None): """ Args: decision_matrix: 决策矩阵 weights: 初始权重 criteria: 准则名称 directions: 准则方向 """ self.decision_matrix = decision_matrix self.base_weights = np.array(weights) / np.sum(weights) self.criteria = criteria self.directions = directions or ['benefit'] * len(criteria) def one_at_a_time(self, variation: float = 0.1) -> Dict[str, Dict]: """ 单因素敏感性分析 (OAT) 每次改变一个准则的权重，观察排名变化 Args: variation: 权重变化幅度（相对于原权重的比例） Returns: 每个准则变化时的排名变化 """ n = len(self.criteria) base_ranking = self._compute_ranking(self.base_weights) results = {} for i, criterion in enumerate(self.criteria): results[criterion] = { 'weight_variations': [], 'rankings': [] } # 增加权重 for delta in np.linspace(-variation, variation, 21): new_weights = self.base_weights.copy() # 调整第i个权重 new_weights[i] = self.base_weights[i] * (1 + delta) # 重新归一化（保持和为1） if new_weights.sum() > 0: new_weights = new_weights / new_weights.sum() else: new_weights = self.base_weights.copy() # 计算新排名 ranking = self._compute_ranking(new_weights) results[criterion]['weight_variations'].append(delta) results[criterion]['rankings'].append(ranking) return results def tornado_analysis(self, variation: float = 0.2) -> Dict: """ 龙卷风图分析 计算每个准则权重变化对最优方案得分的影响 Args: variation: 权重变化幅度 Returns: 龙卷风图数据 """ from topsis import TOPSISAnalyzer base_scores = self._compute_scores(self.base_weights) base_best_score = base_scores.max() impacts = [] for i, criterion in enumerate(self.criteria): impact_high = 0 impact_low = 0 # 增加权重 for sign in [1, -1]: new_weights = self.base_weights.copy() new_weights[i] = self.base_weights[i] * (1 + sign * variation) new_weights = new_weights / new_weights.sum() new_scores = self._compute_scores(new_weights) new_best = new_scores.max() change = new_best - base_best_score if sign > 0: impact_high = change else: impact_low = change impacts.append({ 'criterion': criterion, 'high': impact_high, 'low': impact_low, 'range': impact_high - impact_low }) # 按影响范围排序 impacts.sort(key=lambda x: abs(x['range']), reverse=True) return { 'base_score': base_best_score, 'impacts': impacts } def monte_carlo(self, n_simulations: int = 1000, weight_std: float = 0.1) -> Dict: """ 蒙特卡洛敏感性分析 随机扰动权重，观察结果的统计分布 Args: n_simulations: 模拟次数 weight_std: 权重扰动的标准差 Returns: 统计结果 """ n_alternatives = self.decision_matrix.shape[0] rank_counts = np.zeros(n_alternatives, dtype=int) scores_history = [] for _ in range(n_simulations): # 生成随机权重 random_weights = np.random.normal( self.base_weights, weight_std ) # 确保非负并归一化 random_weights = np.maximum(random_weights, 0) random_weights = random_weights / random_weights.sum() # 计算排名 ranking = self._compute_ranking(random_weights) rank_counts[ranking[0]] += 1 # 统计第一名 # 记录得分 scores = self._compute_scores(random_weights) scores_history.append(scores) # 计算每个方案排第一的概率 probabilities = rank_counts / n_simulations # 得分的统计量 scores_array = np.array(scores_history) scores_stats = { 'mean': scores_array.mean(axis=0), 'std': scores_array.std(axis=0), 'min': scores_array.min(axis=0), 'max': scores_array.max(axis=0), 'percentile_5': np.percentile(scores_array, 5, axis=0), 'percentile_95': np.percentile(scores_array, 95, axis=0) } return { 'rank_probabilities': probabilities, 'scores_stats': scores_stats } def _compute_scores(self, weights: np.ndarray) -> np.ndarray: """计算给定权重下的得分""" # 简单加权求和 normalized = self._normalize_matrix(self.decision_matrix) scores = (normalized * weights).sum(axis=1) return scores def _compute_ranking(self, weights: np.ndarray) -> List[int]: """计算给定权重下的排名""" scores = self._compute_scores(weights) return np.argsort(-scores).tolist() def _normalize_matrix(self, matrix: np.ndarray) -> np.ndarray: """标准化决策矩阵""" normalized = matrix.copy() for j, direction in enumerate(self.directions): col = matrix[:, j] if direction == 'benefit': normalized[:, j] = (col - col.min()) / (col.max() - col.min()) else: # cost normalized[:, j] = (col.max() - col) / (col.max() - col.min()) return normalized # 使用示例 def example_sensitivity(): """敏感性分析示例""" # 决策矩阵 decision_matrix = np.array([ [0.85, 0.72, 0.65], [0.78, 0.85, 0.70], [0.92, 0.68, 0.55], [0.70, 0.90, 0.80], ]) weights = [0.4, 0.35, 0.25] criteria = ['生态价值', '连通性', '可行性'] analyzer = SensitivityAnalyzer( decision_matrix, weights, criteria, directions=['benefit', 'benefit', 'benefit'] ) # 蒙特卡洛分析 mc_results = analyzer.monte_carlo(n_simulations=1000, weight_std=0.1) print("各方案排名第一的概率:") for i, prob in enumerate(mc_results['rank_probabilities']): print(f" 方案{i+1}: {prob:.1%}") # 龙卷风分析 tornado = analyzer.tornado_analysis(variation=0.2) print("\n准则重要性排序:") for impact in tornado['impacts']: print(f" {impact['criterion']}: 影响±{abs(impact['range']):.3f}") return analyzer ``` --- ## 代码示例 ### 空间多准则决策分析 ```python """ 空间多准则决策分析完整示例 """ import numpy as np import matplotlib.pyplot as plt from typing import List, Dict, Optional, Callable from dataclasses import dataclass @dataclass class MCDCriteria: """决策准则""" name: str direction: str # 'benefit' 或 'cost' weight: float raster: Optional[np.ndarray] = None # 栅格数据 normalize_func: Optional[Callable] = None # 标准化函数 class SpatialMCDA: """ 空间多准则决策分析 用于栅格数据的多准则评估 """ def __init__(self, shape: tuple, nodata: float = -9999): """ Args: shape: 栅格形状 (rows, cols) nodata: 无数据值 """ self.shape = shape self.nodata = nodata self.criteria: List[MCDCriteria] = [] self.result: Optional[np.ndarray] = None self.mask: Optional[np.ndarray] = None def add_criterion(self, name: str, direction: str, weight: float, raster: Optional[np.ndarray] = None): """ 添加准则 Args: name: 准则名称 direction: 'benefit'(越大越好) 或 'cost'(越小越好) weight: 权重 raster: 栅格数据 """ criterion = MCDCriteria( name=name, direction=direction, weight=weight, raster=np.asarray(raster) if raster is not None else None ) self.criteria.append(criterion) def set_mask(self, mask: np.ndarray): """ 设置分析掩模 Args: mask: True表示有效区域 """ self.mask = np.asarray(mask, dtype=bool) def normalize_all(self, method: str = 'minmax') -> List[np.ndarray]: """ 标准化所有准则 Args: method: 'minmax', 'zscore', 或 'quantile' Returns: 标准化后的栅格列表 """ normalized = [] for criterion in self.criteria: if criterion.raster is None: raise ValueError(f"准则 {criterion.name} 没有栅格数据") data = criterion.raster.copy() # 处理无数据值 valid_mask = data != self.nodata if method == 'minmax': # 最大-最小标准化 valid_data = data[valid_mask] min_val, max_val = valid_data.min(), valid_data.max() if criterion.direction == 'benefit': data[valid_mask] = (valid_data - min_val) / (max_val - min_val) else: # cost data[valid_mask] = (max_val - valid_data) / (max_val - min_val) elif method == 'zscore': # Z-score标准化 valid_data = data[valid_mask] mean_val, std_val = valid_data.mean(), valid_data.std() data[valid_mask] = (valid_data - mean_val) / std_val # 转换到0-1 data[valid_mask] = (data[valid_mask] - data[valid_mask].min()) / \ (data[valid_mask].max() - data[valid_mask].min()) elif method == 'quantile': # 分位数转换 valid_data = data[valid_mask] ranks = np.argsort(np.argsort(valid_data)) data[valid_mask] = ranks / (len(ranks) - 1) if criterion.direction == 'cost': data[valid_mask] = 1 - data[valid_mask] else: raise ValueError(f"未知方法: {method}") # 应用掩模 if self.mask is not None: data[~self.mask] = self.nodata normalized.append(data) return normalized def compute(self, method: str = 'wlc') -> np.ndarray: """ 计算综合评估结果 Args: method: 'wlc'(加权线性组合) 或 'geometric'(几何平均) Returns: 评估结果栅格 """ # 标准化 normalized = self.normalize_all() # 归一化权重 weights = np.array([c.weight for c in self.criteria]) weights = weights / weights.sum() # 初始化结果 result = np.zeros(self.shape) if method == 'wlc': # 加权线性组合 for i, norm_data in enumerate(normalized): valid_mask = norm_data != self.nodata result[valid_mask] += norm_data[valid_mask] * weights[i] elif method == 'geometric': # 加权几何平均 result.fill(1) for i, norm_data in enumerate(normalized): valid_mask = norm_data != self.nodata result[valid_mask] *= np.power(norm_data[valid_mask], weights[i]) else: raise ValueError(f"未知方法: {method}") # 应用掩模 if self.mask is not None: result[~self.mask] = self.nodata self.result = result return result def get_rankings(self, n_classes: int = 5) -> np.ndarray: """ 将连续结果分级 Args: n_classes: 分类数 Returns: 分级结果 (1=最低, n_classes=最高) """ if self.result is None: self.compute() result = self.result.copy() valid_mask = result != self.nodata # 分位数分级 valid_data = result[valid_mask] thresholds = np.percentile(valid_data, np.linspace(0, 100, n_classes + 1)) rankings = np.zeros_like(result, dtype=int) for i in range(n_classes): mask = (result >= thresholds[i]) & (result < thresholds[i + 1]) rankings[mask] = i + 1 rankings[result >= thresholds[-1]] = n_classes rankings[~valid_mask] = 0 return rankings def sensitivity_analysis(self, criterion_name: str, weight_variation: float = 0.2) -> np.ndarray: """ 单准则权重敏感性分析 Args: criterion_name: 要分析的准则名称 weight_variation: 权重变化范围 (±) Returns: 不同权重下的结果差异 """ base_result = self.compute() variations = [] # 找到准则索引 criterion_idx = None for i, c in enumerate(self.criteria): if c.name == criterion_name: criterion_idx = i break if criterion_idx is None: raise ValueError(f"未找到准则: {criterion_name}") # 变化权重 for delta in np.linspace(-weight_variation, weight_variation, 11): # 修改权重 original_weight = self.criteria[criterion_idx].weight self.criteria[criterion_idx].weight = original_weight * (1 + delta) # 重新计算 new_result = self.compute() variations.append(new_result.copy()) # 恢复权重 self.criteria[criterion_idx].weight = original_weight # 计算标准差作为敏感性度量 variations = np.array(variations) sensitivity = np.std(variations, axis=0) return sensitivity # 使用示例 def example_spatial_mcda(): """空间多准则决策分析示例""" # 创建示例数据 rows, cols = 100, 100 mcda = SpatialMCDA(shape=(rows, cols)) # 创建掩模（例如：排除研究区外的区域） mask = np.zeros((rows, cols), dtype=bool) mask[20:80, 20:80] = True mcda.set_mask(mask) # 添加准则 np.random.seed(42) # 准则1：生境适宜性（效益型，权重0.4） habitat = np.random.rand(rows, cols) * 0.6 + 0.2 habitat[~mask] = mcda.nodata mcda.add_criterion('生境适宜性', 'benefit', 0.4, habitat) # 准则2：连通性（效益型，权重0.3） connectivity = np.random.rand(rows, cols) * 0.7 + 0.15 connectivity[~mask] = mcda.nodata mcda.add_criterion('连通性', 'benefit', 0.3, connectivity) # 准则3：实施成本（成本型，权重0.3） cost = np.random.rand(rows, cols) * 0.8 + 0.1 cost[~mask] = mcda.nodata mcda.add_criterion('实施成本', 'cost', 0.3, cost) # 计算综合评估 result = mcda.compute(method='wlc') print(f"综合评估结果:") print(f" 有效区域均值: {result[mask].mean():.3f}") print(f" 有效区域范围: [{result[mask].min():.3f}, {result[mask].max():.3f}]") # 分级 rankings = mcda.get_rankings(n_classes=5) print(f"\n分级统计:") for i in range(1, 6): count = (rankings == i).sum() print(f" 等级{i}: {count} 个像元") # 敏感性分析 sensitivity = mcda.sensitivity_analysis('生境适宜性', weight_variation=0.3) print(f"\n敏感性分析 (生境适宜性权重±30%):") print(f" 结果最大变化: {sensitivity[mask].max():.3f}") print(f" 结果平均变化: {sensitivity[mask].mean():.3f}") return mcda if __name__ == "__main__": example_spatial_mcda() ``` --- ## 案例分析 ### ENAgent中的生态系统服务评估 ENAgent使用多准则决策分析评估生态系统的综合服务价值： ```python class EcosystemServiceAssessment: """ 生态系统服务评估模块 基于多准则决策分析评估区域生态系统服务价值 """ # 生态系统服务类型 SERVICE_TYPES = { 'provisioning': '供给服务', # 食物、淡水、木材等 'regulating': '调节服务', # 气候调节、洪水调节等 'cultural': '文化服务', # 游憩、美学价值等 'supporting': '支持服务' # 土壤形成、营养循环等 } def __init__(self, study_area_boundary): """ Args: study_area_boundary: 研究区边界(Shapely Polygon) """ self.boundary = study_area_boundary self.services = {} self.weights = {} self.assessment_result = None def add_service_layer(self, service_type: str, value_layer: np.ndarray, weight: float = 1.0): """ 添加生态系统服务图层 Args: service_type: 服务类型 value_layer: 价值评估栅格 weight: 该服务的权重 """ self.services[service_type] = value_layer self.weights[service_type] = weight def assess_habitat_quality(self, land_use_raster, threat_layers: Dict[str, np.ndarray], sensitivity_table: Dict) -> np.ndarray: """ 评估生境质量 (基于InVEST模型思想) Args: land_use_raster: 土地利用栅格 threat_layers: 威胁因子图层 {'name': array} sensitivity_table: 土地类型对威胁的敏感性 Returns: 生境质量指数栅格 """ rows, cols = land_use_raster.shape habitat_quality = np.zeros((rows, cols)) # 计算退化程度 degradation = np.zeros((rows, cols)) for threat_name, threat_layer in threat_layers.items(): # 对每个威胁因子计算影响 # 这里简化处理，实际需要考虑距离衰减等 threat_impact = threat_layer * 0.5 # 简化权重 for land_type, sensitivity in sensitivity_table.items(): mask = (land_use_raster == land_type) # 获取该土地类型对当前威胁的敏感性 sens = sensitivity.get(threat_name, 0.5) degradation[mask] += threat_impact[mask] * sens # 退化程度归一化到0-1 degradation = np.clip(degradation, 0, 1) # 计算生境质量 for land_type in np.unique(land_use_raster): mask = (land_use_raster == land_type) # 生境适宜性 (简化: 林地=高, 建设用地=低) habitat_suitability = self._get_habitat_suitability(land_type) habitat_quality[mask] = habitat_suitability * (1 - degradation[mask]) return habitat_quality def assess_carbon_storage(self, land_use_raster, carbon_table: Dict[int, Dict[str, float]]) -> np.ndarray: """ 评估碳储量 Args: land_use_raster: 土地利用栅格 carbon_table: {土地类型: {'above': x, 'below': y, 'soil': z}} Returns: 碳储量栅格 """ carbon_storage = np.zeros_like(land_use_raster, dtype=float) for land_type, carbon_values in carbon_table.items(): mask = (land_use_raster == land_type) # 总碳储量 = 地上 + 地下 + 土壤 total_carbon = (carbon_values.get('above', 0) + carbon_values.get('below', 0) + carbon_values.get('soil', 0)) carbon_storage[mask] = total_carbon return carbon_storage def multi_service_assessment(self) -> np.ndarray: """ 多服务综合评估 Returns: 综合生态系统服务指数 """ if not self.services: raise ValueError("请先添加服务图层") # 归一化权重 total_weight = sum(self.weights.values()) # 初始化结果 shape = next(iter(self.services.values())).shape result = np.zeros(shape) # 加权求和 for service_type, layer in self.services.items(): weight = self.weights[service_type] / total_weight result += layer * weight self.assessment_result = result return result def identify_priority_areas(self, threshold_percentile: float = 0.75) -> np.ndarray: """ 识别优先保护区域 Args: threshold_percentile: 分位数阈值 Returns: 优先区域布尔栅格 """ if self.assessment_result is None: self.multi_service_assessment() threshold = np.percentile( self.assessment_result[self.assessment_result > 0], threshold_percentile * 100 ) priority_areas = self.assessment_result >= threshold return priority_areas def _get_habitat_suitability(self, land_type) -> float: """获取土地类型的生境适宜性""" # 简化版：根据土地类型返回适宜性 suitability_map = { 1: 1.0, # 森林 2: 0.8, # 灌木 3: 0.6, # 草地 4: 0.4, # 湿地 5: 0.2, # 耕地 6: 0.0, # 建设用地 7: 0.3 # 裸地 } return suitability_map.get(land_type, 0.5) # 使用示例 def example_enagent_assessment(): """ENAgent生态系统服务评估示例""" from shapely.geometry import box # 创建研究区 study_area = box(0, 0, 10000, 10000) # 创建评估器 assessor = EcosystemServiceAssessment(study_area) # 模拟土地利用数据 rows, cols = 100, 100 np.random.seed(42) land_use = np.random.choice([1, 2, 3, 4, 5, 6], size=(rows, cols), p=[0.3, 0.15, 0.2, 0.1, 0.15, 0.1]) # 评估生境质量 threat_layers = { 'roads': np.random.rand(rows, cols) * 0.8, 'urban': np.random.rand(rows, cols) * 0.6 } sensitivity = { 1: {'roads': 0.3, 'urban': 0.8}, # 森林对城市扩张敏感 2: {'roads': 0.5, 'urban': 0.6}, 3: {'roads': 0.7, 'urban': 0.4}, 4: {'roads': 0.8, 'urban': 0.9}, 5: {'roads': 0.2, 'urban': 0.1}, 6: {'roads': 0.0, 'urban': 0.0} } habitat_quality = assessor.assess_habitat_quality( land_use, threat_layers, sensitivity ) # 评估碳储量 carbon_table = { 1: {'above': 150, 'below': 40, 'soil': 100}, 2: {'above': 60, 'below': 20, 'soil': 80}, 3: {'above': 20, 'below': 50, 'soil': 100}, 4: {'above': 80, 'below': 100, 'soil': 150}, 5: {'above': 10, 'below': 10, 'soil': 80}, 6: {'above': 0, 'below': 0, 'soil': 20} } carbon_storage = assessor.assess_carbon_storage(land_use, carbon_table) # 添加服务图层 assessor.add_service_layer('habitat_quality', habitat_quality, weight=0.5) assessor.add_service_layer('carbon_storage', carbon_storage / 300, weight=0.3) # 水文调节服务（模拟） water_regulation = np.random.rand(rows, cols) * 0.5 + 0.3 assessor.add_service_layer('water_regulation', water_regulation, weight=0.2) # 综合评估 result = assessor.multi_service_assessment() print("生态系统服务综合评估:") print(f" 平均服务指数: {result.mean():.3f}") print(f" 高服务区域 (>0.6): {(result > 0.6).sum()} 个像元") # 识别优先区域 priority = assessor.identify_priority_areas(threshold_percentile=0.75) print(f" 优先保护区域: {priority.sum()} 个像元 ({priority.sum()/priority.size*100:.1f}%)") return assessor if __name__ == "__main__": example_enagent_assessment() ``` --- ## 反思与延伸 ### 思考问题 1. **权重来源**：专家判断、文献参考、数据分析，哪种权重确定方式更可靠？ 2. **准则独立性**：当准则之间存在相关性时，MCDA结果会怎样变化？ 3. **不确定性**：除了权重，数据本身的不确定性如何在MCDA中考虑？ 4. **阈值效应**：某些准则是否存在关键阈值？如何处理？ 5. **空间异质性**：不同区域的权重是否应该不同？ ### 延伸阅读 - **"Multi-Criteria Decision Analysis: Methods and Software"** - MCDA方法综述 - **"Spatial Decision Support Systems"** - 空间决策支持系统 - InVEST模型文档 - 生态系统服务评估实践 - **"Decision Analysis"** (Howard) - 决策分析基础理论 --- ## 关键要点 1. **MCDA处理多目标冲突**：通过系统化方法整合多个准则，支持复杂决策 2. **权重是核心**：权重的确定是MCDA的关键，需要专家知识或数据支持 3. **标准化必不可少**：不同量纲的准则必须标准化才能比较和合并 4. **敏感性分析很重要**：评估权重变化对结果的影响，增强决策稳健性 5. **与GIS结合才有意义**：空间MCDA将决策框架落实到地理空间