简介Goff-Gratch算法是气象与气候科学中估算饱和水汽压的经典标准方法由Goff和Gratch于1946年提出适用于湿度计算、露点温度推算及天气模式模拟等场景。压缩包内提供基于MATLAB的Goff_Gratch计算实现共2个文件包括1个.m脚本和1张湿度廓线图整体仅8KB结构精简、即下即用适合气象专业学生、科研人员及工程技术人员参考。MATLAB脚本内置温度修正项并针对冰点以下水汽升华过程进行专门处理可输出任意温度下的饱和水汽压值附带的湿度廓线图则以图形方式直观展示饱和水汽压随温度的变化趋势便于理解温度对空气中最大水汽含量的影响。资源虽小但涵盖了基础方程、温度修正与冰点修正等完整计算流程已有453人学习下载能够帮助读者快速掌握Goff-Gratch公式的实际编程应用为气候建模、农业水分管理或建筑能耗评估提供基础工具。1. 饱和水汽压计算为什么老工程师还在查 Goff-Gratch干过几年探空资料处理的人都有印象相对湿度、露点、比湿整套湿度链都系在饱和水汽压这一个数上。Goff-Gratch 公式虽然写在 1946 年至今仍被 WMO 当作基准方法尤其在对流层顶附近温度落到 -50℃ 以下Magnus 公式的误差会把相对湿度抬高好几个百分点辐射传输参数化跟着出错。这篇文章不是背公式而是把一份 MATLAB 实现拆开从单个Goff_Gratch.m函数扩展到探空数据批量处理和湿度廓线绘制最后加一个露点温度验证工具。适合正在处理 ERA5、探空报批数据或者手写陆面过程代码的人参考。2. Goff-Gratch 公式的结构与 MATLAB 实现2.1 水面方程指数项与对数项怎么配合Goff-Gratch 并不是一个简单多项式。它用三个指数项去逼近 Clausius-Clapeyron 方程在宽温度范围的行为。水面饱和水汽压记为e_w单位 hPa温度 T 用开尔文基准温度T0 273.16e_w 10^( -7.90298*(T0/T - 1) 5.02808*lg(T0/T) - 1.3816e-7*(10^(11.344*(1 - T/T0)) - 1) 8.1328e-3*(10^(-3.49149*(T0/T - 1)) - 1) lg(1013.246) )第一项主导高温区斜率第二项修正温度比值的非线性第三项和第四项分别处理高温端和低温端的气压偏离。1013.246是标准表面气压它把全式锚定到百帕量级。这里lg指常用对数MATLAB 里对应log10不要写成log否则结果会系统性偏小低温层尤其明显。2.2 冰面方程低温升华路径的独立系数当温度低于 0℃ 时水汽对冰面的饱和压低于水面因为相同温度下冰的化学势更低。Goff-Gratch 对冰面单独给出一套拟合系数用于有冰晶、霜的平流层下层或冬季地面观测e_i 10^( -9.09718*(T0/T - 1) - 3.56654*lg(T0/T) 0.876793*(1 - T/T0) lg(6.1071) )常数 6.1071 是冰面方程在 0℃ 附近的基准不是水面的 6.112。由于两套公式的拟合路径不同0℃ 点两侧算出的值会差约 0.005 hPa但这个不连续在气象业务里被广泛接受只要程序能按温度自动切换就足够。下面是公式中几个关键常数的含义写代码时对照着检查不容易抄错系数/常量数值作用T0273.16 K三相点温度基准1013.246hPa水面公式锚定常数6.1071hPa冰面公式基准值7.90298无量纲水面第一项主导蒸发斜率5.02808无量纲水面第二项修正对数偏差8.1328e-3无量纲水面第四项低温区气压修正2.3 Goff_Gratch.m 完整实现与逐段说明下面这个函数接收摄氏温度内部换算成开尔文然后分三段输出。默认把所有小于等于 0℃ 的层按冰面公式处理如果你在研究过冷水云需要把第 4.2 节的混合逻辑插进来替代这一段。function [es, flag] goff_gratch_svp(T) % GOFF_GRATCH_SVP - 饱和水汽压的 Goff-Gratch 实现 % T : 温度单位 ℃支持标量或向量 % es : 饱和水汽压单位 hPa % flag: 1 水面公式2 冰面公式0 超范围 T0 273.16; TK T T0; % 水面公式 term1 -7.90298 .* (T0 ./ TK - 1); term2 5.02808 .* log10(T0 ./ TK); term3 -1.3816e-7 .* (10.^(11.344 .* (1 - TK ./ T0)) - 1); term4 8.1328e-3 .* (10.^(-3.49149 .* (T0 ./ TK - 1)) - 1); lg_ew term1 term2 term3 term4 log10(1013.246); ew 10 .^ lg_ew; % 冰面公式 lg_ei -9.09718 .* (T0 ./ TK - 1) ... - 3.56654 .* log10(T0 ./ TK) ... 0.876793 .* (1 - TK ./ T0) ... log10(6.1071); ei 10 .^ lg_ei; % 分段 es NaN(size(T)); flag zeros(size(T)); mask (TK 173.16) (TK 373.16); flag(mask) 1; iceMask mask (TK T0); flag(iceMask) 2; es(iceMask) ei(iceMask); waterMask mask (TK T0); es(waterMask) ew(waterMask); endT T0对应气象资料里最常见的摄氏温度输入T0取 273.16 而不是 273.15因为公式原文以三相点为基准。10.^(...)的.^保证逐元素运算输入一组探空值时不会被当成矩阵幂。最后用NaN标记超出 -100~100℃ 的层后续循环里配合omitnan或isfinite即可跳过缺测。调用时只要一行[es, flag] goff_gratch_svp(temp)。flag能用于分层统计想单看冰面层的饱和水汽压就筛选flag 2想把水面公式用在过冷水场景需要提前改掉这里的分段策略。代码已经写成数组兼容不需要为每个高度写循环。3. 探空数据批量处理与湿度廓线绘制3.1 读入探空数据并统一单位探空数据常见的格式是每层有高度、温度、气压、相对湿度或露点。下载文件里温度可能是摄氏度也可能直接给开尔文气压可能是 hPa也可能是 Pa。写脚本第一步就把单位归一温度转成 ℃气压转成 hPa。下面的例子按height, temp_C, pres_hPa, RH四列读入 CSV。data readmatrix(sounding_data.csv); height data(:,1); % 高度米 temp data(:,2); % 温度摄氏度 pres data(:,3); % 气压百帕 rh data(:,4); % 相对湿度百分比 es goff_gratch_svp(temp); e rh .* es ./ 100;readmatrix适合有数值表头的简单文件如果第一行是站点名先用readcell读取并手动定位数据起始行。rh .* es ./ 100利用数组点乘把相对湿度换算成实际水汽压e单位和es保持一致。这里的rh必须是 0~100 的量纲如果上游数据是 0~1 小数去掉/100即可。3.2 由饱和水汽压推算比湿与相对湿度有些资料只给温度和露点没有 RH此时先用es goff_gratch_svp(temp)算出饱和水汽压再用露点温度对应的水汽压相除得到 RH。反过来如果有 RH 和温度可以用下面这个转换把湿度场补全e rh .* es ./ 100; % 实际水汽压 q 0.622 .* e ./ (pres - 0.378 .* e); % 比湿 kg/kg rh_check 100 .* q .* pres ./ (0.622 .* es 0.378 .* q .* es);0.622是干空气与水汽分子量之比0.378是同一比例在湿空气修正中的近似。第一行从 RH 得到实际水汽压第二行用气压和实际水汽压推算比湿第三行回代验算能有效检查单位是否正确。对于一般探空分析这个精度已经足够做数值模式输入时再用完整公式替换。转换目标表达式说明实际水汽压 eRH × es / 100es 来自 Goff-Gratch比湿 q0.622 × e / (p − 0.378 × e)百帕量级误差小于 0.1%RH 回验100 × q × p / (0.622 × es 0.378 × q × es)与原始 RH 对比3.3 用 MATLAB 画出湿度廓线湿度廓线.png就是上面结果的可视化横轴水汽压纵轴高度。由于饱和水汽压随温度指数下降图里红色实线通常在高空向左贴和蓝色实际水汽压线之间的距离直观反映了干层与湿层。下面脚本生成常用廓线图figure(Color, w); plot(es, height/1000, r-, LineWidth, 1.5); hold on; plot(e, height/1000, b--, LineWidth, 1.5); set(gca, YDir, reverse); xlabel(水汽压 (hPa)); ylabel(高度 (km)); legend({饱和水汽压 e_s, 实际水汽压 e}, Location, northeast); grid on; saveas(gcf, humidity_profile.png);set(gca,YDir,reverse)翻转 y 轴让高度从地面向上递增。横轴默认线性低层水汽压大值会挤压细节改成set(gca,XScale,log)能放大近地面区域但对数坐标下水汽压为 0 的缺测层会被自动剔除。通常两种图各存一份线性图用于报告对数图用于分析云底高度。4. 极端温度下的数值稳定性与常见误差排查4.1 NaN、缺测与数组传染Goff-Gratch 里的指数项在常规大气温度范围不会溢出真正导致计算失败的是资料里的缺测值。探空文件的 -999、-99 被换算成开尔文后仍然远超出公式范围函数会返回NaN而NaN参与后续 RH、比湿计算时会把一整层结果全部污染。下面的防护代码可以放在读取数据之后valid isfinite(temp) (temp -100) (temp 100); es(~valid) NaN; rh_clean rh .* valid;isfinite同时排除Inf和NaN温度限幅取 -100~100℃ 与第 2 章实现一致。rh_clean是逻辑索引乘以原始值缺测位置会变成 0再配合mean(..., omitnan)做统计时不会把缺测层算进去。4.2 0℃ 附近的相变不连续问题水面公式与冰面公式在 0℃ 处存在约 0.005 hPa 的跳变原因是两套系数分别来自不同的水面和冰面实测数据。对于高空冰晶层直接跳变可以接受但云微物理方案会同时用液态水和冰的饱和水汽压来判断相变。常见做法是加一个混合区间在 -38℃~0℃ 内线性过渡Tmix -38; % 混合相下限 frac max(0, min(1, (temp - Tmix) ./ (-Tmix))); es_mix frac .* es_ice (1 - frac) .* es_water;frac在 0℃ 时为 1在 -38℃ 时为 0。加权后的es_mix让水汽压曲线平滑穿过 0℃避免因不连续导致云底判断在温度小幅波动时来回跳动。这里的es_ice和es_water需要分别从同一温度计算不能直接用第 2 章函数的es单值。4.2.1 不要用单一冰面公式跨全温区有些脚本为了省事把所有低于 0℃ 的层都塞进冰面公式这会低估 -10~0℃ 过冷水云中的饱和水汽压相对湿度被高估 2% 左右云底高度被压得偏低。排查时看廓线在 0℃ 附近是否出现明显折点如果有多半是分段逻辑写在了调用端而不是函数内部。4.3 几个公式的取舍选公式不能只看常温段。下表把常见公式的边界条件放在一起便于根据场景判断公式类型适用温度范围0℃ 参考值主要风险Goff-Gratch 水面基准-70~100℃约 6.11 hPa温度分段不连续Goff-Gratch 冰面基准-100~0℃约 6.107 hPa系数长易抄错Magnus/Tetens简化-40~50℃约 6.11 hPa-20℃ 以下偏差放大Buck1981简化-80~50℃约 6.11 hPa平流层覆盖不足Magnus 公式在中纬度边界层表现不错但 ERA5 的高空层温度经常低于 -40℃这时 Goff-Gratch 仍能保持平滑下降。Buck 公式在 -40℃ 附近比 Magnus 好但系数是为水面设计的冰面场景没有独立处理。因此做探空垂直剖面时Goff-Gratch 依然是相对稳妥的默认选择。4.4 定位误差的一个思路如果输出 RH 出现高于 100% 的异常先看三个量es、e、RH。如果es是NaN查温度单位如果es正常而e异常查 RH 是不是 0~1 小数如果高温层正常、低温层 RH 突变多半是原始露点仪在低温下结冰滞后。把这些量打印在同一行能快速定位是公式问题还是前端数据问题。5. 用露点反演验证 Goff-Gratch 实现5.1 一个反演露点的迭代工具饱和水汽压算得对不对单看数值很难判断把露点温度反算回去是直观的验证给定温度 T 和相对湿度 RH先求实际水汽压 e再迭代反算露点看它与站点报告的露点是否一致。Goff-Gratch 没有显式反函数我用牛顿迭代十次以内通常收敛。function Td dew_point_from_vapor_pressure(e, Tinit) % 牛顿迭代求露点让 goff_gratch_svp(Td) 逼近 e % e : 实际水汽压hPa % Tinit : 初始温度℃ Td Tinit; for k 1:10 es goff_gratch_svp(Td); dT 0.01; df (goff_gratch_svp(TddT) - goff_gratch_svp(Td-dT)) / (2*dT); Td Td - (es - e) / df; if abs(es - e) 1e-6 break; end end end中心差分求导只需要调两次基础函数不用手工推导导数表达式1e-6作为收敛阈值对探空数据足够误差远小于传感器精度。Tinit最好给露点的大致范围比如temp - 5避免在 0℃ 分段处造成振荡。5.2 实测与验证步骤用 20℃、RH50% 的常见工况检查实际水汽压约为 11.69 hPaGoff-Gratch 反算露点应落在 9.3℃ 附近。再试 0℃、RH100%露点应回到 0℃如果结果偏离超过 0.1℃先查输入温度单位。这个验证器放在data_process流程里每次处理完探空资料就自动对一遍站点探空报的露点能在上传数据前发现单位错误和高度错位。更进一步可以把dew_point_from_vapor_pressure用于缺失露点层的插补对同一高度温度场计算temp - Td得到温度露点差再按层结条件判断云区。这样整条湿度链从 Goff-Gratch 单点函数扩展到批量廓线分析不需要额外安装工具箱。把这个函数保存为dew_point_iter.m直接在探空处理脚本里调用即可输出每个标准等压面的露点差再叠加到humidity_profile.png上查看云层与环境湿度结构。本文还有配套的精品资源点击获取 Goff-Gratch饱和水汽压MATLAB实现及探空资料处理
简介Goff-Gratch算法是气象与气候科学中估算饱和水汽压的经典标准方法由Goff和Gratch于1946年提出适用于湿度计算、露点温度推算及天气模式模拟等场景。压缩包内提供基于MATLAB的Goff_Gratch计算实现共2个文件包括1个.m脚本和1张湿度廓线图整体仅8KB结构精简、即下即用适合气象专业学生、科研人员及工程技术人员参考。MATLAB脚本内置温度修正项并针对冰点以下水汽升华过程进行专门处理可输出任意温度下的饱和水汽压值附带的湿度廓线图则以图形方式直观展示饱和水汽压随温度的变化趋势便于理解温度对空气中最大水汽含量的影响。资源虽小但涵盖了基础方程、温度修正与冰点修正等完整计算流程已有453人学习下载能够帮助读者快速掌握Goff-Gratch公式的实际编程应用为气候建模、农业水分管理或建筑能耗评估提供基础工具。1. 饱和水汽压计算为什么老工程师还在查 Goff-Gratch干过几年探空资料处理的人都有印象相对湿度、露点、比湿整套湿度链都系在饱和水汽压这一个数上。Goff-Gratch 公式虽然写在 1946 年至今仍被 WMO 当作基准方法尤其在对流层顶附近温度落到 -50℃ 以下Magnus 公式的误差会把相对湿度抬高好几个百分点辐射传输参数化跟着出错。这篇文章不是背公式而是把一份 MATLAB 实现拆开从单个Goff_Gratch.m函数扩展到探空数据批量处理和湿度廓线绘制最后加一个露点温度验证工具。适合正在处理 ERA5、探空报批数据或者手写陆面过程代码的人参考。2. Goff-Gratch 公式的结构与 MATLAB 实现2.1 水面方程指数项与对数项怎么配合Goff-Gratch 并不是一个简单多项式。它用三个指数项去逼近 Clausius-Clapeyron 方程在宽温度范围的行为。水面饱和水汽压记为e_w单位 hPa温度 T 用开尔文基准温度T0 273.16e_w 10^( -7.90298*(T0/T - 1) 5.02808*lg(T0/T) - 1.3816e-7*(10^(11.344*(1 - T/T0)) - 1) 8.1328e-3*(10^(-3.49149*(T0/T - 1)) - 1) lg(1013.246) )第一项主导高温区斜率第二项修正温度比值的非线性第三项和第四项分别处理高温端和低温端的气压偏离。1013.246是标准表面气压它把全式锚定到百帕量级。这里lg指常用对数MATLAB 里对应log10不要写成log否则结果会系统性偏小低温层尤其明显。2.2 冰面方程低温升华路径的独立系数当温度低于 0℃ 时水汽对冰面的饱和压低于水面因为相同温度下冰的化学势更低。Goff-Gratch 对冰面单独给出一套拟合系数用于有冰晶、霜的平流层下层或冬季地面观测e_i 10^( -9.09718*(T0/T - 1) - 3.56654*lg(T0/T) 0.876793*(1 - T/T0) lg(6.1071) )常数 6.1071 是冰面方程在 0℃ 附近的基准不是水面的 6.112。由于两套公式的拟合路径不同0℃ 点两侧算出的值会差约 0.005 hPa但这个不连续在气象业务里被广泛接受只要程序能按温度自动切换就足够。下面是公式中几个关键常数的含义写代码时对照着检查不容易抄错系数/常量数值作用T0273.16 K三相点温度基准1013.246hPa水面公式锚定常数6.1071hPa冰面公式基准值7.90298无量纲水面第一项主导蒸发斜率5.02808无量纲水面第二项修正对数偏差8.1328e-3无量纲水面第四项低温区气压修正2.3 Goff_Gratch.m 完整实现与逐段说明下面这个函数接收摄氏温度内部换算成开尔文然后分三段输出。默认把所有小于等于 0℃ 的层按冰面公式处理如果你在研究过冷水云需要把第 4.2 节的混合逻辑插进来替代这一段。function [es, flag] goff_gratch_svp(T) % GOFF_GRATCH_SVP - 饱和水汽压的 Goff-Gratch 实现 % T : 温度单位 ℃支持标量或向量 % es : 饱和水汽压单位 hPa % flag: 1 水面公式2 冰面公式0 超范围 T0 273.16; TK T T0; % 水面公式 term1 -7.90298 .* (T0 ./ TK - 1); term2 5.02808 .* log10(T0 ./ TK); term3 -1.3816e-7 .* (10.^(11.344 .* (1 - TK ./ T0)) - 1); term4 8.1328e-3 .* (10.^(-3.49149 .* (T0 ./ TK - 1)) - 1); lg_ew term1 term2 term3 term4 log10(1013.246); ew 10 .^ lg_ew; % 冰面公式 lg_ei -9.09718 .* (T0 ./ TK - 1) ... - 3.56654 .* log10(T0 ./ TK) ... 0.876793 .* (1 - TK ./ T0) ... log10(6.1071); ei 10 .^ lg_ei; % 分段 es NaN(size(T)); flag zeros(size(T)); mask (TK 173.16) (TK 373.16); flag(mask) 1; iceMask mask (TK T0); flag(iceMask) 2; es(iceMask) ei(iceMask); waterMask mask (TK T0); es(waterMask) ew(waterMask); endT T0对应气象资料里最常见的摄氏温度输入T0取 273.16 而不是 273.15因为公式原文以三相点为基准。10.^(...)的.^保证逐元素运算输入一组探空值时不会被当成矩阵幂。最后用NaN标记超出 -100~100℃ 的层后续循环里配合omitnan或isfinite即可跳过缺测。调用时只要一行[es, flag] goff_gratch_svp(temp)。flag能用于分层统计想单看冰面层的饱和水汽压就筛选flag 2想把水面公式用在过冷水场景需要提前改掉这里的分段策略。代码已经写成数组兼容不需要为每个高度写循环。3. 探空数据批量处理与湿度廓线绘制3.1 读入探空数据并统一单位探空数据常见的格式是每层有高度、温度、气压、相对湿度或露点。下载文件里温度可能是摄氏度也可能直接给开尔文气压可能是 hPa也可能是 Pa。写脚本第一步就把单位归一温度转成 ℃气压转成 hPa。下面的例子按height, temp_C, pres_hPa, RH四列读入 CSV。data readmatrix(sounding_data.csv); height data(:,1); % 高度米 temp data(:,2); % 温度摄氏度 pres data(:,3); % 气压百帕 rh data(:,4); % 相对湿度百分比 es goff_gratch_svp(temp); e rh .* es ./ 100;readmatrix适合有数值表头的简单文件如果第一行是站点名先用readcell读取并手动定位数据起始行。rh .* es ./ 100利用数组点乘把相对湿度换算成实际水汽压e单位和es保持一致。这里的rh必须是 0~100 的量纲如果上游数据是 0~1 小数去掉/100即可。3.2 由饱和水汽压推算比湿与相对湿度有些资料只给温度和露点没有 RH此时先用es goff_gratch_svp(temp)算出饱和水汽压再用露点温度对应的水汽压相除得到 RH。反过来如果有 RH 和温度可以用下面这个转换把湿度场补全e rh .* es ./ 100; % 实际水汽压 q 0.622 .* e ./ (pres - 0.378 .* e); % 比湿 kg/kg rh_check 100 .* q .* pres ./ (0.622 .* es 0.378 .* q .* es);0.622是干空气与水汽分子量之比0.378是同一比例在湿空气修正中的近似。第一行从 RH 得到实际水汽压第二行用气压和实际水汽压推算比湿第三行回代验算能有效检查单位是否正确。对于一般探空分析这个精度已经足够做数值模式输入时再用完整公式替换。转换目标表达式说明实际水汽压 eRH × es / 100es 来自 Goff-Gratch比湿 q0.622 × e / (p − 0.378 × e)百帕量级误差小于 0.1%RH 回验100 × q × p / (0.622 × es 0.378 × q × es)与原始 RH 对比3.3 用 MATLAB 画出湿度廓线湿度廓线.png就是上面结果的可视化横轴水汽压纵轴高度。由于饱和水汽压随温度指数下降图里红色实线通常在高空向左贴和蓝色实际水汽压线之间的距离直观反映了干层与湿层。下面脚本生成常用廓线图figure(Color, w); plot(es, height/1000, r-, LineWidth, 1.5); hold on; plot(e, height/1000, b--, LineWidth, 1.5); set(gca, YDir, reverse); xlabel(水汽压 (hPa)); ylabel(高度 (km)); legend({饱和水汽压 e_s, 实际水汽压 e}, Location, northeast); grid on; saveas(gcf, humidity_profile.png);set(gca,YDir,reverse)翻转 y 轴让高度从地面向上递增。横轴默认线性低层水汽压大值会挤压细节改成set(gca,XScale,log)能放大近地面区域但对数坐标下水汽压为 0 的缺测层会被自动剔除。通常两种图各存一份线性图用于报告对数图用于分析云底高度。4. 极端温度下的数值稳定性与常见误差排查4.1 NaN、缺测与数组传染Goff-Gratch 里的指数项在常规大气温度范围不会溢出真正导致计算失败的是资料里的缺测值。探空文件的 -999、-99 被换算成开尔文后仍然远超出公式范围函数会返回NaN而NaN参与后续 RH、比湿计算时会把一整层结果全部污染。下面的防护代码可以放在读取数据之后valid isfinite(temp) (temp -100) (temp 100); es(~valid) NaN; rh_clean rh .* valid;isfinite同时排除Inf和NaN温度限幅取 -100~100℃ 与第 2 章实现一致。rh_clean是逻辑索引乘以原始值缺测位置会变成 0再配合mean(..., omitnan)做统计时不会把缺测层算进去。4.2 0℃ 附近的相变不连续问题水面公式与冰面公式在 0℃ 处存在约 0.005 hPa 的跳变原因是两套系数分别来自不同的水面和冰面实测数据。对于高空冰晶层直接跳变可以接受但云微物理方案会同时用液态水和冰的饱和水汽压来判断相变。常见做法是加一个混合区间在 -38℃~0℃ 内线性过渡Tmix -38; % 混合相下限 frac max(0, min(1, (temp - Tmix) ./ (-Tmix))); es_mix frac .* es_ice (1 - frac) .* es_water;frac在 0℃ 时为 1在 -38℃ 时为 0。加权后的es_mix让水汽压曲线平滑穿过 0℃避免因不连续导致云底判断在温度小幅波动时来回跳动。这里的es_ice和es_water需要分别从同一温度计算不能直接用第 2 章函数的es单值。4.2.1 不要用单一冰面公式跨全温区有些脚本为了省事把所有低于 0℃ 的层都塞进冰面公式这会低估 -10~0℃ 过冷水云中的饱和水汽压相对湿度被高估 2% 左右云底高度被压得偏低。排查时看廓线在 0℃ 附近是否出现明显折点如果有多半是分段逻辑写在了调用端而不是函数内部。4.3 几个公式的取舍选公式不能只看常温段。下表把常见公式的边界条件放在一起便于根据场景判断公式类型适用温度范围0℃ 参考值主要风险Goff-Gratch 水面基准-70~100℃约 6.11 hPa温度分段不连续Goff-Gratch 冰面基准-100~0℃约 6.107 hPa系数长易抄错Magnus/Tetens简化-40~50℃约 6.11 hPa-20℃ 以下偏差放大Buck1981简化-80~50℃约 6.11 hPa平流层覆盖不足Magnus 公式在中纬度边界层表现不错但 ERA5 的高空层温度经常低于 -40℃这时 Goff-Gratch 仍能保持平滑下降。Buck 公式在 -40℃ 附近比 Magnus 好但系数是为水面设计的冰面场景没有独立处理。因此做探空垂直剖面时Goff-Gratch 依然是相对稳妥的默认选择。4.4 定位误差的一个思路如果输出 RH 出现高于 100% 的异常先看三个量es、e、RH。如果es是NaN查温度单位如果es正常而e异常查 RH 是不是 0~1 小数如果高温层正常、低温层 RH 突变多半是原始露点仪在低温下结冰滞后。把这些量打印在同一行能快速定位是公式问题还是前端数据问题。5. 用露点反演验证 Goff-Gratch 实现5.1 一个反演露点的迭代工具饱和水汽压算得对不对单看数值很难判断把露点温度反算回去是直观的验证给定温度 T 和相对湿度 RH先求实际水汽压 e再迭代反算露点看它与站点报告的露点是否一致。Goff-Gratch 没有显式反函数我用牛顿迭代十次以内通常收敛。function Td dew_point_from_vapor_pressure(e, Tinit) % 牛顿迭代求露点让 goff_gratch_svp(Td) 逼近 e % e : 实际水汽压hPa % Tinit : 初始温度℃ Td Tinit; for k 1:10 es goff_gratch_svp(Td); dT 0.01; df (goff_gratch_svp(TddT) - goff_gratch_svp(Td-dT)) / (2*dT); Td Td - (es - e) / df; if abs(es - e) 1e-6 break; end end end中心差分求导只需要调两次基础函数不用手工推导导数表达式1e-6作为收敛阈值对探空数据足够误差远小于传感器精度。Tinit最好给露点的大致范围比如temp - 5避免在 0℃ 分段处造成振荡。5.2 实测与验证步骤用 20℃、RH50% 的常见工况检查实际水汽压约为 11.69 hPaGoff-Gratch 反算露点应落在 9.3℃ 附近。再试 0℃、RH100%露点应回到 0℃如果结果偏离超过 0.1℃先查输入温度单位。这个验证器放在data_process流程里每次处理完探空资料就自动对一遍站点探空报的露点能在上传数据前发现单位错误和高度错位。更进一步可以把dew_point_from_vapor_pressure用于缺失露点层的插补对同一高度温度场计算temp - Td得到温度露点差再按层结条件判断云区。这样整条湿度链从 Goff-Gratch 单点函数扩展到批量廓线分析不需要额外安装工具箱。把这个函数保存为dew_point_iter.m直接在探空处理脚本里调用即可输出每个标准等压面的露点差再叠加到humidity_profile.png上查看云层与环境湿度结构。本文还有配套的精品资源点击获取