🏆 本文已收录于专栏：《滚雪球学数学建模（含历年真题）》

本专栏面向数学建模竞赛学习者，系统覆盖真题解析、建模方法、算法实现、论文写作与 AI 辅助建模等核心环节。无论是建模新手，还是备战华为杯、高教社杯、华数杯、国赛、美赛 MCM/ICM 的参赛者，都能在这里找到清晰、完整、可复用的建模思路，持续更新，长期有效。

🎯 免责声明： 本文题目来源于互联网公开内容，仅供学习交流与建模方法研究，不构成竞赛指导。请遵守相关赛事规则，独立完成竞赛作品，使用本文内容所产生的后果由使用者自行承担。
  
🎉 专栏限时优惠中：一次订阅，永久解锁，后续内容持续更新。 欢迎点击了解 👉 查看专栏详情 👈

2005年C题：雨量预报方法的评价

真题展示

如下为原(真)题，展示如下：

一、前言：为什么这道题值得深入分析？

拿到这道题的第一反应，很多同学会觉得：“不就是对比两组数据误差吗？算个 RMSE 不就完了？”

这种想法恰恰是建模初学者最常见的误区。

这道题看起来是误差计算问题，但核心其实是评价指标体系的设计问题。 它涉及三个层面的挑战：

第一，空间匹配问题。预报数据是网格点上的值（53×47=2491个格点），而实测数据是91个不均匀分布的站点值。两者坐标系不同，你必须先做空间插值或格点映射，才能比较。这一步处理不好，后续所有分析都是错的。

第二，评价指标设计问题。雨量预报的"好坏"不是一个简单的数字能表达的——你要考虑不同量级的雨（小雨和暴雨被预测错的代价是完全不同的）、要考虑无雨/有雨的命中率、要考虑误报和漏报的权衡。

第三，公众感受的数学化问题。题目第二问要求将公众对降雨等级的感知纳入评价体系，这是典型的"主观因素客观化"建模任务，考察的是对权重设计和分类评价的理解。

我在辅导学生参加数学建模竞赛多年，见过很多团队在这道题上"做了但没做好"——代码跑起来了，表格也有了，但模型逻辑站不住脚，论文得分依然不高。

这篇文章的目标，是帮你从零到一地建立起这道题的完整解题框架，让你不仅能写出代码，更能理解每一步建模决策背后的理由。

二、题目背景与现实意义

2.1 背景认知

雨量预报是数值天气预报（NWP）领域的核心任务之一。不同于温度、气压等连续量的预报，降雨预报有几个特殊性：

1. 空间分布极不均匀：一场局地性强降水可能只影响几平方公里，而相邻区域可能滴雨未下。
2. 量级跨度极大：从不足 0.1mm 的痕量降水到超过 100mm 的特大暴雨，跨越三个数量级。
3. 误报与漏报代价不对称：对于防洪、农业灌溉、城市内涝防治，漏报一次暴雨的代价可能远大于误报一次小雨。

2.2 现实意义

题目背景设定在东经120度、北纬32度附近，对应我国江淮地区（苏北、安徽等地），这一地区恰好是梅雨锋系统频繁活跃的区域，6小时短时强降水事件多发，对预报精度要求极高。

气象部门同时研究两种预报方法，需要一个科学、公正的评价框架来确定哪种方法更优秀，或者在什么气象条件下各自的优势更突出，这正是本题的现实驱动。

2.3 建模价值

这道题的建模价值在于：

• 它是一道典型的评价类建模题，覆盖了空间分析、误差度量、等级化评价、权重设计等多个知识模块；
• 它对数据预处理要求非常高，是训练数据处理能力的好题目；
• 第二问引入公众感受，要求将定性判断转化为定量权重，是主观-客观综合赋权的经典场景。

三、题目重述

注意：题目重述不是题目复制。题目重述的目标是用自己的语言提炼题目中的关键信息、已知条件、目标要求，展示你对题目的理解，而不是把原题抄一遍。

3.1 已知条件

类别	内容
研究区域	东经120°、北纬32°附近，53×47等距网格（共2491个格点）
预报时段	每天20点预报，覆盖21点起的4个连续6小时时段
预报方法	两种不同的6小时雨量预报方法（方法一、方法二）
预报数据	41天 × 4时段 × 2方法的网格预报数据，每个文件含2491个值
实测数据	91个站点的实测雨量，站点分布不均匀
时间跨度	41天（具体日期从数据文件名可知）
单位	毫米（mm），< 0.1mm 视为无雨

3.2 待解决问题

问题一：建立数学模型，定量评价两种6小时雨量预报方法的准确性。

问题二：气象部门将6小时降雨量分为6个等级（见下表），若按等级向公众预报，如何将公众感受纳入评价模型？

等级	名称	雨量范围（mm）
1	小雨	0.1 ~ 2.5
2	中雨	2.6 ~ 6.0
3	大雨	6.1 ~ 12.0
4	暴雨	12.1 ~ 25.0
5	大暴雨	25.1 ~ 60.0
6	特大暴雨	> 60.1

3.3 附件数据说明

数据文件夹结构：
├── FORECAST/
│   ├── lon.dat          → 2491个网格点的经度
│   ├── lat.dat          → 2491个网格点的纬度
│   ├── f6181_dis1       → 2002年6月18日，方法一，第1时段预报
│   ├── f6181_dis2       → 2002年6月18日，方法二，第1时段预报
│   ├── f6182_dis1       → 2002年6月18日，方法一，第2时段预报
│   └── ...（命名规则：f[月日][时段]_dis[方法编号]）
│
└── MEASURING/
    ├── 020618.SIX       → 2002年6月18日，91个站点×4时段实测值
    └── ...（命名规则：[年后两位][月][日].SIX）

文件命名解读：

• f6181_dis1：6月18日，第1时段（21:00~03:00），方法1的预报
• f6183_dis2：6月18日，第3时段（09:00~15:00），方法2的预报
• 020618.SIX：2002年6月18日，从21:00开始的连续4个时段，各站点实测值

四、问题分析

4.1 问题一分析

核心挑战：空间数据的对齐

预报值在网格格点上，实测值在91个不规则站点上——这是本题最大的技术障碍。

你不能直接计算误差，因为两者坐标不一致。常见处理方式有两种：

方案A（推荐）：对每个实测站点，在预报网格中找到最近的格点（或做双线性插值），将格点预报值"采样"到站点位置，然后逐站点计算误差。

方案B：将91个站点的实测值通过插值还原成网格数据（如克里金插值），然后与预报网格逐格点比较。

方案A更常用，计算量小，适合竞赛场景；方案B精度更高，但对稀疏不均匀的站点来说插值误差较大。

需要设计的评价指标：

1. 连续量误差指标：MAE、RMSE、偏差
2. 命中率指标：有雨/无雨的正确识别率（TS评分）
3. 综合评价得分

4.2 问题二分析

核心挑战：如何量化"公众感受"

公众对降雨等级的感受不是线性的：

• 把小雨报成中雨，公众感受差异不大；
• 把暴雨漏报为大雨，可能导致防灾不力，公众感受极差；
• 把无雨报成大暴雨，会引发不必要恐慌。

这提示我们：不同等级之间的误判，代价是不对称的。

建模思路：引入等级混淆惩罚矩阵，对不同程度的预报偏差给予不同的惩罚权重，从而在评价中体现公众的差异化感受。

4.3 问题关系分析

具体相关示意图绘制如下，仅供参考：

逻辑关系：问题一是基础框架，问题二是对问题一的扩展——在连续量误差的基础上，叠加等级分类的"感受权重"，形成更贴近实际应用的评价体系。

五、整体建模思路

5.1 建模路线

具体相关示意图绘制如下，仅供参考：

5.2 模型选择依据

模型	选择理由
MAE/RMSE	直观反映预报误差大小，是气象评估的标准指标
TS评分（威胁评分）	国际通用的降水预报检验方法，适合有雨/无雨判断
混淆矩阵	可展示各等级预报与实测的对应关系，直观清晰
等级惩罚矩阵	量化公众感受的不对称性，回答问题二
综合加权评分	将多维指标合并为单一分数，便于两种方法的直接比较

5.3 算法实现思路

1. 最近邻插值：对每个实测站点，找到欧氏距离最近的预报格点；
2. 逐站点逐时段计算各类误差；
3. 等级映射：将连续雨量值映射为1~6等级（或0代表无雨）；
4. 混淆矩阵构建；
5. 惩罚矩阵加权；
6. 综合评分。

5.4 结果验证方法

• 对两种方法分别计算所有指标，做横向对比；
• 对41天数据做日期维度分析，观察哪种方法在哪类天气（强降水/微量降水）下表现更好；
• 对关键参数（如等级惩罚权重）做灵敏度扰动，验证结论稳定性。

六、数据预处理

6.1 数据读取

由于本题原始数据文件来自竞赛附件（C2005Data.rar），以下代码基于真实文件格式设计，可在获取数据后直接运行。

数据文件结构分析

FORECAST文件：每个 f日期i_dis方法 文件包含 2491 个浮点数，对应 53×47 网格的雨量预报值（按行优先排列）。

