2026_DesignAI/officefile/supplements/02-spatial-intelligence/02.3-multi-criteria-decision.md

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()

---

反思与延伸

思考问题

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将决策框架落实到地理空间