2026_DesignAI/officefile/md/supplements/02-spatial-intelligence/02.5-uncertainty-quantification.md

02.5 不确定性量化

核心问题

分析结果有多可信？不确定性的来源有哪些？ 如何在充满未知的世界中做出稳健的决策？

---

概念讲解

什么是不确定性

不确定性是指知识或信息的缺失，在空间分析中无处不在：

┌─────────────────────────────────────────────────────────────┐
│                  不确定性的来源分类                           │
├─────────────────────────────────────────────────────────────┤
│                                                             │
│  1. 数据不确定性 (Data Uncertainty)                         │
│     ┌─────────────────────────────────────────────────┐     │
│     │  - 测量误差: 仪器精度、人为操作                   │     │
│     │  - 采样偏差: 样本不代表总体                       │     │
│     │  - 空间插值: 从点到面的推断误差                   │     │
│     │  - 分类错误: 遥感解译错误                         │     │
│     │  - 过时数据: 数据不能反映当前状况                 │     │
│     └─────────────────────────────────────────────────┘     │
│                                                             │
│  2. 模型不确定性 (Model Uncertainty)                       │
│     ┌─────────────────────────────────────────────────┐     │
│     │  - 结构不确定性: 模型形式选择                     │     │
│     │  - 参数不确定性: 参数估计误差                     │     │
│     │  - 算法近似: 数值计算的近似                       │     │
│     │  - 尺度失配: 模型尺度与过程尺度不一致              │     │
│     └─────────────────────────────────────────────────┘     │
│                                                             │
│  3. 情境不确定性 (Scenario Uncertainty)                    │
│     ┌─────────────────────────────────────────────────┐     │
│     │  - 未来不可预测: 气候变化、政策变动               │     │
│     │  - 行为主体响应: 利益相关者的反应                 │     │
│     │  - 突发事件: 自然灾害、社会事件                   │     │
│     └─────────────────────────────────────────────────┘     │
│                                                             │
└─────────────────────────────────────────────────────────────┘

不确定性的类型

类型	描述	处理方法
随机性 (Aleatory)	系统内在的随机变化，不可减少	概率分布、随机模拟
认知性 (Epistemic)	知识缺失导致的，可通过研究减少	灵敏度分析、更多数据
模糊性 (Vagueness)	概念边界不清晰	模糊集合、模糊逻辑
歧义性 (Ambiguity)	多种解释都合理	情景分析、鲁棒优化

不确定性传播

当多个不确定输入通过模型组合时，不确定性会传播：

输入不确定性 ──→ 模型 ──→ 输出不确定性

                    ┌─────────┐
    x₁ ± Δx₁  ──→│         │
    x₂ ± Δx₂  ──→│  f(x)   │──→  y ± Δy
    x₃ ± Δx₃  ──→│         │
                └─────────┘

传播规则:
- 线性模型: Δy ≈ √(Σ(∂f/∂xi)² × Δxi²)  (误差传播公式)
- 非线性模型: 需要蒙特卡洛模拟
- 相关输入: 需要考虑协方差

不确定性量化的方法谱系

┌─────────────────────────────────────────────────────────────┐
│                  不确定性量化方法                              │
├─────────────────────────────────────────────────────────────┤
│                                                             │
│  1. 确定性敏感性分析 (Deterministic SA)                     │
│     ┌─────────────────────────────────────────────────┐     │
│     │  - OAT (One-at-a-Time): 单因素变化               │     │
│     │  - 局部敏感性: 导数、弹性系数                     │     │
│     │                                                  │     │
│     │  优点: 简单、直观                                 │     │
│     │  缺点: 忽略参数交互                               │     │
│     └─────────────────────────────────────────────────┘     │
│                                                             │
│  2. 全局敏感性分析 (Global SA)                              │
│     ┌─────────────────────────────────────────────────┐     │
│     │  - Sobol指数: 方差分解                           │     │
│     │  - Morris筛选: 定性筛选重要参数                  │     │
│     │  - FAST: 傅里叶幅度敏感性测试                    │     │
│     │                                                  │     │
│     │  优点: 考虑参数空间、交互作用                     │     │
│     │  缺点: 计算成本高                                 │     │
│     └─────────────────────────────────────────────────┘     │
│                                                             │
│  3. 蒙特卡洛方法 (Monte Carlo)                              │
│     ┌─────────────────────────────────────────────────┐     │
│     │  - 随机采样: 从输入分布采样                       │     │
│     │  - LHS: 拉丁超立方采样                            │     │
│     │  - 贝叶斯推断: 更新参数分布                       │     │
│     │                                                  │     │
│     │  优点: 通用、易于实现                             │     │
│     │  缺点: 收敛慢、高维困难                           │     │
│     └─────────────────────────────────────────────────┘     │
│                                                             │
│  4. 场景分析 (Scenario Analysis)                            │
│     ┌─────────────────────────────────────────────────┐     │
│     │  - 定义多个合理情景                               │     │
│     │  - 比较情景结果                                   │     │
│     │  - 识别稳健策略                                   │     │
│     │                                                  │     │
│     │  优点: 直观、易于沟通                             │     │
│     │  缺点: 情景选择主观                               │     │
│     └─────────────────────────────────────────────────┘     │
│                                                             │
└─────────────────────────────────────────────────────────────┘

---

设计原理

蒙特卡洛模拟

蒙特卡洛是最通用的不确定性量化方法：