MEASURING文件：每个 .SIX 文件是纯文本表格，第一行为表头，后续每行为一个站点的7列数据（站号、纬度、经度、4个时段实测值）。

6.2 缺失值处理

实测数据中，部分站点在某些时段可能出现负值（仪器故障）或空缺，处理策略：

• 负值：直接视为无效，置 NaN；
• 无雨阈值：< 0.1mm 统一置为 0；
• 缺失站点：该站点当日该时段数据标记为缺失，不参与误差统计。

6.3 异常值识别

对于实测雨量 > 200mm/6h 的值（历史上极为罕见），应人工检查是否为记录错误。在建模层面，可设置上限阈值进行 Winsorize 处理或标记。

6.4 标准化处理

由于不同时段、不同日期的降水量级差异很大，不建议对雨量做全局标准化（会破坏量级信息）。评价指标的计算应保留原始量纲（mm）。

6.5 可视化分析

在正式建模前，建议画出以下图：

1. 91个站点的地理分布散点图；
2. 预报网格的范围与站点分布叠加图；
3. 41天总降水量的时间序列（实测 vs 方法一 vs 方法二均值）。

七、模型假设

在建立评价模型之前，需要明确以下假设（假设不是越多越好，而是要每条假设都有合理依据）：

假设1：预报数据和实测数据在时间上严格对应，即文件日期标签准确无误，无时区混乱问题。

依据：题目明确说明文件命名规则，可直接按规则解析。

假设2：空间插值采用最近邻方法时，格点间距足够小（约0.1°≈11km），插值引入的误差在可接受范围内。

依据：研究区域范围约5.3°×4.7°，共53×47格点，格间距约0.1°，相对于实测站点的代表性范围合理。

假设3：91个实测站点的测量值具有代表性，可反映其周边区域的实际雨量（忽略微尺度变异）。

假设4：对公众而言，相邻等级之间的误判惩罚低于跨等级误判，惩罚随等级差距递增。

依据：符合认知心理学中关于感知差异的基本规律。

假设5：41天数据量足以对两种预报方法的统计特性做出可靠估计。

假设6：在评价模型中，不同时段（第1～4段）具有同等权重，除非题目另有说明。

八、符号说明

符号	含义	单位
$N$	实测站点总数（= 91）	个
$T$	评价时段总数（41天 × 4时段 = 164）	个
$D$	预报网格点总数（= 53×47 = 2491）	个
${\hat{r}}_{i,t}^{(k)}$	方法 $k$（$k=1,2$）在第 $t$ 时段对站点 $i$ 附近格点的预报雨量	mm
$r_{i,t}$	第 $t$ 时段站点 $i$ 的实测雨量	mm
$e_{i,t}^{(k)}$	方法 $k$ 在站点 $i$、时段 $t$ 的预报误差 $=\hat{r}\ast {i,t}^{(k)}-r\ast i,t$	mm
$L_k$	将雨量 $k$ 映射到降水等级（0~6）的函数	—
${\hat{g}}_{i,t}^{(k)}$	方法 $k$ 对站点 $i$、时段 $t$ 的预报等级	—
$g_{i,t}$	站点 $i$、时段 $t$ 的实测等级	—
$W$	等级惩罚矩阵，$W_{ab}$ 表示预报为等级 $a$、实测为等级 $b$ 时的惩罚	—
${\text{TS}}^{(k)}$	方法 $k$ 的威胁评分（Threat Score）	—
${\text{MAE}}^{(k)}$	方法 $k$ 的平均绝对误差	mm
${\text{RMSE}}^{(k)}$	方法 $k$ 的均方根误差	mm
${\text{BIAS}}^{(k)}$	方法 $k$ 的系统偏差	mm
$S^{(k)}$	方法 $k$ 的综合评价得分	—
$\alpha ,\beta ,\gamma$	综合评价中各分项指标的权重	—

九、模型一：基础误差评价模型

9.1 模型思想

模型一的核心思想是：用统计误差指标直接量化预报值与实测值的差距。这是气象评估中最基础也最直观的方法。

在建立误差模型之前，必须先解决空间匹配问题。我们采用最近邻插值，将每个实测站点映射到最近的预报格点。

9.2 空间匹配的数学表达

设第 $j$ 个站点的坐标为 $({\text{lon}}_j,{\text{lat}}_j)$，预报网格中第 $d$ 个格点的坐标为 $({\text{Lon}}_d,{\text{Lat}}_d)$，则最近格点索引为：

$$d^{\ast }(j)=\arg \underset{d\in 1,\dots ,D}{\min }\sqrt{({\text{lon}}_j-{\text{Lon}}_d{)}^2+({\text{lat}}_j-{\text{Lat}}_d{)}^2}$$

该公式使用欧氏距离（近似，适用于经纬度变化范围较小的区域）找到距离最近的格点，将该格点的预报值作为站点 $j$ 处的预报值。

9.3 平均绝对误差（MAE）

$${\text{MAE}}^{(k)}=\frac{1}{N\cdot T}\sum _{i=1}^N\sum _{t=1}^T|\hat{r}\ast {i,t}^{(k)}-r\ast i,t|$$

• 含义：所有站点、所有时段预报误差绝对值的平均值，单位为 mm；
• 优点：直观、不受极端值过度影响；
• 适用场景：适合作为基础误差指标，越小越好。

9.4 均方根误差（RMSE）

$${\text{RMSE}}^{(k)}=\sqrt{\frac{1}{N\cdot T}\sum _{i=1}^N\sum _{t=1}^T{\left(\hat{r}\ast {i,t}^{(k)}-r\ast i,t\right)}^2}$$

• 含义：误差的均方根，对大误差有更高的惩罚；
• 优点：对大误差（如暴雨漏报）更敏感，在雨量预报评价中非常重要；
• 物理意义：RMSE > MAE 时，说明存在较多大误差事件。

9.5 系统偏差（BIAS）

$${\text{BIAS}}^{(k)}=\frac{1}{N\cdot T}\sum _{i=1}^N\sum _{t=1}^T\left(\hat{r}\ast {i,t}^{(k)}-r\ast i,t\right)$$

• 含义：正值表示系统性高估（预报偏多），负值表示系统性低估（预报偏少）；
• 重要性：BIAS 揭示了预报方法的系统性倾向，MAE/RMSE 无法区分"高估"和"低估"。

9.6 威胁评分（TS Score）

威胁评分是气象学中广泛使用的降水预报检验指标，专门针对"有雨/无雨"的二分类判断：

设阈值 $\theta$（通常取 0.1mm），定义：

• $\text{NA}$：预报有雨（$\hat{r}\ge \theta$）且实测有雨（$r\ge \theta$）→ 命中
• $\text{NB}$：预报有雨但实测无雨 → 空报
• $\text{NC}$：预报无雨但实测有雨 → 漏报

$${\text{TS}}^{(k)}=\frac{\text{NA}}{\text{NA}+\text{NB}+\text{NC}}$$

• 范围：$[0,1]$，值越大表示预报越准；
• 优点：同时考虑空报和漏报，避免只统计命中率时被大量"无雨正确预报"稀释。

9.7 相关系数

$${\rho }^{(k)}=\frac{{\sum }_{i,t}(\hat{r}\ast {i,t}^{(k)}-{\bar{\hat{r}}}^{(k)})(r\ast i,t-\bar{r})}{\sqrt{{\sum }_{i,t}(\hat{r}\ast {i,t}^{(k)}-{\bar{\hat{r}}}^{(k)}{)}^2\cdot \sum \ast i,t(r_{i,t}-\bar{r}{)}^2}}$$

反映预报值与实测值的线性相关性，越接近1越好。

9.8 MATLAB 实现

data_preprocess.m — 数据预处理函数

function [forecast1, forecast2, measured, lon_grid, lat_grid, lon_sta, lat_sta] = data_preprocess(forecast_dir, measure_dir)
% data_preprocess.m
% 功能：读取预报数据和实测数据，返回对齐后的数组
% 输入：
%   forecast_dir - 预报文件夹路径（如 './FORECAST/'）
%   measure_dir  - 实测文件夹路径（如 './MEASURING/'）
% 输出：
%   forecast1 - 方法一预报值 [N_stations × N_time_periods]
%   forecast2 - 方法二预报值 [N_stations × N_time_periods]
%   measured  - 实测值       [N_stations × N_time_periods]
%   lon_grid, lat_grid - 2491个格点的经纬度（列向量）
%   lon_sta, lat_sta   - 91个站点的经纬度（列向量）

%% ===== 第一步：读取格点经纬度 =====
lon_grid = load(fullfile(forecast_dir, 'lon.dat'));  % 2491×1
lat_grid = load(fullfile(forecast_dir, 'lat.dat'));  % 2491×1

%% ===== 第二步：获取所有实测文件列表 =====
six_files = dir(fullfile(measure_dir, '*.SIX'));
n_days = length(six_files);  % 应为41天

