Sobol 全局敏感性

下面给出「Sobol 全局敏感性」完整教学模块，本科零门槛：

• • 学生在前端勾选参数 + 指标 → 平台自动运行 Sobol 序列采样 + 方差分解

• • 输出：

  • • 一阶敏感度 S1（每个参数单独贡献）

  • • 总阶敏感度 ST（参数 + 交互贡献）

  
  

• • 图表：

  • • 分组柱状图（S1 vs ST）

  • • 置信区间（bootstrap）

  
  

• • 无需写代码，点击即运行

✅ 技术路线（后端）

✅ 安装依赖

pip install SALib pandas numpy joblib

✅ 前端页面（templates/sobol.html）

<!DOCTYPE html>
<html lang="zh">
<head>
   <meta charset="UTF-8">
   <title>Sobol 全局敏感性分析</title>
   <script src="https://cdn.jsdelivr.net/npm/chart.js"></script>
   <style>
    body { font-family: Arial; margin: 40px; background: #f9f9f9; }
    .container { max-width: 1000px; margin: auto; background: #fff; padding: 30px; border-radius: 8px; box-shadow: 0 0 10px rgba(0,0,0,0.1); }
    h1 { text-align: center; color: #2e7d32; }
    .row { display: flex; gap: 20px; align-items: center; margin-bottom: 20px; }
    select[multiple] { width: 200px; height: 160px; }
    button { padding: 10px 20px; background: #2e7d32; color: white; border: none; border-radius: 4px; cursor: pointer; }
    button:hover { background: #1b5e20; }
    .chart-container { margin-top: 40px; }
   </style>
</head>
<body>
   <div class="container">
    <h1>🧬 Sobol 全局敏感性分析</h1>
 
    <div class="row">
      <div>
        <label>选择参数（多选）：</label><br>
        <select id="params" multiple>
          <option value="MAXLUE" selected>MAXLUE</option>
          <option value="SPAN" selected>SPAN</option>
          <option value="SLA">SLA</option>
          <option value="AMAXTB">AMAXTB</option>
        </select>
      </div>
      <div>
        <label>选择指标（单选）：</label><br>
        <select id="metrics">
          <option value="LAI" selected>LAI</option>
          <option value="TWSO">产量</option>
          <option value="TAGP">总生物量</option>
          <option value="WUE">水分利用效率</option>
        </select>
      </div>
    </div>
 
    <div>
      <label>采样数（N）：</label>
      <input type="number" id="n_sample" value="512" min="256" max="2048">
      <small>（N=512 约运行 30 秒）</small>
    </div>
 
    <button onclick="runSobol()">开始 Sobol 分析</button>
 
    <div class="chart-container">
      <canvas id="sobolChart"></canvas>
    </div>
   </div>
 
   <script>
    const ctx = document.getElementById('sobolChart').getContext('2d');
    let chart = new Chart(ctx, {
      type: 'bar',
      data: { labels: [], datasets: [] },
      options: {
        responsive: true,
        plugins: {
          title: { display: true, text: 'Sobol 敏感度（S1 vs ST）' },
            tooltip: {
              callbacks: {
                label: function(context) {
                  return context.dataset.label + ': ' + context.parsed.y.toFixed(3);
                }
              }
            }
        },
        scales: {
          y: { title: { display: true, text: '敏感度' }, beginAtZero: true }
        }
      }
    });
 
    async function runSobol() {
      const params = Array.from(document.getElementById('params').selectedOptions).map(o => o.value);
      const metric = document.getElementById('metrics').value;
      const n = parseInt(document.getElementById('n_sample').value);
      if (params.length === 0) { alert('至少选择一个参数'); return; }
 
      const response = await fetch('/sobol', {
        method: 'POST',
        headers: { 'Content-Type': 'application/json' },
        body: JSON.stringify({ params: params, metric: metric, n: n })
      });
 
      const data = await response.json();
 
      chart.data.labels = params;
      chart.data.datasets = [
        {
          label: '一阶敏感度 S1',
          data: data.s1,
          backgroundColor: 'rgba(75, 192, 192, 0.8)'
        },
        {
          label: '总阶敏感度 ST',
          data: data.st,
          backgroundColor: 'rgba(255, 159, 64, 0.8)'
        }
      ];
      chart.update();
    }
   </script>
</body>
</html>

✅ 后端接口（app.py 新增）

import numpy as np
import pandas as pd
import copy
from SALib.sample import saltelli
from SALib.analyze import sobol
from pcse.fileinput import YAMLCropDataProvider
from pcse.base import ParameterProvider
from pcse.db import NASAPowerWeatherDataProvider
from pcse.models import Wofost71_WLP_FD
from pcse.engine import Engine
from flask import Blueprint, request, jsonify
from joblib import Parallel, delayed
 
sobol_bp = Blueprint('sobol', __name__, url_prefix='/')
 
@sobol_bp.route('/sobol')
def sobol_page():
    return render_template('sobol.html')
 
# ---------- 模型运行封装 ----------
def run_wofost_once(crop_dict, metric):
    cropd = YAMLCropDataProvider(crop_dict=crop_dict)
    weather = NASAPowerWeatherDataProvider(latitude=52.0, longitude=5.0)
    params = ParameterProvider(cropdata=cropd, soildata=None, sitedata=None)
    wofost = Wofost71_WLP_FD(params, weather, [])
    eng = Engine(wofost)
    eng.run_till_terminate()
    df = pd.DataFrame(eng.get_output())
    if metric == 'LAI':
        return df['LAI'].max()
    elif metric == 'TWSO':
        return df['TWSO'].iloc[-1] if not df.empty else 0
    elif metric == 'TAGP':
        return df['TAGP'].max()
    elif metric == 'WUE':
        et = df['EVAP'].sum() + df['TRAT'].sum()
        twso = df['TWSO'].iloc[-1] if not df.empty else 0
        return twso / et if et else 0
    else:
        return 0
 
# ---------- Sobol 主接口 ----------
@sobol_bp.route('/sobol', methods=['POST'])
def sobol_analysis():
    data = request.json
    params = data['params']
    metric = data['metric']
    n = data['n']
 
    # 参数边界 [0.5, 1.5] 倍数
    problem = {
        'num_vars': len(params),
        'names': params,
        'bounds': [[0.5, 1.5]] * len(params)
    }
 
    # 生成样本
    param_values = saltelli.sample(problem, n, calc_second_order=False)
    base_crop = {
        'VARIETY_NAME': 'SobolWheat',
        'CropType': 'wheat',
        'MAXLUE': 2.5,
        'SPAN': 30,
        'SLA': 20,
        'AMAXTB': [0, 1.0, 1.5, 1.0],
        'KDIFTB': [0, 0.6, 2, 0.6]
    }
 
    def eval_one(X):
        crop = copy.deepcopy(base_crop)
        for i, p in enumerate(params):
            if isinstance(base_crop[p], list):
                crop[p] = [v * X[i] if isinstance(v, float) else v for v in base_crop[p]]
            else:
                crop[p] *= X[i]
        return run_wofost_once(crop, metric)
 
    # 并行运行
    y = Parallel(n_jobs=-1, verbose=0)(delayed(eval_one)(x) for x in param_values)
 
    # 分析
    si = sobol.analyze(problem, np.array(y), calc_second_order=False)
    s1 = si['S1'].tolist()
    st = si['ST'].tolist()
 
    return jsonify({'s1': s1, 'st': st})

✅ 使用流程（学生/教师）

1. 1. 打开：http://localhost:5000/sobol

2. 2. 左侧 勾选参数（多选）

3. 3. 右侧 选择 1 个指标（单选）

4. 4. 选择采样数（N=512 默认）

5. 5. 点击“开始 Sobol 分析”

6. 6. 等待 30–60 秒 → 看到：

   • • 分组柱状图：S1（浅绿） vs ST（橙色）

   • • 柱子越高 → 参数对指标影响越大

   • • ST > S1 → 存在 交互效应

   
   

✅ 教学提示（可打印）

✅ 后续可拓展

• • ➕ 二阶交互图（heatmap）

• • ➕ 下载 Sobol 报告（CSV + PNG）

• • ➕ 教师后台：全班 Sobol 平均图