import numpy as np
from typing import Callable, List, Dict, Tuple, Optional
import matplotlib.pyplot as plt
from scipy import stats

class MonteCarloSimulator:
    """
    蒙特卡洛模拟器

    核心思想：通过大量随机采样估计输出的概率分布
    """

    def __init__(self,
                 model: Callable,
                 n_simulations: int = 10000,
                 random_seed: Optional[int] = None):
        """
        Args:
            model: 输入→输出的函数
            n_simulations: 模拟次数
            random_seed: 随机种子
        """
        self.model = model
        self.n_simulations = n_simulations
        self.random_seed = random_seed

        self.inputs = None
        self.outputs = None
        self.input_distributions = {}

    def define_input(self, name: str, distribution: str, **params):
        """
        定义输入变量的概率分布

        Args:
            name: 变量名
            distribution: 分布类型 ('normal', 'uniform', 'triangular', 'lognormal', 等)
            **params: 分布参数
        """
        self.input_distributions[name] = {
            'type': distribution,
            'params': params
        }

    def generate_inputs(self, method: str = 'random') -> np.ndarray:
        """
        生成输入样本

        Args:
            method: 'random' 或 'lhs'(拉丁超立方)

        Returns:
            输入样本数组 (n_simulations × n_variables)
        """
        if self.random_seed is not None:
            np.random.seed(self.random_seed)

        n_vars = len(self.input_distributions)
        var_names = list(self.input_distributions.keys())

        if method == 'random':
            samples = np.zeros((self.n_simulations, n_vars))

            for i, name in enumerate(var_names):
                dist_info = self.input_distributions[name]
                samples[:, i] = self._sample_distribution(
                    dist_info['type'],
                    dist_info['params'],
                    self.n_simulations
                )

        elif method == 'lhs':
            # 拉丁超立方采样
            from scipy.stats import qmc

            sampler = qmc.LatinHypercube(d=n_vars, seed=self.random_seed)
            sample_unit = sampler.random(n=self.n_simulations)

            samples = np.zeros_like(sample_unit)

            for i, name in enumerate(var_names):
                dist_info = self.input_distributions[name]
                samples[:, i] = self._transform_unit(
                    sample_unit[:, i],
                    dist_info['type'],
                    dist_info['params']
                )

        self.inputs = samples
        return samples

    def _sample_distribution(self, dist_type: str,
                             params: Dict, size: int) -> np.ndarray:
        """从指定分布采样"""
        if dist_type == 'normal':
            return np.random.normal(params['mean'], params['std'], size)
        elif dist_type == 'uniform':
            return np.random.uniform(params['low'], params['high'], size)
        elif dist_type == 'triangular':
            return np.random.triangular(
                params['left'], params['mode'], params['right'], size
            )
        elif dist_type == 'lognormal':
            return np.random.lognormal(params['mean'], params['sigma'], size)
        elif dist_type == 'beta':
            return np.random.beta(params['a'], params['b'], size)
        else:
            raise ValueError(f"Unknown distribution: {dist_type}")

    def _transform_unit(self, unit_samples: np.ndarray,
                       dist_type: str, params: Dict) -> np.ndarray:
        """将[0,1]均匀分布转换为目标分布"""
        if dist_type == 'normal':
            return stats.norm.ppf(unit_samples, loc=params['mean'], scale=params['std'])
        elif dist_type == 'uniform':
            return params['low'] + unit_samples * (params['high'] - params['low'])
        elif dist_type == 'triangular':
            return stats.triang.ppf(
                unit_samples,
                c=(params['mode']-params['left'])/(params['right']-params['left']),
                loc=params['left'],
                scale=params['right']-params['left']
            )
        elif dist_type == 'lognormal':
            return stats.lognorm.ppf(
                unit_samples,
                s=params['sigma'],
                scale=np.exp(params['mean'])
            )
        else:
            raise ValueError(f"Unknown distribution: {dist_type}")

    def run(self, method: str = 'random') -> np.ndarray:
        """
        运行蒙特卡洛模拟

        Returns:
            输出样本数组
        """
        # 生成输入
        if self.inputs is None:
            self.generate_inputs(method)

        # 运行模型
        var_names = list(self.input_distributions.keys())
        self.outputs = np.zeros(self.n_simulations)

        for i in range(self.n_simulations):
            input_dict = {name: self.inputs[i, j]
                         for j, name in enumerate(var_names)}
            self.outputs[i] = self.model(**input_dict)

        return self.outputs

    def analyze_output(self) -> Dict:
        """
        分析输出分布

        Returns:
            统计摘要
        """
        if self.outputs is None:
            raise RuntimeError("请先运行模拟")

        output = self.outputs[~np.isnan(self.outputs)]  # 移除NaN

        return {
            'mean': np.mean(output),
            'std': np.std(output),
            'median': np.median(output),
            'min': np.min(output),
            'max': np.max(output),
            'percentiles': {
                '5': np.percentile(output, 5),
                '25': np.percentile(output, 25),
                '75': np.percentile(output, 75),
                '95': np.percentile(output, 95)
            },
            'ci_95': (np.percentile(output, 2.5), np.percentile(output, 97.5))
        }

    def plot_output(self, bins: int = 50):
        """绘制输出分布直方图"""
        if self.outputs is None:
            raise RuntimeError("请先运行模拟")

        output = self.outputs[~np.isnan(self.outputs)]

        fig, axes = plt.subplots(1, 2, figsize=(12, 4))

        # 直方图
        axes[0].hist(output, bins=bins, density=True, alpha=0.7, edgecolor='black')
        axes[0].axvline(np.mean(output), color='red', linestyle='--', label='Mean')
        axes[0].axvline(np.median(output), color='green', linestyle='--', label='Median')
        axes[0].set_xlabel('Output Value')
        axes[0].set_ylabel('Probability Density')
        axes[0].set_title('Output Distribution')
        axes[0].legend()

        # 累积分布
        sorted_output = np.sort(output)
        cumulative = np.arange(1, len(sorted_output) + 1) / len(sorted_output)

        axes[1].plot(sorted_output, cumulative, linewidth=2)
        axes[1].axvline(np.percentile(output, 5), color='orange', linestyle='--',
                       label='5th percentile')
        axes[1].axvline(np.percentile(output, 95), color='orange', linestyle='--',
                       label='95th percentile')
        axes[1].set_xlabel('Output Value')
        axes[1].set_ylabel('Cumulative Probability')
        axes[1].set_title('Cumulative Distribution Function')
        axes[1].legend()

        plt.tight_layout()
        return fig

    def sensitivity_analysis(self) -> Dict[str, float]:
        """
        计算敏感性指标（基于相关性）

        Returns:
            各输入的敏感性系数
        """
        if self.inputs is None or self.outputs is None:
            raise RuntimeError("请先运行模拟")

        var_names = list(self.input_distributions.keys())
        sensitivity = {}

        for i, name in enumerate(var_names):
            # Spearman秩相关（对单调关系更稳健）
            corr, _ = stats.spearmanr(self.inputs[:, i], self.outputs)
            sensitivity[name] = abs(corr)

        # 归一化
        total = sum(sensitivity.values())
        if total > 0:
            sensitivity = {k: v/total for k, v in sensitivity.items()}

        return sensitivity