% 初始化输出变量（先从第一个文件获取站点数）
first_file = fullfile(measure_dir, six_files(1).name);
raw = readmatrix(first_file, 'NumHeaderLines', 1);
N = size(raw, 1);  % 站点数（91）
lon_sta = raw(:, 3);  % 经度（第3列）
lat_sta = raw(:, 2);  % 纬度（第2列）

n_periods = 4;  % 每天4个时段
T = n_days * n_periods;  % 总时段数

forecast1 = zeros(N, T);  % 方法一预报，行=站点，列=时段
forecast2 = zeros(N, T);  % 方法二预报
measured  = zeros(N, T);  % 实测值

%% ===== 第三步：空间匹配（最近邻插值）=====
% 找到每个站点对应的最近预报格点索引
nearest_idx = zeros(N, 1);
for j = 1:N
    dist = sqrt((lon_grid - lon_sta(j)).^2 + (lat_grid - lat_sta(j)).^2);
    [~, nearest_idx(j)] = min(dist);
end

%% ===== 第四步：逐天逐时段读取数据 =====
for d = 1:n_days
    fname = six_files(d).name;
    % 解析日期：文件名格式为 YYMMDD.SIX，例如 020618.SIX
    date_str = fname(1:6);  % '020618'
    month_str = date_str(3:4);  % '06'
    day_str   = date_str(5:6);  % '18'
    
    % 读取实测数据
    meas_file = fullfile(measure_dir, fname);
    meas_data = readmatrix(meas_file, 'NumHeaderLines', 1);
    % meas_data 格式：[站号, 纬度, 经度, 第1段, 第2段, 第3段, 第4段]
    
    for p = 1:n_periods  % p = 时段编号 1~4
        col_idx = (d-1)*n_periods + p;  % 对应全局时段列索引
        
        % 读取实测值（第4~7列对应4个时段）
        meas_values = meas_data(:, 3 + p);
        meas_values(meas_values < 0.1) = 0;  % 低于阈值视为无雨
        measured(:, col_idx) = meas_values;
        
        % 构造预报文件名
        % 格式：f[月日][时段]_dis[方法]
        forecast_basename1 = sprintf('f%s%s%d_dis1', month_str, day_str, p);
        forecast_basename2 = sprintf('f%s%s%d_dis2', month_str, day_str, p);
        
        fpath1 = fullfile(forecast_dir, forecast_basename1);
        fpath2 = fullfile(forecast_dir, forecast_basename2);
        
        % 读取预报数据（2491个值）
        if exist(fpath1, 'file')
            fc1_all = load(fpath1);  % 2491×1 向量
            fc1_all(fc1_all < 0) = 0;
            forecast1(:, col_idx) = fc1_all(nearest_idx);  % 按最近格点取值
        else
            forecast1(:, col_idx) = NaN;  % 文件不存在时置NaN
            warning('预报文件不存在：%s', fpath1);
        end
        
        if exist(fpath2, 'file')
            fc2_all = load(fpath2);
            fc2_all(fc2_all < 0) = 0;
            forecast2(:, col_idx) = fc2_all(nearest_idx);
        else
            forecast2(:, col_idx) = NaN;
            warning('预报文件不存在：%s', fpath2);
        end
        
    end
end

fprintf('数据预处理完成。站点数=%d，总时段数=%d\n', N, T);
end

代码解析：

1. 这段代码解决什么问题：完成了从原始文件到结构化矩阵的全部数据读取和空间匹配工作；
2. 为什么这样写：将预处理封装为独立函数，主程序调用后得到干净的数据矩阵，便于后续模型函数复用；
3. 和数学模型的对应：nearest_idx 实现了 $d^{\ast }(j)$ 的计算；forecast1(:, col_idx) = fc1_all(nearest_idx) 实现了 ${\hat{r}}_{i,t}^{(1)}$ 的提取；
4. 初学者注意：readmatrix 的 NumHeaderLines 参数必须正确，否则会把表头读成数据；文件名中的月日必须是两位数（6月→’06’，18日→’18’），注意用 %02s 格式；
5. 竞赛改进建议：可以增加文件存在性检查的统计，输出"共找到XXX个有效文件对"，便于调试。

build_model.m — 误差指标计算函数

function metrics = build_model(forecast, measured)
% build_model.m
% 功能：计算给定预报矩阵与实测矩阵之间的各类误差指标
% 输入：
%   forecast - 预报值矩阵 [N × T]，N为站点数，T为总时段数
%   measured - 实测值矩阵 [N × T]
% 输出：
%   metrics - 包含各指标的结构体

% 去除NaN行（有缺失数据的时段）
valid_mask = ~isnan(forecast) & ~isnan(measured);

fc_valid = forecast(valid_mask);  % 展平为列向量
ms_valid = measured(valid_mask);

n_valid = length(fc_valid);

%% ===== 平均绝对误差 MAE =====
metrics.MAE = mean(abs(fc_valid - ms_valid));

%% ===== 均方根误差 RMSE =====
metrics.RMSE = sqrt(mean((fc_valid - ms_valid).^2));

%% ===== 系统偏差 BIAS =====
metrics.BIAS = mean(fc_valid - ms_valid);

%% ===== 相关系数 =====
metrics.CORR = corr(fc_valid, ms_valid);

%% ===== 威胁评分 TS =====
threshold = 0.1;  % 有雨阈值（mm）
fc_rain = fc_valid >= threshold;  % 预报有雨
ms_rain = ms_valid >= threshold;  % 实测有雨

NA = sum(fc_rain & ms_rain);    % 命中（Hit）
NB = sum(fc_rain & ~ms_rain);   % 空报（False Alarm）
NC = sum(~fc_rain & ms_rain);   % 漏报（Miss）
ND = sum(~fc_rain & ~ms_rain);  % 正确无雨（Correct Null）

metrics.NA = NA;
metrics.NB = NB;
metrics.NC = NC;
metrics.ND = ND;

if (NA + NB + NC) > 0
    metrics.TS = NA / (NA + NB + NC);
else
    metrics.TS = NaN;
end

% 命中率 POD（Probability of Detection）
if (NA + NC) > 0
    metrics.POD = NA / (NA + NC);
else
    metrics.POD = NaN;
end

% 空报率 FAR（False Alarm Ratio）
if (NA + NB) > 0
    metrics.FAR = NB / (NA + NB);
else
    metrics.FAR = NaN;
end

% 偏频 BIAS_freq = 频率偏差（预报有雨频次/实测有雨频次）
if (NA + NC) > 0
    metrics.BIAS_freq = (NA + NB) / (NA + NC);
else
    metrics.BIAS_freq = NaN;
end

fprintf('MAE=%.4f mm, RMSE=%.4f mm, BIAS=%.4f mm\n', ...
    metrics.MAE, metrics.RMSE, metrics.BIAS);
fprintf('TS=%.4f, POD=%.4f, FAR=%.4f\n', ...
    metrics.TS, metrics.POD, metrics.FAR);
end

代码解析：

这段代码实现了模型一中全部核心误差指标的计算。valid_mask 的处理非常重要——竞赛中经常因为忽略 NaN 值导致计算结果为 NaN，务必在每个指标计算前过滤无效值。TS 评分的三个分量（NA、NB、NC）建议同时输出到结果结构体中，便于后续分析空报/漏报的来源。

十、模型二：基于等级分类的改进评价模型

10.1 基础模型的不足

模型一（连续量误差）存在以下问题：

1. 对量级不敏感：预报0.5mm实测0mm，误差0.5mm；预报15mm实测14.5mm，误差也是0.5mm。但前者是"无雨→小雨"的误判，后者是同一量级内的精确偏差，现实危害完全不同；
2. RMSE对大雨过于敏感：极端强降水事件（暴雨/大暴雨）的误差会主导整体RMSE，可能掩盖对小雨的评价；
3. 未区分误判方向：高估和低估的RMSE贡献相同，但实际上漏报暴雨比误报暴雨危害更大。

10.2 改进思路：等级化评价

将连续雨量映射到离散等级，在等级空间中构建混淆矩阵和惩罚加权评分。

10.3 雨量等级映射函数

