# 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) │ │ ┌─────────────────────────────────────────────────┐ │ │ │ - 定义多个合理情景 │ │ │ │ - 比较情景结果 │ │ │ │ - 识别稳健策略 │ │ │ │ │ │ │ │ 优点: 直观、易于沟通 │ │ │ │ 缺点: 情景选择主观 │ │ │ └─────────────────────────────────────────────────┘ │ │ │ └─────────────────────────────────────────────────────────────┘ ``` --- ## 设计原理 ### 蒙特卡洛模拟 蒙特卡洛是最通用的不确定性量化方法： ```python 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指数基于方差分解，是全局敏感性分析的金标准： ```python 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 ``` ### 不确定性可视化 空间不确定性需要特别的可视化方法： ```python 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 ``` --- ## 代码示例 ### 阻力面不确定性分析 ```python """ 阻力面不确定性量化示例 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() ``` ### 情景分析框架 ```python """ 情景分析框架 处理情境不确定性 """ 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在生态网络分析中系统性地应用不确定性量化： ```python 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. **稳健性比精确性更重要**：在不确定条件下寻找稳健的解决方案