Sobol全局敏感性分析

Sobol指数基于方差分解，是全局敏感性分析的金标准：

import numpy as np
from typing import Callable, Dict, List
from scipy import stats

class SobolAnalyzer:
    """
    Sobol全局敏感性分析

    基于方差分解，计算一阶和高阶敏感性指数
    """

    def __init__(self,
                 model: Callable,
                 n_vars: int,
                 bounds: List[Tuple[float, float]]):
        """
        Args:
            model: 输入→输出的函数
            n_vars: 输入变量数量
            bounds: 每个变量的边界 [(min, max), ...]
        """
        self.model = model
        self.n_vars = n_vars
        self.bounds = bounds

        self.S1 = None  # 一阶效应
        self.ST = None  # 总效应

    def generate_samples(self, N: int) -> Dict[str, np.ndarray]:
        """
        生成Sobol序列样本

        使用Saltelli采样方案

        Args:
            N: 基础样本数

        Returns:
            包含A、B、AB矩阵的字典
        """
        # 生成两个基础样本矩阵
        A = np.zeros((N, self.n_vars))
        B = np.zeros((N, self.n_vars))

        for i in range(self.n_vars):
            A[:, i] = np.random.uniform(self.bounds[i][0], self.bounds[i][1], N)
            B[:, i] = np.random.uniform(self.bounds[i][0], self.bounds[i][1], N)

        # 生成AB矩阵（每次替换一列）
        AB = np.zeros((self.n_vars, N, self.n_vars))
        for i in range(self.n_vars):
            AB[i] = A.copy()
            AB[i][:, i] = B[:, i]

        return {'A': A, 'B': B, 'AB': AB}

    def compute(self, N: int) -> Dict[str, np.ndarray]:
        """
        计算Sobol指数

        Args:
            N: 基础样本数

        Returns:
            包含S1和ST的字典
        """
        samples = self.generate_samples(N)

        # 计算模型输出
        fA = np.array([self.model(*x) for x in samples['A']])
        fB = np.array([self.model(*x) for x in samples['B']])

        fAB = np.zeros((self.n_vars, N))
        for i in range(self.n_vars):
            fAB[i] = np.array([self.model(*x) for x in samples['AB'][i]])

        # 计算总方差
        all_outputs = np.concatenate([fA, fB])
        V = np.var(all_outputs)

        if V == 0:
            raise ValueError("模型输出方差为0，无法计算敏感性")

        # 一阶敏感性 (S1)
        S1 = np.zeros(self.n_vars)
        for i in range(self.n_vars):
            numerator = np.mean(fA * fAB[i]) - np.mean(fA) ** 2
            S1[i] = numerator / V

        # 总效应 (ST)
        ST = np.zeros(self.n_vars)
        for i in range(self.n_vars):
            numerator = np.mean((fA - fAB[i]) ** 2)
            ST[i] = numerator / (2 * V)

        self.S1 = np.maximum(S1, 0)  # 确保非负
        self.ST = np.maximum(ST, 0)

        return {'S1': self.S1, 'ST': self.ST}

    def plot_sensitivity(self, var_names: List[str] = None):
        """绘制敏感性指数图"""
        if self.S1 is None or self.ST is None:
            raise RuntimeError("请先运行计算")

        if var_names is None:
            var_names = [f'X{i+1}' for i in range(self.n_vars)]

        fig, ax = plt.subplots(figsize=(10, 6))

        x = np.arange(self.n_vars)
        width = 0.35

        ax.bar(x - width/2, self.S1, width, label='First Order (S1)', alpha=0.8)
        ax.bar(x + width/2, self.ST, width, label='Total Effect (ST)', alpha=0.8)

        ax.set_xlabel('Input Variables')
        ax.set_ylabel('Sensitivity Index')
        ax.set_title('Sobol Sensitivity Indices')
        ax.set_xticks(x)
        ax.set_xticklabels(var_names, rotation=45)
        ax.legend()

        plt.tight_layout()
        return fig