$$L(r)=\{\begin{array}{llcccccc}0& r<0.11& 0.1\le r\le 2.52& 2.6\le r\le 6.03& 6.1\le r\le 12.04& 12.1\le r\le 25.05& 25.1\le r\le 60.06& r>60.0\end{array}$$

其中 $L(r)$ 将雨量 $r$ 映射到等级 $0,1,2,3,4,5,6$，0代表无雨，1~6对应六个雨量等级。

10.4 混淆矩阵

定义 $7\times 7$ 的混淆矩阵 $C$，其中 $C_{ab}$ 表示实测为等级 $b$、预报为等级 $a$ 的样本数：

$$C_{ab}=\sum _{i=1}^N\sum _{t=1}^{T}1\left[L(\hat{r}\ast {i,t}^{(k)})=a\right]\cdot 1\left[L(r\ast i,t)=b\right]$$

混淆矩阵的对角元素代表预报正确的样本，非对角元素代表各类误判。

10.5 等级准确率

KaTeX parse error: Expected 'EOF', got '_' at position 10: \text{Acc_̲grade}^{(k)} = …

这是最简单的等级评价指标，但它对所有等级同等对待（把"无雨被预报为特大暴雨"和"小雨被预报为中雨"惩罚相同），这不符合公众感受。

10.6 MATLAB 实现

function [C, grade_acc, ts_by_grade] = grade_evaluation(forecast, measured)
% grade_evaluation.m
% 功能：计算等级混淆矩阵和分等级TS评分
% 输入：forecast, measured 均为 [N × T] 矩阵（单位mm）
% 输出：
%   C            - 7×7 混淆矩阵（等级0~6）
%   grade_acc    - 等级总准确率
%   ts_by_grade  - 7×1向量，各等级的TS评分

%% ===== 等级映射函数（嵌套函数）=====
grade = @(r) (r < 0.1)*0 + (r >= 0.1 & r <= 2.5)*1 + ...
             (r > 2.5  & r <= 6.0)*2 + ...
             (r > 6.0  & r <= 12.0)*3 + ...
             (r > 12.0 & r <= 25.0)*4 + ...
             (r > 25.0 & r <= 60.0)*5 + ...
             (r > 60.0)*6;

%% ===== 将所有数据转换为等级 =====
fc_grade = grade(forecast);   % N×T 等级矩阵（预报）
ms_grade = grade(measured);   % N×T 等级矩阵（实测）

% 去除NaN
valid_mask = ~isnan(forecast) & ~isnan(measured);
fc_g = fc_grade(valid_mask);
ms_g = ms_grade(valid_mask);

%% ===== 构建 7×7 混淆矩阵 =====
C = zeros(7, 7);  % 行=预报等级, 列=实测等级
for a = 0:6
    for b = 0:6
        C(a+1, b+1) = sum(fc_g == a & ms_g == b);
    end
end

%% ===== 等级总准确率 =====
grade_acc = sum(diag(C)) / sum(C(:));

%% ===== 分等级 TS 评分 =====
ts_by_grade = zeros(7, 1);
for g = 0:6
    % 对等级g的二元分类：是/否属于该等级
    hit    = C(g+1, g+1);                          % NA
    f_alarm = sum(C(g+1, :)) - C(g+1, g+1);       % NB（预报g但实测非g）
    miss   = sum(C(:, g+1)) - C(g+1, g+1);        % NC（实测g但预报非g）
    
    denom = hit + f_alarm + miss;
    if denom > 0
        ts_by_grade(g+1) = hit / denom;
    else
        ts_by_grade(g+1) = NaN;
    end
end

fprintf('等级准确率: %.2f%%\n', grade_acc * 100);
fprintf('分等级TS评分：');
fprintf(' %.3f', ts_by_grade);
fprintf('\n');
end

代码解析：

等级映射函数使用了 MATLAB 的逻辑运算符链式乘法，这种写法在矩阵输入时可以向量化运行，速度远快于 for 循环遍历。混淆矩阵用双重 for 循环构建，逻辑清晰但效率一般——在 N×T 较大时（本题约 91×164 ≈ 15000 个样本），速度可接受。分等级 TS 评分将每个等级视为二分类问题，这是气象领域标准做法，建议在论文表格中展示每个等级的 TS 值，便于找出两种方法各自的优势区间。

10.7 对比分析

模型一与模型二的对比：

维度	模型一（连续量误差）	模型二（等级化评价）
信息粒度	精细（mm级）	粗粒（6个等级）
对极端值敏感性	RMSE极高	一个等级算一次误判
符合气象实践	适合科研	适合业务预报
符合公众感知	不直接符合	较好符合
是否能区分误判方向	BIAS可以	混淆矩阵可以

十一、模型三：基于公众感受的加权综合评价模型

11.1 综合建模目标

回答问题二的核心：如何将公众的感受量化为评价权重？

公众对降水等级误判的感受遵循以下规律：

1. 误判的等级差越大，感受越差；
2. 漏报高等级降水（如将暴雨报为中雨）比空报高等级降水危害更大；
3. 相邻等级的误判（如将小雨报为中雨）感受差异不大。

11.2 等级惩罚矩阵设计

定义 $7\times 7$ 的惩罚矩阵 $W$，$W_{ab}$ 表示实测为等级 $b$、预报为等级 $a$ 时的惩罚分数：

基础对称惩罚（只考虑等级差距）：

$$W_{ab}^{(0)}=|a-b|$$

非对称惩罚（考虑漏报高等级比空报惩罚更重）：

$$W_{ab}=\{\begin{array}{llc}|a-b{|}^{\alpha }& a\le b\text{（漏报：预报值偏低）}\beta \cdot |a-b{|}^{\alpha }& a>b\text{（空报：预报值偏高）}\end{array}$$

其中 $\alpha >1$ 使惩罚随等级差非线性增长，$\beta <1$ 使空报惩罚轻于漏报（从防灾角度，空报可接受，漏报不可接受）。

典型参数选择：$\alpha =1.5$，$\beta =0.7$（可通过专家咨询或灵敏度分析调整）。

11.3 加权评价得分

$$S_{\text{grade}}^{(k)}=1-\frac{{\sum }_{a=0}^6{\sum }_{b=0}^6{W}_{ab}\cdot C_{ab}^{(k)}}{{\sum }_{a=0}^6{\sum }_{b=0}^6{C}_{ab}^{(k)}\cdot W_{ab}^{\max }}$$

其中 $W_{ab}^{\max }$ 是归一化因子（最大可能惩罚），确保 $S_{\text{grade}}^{(k)}\in [0,1]$，值越大表示预报越好。

11.4 综合评价得分

融合模型一（连续误差）和模型三（等级加权）：

$$S^{(k)}=w_1\cdot S_{\text{cont}}^{(k)}+w_2\cdot S_{\text{grade}}^{(k)}$$

其中 $S_{\text{cont}}^{(k)}$ 为连续量评分（需归一化），$w_1+w_2=1$，可根据评价目的调整权重（科研导向则 $w_1$ 大，业务导向则 $w_2$ 大）。

对连续量误差的归一化评分：

$$S_{\text{cont}}^{(k)}=1-\frac{{\text{RMSE}}^{(k)}}{{\text{RMSE}}_{\text{climatology}}}$$

其中 ${\text{RMSE}}_{\text{climatology}}$ 是以历史气候均值作为预报值时的 RMSE，作为基准线。

11.5 MATLAB 实现

function [S_weighted, W_penalty] = weighted_grade_score(forecast, measured, alpha, beta)
% weighted_grade_score.m
% 功能：计算基于公众感受惩罚矩阵的加权等级评分
% 输入：
%   forecast, measured - 数据矩阵 [N × T]
%   alpha - 惩罚非线性系数（默认1.5）
%   beta  - 空报相对漏报的惩罚折扣（默认0.7）
% 输出：
%   S_weighted - 综合加权评分（0~1，越大越好）
%   W_penalty  - 7×7 惩罚矩阵

if nargin < 3, alpha = 1.5; end
if nargin < 4, beta  = 0.7; end

%% ===== 构建惩罚矩阵 =====
n_levels = 7;  % 等级0~6，共7级
W_penalty = zeros(n_levels, n_levels);
for a = 0:6
    for b = 0:6
        level_diff = abs(a - b);
        if a <= b
            % 漏报（预报偏低），惩罚较重
            W_penalty(a+1, b+1) = level_diff^alpha;
        else
            % 空报（预报偏高），惩罚折扣
            W_penalty(a+1, b+1) = beta * level_diff^alpha;
        end
    end
end
% 对角线（预报正确）惩罚为0，已自动满足

%% ===== 获取混淆矩阵 =====
[C, ~, ~] = grade_evaluation(forecast, measured);

%% ===== 计算总加权惩罚 =====
total_penalty = sum(sum(W_penalty .* C));
total_samples = sum(C(:));

% 计算理论最大惩罚（如果每个样本都被预报为最远等级）
% 即惩罚矩阵中每列最大惩罚乘以该实测等级的样本数
max_penalty = 0;
for b = 0:6
    col_samples = sum(C(:, b+1));  % 该实测等级的总样本数
    col_max_w = max(W_penalty(:, b+1));  % 该等级被最差预报的惩罚
    max_penalty = max_penalty + col_samples * col_max_w;
end

%% ===== 计算加权评分 =====
if max_penalty > 0
    S_weighted = 1 - total_penalty / max_penalty;
else
    S_weighted = 1;
end

fprintf('惩罚矩阵参数：alpha=%.2f, beta=%.2f\n', alpha, beta);
fprintf('加权等级评分：S = %.4f\n', S_weighted);

%% ===== 可视化惩罚矩阵 =====
figure;
heatmap(0:6, 0:6, W_penalty, ...
    'Title', 'Grade Penalty Matrix W', ...
    'XLabel', 'Forecast Grade', ...
    'YLabel', 'Observed Grade', ...
    'Colormap', hot, ...
    'ColorbarVisible', 'on');
end

代码解析：

惩罚矩阵的设计是本模型的核心创新点。a <= b 对应"漏报"（预报偏低，如把暴雨报成大雨），a > b 对应"空报"（预报偏高，如把无雨报成小雨）。在防洪、应急场景下，漏报高等级降水的后果更严重，所以不打折扣；空报会引发公众警觉疲劳，但后果相对可控，所以乘以折扣 $\beta <1$。alpha > 1 确保等级差越大惩罚增速越快（非线性），符合公众感受的心理特征。

十二、算法流程设计

具体相关示意图绘制如下，仅供参考：

流程说明：

• 步骤 B~H 为数据层，重点在空间匹配；
• 步骤 J~L 为模型层，从简到繁递进构建；
• 步骤 M~N 为分析层，横向对比与鲁棒性验证；
• 步骤 O 为呈现层，结果可视化。

十三、MATLAB 完整代码

13.1 主程序 main.m

%% main.m
% 2005 CUMCM C题：雨量预报方法评价 主程序
% 运行前请确保以下文件存在于 MATLAB 路径中：
%   data_preprocess.m, build_model.m, grade_evaluation.m,
%   weighted_grade_score.m, plot_results.m, sensitivity_analysis.m

clc;
clear;
close all;

%% ===== 路径设置 =====
forecast_dir = './FORECAST/';   % 预报文件夹路径
measure_dir  = './MEASURING/';  % 实测文件夹路径

%% ===== 第一步：数据预处理 =====
fprintf('===== 数据预处理中... =====\n');
[forecast1, forecast2, measured, lon_grid, lat_grid, lon_sta, lat_sta] = ...
    data_preprocess(forecast_dir, measure_dir);

%% ===== 第二步：模型一 - 连续量误差评价 =====
fprintf('\n===== 模型一：连续量误差评价 =====\n');
fprintf('-- 方法一 --\n');
metrics1 = build_model(forecast1, measured);

fprintf('-- 方法二 --\n');
metrics2 = build_model(forecast2, measured);

%% ===== 第三步：模型二 - 等级化评价 =====
fprintf('\n===== 模型二：等级分类评价 =====\n');
fprintf('-- 方法一 --\n');
[C1, ga1, ts1] = grade_evaluation(forecast1, measured);

fprintf('-- 方法二 --\n');
[C2, ga2, ts2] = grade_evaluation(forecast2, measured);

%% ===== 第四步：模型三 - 公众感受加权评价 =====
fprintf('\n===== 模型三：公众感受加权评价 =====\n');
alpha_param = 1.5;
beta_param  = 0.7;

fprintf('-- 方法一 --\n');
[S1, W_mat] = weighted_grade_score(forecast1, measured, alpha_param, beta_param);

fprintf('-- 方法二 --\n');
[S2, ~] = weighted_grade_score(forecast2, measured, alpha_param, beta_param);

%% ===== 第五步：汇总结果表格 =====
fprintf('\n===== 综合评价结果汇总 =====\n');
fprintf('%-20s %-12s %-12s\n', '指标', '方法一', '方法二');
fprintf('%s\n', repmat('-', 1, 46));
fprintf('%-20s %-12.4f %-12.4f\n', 'MAE (mm)',   metrics1.MAE,   metrics2.MAE);
fprintf('%-20s %-12.4f %-12.4f\n', 'RMSE (mm)',  metrics1.RMSE,  metrics2.RMSE);
fprintf('%-20s %-12.4f %-12.4f\n', 'BIAS (mm)',  metrics1.BIAS,  metrics2.BIAS);
fprintf('%-20s %-12.4f %-12.4f\n', 'CORR',       metrics1.CORR,  metrics2.CORR);
fprintf('%-20s %-12.4f %-12.4f\n', 'TS (全部)',  metrics1.TS,    metrics2.TS);
fprintf('%-20s %-12.4f %-12.4f\n', 'POD',        metrics1.POD,   metrics2.POD);
fprintf('%-20s %-12.4f %-12.4f\n', 'FAR',        metrics1.FAR,   metrics2.FAR);
fprintf('%-20s %-12.2f %-12.2f\n', '等级准确率%', ga1*100, ga2*100);
fprintf('%-20s %-12.4f %-12.4f\n', '加权等级评分', S1, S2);

%% ===== 第六步：可视化 =====
plot_results(metrics1, metrics2, C1, C2, ts1, ts2, W_mat, ...
             forecast1, forecast2, measured, lon_sta, lat_sta);

%% ===== 第七步：灵敏度分析 =====
sensitivity_analysis(forecast1, forecast2, measured);

13.2 数据预处理函数

（已在第九节完整给出，见 data_preprocess.m）

13.3 模型求解函数

（已在第九节和第十节完整给出，见 build_model.m、grade_evaluation.m、weighted_grade_score.m）

13.4 结果可视化函数

function plot_results(m1, m2, C1, C2, ts1, ts2, W_mat, fc1, fc2, ms, lon_sta, lat_sta)
% plot_results.m
% 功能：绘制所有评价结果图

grade_labels = {'No Rain','Light','Moderate','Heavy','Storm','Heavy Storm','Extreme'};

%% ===== 图1：各指标对比柱状图 =====
figure('Name', 'Metric Comparison', 'Position', [100 100 900 500]);
metrics_names = {'MAE', 'RMSE', 'TS', 'POD', 'CORR'};
vals1 = [m1.MAE, m1.RMSE, m1.TS, m1.POD, m1.CORR];
vals2 = [m2.MAE, m2.RMSE, m2.TS, m2.POD, m2.CORR];

subplot(1,2,1);
bar([vals1(1:2); vals2(1:2)]', 'grouped');
set(gca, 'XTickLabel', {'MAE', 'RMSE'});
legend({'Method 1', 'Method 2'}, 'Location', 'northwest');
ylabel('Error (mm)');
title('Continuous Error Metrics');
grid on;

subplot(1,2,2);
bar([vals1(3:5); vals2(3:5)]', 'grouped');
set(gca, 'XTickLabel', {'TS', 'POD', 'CORR'});
legend({'Method 1', 'Method 2'}, 'Location', 'southwest');
ylabel('Score (0-1)');
title('Skill Scores');
grid on;

%% ===== 图2：混淆矩阵热力图 =====
figure('Name', 'Confusion Matrix', 'Position', [100 100 1200 500]);

subplot(1,2,1);
% 将混淆矩阵转为百分比（按列归一化，即已知实测等级下的分布）
C1_norm = C1 ./ (sum(C1, 1) + eps);  % 避免除以0
heatmap(grade_labels, grade_labels, round(C1_norm*100)/100, ...
    'Title', 'Confusion Matrix - Method 1 (%)', ...
    'XLabel', 'Observed Grade', ...
    'YLabel', 'Forecast Grade', ...
    'Colormap', parula);

subplot(1,2,2);
C2_norm = C2 ./ (sum(C2, 1) + eps);
heatmap(grade_labels, grade_labels, round(C2_norm*100)/100, ...
    'Title', 'Confusion Matrix - Method 2 (%)', ...
    'XLabel', 'Observed Grade', ...
    'YLabel', 'Forecast Grade', ...
    'Colormap', parula);

%% ===== 图3：分等级TS评分对比 =====
figure('Name', 'TS by Grade', 'Position', [100 100 700 450]);
bar_data = [ts1, ts2];
% 替换NaN为0用于绘图
bar_data(isnan(bar_data)) = 0;
bar(bar_data, 'grouped');
set(gca, 'XTickLabel', grade_labels);
xtickangle(30);
legend({'Method 1', 'Method 2'});
ylabel('Threat Score (TS)');
title('TS Score by Precipitation Grade');
grid on;
ylim([0, 1]);

%% ===== 图4：惩罚矩阵可视化 =====
figure('Name', 'Penalty Matrix', 'Position', [100 100 600 500]);
heatmap(grade_labels, grade_labels, W_mat, ...
    'Title', 'Public Perception Penalty Matrix', ...
    'XLabel', 'Forecast Grade', ...
    'YLabel', 'Observed Grade', ...
    'Colormap', hot);

%% ===== 图5：预报值 vs 实测值散点图（随机抽样2000个点）=====
figure('Name', 'Scatter Plot', 'Position', [100 100 900 400]);
% 展平数据并随机抽取
fc1_flat = fc1(:);
fc2_flat = fc2(:);
ms_flat  = ms(:);
valid_idx = find(~isnan(fc1_flat) & ~isnan(fc2_flat));
sample_idx = valid_idx(randperm(length(valid_idx), min(2000, length(valid_idx))));

subplot(1,2,1);
scatter(ms_flat(sample_idx), fc1_flat(sample_idx), 5, 'b', 'filled', 'MarkerFaceAlpha', 0.3);
hold on;
max_val = max([ms_flat(sample_idx); fc1_flat(sample_idx)]);
plot([0 max_val], [0 max_val], 'r--', 'LineWidth', 1.5);
xlabel('Observed Rainfall (mm)');
ylabel('Forecast Rainfall (mm)');
title('Method 1: Forecast vs Observed');
grid on;

subplot(1,2,2);
scatter(ms_flat(sample_idx), fc2_flat(sample_idx), 5, 'g', 'filled', 'MarkerFaceAlpha', 0.3);
hold on;
plot([0 max_val], [0 max_val], 'r--', 'LineWidth', 1.5);
xlabel('Observed Rainfall (mm)');
ylabel('Forecast Rainfall (mm)');
title('Method 2: Forecast vs Observed');
grid on;

%% ===== 图6：站点地理分布 =====
figure('Name', 'Station Distribution', 'Position', [100 100 600 500]);
scatter(lon_sta, lat_sta, 50, 'r', 'filled');
xlabel('Longitude (°E)');
ylabel('Latitude (°N)');
title('Observation Station Distribution (91 Stations)');
grid on;
axis equal;

fprintf('所有图表已生成。\n');
end

13.5 灵敏度分析函数

function sensitivity_analysis(forecast1, forecast2, measured)
% sensitivity_analysis.m
% 功能：对惩罚矩阵参数 alpha 和 beta 做灵敏度分析

alpha_range = 1.0 : 0.2 : 2.5;  % alpha 的扰动范围
beta_range  = 0.4 : 0.1 : 1.0;  % beta 的扰动范围

n_alpha = length(alpha_range);
n_beta  = length(beta_range);

S1_matrix = zeros(n_alpha, n_beta);  % 方法一得分矩阵
S2_matrix = zeros(n_alpha, n_beta);  % 方法二得分矩阵
diff_matrix = zeros(n_alpha, n_beta); % 两方法得分之差

fprintf('灵敏度分析进行中，共 %d 个参数组合...\n', n_alpha * n_beta);

for i = 1:n_alpha
    for j = 1:n_beta
        a = alpha_range(i);
        b = beta_range(j);
        % 关闭内部输出
        S1 = weighted_grade_score_silent(forecast1, measured, a, b);
        S2 = weighted_grade_score_silent(forecast2, measured, a, b);
        S1_matrix(i,j) = S1;
        S2_matrix(i,j) = S2;
        diff_matrix(i,j) = S1 - S2;  % 正值表示方法一更优
    end
end

%% ===== 绘制参数热力图 =====
figure('Name', 'Sensitivity Analysis', 'Position', [100 100 1100 400]);

subplot(1,3,1);
heatmap(beta_range, alpha_range, S1_matrix, ...
    'Title', 'Method 1 Weighted Score', ...
    'XLabel', 'beta', 'YLabel', 'alpha', ...
    'Colormap', cool);

subplot(1,3,2);
heatmap(beta_range, alpha_range, S2_matrix, ...
    'Title', 'Method 2 Weighted Score', ...
    'XLabel', 'beta', 'YLabel', 'alpha', ...
    'Colormap', cool);

subplot(1,3,3);
heatmap(beta_range, alpha_range, diff_matrix, ...
    'Title', 'Score Diff (M1 - M2)', ...
    'XLabel', 'beta', 'YLabel', 'alpha', ...
    'Colormap', rdbu);  % 正值蓝色（方法一更优），负值红色

%% ===== 统计结论稳定性 =====
n_m1_better = sum(diff_matrix(:) > 0);
n_m2_better = sum(diff_matrix(:) < 0);
n_tie       = sum(diff_matrix(:) == 0);

fprintf('\n===== 灵敏度分析结论 =====\n');
fprintf('参数组合总数：%d\n', n_alpha * n_beta);
fprintf('方法一更优的参数比例：%.1f%%\n', n_m1_better / (n_alpha*n_beta) * 100);
fprintf('方法二更优的参数比例：%.1f%%\n', n_m2_better / (n_alpha*n_beta) * 100);
fprintf('得分差均值：%.4f，标准差：%.4f\n', mean(diff_matrix(:)), std(diff_matrix(:)));

if n_m1_better > n_m2_better
    fprintf('结论稳健：在%.1f%%的参数组合下，方法一优于方法二。\n', ...
        n_m1_better/(n_alpha*n_beta)*100);
else
    fprintf('结论稳健：在%.1f%%的参数组合下，方法二优于方法一。\n', ...
        n_m2_better/(n_alpha*n_beta)*100);
end
end

%% ===== 辅助函数：静默版加权评分（不打印中间输出）=====
function S = weighted_grade_score_silent(forecast, measured, alpha, beta)
% 与 weighted_grade_score 相同，但不输出中间信息，专用于循环调用

grade = @(r) (r < 0.1)*0 + (r >= 0.1 & r <= 2.5)*1 + ...
             (r > 2.5  & r <= 6.0)*2 + ...
             (r > 6.0  & r <= 12.0)*3 + ...
             (r > 12.0 & r <= 25.0)*4 + ...
             (r > 25.0 & r <= 60.0)*5 + (r > 60.0)*6;

fc_g = grade(forecast(:));
ms_g = grade(measured(:));
valid = ~isnan(forecast(:)) & ~isnan(measured(:));
fc_g = fc_g(valid);
ms_g = ms_g(valid);

W = zeros(7,7);
for a = 0:6
    for b = 0:6
        d = abs(a-b);
        if a <= b
            W(a+1,b+1) = d^alpha;
        else
            W(a+1,b+1) = beta * d^alpha;
        end
    end
end

C = zeros(7,7);
for a = 0:6
    for b = 0:6
        C(a+1,b+1) = sum(fc_g==a & ms_g==b);
    end
end

total_penalty = sum(sum(W .* C));
max_penalty = 0;
for b = 0:6
    col_samples = sum(C(:,b+1));
    col_max_w = max(W(:,b+1));
    max_penalty = max_penalty + col_samples * col_max_w;
end

if max_penalty > 0
    S = 1 - total_penalty / max_penalty;
else
    S = 1;
end
end

代码解析：

灵敏度分析是竞赛论文中非常重要但经常被忽视的环节。这里通过双层 for 循环遍历 alpha 和 beta 的参数空间，记录每种参数组合下两种预报方法的得分差，最终以热力图形式展示。关键结论：如果在几乎所有参数组合下，某一方法都更优，说明结论是鲁棒的；如果结论随参数变化而反转，则需要额外分析原因（可能两种方法各有擅长的量级区间）。注意提供了 weighted_grade_score_silent 的静默版本，避免循环时控制台被大量打印信息淹没。

十四、结果展示与分析

14.1 模拟结果示例（明确标注）

⚠️ 说明：由于原始竞赛数据需从竞赛网站下载，以下为基于题目背景的模拟分析示例，仅用于说明结果呈现方式。真实结果需运行上述代码获取。

连续误差指标示例表（模拟）

指标	方法一	方法二	更优
MAE (mm)	2.134	1.872	方法二
RMSE (mm)	5.871	6.423	方法一
BIAS (mm)	+0.312	-0.089	方法二更中性
CORR	0.683	0.721	方法二
TS	0.452	0.438	方法一
POD	0.613	0.578	方法一
FAR	0.237	0.215	方法二

14.2 结果解读

这个模拟结果揭示了一个非常值得深入分析的现象：两种方法在不同指标下各有优势，没有绝对的胜者。

• 方法一 RMSE 更低，说明它对极端降水的预报偏差更小（RMSE 对大误差敏感）；
• 方法二 MAE 更低，说明它在总体平均误差上表现更好；
• 方法一 TS/POD 更高，说明它对有雨/无雨的识别能力更强；
• 方法二 FAR 更低，说明它空报率更少。

这一发现具有重要现实意义：没有一种预报方法在所有指标上都占优势，这是气象预报领域普遍存在的现象，也是气象部门需要建立综合评价模型的根本原因。

14.3 分等级TS评分分析（关键结论）

通常，对于无雨和低等级降水（小雨、中雨），样本量大，两种方法的 TS 评分都较高；而对于高等级降水（暴雨、大暴雨），样本量极少，TS 评分差异会显著扩大——这恰恰是最有价值的对比区间。

14.4 公众感受加权评分的意义

当采用非对称惩罚矩阵时，结论可能发生变化：某方法在连续误差上表现较差，但因为它很少漏报高等级降水（即便偶有空报），其公众感受评分反而更高。这一现象具有强烈的实践指导意义：对于防灾减灾应用，宁可空报也不能漏报。

十五、模型检验

15.1 误差分析

问：为什么不同时段的误差会有差异？

第 1 时段（21:00~03:00）是预报发出后 1~7 小时内的降水，预报时效最短，误差通常最小；第 4 时段（15:00~21:00）是发出预报后 19~25 小时的降水，时效最长，误差通常最大。建议对4个时段分别计算误差，并绘制"时效-误差"关系图，分析两种方法随预报时效增加而误差增大的速率。

误差分解：

KaTeX parse error: Expected 'EOF', got '_' at position 51: …2} + \text{RMSE_̲centered}^{(k)2…

其中均方误差可以分解为系统误差（偏差的平方）和随机误差两部分，这是评估预报方法的深层指标。

15.2 灵敏度分析

对评价模型中的关键参数进行扰动：

1. 最近邻 vs 双线性插值：改变空间匹配方式，对比两种插值方法对最终评分结论的影响；
2. 有雨阈值敏感性：将 0.1mm 阈值改为 0.05mm 或 0.2mm，观察 TS/POD/FAR 的变化；
3. 惩罚矩阵参数：已在 sensitivity_analysis.m 中实现；
4. 时段权重：将4个时段赋予不同权重（如晚间权重更高），分析对总分的影响。

15.3 稳定性分析

对41天数据进行Bootstrap抽样检验：

% 稳定性分析 - Bootstrap
n_bootstrap = 500;
T_total = size(measured, 2);  % 总时段数（164）
mae1_boot = zeros(n_bootstrap, 1);
mae2_boot = zeros(n_bootstrap, 1);

for b = 1:n_bootstrap
    % 有放回地抽取 T_total 个时段
    sample_idx = randsample(T_total, T_total, true);
    fc1_sample = forecast1(:, sample_idx);
    fc2_sample = forecast2(:, sample_idx);
    ms_sample  = measured(:, sample_idx);
    
    m1_b = build_model(fc1_sample, ms_sample);
    m2_b = build_model(fc2_sample, ms_sample);
    mae1_boot(b) = m1_b.MAE;
    mae2_boot(b) = m2_b.MAE;
end

% 95% 置信区间
ci1 = prctile(mae1_boot, [2.5, 97.5]);
ci2 = prctile(mae2_boot, [2.5, 97.5]);
fprintf('方法一 MAE 95%%CI: [%.3f, %.3f]\n', ci1(1), ci1(2));
fprintf('方法二 MAE 95%%CI: [%.3f, %.3f]\n', ci2(1), ci2(2));

如果两种方法的置信区间不重叠，说明差异显著；若重叠，则需要更多数据才能得出确定结论。

15.4 鲁棒性分析

对91个站点做Leave-One-Out测试：每次去掉一个站点，计算剩余90个站点的评分，观察结论是否随站点变化而发生反转。这一方法可以识别是否存在"关键站点"——某个站点的数据对评价结论有不成比例的影响（这在站点分布不均匀时很常见）。

十六、模型优缺点

16.1 模型一（连续量误差评价）

优点：

• 直观、有物理量纲（mm），便于解释；
• MAE/RMSE 是国际通用标准，与其他研究有可比性；
• 计算简单，不依赖主观参数。

缺点：

• 对量级不敏感（0.5mm误差和15mm误差同样计算）；
• RMSE 被少数强降水事件主导，可能掩盖对常见降水的评价；
• 不能直接反映公众感受。

改进方向：引入对数变换（$\log (r+1)$）再计算误差，削弱强降水的主导效应；或分层计算（分别对小雨、中雨、大雨各层计算误差）。

16.2 模型二（等级化评价）

优点：

• 与气象业务预报直接对应（公众接收的就是等级信息）；
• 混淆矩阵提供完整的误判分布信息；
• 分等级 TS 评分揭示各量级的预报能力。

缺点：

• 等级边界处的误差被放大（如 6.0mm 和 6.1mm 判断不同等级，但实际误差极小）；
• 等级映射丢失了连续量的精度信息；
• 高等级样本量极少，分等级TS评分统计稳定性差。

改进方向：采用模糊等级边界（如用隶属函数代替硬分界）；对高等级评分引入置信区间。

16.3 模型三（公众感受加权评价）

优点：

• 将主观感受通过参数化方法客观化；
• 非对称惩罚符合防灾减灾的实际需要；
• 参数化设计便于针对不同应用场景调整。

缺点：

• 惩罚矩阵参数（$\alpha$, $\beta$）的选取带有主观性，不同参数可能导致不同结论；
• 缺乏实际公众调查数据支撑，参数难以精确标定；
• 不同区域、不同人群对降雨的感受存在差异，统一参数可能不合适。

改进方向：结合公众问卷调查或专家打分确定惩罚权重；引入层次分析法（AHP）建立更系统的权重体系。

十七、论文写作建议

17.1 摘要写法

竞赛摘要是论文最重要的部分（评委很多时候只看摘要），必须包含四要素：背景-方法-结果-结论。

• 背景（1~2句）：点明题目核心问题；
• 方法（3~5句）：列出建立了哪些模型、用了哪些核心指标；
• 结果（2~3句）：给出核心数值结论（不能回避具体数字）；
• 结论（1~2句）：两种方法的比较结论及适用建议。

摘要中必须出现具体数字，例如"方法一的MAE为X.XX mm，RMSE为X.XX mm，TS评分为0.XXX"，不能只写"计算了误差指标"。

17.2 关键词选择

建议关键词：雨量预报评价、威胁评分（TS）、混淆矩阵、等级惩罚矩阵、最近邻插值、综合评价模型

17.3 模型假设写法

模型假设每一条要有假设内容和假设依据，格式如下：

假设1：空间插值采用最近邻方法，认为格点间距（约0.1°）足够小，插值误差可忽略。
依据：研究区经纬度跨度约5°×5°，格间距约11km，与气象尺度的中尺度系统（空间尺度>10km）相匹配。

避免写成"假设数据是可靠的"这种废话式假设。

17.4 符号说明写法

符号表中不能只写符号和名称，还要注明：单位、取值范围（如果有）、在哪个公式中出现。竞赛中有10~15个符号是合理的，太少说明建模深度不够，太多则臃肿难读。

17.5 模型建立写法

每个模型的呈现顺序建议：

1. 为什么选这个模型（1段）
2. 数学定义/公式（带编号的公式）
3. 公式中各符号的含义（用文字解释，不要只靠符号表）
4. 该模型解决了题目中的哪个子问题

常见错误：把 MAE 的公式放出来，后面就直接讲结果了，中间没有说明为什么选 MAE、MAE 在本题中的适用性在哪里。评委会觉得"这个公式是照搬的，作者不理解"。

17.6 结果分析写法

结果分析不能只写"如表3所示，方法一的RMSE为5.87mm，方法二为6.42mm，说明方法一更好"。要进一步写：

• “RMSE较低说明方法一对大误差事件（强降水）的预报偏差更小，这在防汛场景中尤为重要”；
• “但方法二的FAR（空报率）更低，说明当方法二预报有雨时，更可能确实发生降水，这对需要减少不必要防灾响应的部门更有价值”；
• “两种方法各有优势这一现象，反映了降水预报本质上是一个准确性与保守性之间权衡的问题（recall-precision tradeoff）”。

17.7 参考文献写法

竞赛参考文献至少5~8篇，建议包含：

1. 气象领域的降水预报评估经典论文（证明模型选择有依据）；
2. 数学建模教材或综合评价方法参考书；
3. 中国气象局相关技术规范。

十八、数学建模论文摘要示例

摘要

本文针对2005年全国大学生数学建模竞赛C题，建立了一套用于评价6小时雨量预报方法准确性的数学模型体系。

针对问题一，我们首先利用最近邻插值方法，将53×47预报网格中的预报值映射至91个不均匀分布的实测站点，实现了预报数据与实测数据的空间对齐。在此基础上，从连续量误差和分类能力两个维度建立评价模型：连续量误差模型计算了MAE、RMSE、系统偏差和相关系数；分类评价模型引入气象领域通用的威胁评分（TS）、命中率（POD）和空报率（FAR），并构建7×7等级混淆矩阵，分析两种方法在各雨量等级上的预报能力差异。

针对问题二，我们提出了基于公众感受的非对称惩罚矩阵模型。考虑到漏报高等级降水（如将暴雨预报为大雨）对公众安全的危害大于空报，设计了惩罚权重 $W_{ab}=|a-b{|}^{1.5}$（漏报）和 $W_{ab}=0.7\times |a-b{|}^{1.5}$（空报）的不对称惩罚机制，并将其与混淆矩阵结合，计算两种方法的公众感受加权评分。灵敏度分析表明，在惩罚参数 $\alpha \in [1.0,2.5]$、$\beta \in [0.4,1.0]$ 的范围内，评价结论保持稳健。

综合分析结果显示（基于模拟示例），方法一在威胁评分（TS=0.452）和均方根误差（RMSE=5.871mm）上表现更优，适合对强降水事件识别要求较高的防灾场景；方法二在空报率（FAR=0.215）和公众感受加权评分上表现更优，适合面向公众发布的业务预报。本文建立的综合评价框架可为气象部门在不同应用场景下选择预报方法提供定量依据。

关键词：雨量预报评价；威胁评分（TS）；等级混淆矩阵；非对称惩罚矩阵；最近邻插值；综合评价模型

十九、常见问题与踩坑总结

Q1：拿到数学建模题目后为什么不能马上写代码？

因为写代码之前你需要明确：要解决的目标是什么、输入输出是什么、数学模型是什么。很多同学写了两小时代码后才发现——自己算的东西根本没有回答题目的问题。正确顺序是：读题→理解问题→画变量关系图→确定模型→再写代码。本题中，如果你不先想清楚"预报网格和实测站点如何对齐"，写代码就是在无效劳动。

Q2：问题重述和题目复述有什么区别？

题目复述是把原题照抄或稍加改动，这没有任何建模价值，而且通常是丢分项。问题重述是用建模的语言重新表达题目：哪些是已知量？哪些是目标量？哪些是决策变量？有哪些约束？本题的问题重述应该明确指出"需要建立从预报数据到评价得分的映射函数"，而不是仅仅说"评价两种预报方法"。

Q3：模型假设是不是越多越好？

不是。假设的数量不是越多越好，质量才是关键。每条假设应该：（1）有存在的必要（不假设这条，模型无法建立）；（2）有现实依据（不是凭空杜撰）；（3）可以事后检验或说明其合理性。本题中"假设站点分布均匀"就是一条不应该出现的假设——因为题目明确说了站点分布不均匀，这正是本题的难点，不能假设掉。

Q4：为什么公式很多但论文依然得分不高？

因为公式是手段而非目的。如果你有20个公式，但没有一个公式后面解释"这个公式反映了什么现实关系、为什么这样建立、如何和数据对应"，评委会认为这些公式是从网上复制的，你自己并不理解。公式少但理解深，远比公式多但理解浅得分高。

Q5：MATLAB 代码结果如何对应论文表格？

在代码中，建议将所有关键结果存入结构体（如 results.MAE1 = metrics1.MAE），最后统一格式化输出，并将输出结果直接复制到论文表格中。一定要保持代码输出和论文表格中的数字完全一致，评委如果发现两者不一致会严重质疑论文的可信度。

Q6：没有附件数据时如何构建合理分析框架？

本题数据来自竞赛附件，如果暂时没有数据，可以：（1）用模拟数据（正态分布或历史降水统计特征生成）验证代码逻辑；（2）明确标注"以下为基于模拟数据的示例"；（3）描述当真实数据可用时，结果应呈现的形态和解读方向。绝对不能把模拟结果当真实结果呈现。

Q7：预测模型如何选择误差指标？

选择误差指标要根据应用场景：MAE 不受极端值影响，适合总体表现评价；RMSE 对大误差敏感，适合强调大偏差代价的场景（如本题的暴雨预报）；MAPE 适合相对误差的场景（但本题雨量可以为0，不适合用MAPE）；$R^2$ 反映方差解释比例，适合和回归模型配合使用。本题中，RMSE 和 TS 是最核心的两个指标。

Q8：评价模型中权重如何确定？

权重确定方法主要有三类：主观赋权（专家打分、层次分析法AHP）、客观赋权（熵权法、PCA方差贡献率）、组合赋权（两者结合）。本题中，公众感受惩罚矩阵的参数（$\alpha$, $\beta$）本质上就是一种权重设定，可以通过文献调研（气象学中已有相关研究）或专家咨询确定。无论选哪种方法，都要在论文中解释选择依据，并通过灵敏度分析验证结论对权重变化的鲁棒性。

Q9：优化模型如何确定目标函数和约束条件？

本题是评价模型而非优化模型，但对于一般优化题：目标函数来自题目中"最大化/最小化"的目标（要量化为具体表达式）；约束条件来自题目中"满足/不超过/不低于"的限制条件（要写成不等式/等式形式）。如果题目没有明确说"最优化"，也可以把评价转化为"寻找使综合得分最高的方案"，从而引入优化框架。

Q10：国赛论文和美赛论文写法有什么区别？

国赛论文：中文写作，一般10000_{20000字，格式规范（摘要、正文、参考文献、附录），评价侧重模型的数学严谨性和创新性，摘要单独装订评分。美赛论文：英文写作，一般20}25页，Summary（摘要）极为重要（评委可能只看Summary），行文更重视"故事性"（narrative）和实际意义，需要有清晰的Executive Summary，可读性要求高于国赛。

Q11：如何避免论文像代码说明书？

论文的主角是"模型"而非"代码"。代码说明书的特征是：大段代码后面接"该代码实现了XXX功能"，然后直接给结果。建模论文的正确写法是：先建立数学模型，然后说明"用MATLAB实现了该模型，代码见附录"，正文中只需要说明算法的关键步骤，不需要把代码大段贴在正文里。核心代码可以放附录，正文聚焦于数学表达和结果分析。

Q12：如何写出高质量摘要？

高质量摘要的黄金模板：（1）“针对[题目核心问题]，本文建立了[核心模型名称]。”（2）“[针对问题一，做了什么，用了什么方法]；[针对问题二，做了什么，用了什么方法]。”（3）“结果表明，[方法一的核心指标具体数值]，[方法二的核心指标具体数值]。”（4）"综合分析，[结论]，具有[现实意义]。“切记：摘要不是大纲，不能只写"分析了问题一”“建立了模型”，必须有具体内容。

Q13：如何自然地提出模型改进？

改进要从"当前模型的缺陷"出发，自然过渡到改进方向。例如：在模型一的结果分析部分写"我们注意到RMSE受极端降水事件主导，为了更公平地评价各量级的预报质量，我们进一步建立了等级化评价模型（模型二）"。这样的改进是有逻辑驱动的，而不是"除此之外，我们还建立了模型二"这种无意义的并列。

Q14：模型优缺点如何写得具体？

优缺点不能写"本模型计算简单，结果直观"（这是废话），要写：（1）优点：为什么这个模型适合本题，解决了哪个具体问题；（2）缺点：这个模型在什么条件下会失效，忽略了哪些现实因素；（3）改进：如果有更多时间/数据，可以如何优化。例如本题最近邻插值的缺点是"站点密集区域格点被多个站点共享，而稀疏区域的格点无对应站点，评价区域存在空间偏倚"，改进方向是"采用克里金插值或自然邻域插值提高精度"。

Q15：附录代码应该如何整理？

附录代码建议按以下顺序排列：（1）主程序 main.m（含注释说明调用关系）；（2）数据预处理函数；（3）模型一求解函数；（4）模型二求解函数；（5）可视化函数；（6）灵敏度分析函数。每个函数文件开头要有标准的函数说明注释（功能、输入、输出），便于评委快速理解。不要把所有代码合并成一个文件——这会让代码看起来像"一团乱麻"，给人代码能力差的印象。

二十、总结

这道 2005 年国赛 C 题在表面上是一道"计算误差"的题目，但真正深入分析后你会发现它涵盖了空间数据处理、多指标评价设计、主观感受量化三大核心建模能力，是评价类建模题中难度和综合性都属上乘的经典案例。

回顾全文，我们做了以下工作：

数据层：通过最近邻插值解决了网格预报数据与站点实测数据的空间匹配问题，建立了 [N站点 × T时段] 的统一数据矩阵。

模型一（基础层）：计算 MAE、RMSE、BIAS、相关系数、TS/POD/FAR 等连续量误差和分类能力指标，建立了科学可比较的基础评价体系。

模型二（分类层）：将连续雨量映射到7个离散等级，构建混淆矩阵和分等级 TS 评分，从等级分类维度揭示两种方法的差异化表现区间。

模型三（综合层）：引入非对称惩罚矩阵量化公众感受，将漏报高等级降水的代价建模为高于空报，形成加权综合评分，直接回答了题目第二问。

检验层：通过 Bootstrap 稳定性分析、惩罚参数灵敏度分析、LOO 鲁棒性分析，验证了模型结论的可靠性。

最后，想对所有正在学习数学建模的同学说：

建模的本质是"翻译"——把真实世界的问题翻译成数学语言，再把数学结论翻译回现实意义。 代码只是工具，公式只是符号，真正重要的是：你理解了这个问题，你能解释为什么这样建模，你能说清楚结论意味着什么。

做到这三点，哪怕代码只跑通了一半，哪怕模型不够精妙，你的论文也能给评委留下"这个团队真的懂建模"的深刻印象。

加油！建模之路虽然充满挑战，但每解决一道真题，你对数学与现实世界联系的理解就深一分。

声明：以上内容部分基于人工智能辅助生成，仅供参考交流，不构成任何专业建议。模型输出可能存在偏差，使用前请自行核实，后果自负。欢迎理性讨论。

若需原题 PDF、附件或历年高教社杯真题，关注技术号 「猿圈奇妙屋」，回复【高教社杯】即可获取。

🎁 文末福利

本专栏内容源自实际建模经验、竞赛题目及读者需求。如涉及版权问题，请告知，将立即处理。部分解法思路参考了网络优秀文章，若未能完全契合你的场景，欢迎在评论区分享更优解法，共同探讨、共同进步！

更多建模方法、工具与竞赛题解，欢迎访问专栏 👉 《《滚雪球学数学建模（含历年真题）》

如果本文对你有帮助，欢迎点赞、收藏、关注，你的支持是我持续创作的动力！

同时推荐关注技术号 「猿圈奇妙屋」，获取建模干货、竞赛真题解析、4000G 技术资料、简历模板等海量内容，助你快速突破瓶颈。

🫵 关于作者

我是 bug菌，数学建模竞赛指导教师，曾指导学生斩获国赛一等奖、美赛 M 奖等，擅长运动学建模、优化模型、评价模型等方向。

活跃于 CSDN · 掘金 · InfoQ · 51CTO · 华为云 · 阿里云 · 腾讯云 · 开源中国 · 博客园 · 墨天轮 等平台

🏅 CSDN 博客之星 Top30 · 华为云十佳博主 · 掘金人气作者 Top40 · 多平台签约优质作者 · 全网粉丝 30w+

更多优质内容与成长资料 👉 点击查看 👈

欢迎加入硬核技术号 「猿圈奇妙屋」，一起进阶打怪！

- End -