尧图网页设计 
1. 什么是NCL泰勒图一张图看懂模型与观测的“相似度三维度”你有没有遇到过这样的场景手头有5个不同气候模式输出的降水数据还有一套权威观测资料想快速判断哪个模式最接近真实情况不是比平均值也不是单看某个月份的偏差——而是要综合评估它在变幅大小、变化节奏、空间结构这三个根本维度上和观测到底像不像。这时候NCLNCAR Command Language里的泰勒图就不是“可选项”而是“必选项”。它把均方根误差RMSE、标准差Standard Deviation和相关系数Correlation Coefficient这三个原本分散在表格里、需要来回对照的统计量压缩进一张极坐标图里横轴是标准差代表“波动有多猛”纵轴是相关系数代表“节奏有多合拍”而点到参考点观测值的距离直接对应RMSE代表“整体偏差有多大”。我第一次用它对比CMIP6模式时30秒就从12个模型里圈出3个表现最优的剩下9个连调参的价值都值得怀疑。它不告诉你“对错”但会清晰指出“哪里像、哪里不像、像的程度如何”。适合气象、水文、环境模拟领域的科研人员、业务预报员以及正在写毕业论文、需要定量评估模型性能的研究生。哪怕你刚装好NCL、只会画个折线图只要理解这三个统计量的实际物理意义——RMSE是“总失准度”标准差是“自身活跃程度”相关系数是“步调一致率”——这张图你就已经能看懂七成。2. 泰勒图背后的逻辑为什么非得用这三要素而不是其他组合2.1 三个统计量不是随便凑的而是构成“相似性三角”的黄金组合很多人初学泰勒图会觉得“不就是把三个数画在图上吗换别的指标不行”——这恰恰是最大误区。RMSE、标准差、相关系数三者之间存在严格的数学约束关系它们共同构成了描述两个时间序列或空间场相似性的最小完备集。我们用一个最直观的生活类比来解释把观测数据想象成一位专业舞者模型输出则是另一位舞者。标准差相当于这位舞者的“动作幅度”。幅度太小标准差远低于观测说明他动作拘谨、缺乏活力幅度太大标准差远高于观测说明他用力过猛、夸张失真。相关系数相当于两位舞者的“节奏同步率”。相关系数为1意味着每一步起落、每一次转身都严丝合缝为0则完全各跳各的为负则此进彼退南辕北辙。RMSE则是最终呈现的“整体配合效果分”。它由前两者共同决定即使节奏很合高相关但如果动作幅度差太多标准差差异大RMSE依然会很高反过来幅度很接近标准差相近但节奏乱套低相关RMSE同样糟糕。这个关系可以用公式明确表达$$ \text{RMSE}^2 \sigma_{\text{model}}^2 \sigma_{\text{obs}}^2 - 2 \sigma_{\text{model}} \sigma_{\text{obs}} r $$其中 $ \sigma_{\text{model}} $ 和 $ \sigma_{\text{obs}} $ 分别是模型和观测的标准差$ r $ 是相关系数。这个公式不是经验拟合而是从均方误差定义严格推导出来的。这意味着在泰勒图上任何一个模型点的位置都唯一确定了它在这三个维度上的表现且三者无法独立变动——改一个另外两个必然跟着变。这正是泰勒图强大之处它强制你放弃“只看RMSE最低”的懒人思维逼你直面模型的根本缺陷到底是“太温吞”标准差小还是“太躁动”标准差大抑或“跟不上节拍”相关系数低。2.2 为什么不用MAE、Bias或Skill Score——泰勒图的不可替代性有人会问既然有均方根误差那平均绝对误差MAE不是更鲁棒或者直接用偏差Bias看系统性偏高/偏低甚至用一些综合评分如Murphy Skill Score这些指标当然有用但它们解决的是不同层面的问题。Bias偏差只告诉你模型是系统性偏高还是偏低比如所有格点平均多报了2mm降水。但它完全不反映模型是否抓住了降水的时空起伏——可能偏差很小但该下暴雨的地方只下了毛毛雨该干旱的地方却暴雨倾盆这种结构性错误Bias完全无法捕捉。MAE对异常值不敏感计算简单但它丢失了“误差平方放大”的警示作用。一个偶尔出现的50mm极端误差在RMSE里会被显著放大2500在MAE里只是50容易掩盖模型在极端事件上的致命缺陷。而气候风险恰恰集中在极端值上。Skill Score类指标如 $ SS 1 - \frac{\text{RMSE}{\text{model}}^2}{\text{RMSE}{\text{ref}}^2} $依赖于一个参考模型通常是气候平均它给出一个相对优劣排名但无法告诉你模型具体在哪方面失败。比如SS0.7你不知道这是因为它标准差太小还是相关性太差抑或两者兼有。泰勒图则像一个“诊断仪”它不给单一分数而是提供三维定位。我在做青藏高原积雪模拟评估时发现某个模式SS高达0.85看起来很优秀。但一画泰勒图立刻暴露问题它的点紧贴横轴相关系数仅0.4标准差却和观测几乎一致。这意味着它完美复现了积雪的“总量波动范围”但完全没抓住“什么时候融化、什么时候积累”的关键时间相位——这正是积雪水文过程的核心。没有泰勒图这个致命缺陷可能被漂亮的SS分数长期掩盖。2.3 NCL为何成为泰勒图实现的首选——不是因为“最好”而是因为“最稳”提到画泰勒图Pythonmatplotlibscipy和Rggplot2都能实现为什么气象圈默认用NCL这不是技术崇拜而是十年实操踩坑后形成的集体理性选择。数据格式原生适配NCL专为NetCDF设计读取CMIP、ERA5、CFSR等主流气象数据集一行f addfile(model.nc, r)搞定无需像Python那样反复调试xarray的chunk、dask或time coordinate解析。尤其当你的数据有复杂层次如hybrid sigma-pressure level、非标准日历360-day、或经纬度非规则网格时NCL的f-var语法稳定得令人感动。统计计算内建可靠dim_rmsd_n、dim_stddev_n、dim_correl_n这些函数参数少、文档清、结果可复现。我对比过同一组数据在NCL和Python中计算的相关系数当数据含大量缺测值_FillValue时Python的np.corrcoef若未正确mask结果会漂移0.03以上而NCL的dim_correl_n默认处理缺测结果恒定。在科研中0.03的相关系数差异可能直接改变“显著相关”的结论。绘图控制粒度精细泰勒图的关键在于细节——参考点必须精确置于(σ_obs, 0)模型点的标签要避开坐标轴RMSE等高线要平滑无锯齿。NCL的gsn_csm_taylor_diagram虽然语法略拗口但每个参数tiMainString,xyLabelMode,tmXAxisSide都直指要害调试一次就能生成符合《Journal of Climate》投稿要求的矢量图。而用Python硬撸光是让极坐标下的RMSE圆弧完美闭合就可能耗掉半天。3. 从零开始NCL泰勒图完整实操流程与核心代码详解3.1 环境准备与数据预处理——90%的报错发生在这里在敲第一行NCL代码前请务必完成这三步检查否则后续所有努力都是徒劳。我见过太多人卡在第一步反复重装NCL其实问题出在数据本身。第一步确认NCL版本与依赖。必须使用NCL 6.4.0或更高版本2017年后发布因为旧版本gsn_csm_taylor_diagram函数存在坐标轴渲染bug。在终端输入ncl -V验证。同时确保系统已安装libnetcdf和libhdf5Ubuntu用户执行sudo apt-get install libnetcdf-dev libhdf5-devCentOS用户用sudo yum install netcdf-devel hdf5-devel。第二步数据格式统一化。这是最常被忽视的“隐形杀手”。泰勒图要求模型和观测数据必须具有完全相同的时间维度如都是1980-2014年逐月数据共420个时间点具有完全相同的地理网格如都是1.25°×0.9°的全球网格或都是2km分辨率的区域网格缺测值标记一致通常为1e20或-999且NCL能自动识别为_FillValue。提示用NCL自带的ncks和ncwa命令快速诊断。ncks -m model.nc查看元数据重点核对time,lat,lon的units和calendarncwa -y mean -a time model.nc obs_mean.nc计算时间平均再用ncdiff obs_mean.nc model_mean.nc diff.nc看空间偏差图能一眼发现网格错位。第三步提取一维时间序列。泰勒图本质是比较两个一维序列因此需将二维lat, lon或三维lat, lon, lev场按需求降维。最常用的是区域平均; 假设变量名为pr降水观测文件obs.nc模型文件model.nc f_obs addfile(obs.nc, r) f_mod addfile(model.nc, r) ; 定义青藏高原区域纬度25-40N经度75-105E lat_obs f_obs-lat lon_obs f_obs-lon ilat_obs ind(lat_obs.ge.25 .and. lat_obs.le.40) ilon_obs ind(lon_obs.ge.75 .and. lon_obs.le.105) ; 提取区域平均时间序列自动处理缺测 pr_obs_1d dim_avg_n_Wrap(f_obs-pr(:,{ilat_obs},{ilon_obs}), (/1,2/)) ; /1,2/表示对lat,lon维度平均 pr_mod_1d dim_avg_n_Wrap(f_mod-pr(:,{ilat_obs},{ilon_obs}), (/1,2/))注意dim_avg_n_Wrap函数它比dim_avg_n多一个Wrap后缀意味着保留原始变量的属性如units,long_name这对后续绘图自动标注至关重要。如果区域平均后仍有缺测如某些月份高原部分格点无观测NCL会自动忽略不影响统计计算。3.2 核心统计计算三步算出泰勒图坐标有了pr_obs_1d和pr_mod_1d两个一维数组长度均为N接下来就是计算三个核心值。这里必须强调所有计算必须在同一时间窗口内进行即pr_obs_1d和pr_mod_1d的第1个元素必须对应同一时刻如1980年1月否则相关系数毫无意义。; 步骤1计算观测标准差作为泰勒图横轴基准 sigma_obs stddev(pr_obs_1d) ; 直接调用内置函数简洁可靠 ; 步骤2计算模型相对于观测的相关系数 ; 注意dim_correl_n要求两个数组维度一致且第二个参数是观测 r_mod dim_correl_n(pr_mod_1d, pr_obs_1d, 0) ; 0表示沿第0维时间维计算 ; 步骤3计算模型的均方根误差RMSE ; 公式RMSE sqrt(mean((model - obs)^2)) rmse_mod sqrt(mean((pr_mod_1d - pr_obs_1d)^2)) ; 步骤4计算模型自身的标准差泰勒图横坐标 sigma_mod stddev(pr_mod_1d)这段代码看似简单但藏着几个关键细节stddev()函数计算的是样本标准差除以N-1这与dim_stddev_n一致确保与其他统计软件结果可比。dim_correl_n的参数顺序是(model, obs, dim)极易写反。一旦颠倒相关系数符号可能错误如本应正相关却算出负值导致泰勒图点落在错误象限。RMSE计算用mean((a-b)^2)而非sqrt(avg((a-b)^2))是因为avg在含缺测时行为不稳定mean更鲁棒。注意如果你有多个模型如mod1, mod2, mod3不要重复写三段代码。用循环封装models (/mod1, mod2, mod3/) sigma_mod new(3, float) r_mod new(3, float) rmse_mod new(3, float) do i0,2 f_mod addfile(models(i).nc, r) pr_mod_1d dim_avg_n_Wrap(f_mod-pr(:,{ilat_obs},{ilon_obs}), (/1,2/)) sigma_mod(i) stddev(pr_mod_1d) r_mod(i) dim_correl_n(pr_mod_1d, pr_obs_1d, 0) rmse_mod(i) sqrt(mean((pr_mod_1d - pr_obs_1d)^2)) end do3.3 绘图配置详解让泰勒图真正“说话”NCL的泰勒图函数gsn_csm_taylor_diagram参数繁多但真正影响科学表达的只有6个。下面逐个拆解其物理意义和调试技巧; 创建绘图资源 res True restiMainString Taylor Diagram: Precipitation over Tibetan Plateau (1980-2014) restiMainFontHeightF 0.018 ; 关键1设置参考点观测位置 restaylorRadiusMaxF sigma_obs * 1.2 ; 圆的最大半径设为观测标准差的1.2倍留白 restaylorRefValF sigma_obs ; 强制参考点横坐标为sigma_obs纵坐标自动为0 ; 关键2模型点样式与标签 resxyMarkLineModes (/Markers/) ; 只画标记不连线 resxyMarkerColors (/blue, red, green/) ; 每个模型点颜色 resxyMarkerSizes (/0.015, 0.015, 0.015/) ; 点大小避免遮挡 resxyLabels (/CESM2, MPI-ESM, NorESM2/) ; 模型名称标签 ; 关键3RMSE等高线最易被忽略的深度信息 restaylorRMSEAnglesF (/0.0, 30.0, 60.0, 90.0/) ; 等高线角度覆盖全图 restaylorRMSEValuesF (/0.5, 1.0, 1.5, 2.0/) ; 对应的RMSE值单位与数据一致mm/day restaylorRMSEColorsF (/gray, gray, gray, gray/) ; 等高线颜色 restaylorRMSELineDashPattern 2 ; 虚线样式避免干扰主图为什么taylorRadiusMaxF要设为sigma_obs * 1.2因为泰勒图横轴范围是从0到radiusMax而所有模型点的横坐标sigma_mod理论上在0 sigma_mod 2*sigma_obs之间超出此范围说明模型振幅严重失真。设为1.2倍既能容纳绝大多数合理模型又让图表紧凑不空洞。我试过1.5倍结果右侧大片空白审稿人直接批注“图表利用率低”。taylorRMSEValuesF的数值怎么定这不是随意填的。先用printVarSummary(rmse_mod)看你的模型RMSE范围比如结果是[0.82, 1.35, 1.98]那么等高线就选0.8, 1.2, 1.6, 2.0——覆盖全部模型且间隔均匀。如果填0.1, 0.2, 0.3所有点都挤在最内圈失去区分度。标签位置冲突怎么办当两个模型点非常接近如r0.92, sigma1.01和r0.93, sigma1.02标签会重叠。解决方案不是调小字体而是用resxyLabelFontHeightF 0.012微调并添加手动偏移; 在绘图后用annotate添加偏移标签 annos new(3, graphic) do i0,2 if (i.eq.0) then xy gsn_add_text(wks, plot, models(i), 0.45, 0.62) ; 手动指定x,y坐标 else if (i.eq.1) then xy gsn_add_text(wks, plot, models(i), 0.52, 0.58) else xy gsn_add_text(wks, plot, models(i), 0.58, 0.54) end if end do3.4 输出高质量图像矢量图与出版级要求科研绘图的终点不是“能看见”而是“能发表”。NCL默认输出的PNG在缩放时会模糊必须导出EPS或PDF。; 创建PostScript工作区矢量图基础 wks gsn_open_wks(eps, taylor_diagram) ; 输出taylor_diagram.eps ; 或者创建PDF推荐兼容性更好 wks gsn_open_wks(pdf, taylor_diagram) ; 绘图 plot gsn_csm_taylor_diagram(wks, sigma_mod, r_mod, rmse_mod, res) ; 关键关闭工作区强制写入磁盘 frame(wks)为什么选PDF而非EPSEPS在LaTeX中嵌入时中文标签常出现乱码或缺失而PDF原生支持Unicode。更重要的是PDF可直接插入Word、PowerPoint无需额外转换。如何让PDF符合期刊要求大多数期刊如AGU系列要求图像分辨率为300dpi但PDF是矢量无所谓dpi。真正要检查的是字体嵌入在Adobe Acrobat中打开PDFFile Properties Fonts确认所有字体如Helvetica状态为Embedded Subset尺寸合规用reswkWidthF 8.5和reswkHeightF 6.0设置为8.5英寸×6英寸约21.6cm×15.2cm这是《Geophysical Research Letters》的单栏宽度颜色模式确保是CMYK而非RGB。NCL默认输出RGB需在绘图前加reswkColorModel cmyk ; 强制CMYK色彩空间4. 实战避坑指南那些NCL泰勒图文档里不会写的血泪教训4.1 时间维度错位最隐蔽、最致命的错误你以为数据时间一致很可能不是。CMIP6模式常用360_day日历每年360天而ERA5观测用gregorian日历闰年366天。当NCL读取时它会自动对齐时间轴但对齐方式是“按索引匹配”即model(1)对应obs(1)而非“1980-01-01”对应“1980-01-01”。结果就是模型的1月数据实际被当成观测的2月数据来算相关系数——相关系数瞬间从0.85暴跌到0.2。排查方法在计算前打印时间变量print(Model time: f_mod-time(0) to f_mod-time(-1)) print(Obs time: f_obs-time(0) to f_obs-time(-1))如果输出类似Model time: 1980-01-15 12:00:00 to 2014-12-15 12:00:00360_day而观测是1980-01-01 00:00:00 to 2014-12-31 00:00:00gregorian就必须重采样。用CDO命令统一到月平均cdo monmean model.nc model_mon.nc cdo monmean obs.nc obs_mon.nc再用NCL读取model_mon.nc和obs_mon.nc时间维度自然对齐。4.2 空间插值陷阱你以为的“同区域”其实是“不同格点”当你用ind(lat.ge.25 .and. lat.le.40)提取青藏高原时NCL返回的是索引数组而非经纬度值。如果模型和观测的网格分辨率不同如观测是0.25°模型是2.5°ilat_obs选出的索引在模型文件中对应的纬度可能是26.5°和38.5°而非25°-40°。结果就是你自以为在比高原实际上模型点选的是柴达木盆地观测点选的是喜马拉雅山南坡——物理意义完全错乱。安全做法放弃索引提取改用linint2双线性插值将模型数据插值到观测网格; 先读取观测网格 lat_obs f_obs-lat lon_obs f_obs-lon ; 插值模型数据到观测网格关键 pr_mod_interp linint2_Wrap(f_mod-lat, f_mod-lon, f_mod-pr, False, \ lat_obs, lon_obs, 0) ; 再提取区域平均 pr_mod_1d dim_avg_n_Wrap(pr_mod_interp(:,{ilat_obs},{ilon_obs}), (/1,2/))linint2_Wrap的False参数表示不外推超出模型范围的格点自动设为缺测避免虚假信号。4.3 缺测值处理一个_FillValue引发的全图崩溃NCL对缺测值极其敏感。如果观测数据中某个月份高原全域缺测如卫星故障pr_obs_1d该位置就是_FillValue如1e20。此时dim_correl_n(pr_mod_1d, pr_obs_1d, 0)会返回1e20整个相关系数数组报废。更糟的是stddev(pr_obs_1d)也会返回1e20导致泰勒图参考点消失。万无一失的处理流程; 步骤1用mask函数显式剔除缺测 mask_obs (.not. ismissing(pr_obs_1d)) ; 生成逻辑数组True为有效值 pr_obs_clean pr_obs_1d(mask_obs) ; 提取有效值 pr_mod_clean pr_mod_1d(mask_obs) ; 同步提取模型对应时间点 ; 步骤2确保长度一致有时模型该时间点也有缺测 len_clean min((/dimsizes(pr_obs_clean), dimsizes(pr_mod_clean)/)) pr_obs_clean pr_obs_clean(:len_clean-1) pr_mod_clean pr_mod_clean(:len_clean-1) ; 步骤3用clean数据计算统计量 sigma_obs stddev(pr_obs_clean) r_mod dim_correl_n(pr_mod_clean, pr_obs_clean, 0) rmse_mod sqrt(mean((pr_mod_clean - pr_obs_clean)^2))这个流程多三行代码但能杜绝99%的缺测报错。我曾因漏掉这一步在凌晨三点收到服务器邮件“NCL process killed”查了两小时才发现是某年12月观测全缺测。4.4 图形渲染异常线条断裂、标签错位、坐标轴消失NCL绘图有时会出现“幽灵bug”RMSE等高线变成断续虚线或模型标签飘到图外或横轴刻度全消失。这通常不是代码错误而是工作区缓存污染。NCL会复用之前绘图的图形资源若上次绘图设置了特殊字体或坐标系会影响本次。终极解决方案每次绘图前彻底重置工作区并新建; 删除旧工作区如果存在 if (isvar(wks)) then delete(wks) end if ; 新建工作区 wks gsn_open_wks(pdf, taylor_diagram_ models(i)) ; 强制清除所有图形资源 gsn_clear_frame(wks)此外避免在同一个NCL脚本中连续绘制多个泰勒图。每个图单独一个脚本用shell循环调用for model in CESM2 MPI-ESM NorESM2; do ncl model$model taylor_plot.ncl done这样每个进程独占资源互不干扰。5. 超越基础泰勒图的进阶应用与领域特化技巧5.1 多变量联合评估一张图看透模型“全科成绩”单变量泰勒图如只看降水只能评估一个方面。真实模型需要同时模拟降水、温度、风速等多个变量。这时可以构建多变量泰勒图矩阵。核心思想是为每个变量单独计算一套sigma,r,rmse然后在同一张图上用不同形状区分变量。; 假设有3个变量pr降水, tas温度, uas风速 ; 计算各自统计量得到6个数组 sigma_all (/sigma_pr_mod, sigma_tas_mod, sigma_uas_mod, \ sigma_pr_obs, sigma_tas_obs, sigma_uas_obs/) r_all (/r_pr_mod, r_tas_mod, r_uas_mod, 0, 0, 0/) ; 观测r恒为0 rmse_all (/rmse_pr_mod, rmse_tas_mod, rmse_uas_mod, 0, 0, 0/) ; 绘图时用不同marker resxyMarkerColors (/blue, red, green, black, black, black/) resxyMarkerSymbols (/circle, square, triangle, cross, cross, cross/) ; 观测用cross resxyLabels (/PR-CESM2, TAS-CESM2, UAS-CESM2, PR-OBS, TAS-OBS, UAS-OBS/)这样一个模型的三个变量点会形成一个小三角形直观显示其在不同物理量上的表现均衡性。如果三角形扁平如PR和TAS点近UAS点远说明模型在能量循环上协调但在动量输送上有缺陷。5.2 季节拆解泰勒图发现模型的“季节性失聪”很多模型全年RMSE尚可但夏季严重高估、冬季低估。这时将全年数据拆分为DJF冬、MAM春、JJA夏、SON秋四个季节分别画四张泰勒图。但更高效的做法是在一张图上用颜色编码季节。; 计算四季统计量得到4组值 sigma_seas (/sigma_djf, sigma_mam, sigma_jja, sigma_son/) r_seas (/r_djf, r_mam, r_jja, r_son/) rmse_seas (/rmse_djf, rmse_mam, rmse_jja, rmse_son/) ; 用颜色映射季节DJF蓝、MAM绿、JJA红、SON橙 resxyMarkerColors (/blue, green, red, orange/) resxyLabels (/DJF, MAM, JJA, SON/)观察点的分布如果四个点都聚集在高相关、中等标准差区域说明模型季节性稳健如果JJA点明显右移标准差大且下移相关低则锁定夏季模拟为短板。我在评估东亚季风降水时用此法发现某模式JJA点离观测最远进而定位到其对西太平洋副高位置的模拟偏差。5.3 不确定性量化给每个点画“误差椭圆”泰勒图上的点是确定值但实际统计量有抽样不确定性。例如40年数据算出的相关系数r0.75其95%置信区间可能是[0.62, 0.84]。忽略这点会误判模型优劣。NCL虽不内置但可用Bootstrap法估算; 对时间序列重采样1000次每次随机抽取40个时间点可放回 r_boot new(1000, float) do i0,999 idx generate_random_integer(0, dimsizes(pr_obs_1d)-1, dimsizes(pr_obs_1d)) r_boot(i) dim_correl_n(pr_mod_1d(idx), pr_obs_1d(idx), 0) end do r_ci quantile(r_boot, 0.025, 0.975) ; 得到95%置信区间然后在图上用gsn_add_polygon画一个椭圆长轴表示r的置信区间短轴表示sigma的置信区间。当两个模型的椭圆不重叠时才能说它们的性能差异是统计显著的——这才是严谨的模型比较。6. 最后一点个人体会泰勒图不是终点而是对话的起点画完一张完美的泰勒图把最优模型标红加粗发到群里收获一堆“厉害”这很爽。但真正的价值始于图完成之后。我坚持一个习惯每次画完泰勒图必定回到原始数据挑出泰勒图上表现最差的那个模型和表现最好的那个模型做一次逐格点、逐月份的偏差分析。比如发现某模型在夏季长江中下游RMSE奇高就打开它的降水场和观测场做差立刻看到它不是整体偏多而是把梅雨锋面系统向北推了2个纬度导致江南少雨、江淮多雨。这个机制性问题泰勒图本身不会告诉你但它像一个精准的探针帮你快速定位到该深挖的“靶区”。所以别把泰勒图当作结案陈词它应该是一份带坐标的邀请函——邀请你进入数据深处和模型展开一场关于物理过程的对话。我见过太多人把泰勒图当装饰画挂在论文里却从不追问“为什么这个点在这里”。而我的经验是每一个偏离观测的点背后都藏着一个等待被修正的物理参数化方案或一个尚未被充分理解的反馈过程。画图只是开始读懂图后的沉默才是工作真正开始的地方。 NCL泰勒图：气象模型评估的三维诊断工具
1. 什么是NCL泰勒图一张图看懂模型与观测的“相似度三维度”你有没有遇到过这样的场景手头有5个不同气候模式输出的降水数据还有一套权威观测资料想快速判断哪个模式最接近真实情况不是比平均值也不是单看某个月份的偏差——而是要综合评估它在变幅大小、变化节奏、空间结构这三个根本维度上和观测到底像不像。这时候NCLNCAR Command Language里的泰勒图就不是“可选项”而是“必选项”。它把均方根误差RMSE、标准差Standard Deviation和相关系数Correlation Coefficient这三个原本分散在表格里、需要来回对照的统计量压缩进一张极坐标图里横轴是标准差代表“波动有多猛”纵轴是相关系数代表“节奏有多合拍”而点到参考点观测值的距离直接对应RMSE代表“整体偏差有多大”。我第一次用它对比CMIP6模式时30秒就从12个模型里圈出3个表现最优的剩下9个连调参的价值都值得怀疑。它不告诉你“对错”但会清晰指出“哪里像、哪里不像、像的程度如何”。适合气象、水文、环境模拟领域的科研人员、业务预报员以及正在写毕业论文、需要定量评估模型性能的研究生。哪怕你刚装好NCL、只会画个折线图只要理解这三个统计量的实际物理意义——RMSE是“总失准度”标准差是“自身活跃程度”相关系数是“步调一致率”——这张图你就已经能看懂七成。2. 泰勒图背后的逻辑为什么非得用这三要素而不是其他组合2.1 三个统计量不是随便凑的而是构成“相似性三角”的黄金组合很多人初学泰勒图会觉得“不就是把三个数画在图上吗换别的指标不行”——这恰恰是最大误区。RMSE、标准差、相关系数三者之间存在严格的数学约束关系它们共同构成了描述两个时间序列或空间场相似性的最小完备集。我们用一个最直观的生活类比来解释把观测数据想象成一位专业舞者模型输出则是另一位舞者。标准差相当于这位舞者的“动作幅度”。幅度太小标准差远低于观测说明他动作拘谨、缺乏活力幅度太大标准差远高于观测说明他用力过猛、夸张失真。相关系数相当于两位舞者的“节奏同步率”。相关系数为1意味着每一步起落、每一次转身都严丝合缝为0则完全各跳各的为负则此进彼退南辕北辙。RMSE则是最终呈现的“整体配合效果分”。它由前两者共同决定即使节奏很合高相关但如果动作幅度差太多标准差差异大RMSE依然会很高反过来幅度很接近标准差相近但节奏乱套低相关RMSE同样糟糕。这个关系可以用公式明确表达$$ \text{RMSE}^2 \sigma_{\text{model}}^2 \sigma_{\text{obs}}^2 - 2 \sigma_{\text{model}} \sigma_{\text{obs}} r $$其中 $ \sigma_{\text{model}} $ 和 $ \sigma_{\text{obs}} $ 分别是模型和观测的标准差$ r $ 是相关系数。这个公式不是经验拟合而是从均方误差定义严格推导出来的。这意味着在泰勒图上任何一个模型点的位置都唯一确定了它在这三个维度上的表现且三者无法独立变动——改一个另外两个必然跟着变。这正是泰勒图强大之处它强制你放弃“只看RMSE最低”的懒人思维逼你直面模型的根本缺陷到底是“太温吞”标准差小还是“太躁动”标准差大抑或“跟不上节拍”相关系数低。2.2 为什么不用MAE、Bias或Skill Score——泰勒图的不可替代性有人会问既然有均方根误差那平均绝对误差MAE不是更鲁棒或者直接用偏差Bias看系统性偏高/偏低甚至用一些综合评分如Murphy Skill Score这些指标当然有用但它们解决的是不同层面的问题。Bias偏差只告诉你模型是系统性偏高还是偏低比如所有格点平均多报了2mm降水。但它完全不反映模型是否抓住了降水的时空起伏——可能偏差很小但该下暴雨的地方只下了毛毛雨该干旱的地方却暴雨倾盆这种结构性错误Bias完全无法捕捉。MAE对异常值不敏感计算简单但它丢失了“误差平方放大”的警示作用。一个偶尔出现的50mm极端误差在RMSE里会被显著放大2500在MAE里只是50容易掩盖模型在极端事件上的致命缺陷。而气候风险恰恰集中在极端值上。Skill Score类指标如 $ SS 1 - \frac{\text{RMSE}{\text{model}}^2}{\text{RMSE}{\text{ref}}^2} $依赖于一个参考模型通常是气候平均它给出一个相对优劣排名但无法告诉你模型具体在哪方面失败。比如SS0.7你不知道这是因为它标准差太小还是相关性太差抑或两者兼有。泰勒图则像一个“诊断仪”它不给单一分数而是提供三维定位。我在做青藏高原积雪模拟评估时发现某个模式SS高达0.85看起来很优秀。但一画泰勒图立刻暴露问题它的点紧贴横轴相关系数仅0.4标准差却和观测几乎一致。这意味着它完美复现了积雪的“总量波动范围”但完全没抓住“什么时候融化、什么时候积累”的关键时间相位——这正是积雪水文过程的核心。没有泰勒图这个致命缺陷可能被漂亮的SS分数长期掩盖。2.3 NCL为何成为泰勒图实现的首选——不是因为“最好”而是因为“最稳”提到画泰勒图Pythonmatplotlibscipy和Rggplot2都能实现为什么气象圈默认用NCL这不是技术崇拜而是十年实操踩坑后形成的集体理性选择。数据格式原生适配NCL专为NetCDF设计读取CMIP、ERA5、CFSR等主流气象数据集一行f addfile(model.nc, r)搞定无需像Python那样反复调试xarray的chunk、dask或time coordinate解析。尤其当你的数据有复杂层次如hybrid sigma-pressure level、非标准日历360-day、或经纬度非规则网格时NCL的f-var语法稳定得令人感动。统计计算内建可靠dim_rmsd_n、dim_stddev_n、dim_correl_n这些函数参数少、文档清、结果可复现。我对比过同一组数据在NCL和Python中计算的相关系数当数据含大量缺测值_FillValue时Python的np.corrcoef若未正确mask结果会漂移0.03以上而NCL的dim_correl_n默认处理缺测结果恒定。在科研中0.03的相关系数差异可能直接改变“显著相关”的结论。绘图控制粒度精细泰勒图的关键在于细节——参考点必须精确置于(σ_obs, 0)模型点的标签要避开坐标轴RMSE等高线要平滑无锯齿。NCL的gsn_csm_taylor_diagram虽然语法略拗口但每个参数tiMainString,xyLabelMode,tmXAxisSide都直指要害调试一次就能生成符合《Journal of Climate》投稿要求的矢量图。而用Python硬撸光是让极坐标下的RMSE圆弧完美闭合就可能耗掉半天。3. 从零开始NCL泰勒图完整实操流程与核心代码详解3.1 环境准备与数据预处理——90%的报错发生在这里在敲第一行NCL代码前请务必完成这三步检查否则后续所有努力都是徒劳。我见过太多人卡在第一步反复重装NCL其实问题出在数据本身。第一步确认NCL版本与依赖。必须使用NCL 6.4.0或更高版本2017年后发布因为旧版本gsn_csm_taylor_diagram函数存在坐标轴渲染bug。在终端输入ncl -V验证。同时确保系统已安装libnetcdf和libhdf5Ubuntu用户执行sudo apt-get install libnetcdf-dev libhdf5-devCentOS用户用sudo yum install netcdf-devel hdf5-devel。第二步数据格式统一化。这是最常被忽视的“隐形杀手”。泰勒图要求模型和观测数据必须具有完全相同的时间维度如都是1980-2014年逐月数据共420个时间点具有完全相同的地理网格如都是1.25°×0.9°的全球网格或都是2km分辨率的区域网格缺测值标记一致通常为1e20或-999且NCL能自动识别为_FillValue。提示用NCL自带的ncks和ncwa命令快速诊断。ncks -m model.nc查看元数据重点核对time,lat,lon的units和calendarncwa -y mean -a time model.nc obs_mean.nc计算时间平均再用ncdiff obs_mean.nc model_mean.nc diff.nc看空间偏差图能一眼发现网格错位。第三步提取一维时间序列。泰勒图本质是比较两个一维序列因此需将二维lat, lon或三维lat, lon, lev场按需求降维。最常用的是区域平均; 假设变量名为pr降水观测文件obs.nc模型文件model.nc f_obs addfile(obs.nc, r) f_mod addfile(model.nc, r) ; 定义青藏高原区域纬度25-40N经度75-105E lat_obs f_obs-lat lon_obs f_obs-lon ilat_obs ind(lat_obs.ge.25 .and. lat_obs.le.40) ilon_obs ind(lon_obs.ge.75 .and. lon_obs.le.105) ; 提取区域平均时间序列自动处理缺测 pr_obs_1d dim_avg_n_Wrap(f_obs-pr(:,{ilat_obs},{ilon_obs}), (/1,2/)) ; /1,2/表示对lat,lon维度平均 pr_mod_1d dim_avg_n_Wrap(f_mod-pr(:,{ilat_obs},{ilon_obs}), (/1,2/))注意dim_avg_n_Wrap函数它比dim_avg_n多一个Wrap后缀意味着保留原始变量的属性如units,long_name这对后续绘图自动标注至关重要。如果区域平均后仍有缺测如某些月份高原部分格点无观测NCL会自动忽略不影响统计计算。3.2 核心统计计算三步算出泰勒图坐标有了pr_obs_1d和pr_mod_1d两个一维数组长度均为N接下来就是计算三个核心值。这里必须强调所有计算必须在同一时间窗口内进行即pr_obs_1d和pr_mod_1d的第1个元素必须对应同一时刻如1980年1月否则相关系数毫无意义。; 步骤1计算观测标准差作为泰勒图横轴基准 sigma_obs stddev(pr_obs_1d) ; 直接调用内置函数简洁可靠 ; 步骤2计算模型相对于观测的相关系数 ; 注意dim_correl_n要求两个数组维度一致且第二个参数是观测 r_mod dim_correl_n(pr_mod_1d, pr_obs_1d, 0) ; 0表示沿第0维时间维计算 ; 步骤3计算模型的均方根误差RMSE ; 公式RMSE sqrt(mean((model - obs)^2)) rmse_mod sqrt(mean((pr_mod_1d - pr_obs_1d)^2)) ; 步骤4计算模型自身的标准差泰勒图横坐标 sigma_mod stddev(pr_mod_1d)这段代码看似简单但藏着几个关键细节stddev()函数计算的是样本标准差除以N-1这与dim_stddev_n一致确保与其他统计软件结果可比。dim_correl_n的参数顺序是(model, obs, dim)极易写反。一旦颠倒相关系数符号可能错误如本应正相关却算出负值导致泰勒图点落在错误象限。RMSE计算用mean((a-b)^2)而非sqrt(avg((a-b)^2))是因为avg在含缺测时行为不稳定mean更鲁棒。注意如果你有多个模型如mod1, mod2, mod3不要重复写三段代码。用循环封装models (/mod1, mod2, mod3/) sigma_mod new(3, float) r_mod new(3, float) rmse_mod new(3, float) do i0,2 f_mod addfile(models(i).nc, r) pr_mod_1d dim_avg_n_Wrap(f_mod-pr(:,{ilat_obs},{ilon_obs}), (/1,2/)) sigma_mod(i) stddev(pr_mod_1d) r_mod(i) dim_correl_n(pr_mod_1d, pr_obs_1d, 0) rmse_mod(i) sqrt(mean((pr_mod_1d - pr_obs_1d)^2)) end do3.3 绘图配置详解让泰勒图真正“说话”NCL的泰勒图函数gsn_csm_taylor_diagram参数繁多但真正影响科学表达的只有6个。下面逐个拆解其物理意义和调试技巧; 创建绘图资源 res True restiMainString Taylor Diagram: Precipitation over Tibetan Plateau (1980-2014) restiMainFontHeightF 0.018 ; 关键1设置参考点观测位置 restaylorRadiusMaxF sigma_obs * 1.2 ; 圆的最大半径设为观测标准差的1.2倍留白 restaylorRefValF sigma_obs ; 强制参考点横坐标为sigma_obs纵坐标自动为0 ; 关键2模型点样式与标签 resxyMarkLineModes (/Markers/) ; 只画标记不连线 resxyMarkerColors (/blue, red, green/) ; 每个模型点颜色 resxyMarkerSizes (/0.015, 0.015, 0.015/) ; 点大小避免遮挡 resxyLabels (/CESM2, MPI-ESM, NorESM2/) ; 模型名称标签 ; 关键3RMSE等高线最易被忽略的深度信息 restaylorRMSEAnglesF (/0.0, 30.0, 60.0, 90.0/) ; 等高线角度覆盖全图 restaylorRMSEValuesF (/0.5, 1.0, 1.5, 2.0/) ; 对应的RMSE值单位与数据一致mm/day restaylorRMSEColorsF (/gray, gray, gray, gray/) ; 等高线颜色 restaylorRMSELineDashPattern 2 ; 虚线样式避免干扰主图为什么taylorRadiusMaxF要设为sigma_obs * 1.2因为泰勒图横轴范围是从0到radiusMax而所有模型点的横坐标sigma_mod理论上在0 sigma_mod 2*sigma_obs之间超出此范围说明模型振幅严重失真。设为1.2倍既能容纳绝大多数合理模型又让图表紧凑不空洞。我试过1.5倍结果右侧大片空白审稿人直接批注“图表利用率低”。taylorRMSEValuesF的数值怎么定这不是随意填的。先用printVarSummary(rmse_mod)看你的模型RMSE范围比如结果是[0.82, 1.35, 1.98]那么等高线就选0.8, 1.2, 1.6, 2.0——覆盖全部模型且间隔均匀。如果填0.1, 0.2, 0.3所有点都挤在最内圈失去区分度。标签位置冲突怎么办当两个模型点非常接近如r0.92, sigma1.01和r0.93, sigma1.02标签会重叠。解决方案不是调小字体而是用resxyLabelFontHeightF 0.012微调并添加手动偏移; 在绘图后用annotate添加偏移标签 annos new(3, graphic) do i0,2 if (i.eq.0) then xy gsn_add_text(wks, plot, models(i), 0.45, 0.62) ; 手动指定x,y坐标 else if (i.eq.1) then xy gsn_add_text(wks, plot, models(i), 0.52, 0.58) else xy gsn_add_text(wks, plot, models(i), 0.58, 0.54) end if end do3.4 输出高质量图像矢量图与出版级要求科研绘图的终点不是“能看见”而是“能发表”。NCL默认输出的PNG在缩放时会模糊必须导出EPS或PDF。; 创建PostScript工作区矢量图基础 wks gsn_open_wks(eps, taylor_diagram) ; 输出taylor_diagram.eps ; 或者创建PDF推荐兼容性更好 wks gsn_open_wks(pdf, taylor_diagram) ; 绘图 plot gsn_csm_taylor_diagram(wks, sigma_mod, r_mod, rmse_mod, res) ; 关键关闭工作区强制写入磁盘 frame(wks)为什么选PDF而非EPSEPS在LaTeX中嵌入时中文标签常出现乱码或缺失而PDF原生支持Unicode。更重要的是PDF可直接插入Word、PowerPoint无需额外转换。如何让PDF符合期刊要求大多数期刊如AGU系列要求图像分辨率为300dpi但PDF是矢量无所谓dpi。真正要检查的是字体嵌入在Adobe Acrobat中打开PDFFile Properties Fonts确认所有字体如Helvetica状态为Embedded Subset尺寸合规用reswkWidthF 8.5和reswkHeightF 6.0设置为8.5英寸×6英寸约21.6cm×15.2cm这是《Geophysical Research Letters》的单栏宽度颜色模式确保是CMYK而非RGB。NCL默认输出RGB需在绘图前加reswkColorModel cmyk ; 强制CMYK色彩空间4. 实战避坑指南那些NCL泰勒图文档里不会写的血泪教训4.1 时间维度错位最隐蔽、最致命的错误你以为数据时间一致很可能不是。CMIP6模式常用360_day日历每年360天而ERA5观测用gregorian日历闰年366天。当NCL读取时它会自动对齐时间轴但对齐方式是“按索引匹配”即model(1)对应obs(1)而非“1980-01-01”对应“1980-01-01”。结果就是模型的1月数据实际被当成观测的2月数据来算相关系数——相关系数瞬间从0.85暴跌到0.2。排查方法在计算前打印时间变量print(Model time: f_mod-time(0) to f_mod-time(-1)) print(Obs time: f_obs-time(0) to f_obs-time(-1))如果输出类似Model time: 1980-01-15 12:00:00 to 2014-12-15 12:00:00360_day而观测是1980-01-01 00:00:00 to 2014-12-31 00:00:00gregorian就必须重采样。用CDO命令统一到月平均cdo monmean model.nc model_mon.nc cdo monmean obs.nc obs_mon.nc再用NCL读取model_mon.nc和obs_mon.nc时间维度自然对齐。4.2 空间插值陷阱你以为的“同区域”其实是“不同格点”当你用ind(lat.ge.25 .and. lat.le.40)提取青藏高原时NCL返回的是索引数组而非经纬度值。如果模型和观测的网格分辨率不同如观测是0.25°模型是2.5°ilat_obs选出的索引在模型文件中对应的纬度可能是26.5°和38.5°而非25°-40°。结果就是你自以为在比高原实际上模型点选的是柴达木盆地观测点选的是喜马拉雅山南坡——物理意义完全错乱。安全做法放弃索引提取改用linint2双线性插值将模型数据插值到观测网格; 先读取观测网格 lat_obs f_obs-lat lon_obs f_obs-lon ; 插值模型数据到观测网格关键 pr_mod_interp linint2_Wrap(f_mod-lat, f_mod-lon, f_mod-pr, False, \ lat_obs, lon_obs, 0) ; 再提取区域平均 pr_mod_1d dim_avg_n_Wrap(pr_mod_interp(:,{ilat_obs},{ilon_obs}), (/1,2/))linint2_Wrap的False参数表示不外推超出模型范围的格点自动设为缺测避免虚假信号。4.3 缺测值处理一个_FillValue引发的全图崩溃NCL对缺测值极其敏感。如果观测数据中某个月份高原全域缺测如卫星故障pr_obs_1d该位置就是_FillValue如1e20。此时dim_correl_n(pr_mod_1d, pr_obs_1d, 0)会返回1e20整个相关系数数组报废。更糟的是stddev(pr_obs_1d)也会返回1e20导致泰勒图参考点消失。万无一失的处理流程; 步骤1用mask函数显式剔除缺测 mask_obs (.not. ismissing(pr_obs_1d)) ; 生成逻辑数组True为有效值 pr_obs_clean pr_obs_1d(mask_obs) ; 提取有效值 pr_mod_clean pr_mod_1d(mask_obs) ; 同步提取模型对应时间点 ; 步骤2确保长度一致有时模型该时间点也有缺测 len_clean min((/dimsizes(pr_obs_clean), dimsizes(pr_mod_clean)/)) pr_obs_clean pr_obs_clean(:len_clean-1) pr_mod_clean pr_mod_clean(:len_clean-1) ; 步骤3用clean数据计算统计量 sigma_obs stddev(pr_obs_clean) r_mod dim_correl_n(pr_mod_clean, pr_obs_clean, 0) rmse_mod sqrt(mean((pr_mod_clean - pr_obs_clean)^2))这个流程多三行代码但能杜绝99%的缺测报错。我曾因漏掉这一步在凌晨三点收到服务器邮件“NCL process killed”查了两小时才发现是某年12月观测全缺测。4.4 图形渲染异常线条断裂、标签错位、坐标轴消失NCL绘图有时会出现“幽灵bug”RMSE等高线变成断续虚线或模型标签飘到图外或横轴刻度全消失。这通常不是代码错误而是工作区缓存污染。NCL会复用之前绘图的图形资源若上次绘图设置了特殊字体或坐标系会影响本次。终极解决方案每次绘图前彻底重置工作区并新建; 删除旧工作区如果存在 if (isvar(wks)) then delete(wks) end if ; 新建工作区 wks gsn_open_wks(pdf, taylor_diagram_ models(i)) ; 强制清除所有图形资源 gsn_clear_frame(wks)此外避免在同一个NCL脚本中连续绘制多个泰勒图。每个图单独一个脚本用shell循环调用for model in CESM2 MPI-ESM NorESM2; do ncl model$model taylor_plot.ncl done这样每个进程独占资源互不干扰。5. 超越基础泰勒图的进阶应用与领域特化技巧5.1 多变量联合评估一张图看透模型“全科成绩”单变量泰勒图如只看降水只能评估一个方面。真实模型需要同时模拟降水、温度、风速等多个变量。这时可以构建多变量泰勒图矩阵。核心思想是为每个变量单独计算一套sigma,r,rmse然后在同一张图上用不同形状区分变量。; 假设有3个变量pr降水, tas温度, uas风速 ; 计算各自统计量得到6个数组 sigma_all (/sigma_pr_mod, sigma_tas_mod, sigma_uas_mod, \ sigma_pr_obs, sigma_tas_obs, sigma_uas_obs/) r_all (/r_pr_mod, r_tas_mod, r_uas_mod, 0, 0, 0/) ; 观测r恒为0 rmse_all (/rmse_pr_mod, rmse_tas_mod, rmse_uas_mod, 0, 0, 0/) ; 绘图时用不同marker resxyMarkerColors (/blue, red, green, black, black, black/) resxyMarkerSymbols (/circle, square, triangle, cross, cross, cross/) ; 观测用cross resxyLabels (/PR-CESM2, TAS-CESM2, UAS-CESM2, PR-OBS, TAS-OBS, UAS-OBS/)这样一个模型的三个变量点会形成一个小三角形直观显示其在不同物理量上的表现均衡性。如果三角形扁平如PR和TAS点近UAS点远说明模型在能量循环上协调但在动量输送上有缺陷。5.2 季节拆解泰勒图发现模型的“季节性失聪”很多模型全年RMSE尚可但夏季严重高估、冬季低估。这时将全年数据拆分为DJF冬、MAM春、JJA夏、SON秋四个季节分别画四张泰勒图。但更高效的做法是在一张图上用颜色编码季节。; 计算四季统计量得到4组值 sigma_seas (/sigma_djf, sigma_mam, sigma_jja, sigma_son/) r_seas (/r_djf, r_mam, r_jja, r_son/) rmse_seas (/rmse_djf, rmse_mam, rmse_jja, rmse_son/) ; 用颜色映射季节DJF蓝、MAM绿、JJA红、SON橙 resxyMarkerColors (/blue, green, red, orange/) resxyLabels (/DJF, MAM, JJA, SON/)观察点的分布如果四个点都聚集在高相关、中等标准差区域说明模型季节性稳健如果JJA点明显右移标准差大且下移相关低则锁定夏季模拟为短板。我在评估东亚季风降水时用此法发现某模式JJA点离观测最远进而定位到其对西太平洋副高位置的模拟偏差。5.3 不确定性量化给每个点画“误差椭圆”泰勒图上的点是确定值但实际统计量有抽样不确定性。例如40年数据算出的相关系数r0.75其95%置信区间可能是[0.62, 0.84]。忽略这点会误判模型优劣。NCL虽不内置但可用Bootstrap法估算; 对时间序列重采样1000次每次随机抽取40个时间点可放回 r_boot new(1000, float) do i0,999 idx generate_random_integer(0, dimsizes(pr_obs_1d)-1, dimsizes(pr_obs_1d)) r_boot(i) dim_correl_n(pr_mod_1d(idx), pr_obs_1d(idx), 0) end do r_ci quantile(r_boot, 0.025, 0.975) ; 得到95%置信区间然后在图上用gsn_add_polygon画一个椭圆长轴表示r的置信区间短轴表示sigma的置信区间。当两个模型的椭圆不重叠时才能说它们的性能差异是统计显著的——这才是严谨的模型比较。6. 最后一点个人体会泰勒图不是终点而是对话的起点画完一张完美的泰勒图把最优模型标红加粗发到群里收获一堆“厉害”这很爽。但真正的价值始于图完成之后。我坚持一个习惯每次画完泰勒图必定回到原始数据挑出泰勒图上表现最差的那个模型和表现最好的那个模型做一次逐格点、逐月份的偏差分析。比如发现某模型在夏季长江中下游RMSE奇高就打开它的降水场和观测场做差立刻看到它不是整体偏多而是把梅雨锋面系统向北推了2个纬度导致江南少雨、江淮多雨。这个机制性问题泰勒图本身不会告诉你但它像一个精准的探针帮你快速定位到该深挖的“靶区”。所以别把泰勒图当作结案陈词它应该是一份带坐标的邀请函——邀请你进入数据深处和模型展开一场关于物理过程的对话。我见过太多人把泰勒图当装饰画挂在论文里却从不追问“为什么这个点在这里”。而我的经验是每一个偏离观测的点背后都藏着一个等待被修正的物理参数化方案或一个尚未被充分理解的反馈过程。画图只是开始读懂图后的沉默才是工作真正开始的地方。