不确定性可视化

空间不确定性需要特别的可视化方法：

import numpy as np
import matplotlib.pyplot as plt
from typing import List, Dict, Optional

class UncertaintyVisualizer:
    """
    空间不确定性可视化
    """

    @staticmethod
    def probability_map(mean: np.ndarray,
                       threshold: float,
                       direction: str = 'above') -> np.ndarray:
        """
        创建概率地图

        Args:
            mean: 均值栅格
            threshold: 阈值
            direction: 'above' 或 'below'

        Returns:
            超过/低于阈值的概率（需要配合std）
        """
        # 这里简化处理，实际需要蒙特卡洛结果
        if direction == 'above':
            # 假设正态分布
            # 实际应使用MC结果的累积分布
            pass
        return mean > threshold

    @staticmethod
    def confidence_interval(mean: np.ndarray,
                           std: np.ndarray,
                           confidence: float = 0.95) -> Dict[str, np.ndarray]:
        """
        计算置信区间

        Args:
            mean: 均值栅格
            std: 标准差栅格
            confidence: 置信水平

        Returns:
            包含下界和上界的字典
        """
        from scipy import stats

        z = stats.norm.ppf(1 - (1 - confidence) / 2)

        return {
            'lower': mean - z * std,
            'upper': mean + z * std,
            'margin': z * std
        }

    @staticmethod
    def uncertainty_classification(mean: np.ndarray,
                                  std: np.ndarray,
                                  n_classes: int = 5) -> np.ndarray:
        """
        基于均值和不确定性的分类

        结合期望值和不确定性进行决策分类

        Args:
            mean: 均值栅格
            std: 标准差栅格
            n_classes: 分类数

        Returns:
            分类栅格
        """
        # 标准化
        mean_norm = (mean - mean.min()) / (mean.max() - mean.min())
        std_norm = (std - std.min()) / (std.max() - std.min() + 1e-10)

        # 决策分类
        # 高值+低不确定性 = 高优先级
        # 高值+高不确定性 = 需要更多信息
        # 低值+低不确定性 = 低优先级
        # 低值+高不确定性 = 不确定

        classification = np.zeros_like(mean, dtype=int)

        # 定义阈值
        high_value = mean_norm > 0.6
        low_uncertainty = std_norm < 0.4

        classification[high_value & low_uncertainty] = 5  # 高值，确定
        classification[high_value & ~low_uncertainty] = 4  # 高值，不确定
        classification[~high_value & low_uncertainty] = 2  # 低值，确定
        classification[~high_value & ~low_uncertainty] = 3  # 低值，不确定

        return classification

    @staticmethod
    def plot_with_uncertainty(mean: np.ndarray,
                             std: np.ndarray,
                             title: str = 'Value with Uncertainty'):
        """
        绘制带不确定性的地图

        使用颜色表示值，透明度表示不确定性
        """
        fig, axes = plt.subplots(1, 3, figsize=(15, 5))

        # 均值图
        im1 = axes[0].imshow(mean, cmap='RdYlGn')
        axes[0].set_title('Mean Value')
        plt.colorbar(im1, ax=axes[0])

        # 不确定性图
        im2 = axes[1].imshow(std, cmap='Oranges')
        axes[1].set_title('Uncertainty (Std)')
        plt.colorbar(im2, ax=axes[1])

        # 组合图（值+不确定性）
        # 归一化
        mean_norm = (mean - mean.min()) / (mean.max() - mean.min())
        std_norm = (std - std.min()) / (std.max() - std.min())

        # 创建RGBA图像
        cmap = plt.cm.RdYlGn
        rgba = cmap(mean_norm)

        # 用alpha通道表示不确定性（高不确定性=低透明度）
        rgba[:, :, 3] = 1 - std_norm * 0.7

        axes[2].imshow(rgba)
        axes[2].set_title('Value (color) + Uncertainty (alpha)')

        plt.suptitle(title)
        plt.tight_layout()
        return fig

---

代码示例

阻力面不确定性分析

"""
阻力面不确定性量化示例

ENAgent中阻力面的参数不确定性分析
"""
import numpy as np
import matplotlib.pyplot as plt
from typing import Dict, List, Tuple, Callable

