pengxiao
a90f7adfa1
将 Markdown 源文件移入 md/，LaTeX 工作目录保留在 latex/， Word 导出移入 word/；删除临时脚本、调试截图和空 stub。 Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
48 KiB
48 KiB
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通过成对比较确定权重，具有结构化和一致性检验的优点：
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通过计算与理想解的距离进行排序：
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的重要组成部分，评估权重变化对决策结果的影响：
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
代码示例
空间多准则决策分析
"""
空间多准则决策分析完整示例
"""
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使用多准则决策分析评估生态系统的综合服务价值：
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()
反思与延伸
思考问题
-
权重来源：专家判断、文献参考、数据分析，哪种权重确定方式更可靠？
-
准则独立性：当准则之间存在相关性时，MCDA结果会怎样变化？
-
不确定性：除了权重，数据本身的不确定性如何在MCDA中考虑？
-
阈值效应：某些准则是否存在关键阈值？如何处理？
-
空间异质性：不同区域的权重是否应该不同？
延伸阅读
- "Multi-Criteria Decision Analysis: Methods and Software" - MCDA方法综述
- "Spatial Decision Support Systems" - 空间决策支持系统
- InVEST模型文档 - 生态系统服务评估实践
- "Decision Analysis" (Howard) - 决策分析基础理论
关键要点
-
MCDA处理多目标冲突：通过系统化方法整合多个准则，支持复杂决策
-
权重是核心：权重的确定是MCDA的关键，需要专家知识或数据支持
-
标准化必不可少：不同量纲的准则必须标准化才能比较和合并
-
敏感性分析很重要：评估权重变化对结果的影响，增强决策稳健性
-
与GIS结合才有意义：空间MCDA将决策框架落实到地理空间