class ResistanceUncertaintyAnalyzer:
    """
    阻力面不确定性分析器

    分析不同土地类型阻力值的不确定性对结果的影响
    """

    def __init__(self,
                 land_use_raster: np.ndarray,
                 base_resistance_dict: Dict[int, float]):
        """
        Args:
            land_use_raster: 土地利用栅格
            base_resistance_dict: 基准阻力值 {土地类型: 阻力值}
        """
        self.land_use = land_use_raster
        self.base_resistance = base_resistance_dict
        self.land_types = list(base_resistance_dict.keys())

        # 阻力面
        self.base_surface = self._create_surface(base_resistance_dict)

        # 定义阻力值的先验分布
        self.resistance_distributions = self._define_distributions()

    def _create_surface(self, resistance_dict: Dict[int, float]) -> np.ndarray:
        """根据阻力字典创建阻力面"""
        surface = np.zeros_like(self.land_use, dtype=float)
        for land_type, resistance in resistance_dict.items():
            surface[self.land_use == land_type] = resistance
        return surface

    def _define_distributions(self) -> Dict[int, Dict]:
        """
        定义阻力值的概率分布

        假设阻力值服从对数正态分布
        """
        distributions = {}
        for land_type, base_value in self.base_resistance.items():
            # 变异系数随阻力值增加
            cv = 0.3 if base_value > 50 else 0.2

            # 对数正态分布参数
            # 如果X ~ Lognormal(μ, σ), 则 E[X] = exp(μ + σ²/2)
            # Var[X] = (exp(σ²) - 1) * exp(2μ + σ²)
            sigma2 = np.log(1 + cv**2)
            mu = np.log(base_value) - sigma2 / 2

            distributions[land_type] = {
                'type': 'lognormal',
                'mean': mu,
                'sigma': np.sqrt(sigma2),
                'base': base_value
            }

        return distributions

    def sample_resistance_surface(self, n_samples: int = 1) -> List[np.ndarray]:
        """
        采样阻力面

        Args:
            n_samples: 采样次数

        Returns:
            采样得到的阻力面列表
        """
        surfaces = []

        for _ in range(n_samples):
            resistance_dict = {}
            for land_type, dist in self.resistance_distributions.items():
                # 从对数正态分布采样
                value = np.random.lognormal(dist['mean'], dist['sigma'])
                resistance_dict[land_type] = value

            surfaces.append(self._create_surface(resistance_dict))

        return surfaces

    def monte_carlo_connectivity(self,
                                 sources: List[Tuple[int, int]],
                                 n_samples: int = 100) -> Dict:
        """
        蒙特卡洛连通性分析

        Args:
            sources: 源地列表
            n_samples: 采样次数

        Returns:
            统计结果
        """
        from connectivity import compute_cost_distance

        all_cost_distances = []
        all_paths = []

        for i in range(n_samples):
            # 采样阻力面
            surface = self.sample_resistance_surface(1)[0]

            # 计算成本距离
            cost_distance = compute_cost_distance(surface, sources)

            all_cost_distances.append(cost_distance)

        # 转换为数组
        all_cost_distances = np.array(all_cost_distances)

        # 统计分析
        mean_cost = np.mean(all_cost_distances, axis=0)
        std_cost = np.std(all_cost_distances, axis=0)
        cv_cost = std_cost / (mean_cost + 1e-10)  # 变异系数

        # 找出不确定性高的区域
        high_uncertainty_mask = cv_cost > cv_cost.mean() + cv_cost.std()

        return {
            'mean': mean_cost,
            'std': std_cost,
            'cv': cv_cost,
            'high_uncertainty': high_uncertainty_mask,
            'all_samples': all_cost_distances
        }

    def sensitivity_to_resistance(self,
                                  sources: List[Tuple[int, int]],
                                  variation: float = 0.3) -> Dict:
        """
        阻力值敏感性分析 (OAT方法)

        每次改变一个土地类型的阻力值

        Args:
            sources: 源地列表
            variation: 变化幅度 (±30%)

        Returns:
            敏感性结果
        """
        from connectivity import compute_cost_distance

        # 基准结果
        base_cost = compute_cost_distance(self.base_surface, sources)

        sensitivity = {}

        for land_type in self.land_types:
            base_value = self.base_resistance[land_type]

            # 测试不同阻力值
            test_values = [
                base_value * (1 - variation),  # 减少
                base_value * (1 + variation)   # 增加
            ]

            results = []
            for test_value in test_values:
                test_dict = self.base_resistance.copy()
                test_dict[land_type] = test_value
                test_surface = self._create_surface(test_dict)
                test_cost = compute_cost_distance(test_surface, sources)

                # 计算与基准的差异
                diff = np.abs(test_cost - base_cost).mean()
                results.append(diff)

            # 敏感性指标 = 平均绝对变化
            sensitivity[land_type] = np.mean(results)

        # 归一化
        total = sum(sensitivity.values())
        if total > 0:
            sensitivity = {k: v/total for k, v in sensitivity.items()}

        return sensitivity

    def robust_corridor_selection(self,
                                  sources: List[Tuple[int, int]],
                                  cost_threshold: float,
                                  n_samples: int = 100,
                                  reliability: float = 0.8) -> np.ndarray:
        """
        稳健廊道选择

        只选择在多数采样中都满足阈值的像元

        Args:
            sources: 源地列表
            cost_threshold: 成本阈值
            n_samples: 采样次数
            reliability: 可靠性要求 (80%的采样满足)

        Returns:
            稳健廊道掩模
        """
        mc_results = self.monte_carlo_connectivity(sources, n_samples)

        # 计算每个像元在多少比例的采样中满足阈值
        satisfy_count = np.zeros_like(mc_results['mean'])

        for sample_cost in mc_results['all_samples']:
            satisfy_count += (sample_cost < cost_threshold).astype(int)

        satisfy_ratio = satisfy_count / n_samples

        # 只选择满足可靠性要求的像元
        robust_corridor = satisfy_ratio >= reliability

        return robust_corridor

# 使用示例
def example_resistance_uncertainty():
    """阻力面不确定性分析示例"""
    # 创建示例土地利用数据
    np.random.seed(42)
    rows, cols = 50, 50

    # 土地类型: 1=森林, 2=灌木, 3=草地, 4=农田, 5=建设用地, 6=水体
    land_use = np.random.choice(
        [1, 2, 3, 4, 5, 6],
        size=(rows, cols),
        p=[0.25, 0.15, 0.20, 0.25, 0.10, 0.05]
    )

    # 基准阻力值
    base_resistance = {
        1: 10,   # 森林 - 低阻力
        2: 30,   # 灌木
        3: 50,   # 草地
        4: 80,   # 农田 - 高阻力
        5: 100,  # 建设用地 - 最高阻力
        6: 150   # 水体 - 障碍
    }

    # 创建分析器
    analyzer = ResistanceUncertaintyAnalyzer(land_use, base_resistance)

    # 源地
    sources = [(10, 10), (40, 40)]

    # 蒙特卡洛分析
    print("Running Monte Carlo analysis...")
    mc_results = analyzer.monte_carlo_connectivity(sources, n_samples=50)

    print(f"Mean cost distance: {mc_results['mean'][sources[0]]:.2f}")
    print(f"Std of cost distance: {mc_results['std'][sources[0]]:.2f}")
    print(f"High uncertainty pixels: {mc_results['high_uncertainty'].sum()}")

    # 敏感性分析
    print("\nRunning sensitivity analysis...")
    sensitivity = analyzer.sensitivity_to_resistance(sources, variation=0.3)

    print("Sensitivity to resistance values:")
    land_type_names = {1: '森林', 2: '灌木', 3: '草地', 4: '农田', 5: '建设用地', 6: '水体'}
    for land_type, sens in sorted(sensitivity.items(), key=lambda x: -x[1]):
        print(f"  {land_type_names[land_type]}: {sens:.3f}")

    # 稳健廊道选择
    print("\nSelecting robust corridors...")
    threshold = np.percentile(mc_results['mean'], 60)
    robust_corridor = analyzer.robust_corridor_selection(
        sources, threshold, n_samples=50, reliability=0.7
    )

    print(f"Robust corridor pixels: {robust_corridor.sum()}")

    return analyzer

if __name__ == "__main__":
    example_resistance_uncertainty()

情景分析框架

"""
情景分析框架

处理情境不确定性
"""
import numpy as np
from typing import Dict, List, Callable, Any
from dataclasses import dataclass
from enum import Enum

class ScenarioType(Enum):
    """情景类型"""
    OPTIMISTIC = "乐观"
    PESSIMISTIC = "悲观"
    BUSINESS_AS_USUAL = "照常"
    SUSTAINABLE = "可持续"

@dataclass
class Scenario:
    """情景定义"""
    name: str
    description: str
    parameters: Dict[str, Any]
    probability: float = 1.0  # 情景发生的概率

class ScenarioAnalyzer:
    """
    情景分析器

    通过定义多个合理情景来处理情境不确定性
    """

    def __init__(self, model: Callable):
        """
        Args:
            model: 接受参数字典并返回结果的函数
        """
        self.model = model
        self.scenarios: Dict[str, Scenario] = {}
        self.results: Dict[str, Any] = {}

    def add_scenario(self, scenario_id: str, scenario: Scenario):
        """添加情景"""
        self.scenarios[scenario_id] = scenario

    def define_land_use_scenarios(self,
                                  base_year: int,
                                  target_year: int) -> Dict[str, Scenario]:
        """
        定义土地利用变化情景

        Args:
            base_year: 基准年
            target_year: 目标年

        Returns:
            定义的情景字典
        """
        scenarios = {
            'bau': Scenario(
                name='Business as Usual',
                description='延续当前发展趋势',
                parameters={
                    'urban_expansion_rate': 0.02,
                    'forest_loss_rate': 0.01,
                    'agriculture_intensity': 1.0
                },
                probability=0.5
            ),
            'optimistic': Scenario(
                name='Optimistic',
                description='生态保护加强，可持续发展',
                parameters={
                    'urban_expansion_rate': 0.01,
                    'forest_loss_rate': -0.005,  # 森林恢复
                    'agriculture_intensity': 1.2  # 精准农业
                },
                probability=0.2
            ),
            'pessimistic': Scenario(
                name='Pessimistic',
                description='快速城市化，生态退化',
                parameters={
                    'urban_expansion_rate': 0.04,
                    'forest_loss_rate': 0.03,
                    'agriculture_intensity': 0.8
                },
                probability=0.3
            )
        }

        for scenario_id, scenario in scenarios.items():
            self.add_scenario(scenario_id, scenario)

        return scenarios

    def run_scenarios(self) -> Dict[str, Any]:
        """运行所有情景"""
        self.results = {}

        for scenario_id, scenario in self.scenarios.items():
            try:
                result = self.model(scenario.parameters)
                self.results[scenario_id] = {
                    'result': result,
                    'scenario': scenario
                }
            except Exception as e:
                self.results[scenario_id] = {
                    'error': str(e),
                    'scenario': scenario
                }

        return self.results

    def compare_results(self, metric_extractor: Callable = None) -> Dict:
        """
        比较情景结果

        Args:
            metric_extractor: 从结果中提取比较指标的函数

        Returns:
            比较结果
        """
        comparison = {}

        for scenario_id, result_data in self.results.items():
            if 'error' in result_data:
                comparison[scenario_id] = {'error': result_data['error']}
            else:
                result = result_data['result']
                scenario = result_data['scenario']

                if metric_extractor:
                    metrics = metric_extractor(result)
                else:
                    # 默认: 直接使用结果
                    metrics = {'value': result}

                comparison[scenario_id] = {
                    'metrics': metrics,
                    'probability': scenario.probability,
                    'description': scenario.description
                }

        return comparison

    def weighted_outcome(self, metric_name: str = 'value') -> float:
        """
        计算加权期望结果

        Args:
            metric_name: 要加权平均的指标名称

        Returns:
            期望值
        """
        total = 0
        total_prob = 0

        for scenario_id, result_data in self.results.items():
            if 'error' not in result_data:
                scenario = result_data['scenario']
                result = result_data['result']

                if isinstance(result, dict):
                    value = result.get(metric_name, 0)
                else:
                    value = result

                total += value * scenario.probability
                total_prob += scenario.probability

        return total / total_prob if total_prob > 0 else 0

    def identify_robust_strategy(self,
                                 strategies: Dict[str, Dict],
                                 criterion: str = 'maximin') -> str:
        """
        识别稳健策略

        Args:
            strategies: 策略字典 {strategy_name: parameters}
            criterion: 'maximin' (最大化最小收益) 或 'maximize_expected'

        Returns:
            最稳健的策略名称
        """
        strategy_outcomes = {}

        # 评估每个策略在各情景下的表现
        for strategy_name, strategy_params in strategies.items():
            outcomes = []

            for scenario_id, scenario in self.scenarios.items():
                # 合并策略和情景参数
                combined_params = {**strategy_params, **scenario.parameters}

                try:
                    result = self.model(combined_params)

                    if isinstance(result, dict):
                        value = result.get('value', result.get('score', 0))
                    else:
                        value = result

                    outcomes.append(value)
                except:
                    outcomes.append(-float('inf'))

            strategy_outcomes[strategy_name] = outcomes

        # 根据准则选择
        if criterion == 'maximin':
            # 选择最差情景下表现最好的策略
            best_strategy = max(
                strategy_outcomes.keys(),
                key=lambda s: min(strategy_outcomes[s])
            )
        elif criterion == 'maximize_expected':
            # 选择期望值最高的策略
            best_strategy = max(
                strategy_outcomes.keys(),
                key=lambda s: np.mean(strategy_outcomes[s])
            )
        else:
            raise ValueError(f"Unknown criterion: {criterion}")

        return best_strategy

---

案例分析

ENAgent中的不确定性量化实践

ENAgent在生态网络分析中系统性地应用不确定性量化：

class ENAgentUncertaintyModule:
    """
    ENAgent不确定性量化模块

    整合多种方法处理生态网络分析中的不确定性
    """

    def __init__(self, enagent_core):
        """
        Args:
            enagent_core: ENAgent核心实例
        """
        self.core = enagent_core
        self.uncertainty_results = {}

    def full_uncertainty_analysis(self,
                                   sources: List[Dict],
                                   n_samples: int = 100) -> Dict:
        """
        完整的不确定性分析

        包括:
        1. 数据不确定性 (阻力值变化)
        2. 参数不确定性 (源地质量)
        3. 情景不确定性 (不同发展情景)

        Args:
            sources: 源地列表
            n_samples: 蒙特卡洛采样次数

        Returns:
            不确定性分析结果
        """
        results = {}

        # 1. 数据不确定性 - 阻力面
        print("Analyzing data uncertainty (resistance surface)...")
        results['data'] = self._resistance_uncertainty(sources, n_samples)

        # 2. 参数不确定性 - 源地质量
        print("Analyzing parameter uncertainty (source quality)...")
        results['parameter'] = self._source_quality_uncertainty(sources, n_samples)

        # 3. 情景不确定性
        print("Analyzing scenario uncertainty...")
        results['scenario'] = self._scenario_uncertainty(sources)

        # 4. 综合分析
        results['summary'] = self._synthesize_results(results)

        self.uncertainty_results = results
        return results

    def _resistance_uncertainty(self, sources, n_samples):
        """阻力面数据不确定性分析"""
        # 获取基准阻力
        base_resistance = self.core.resistance_surface

        # 定义阻力值的变异系数
        resistance_cv = {
            'forest': 0.2,
            'grassland': 0.3,
            'agriculture': 0.25,
            'urban': 0.15,
            'water': 0.1
        }

        # 蒙特卡洛采样
        connectivity_results = []

        for _ in range(n_samples):
            # 扰动阻力值
            perturbed_resistance = self._perturb_resistance(
                base_resistance, resistance_cv
            )

            # 计算连通性
            connectivity = self.core.compute_connectivity(
                sources, perturbed_resistance
            )
            connectivity_results.append(connectivity)

        # 统计分析
        connectivity_array = np.array(connectivity_results)

        return {
            'mean': connectivity_array.mean(axis=0),
            'std': connectivity_array.std(axis=0),
            'percentiles': {
                '5': np.percentile(connectivity_array, 5, axis=0),
                '25': np.percentile(connectivity_array, 25, axis=0),
                '75': np.percentile(connectivity_array, 75, axis=0),
                '95': np.percentile(connectivity_array, 95, axis=0)
            }
        }

    def _source_quality_uncertainty(self, sources, n_samples):
        """源地质量参数不确定性分析"""
        results = []

        for _ in range(n_samples):
            # 扰动源地质量 (假设±20%)
            perturbed_sources = []
            for source in sources:
                perturbed = source.copy()
                perturbed['quality'] = source['quality'] * np.random.uniform(0.8, 1.2)
                perturbed_sources.append(perturbed)

            # 计算网络指标
            metrics = self.core.compute_network_metrics(perturbed_sources)
            results.append(metrics)

        # 统计
        return self._summarize_metrics(results)

    def _scenario_uncertainty(self, sources):
        """情景不确定性分析"""
        scenarios = {
            'current': {
                'description': '当前状况',
                'urban_expansion': 0,
                'restoration': 0
            },
            'urban_growth': {
                'description': '城市扩张情景',
                'urban_expansion': 0.5,  # 50%扩张
                'restoration': 0
            },
            'restoration': {
                'description': '生态修复情景',
                'urban_expansion': 0,
                'restoration': 0.3  # 30%修复
            }
        }

        results = {}
        for scenario_name, params in scenarios.items():
            # 应用情景参数
            modified_resistance = self.core.apply_scenario(params)

            # 计算结果
            connectivity = self.core.compute_connectivity(
                sources, modified_resistance
            )
            metrics = self.core.compute_network_metrics(sources)

            results[scenario_name] = {
                'connectivity': connectivity,
                'metrics': metrics,
                'description': params['description']
            }

        return results

    def _synthesize_results(self, results):
        """综合各种不确定性分析结果"""
        synthesis = {
            'recommendations': [],
            'confidence_levels': {}
        }

        # 分析不同来源的不确定性
        data_cv = results['data']['std'].mean() / results['data']['mean'].mean()
        parameter_cv = results['parameter']['std'].mean() / results['parameter']['mean'].mean()

        # 根据不确定性的相对大小给出建议
        if data_cv > parameter_cv * 1.5:
            synthesis['recommendations'].append(
                "数据不确定性是主要来源，建议提高阻力面数据质量"
            )
        elif parameter_cv > data_cv * 1.5:
            synthesis['recommendations'].append(
                "参数不确定性是主要来源，建议更精确地评估源地质量"
            )

        # 情景比较
        scenario_results = results['scenario']
        if 'restoration' in scenario_results and 'urban_growth' in scenario_results:
            restoration_connectivity = scenario_results['restoration']['connectivity'].mean()
            urban_connectivity = scenario_results['urban_growth']['connectivity'].mean()

            if restoration_connectivity > urban_connectivity * 1.2:
                synthesis['recommendations'].append(
                    "生态修复情景显著改善连通性，建议优先考虑生态修复措施"
                )

        return synthesis

    def generate_uncertainty_report(self) -> str:
        """生成不确定性分析报告"""
        if not self.uncertainty_results:
            return "请先运行不确定性分析"

        report = []
        report.append("=" * 60)
        report.append("生态网络不确定性分析报告")
        report.append("=" * 60)
        report.append("")

        # 数据不确定性
        data_results = self.uncertainty_results['data']
        report.append("1. 数据不确定性 (阻力面)")
        report.append(f"   平均变异系数: {data_results['std'].mean() / data_results['mean'].mean():.2%}")
        report.append("")

        # 参数不确定性
        param_results = self.uncertainty_results['parameter']
        report.append("2. 参数不确定性 (源地质量)")
        report.append(f"   平均变异系数: {param_results['std'].mean() / param_results['mean'].mean():.2%}")
        report.append("")

        # 情景不确定性
        scenario_results = self.uncertainty_results['scenario']
        report.append("3. 情景不确定性")
        for name, result in scenario_results.items():
            report.append(f"   {name}: {result['description']}")
            report.append(f"     平均连通性: {result['connectivity'].mean():.2f}")
        report.append("")

        # 建议
        summary = self.uncertainty_results['summary']
        report.append("4. 建议")
        for i, rec in enumerate(summary['recommendations'], 1):
            report.append(f"   {i}. {rec}")

        report.append("")
        report.append("=" * 60)

        return "\n".join(report)

---

反思与延伸

思考问题

1. 不确定性与风险的区分：什么情况下是真正的风险，什么情况只是不确定性？

2. 可接受的不确定性水平：在实践中，什么样的不确定性水平是可以接受的？

3. 不确定性的传播：多个不确定性因素组合时，是相互放大还是相互抵消？

4. 减少不确定性的成本：何时值得投入资源获取更精确的数据或更复杂的模型？

5. 沟通挑战：如何向决策者有效传达分析结果的不确定性？

延伸阅读

• "Uncertainty Quantification in Predictive Modeling" - 不确定性量化综述
• "Risk Assessment and Decision Analysis" - 风险评估与决策分析
• "Spatial Uncertainty" ( Zhang & Goodchild) - 空间不确定性专门著作
• IPCC不确定性指南 - 气候变化中的不确定性处理实践

---

关键要点

1. 不确定性无处不在：数据、模型、情境都可能引入不确定性

2. 区分不确定性类型：随机性 vs 认知性，需要不同的处理方式

3. 蒙特卡洛是通用工具：通过随机采样估计输出分布

4. 敏感性分析识别关键因素：找出对结果影响最大的输入

5. 稳健性比精确性更重要：在不确定条件下寻找稳健的解决方案