1 引言

潜热通量改变影响大气过程，准确量化地表水热通量对于认识与理解水循环演变、气候与水资源变化过程具有重要意义。基于Jarvis冠层阻力方案的Penman-Monteith（P-M）模型估计的潜热通量结果（Jarvis, 1976）依赖于最小冠层阻力参数（Kumar et al., 2014），它与冠层内部特性密切相关，在常规方法中将其设为定值的做法导致了模拟的不确定性（Chen and Dudhia, 2001; Douglas et al., 2009）。刘春伟等（2016）模拟冬小麦潜热通量过程中发现在生育后期潜热模拟值明显高估，原因在于该方案忽略了最小冠层阻力的季节变化特征（Bie et al., 2015）；Alfieri et al.（2008）通过对最小冠层阻力的估算指出其具有日变化及更长时间尺度的变化特征，不能简单设为定值（Kumar et al., 2014）。同时，敏感性研究表明定值方案的不足在陆地水循环模拟研究中将对气候预估产生极大偏差（Alapaty et al., 1997; Cooter and Schwede, 2000）。为了规避最小冠层阻力为定值的不足，通过分时段修正最小冠层阻力值改善模拟性能，并在不同生育期小麦蒸散估计应用中取得了较好效果（Fang et al., 2014；陈志忠，2016），然而这种修正方案使最小冠层阻力值在生育过渡期不连续，同时也没有合理刻画其在全生育期中的季节分布特征。

本文利用多年高质量通量、气象（风速、温文、湿度和压强）及MODIS卫星的地表叶面积指数数据，反向利用P-M模型，采用整体和分段相结合的最小冠层阻力拟合迭代算法，将最小冠层阻力分布扩展至全年以表现其在全生育期的季节分布特征，弥补常规Jarvis阻力参数化中最小冠层阻力不随生育期变化的缺陷，进而改善潜热通量模拟，并以宁乡通量观测站为例进行验证。同时，值得关注的是受湍流瞬时变化、大气层结状态、降水、电力故障和人为等因素潜热数据连续性和数据质量也存在不确定性（Jia et al., 2015），针对该站点涡度相关通量缺测问题，应用本文改进后的冠层阻力方案以模拟值插补缺测数据，为涡度相关潜热通量数据完整性提供新的方式。同时采用P-M模型模拟潜热通量还依赖其他通量及环境数据，为了定量评估数据输入误差对模拟精度的影响，本文对模型采用微分误差分析方法，以确定不同变量对潜热通量模拟的误差贡献。

2 研究方法 2.1 通量估计方法

P-M模型：

$ \lambda E = \frac{{\Delta ({R_{\rm{n}}} - {G_{\rm{5}}} - S) + \rho \cdot {c_p} \cdot ({e_{\rm{s}}} - e)/{r_{\rm{a}}}}}{{\Delta + \gamma (1 + {r_{\rm{s}}}/{r_{\rm{a}}})}}, $

(1)

其中，λE为潜热通量（单位：W m⁻²）；ρ为空气密度（单位：kg m⁻³）；c_p为空气比热（单位：J kg⁻¹ ℃⁻¹）；e_s为饱和水汽压（单位：Pa）；e为实际水汽压（单位：Pa）；R_n为净辐射（单位：W m⁻²）；Δ为饱和水汽压随温度变化率（单位：Pa ℃⁻¹）；γ为干湿球常数；G₅为5 cm土壤热通量（单位：W m⁻²）；储热项S包含0~5 cm的土壤储热S_soil（单位：W m⁻²）、潜热储存S_LE（单位：W m⁻²）、显热储存S_H（单位：W m⁻²）和冠层储热S_canopy（单位：W m⁻²），即S＝S_soil+S_LE+S_H+S_canopy，储热项详细计算参见文献（Meesters and Vugts, 1996）；r_s、r_a分别为冠层阻力和空气动力学阻力（单位：s m⁻¹）。参数敏感性分析中已有研究表明P-M模型中冠层阻力对模型模拟敏感性高于边界层阻力（Ortega-Farias et al., 2007; 李俊等，2014），这也表明了本文通过改善冠层阻力方案提高模拟能力具有合理性，本文不做过多赘述。通量估计流程如图 1所示。

在通量估计流程中本文设计了整体与分段拟合相结合的迭代算法，发展了基于最小冠层阻力分布的阻力方案，从而改进了原始冠层阻力方案。这是本文有别于同类P-M模型最主要的方面，敏感性分析旨在表明冠层阻力对潜热模拟的重要性，参数率定采用最小二乘法以得到最优参数（Kişi，2006）。

（1）最小冠层阻力分布拟合与参数估计

冠层阻力方案采用Jarvis冠层阻力模型[公式（2）]，其表述为最小冠层阻力与环境因子胁迫函数的乘积形式；叶面积、净辐射、水汽压差、空气温度和土壤湿度胁迫的参数化方案见公式（3）至（7），本文为了兼顾参数连续性和模拟效果，本文采用了新的最小冠层阻力拟合迭代算法，基于整体参数拟合以保持与净辐射、水汽压差、空气温度和土壤湿度相关的参数在全生育期的连续性；并借助分段参数拟合以限定最小冠层阻力的分布并且降低不同生育期的模拟误差。分段参数拟合以两个月为间隔，每个间隔期采用最小二乘非线性同时拟合得出每个间隔期的气象胁迫参数（叶面积、土壤湿度、水汽压差、温度）。根据公式（8），最小冠层阻力由分段和整体拟合相结合共同得到分布，它具体迭代流程见图 2；并重复上述流程经过迭代以得到更好模拟效果，当前后两次迭代潜热模拟值变化小于阈值，或者模拟精度降低时迭代终止。

$ {r_{\rm{s}}} = {r_{\rm{0}}}\frac{1}{{f({\rm{LAI}}) \cdot f({R_{\rm{n}}}) \cdot f({\rm{VPD}}) \cdot f(T) \cdot f(\theta )}}, $

(2)

$ f({\rm{LAI}}) = \frac{{{\rm{LAI}}}}{{1 + 0.5{\rm{LAI}}}}, $

(3)

$ f({R_{\rm{n}}}) = 1 - {e^{ - {R_n}/{a_1}}}, $

(4)

$ f({p_{{\rm{VPD}}}}) = 1 - {a_2} \cdot {p_{{\rm{VPD}}}}, $

(5)

$ f(T) = 1 - {a_3} \cdot {(28 - T)^2}, $

(6)

$ f(\theta ) = {a_4}\frac{{\theta - {\theta _{\rm{l}}}}}{{{\theta _{\rm{u}}} - {\theta _{\rm{l}}}}}. $

(7)

其中，r₀为最小冠层阻力（单位：s m⁻¹），其表征的是下垫面未受环境强迫状况时的冠层阻力；LAI为叶面积指数；R_n为净辐射（单位：W m⁻²）；p_VPD为水汽压差（单位：Pa）；T为空气温度（单位：℃）；θ为土壤湿度；θ_u、θ_l分别为土壤湿度上下限；a₁至a₄为净辐射、水汽压差、温度和土壤湿度的胁迫参数。

$ f({r_0}) = \frac{{{r_{{\rm{0}},i}}f({\rm{LAI}})f({R_{\rm{n}}})f({p_{{\rm{VPD}}}})f(T)f(\theta )}}{{{f_i}({\rm{LAI}}){f_i}({R_{\rm{n}}}){f_i}({p_{{\rm{VPD}}}}){f_i}(T){f_i}(\theta )}},\;\;i = 1,2 \cdot \cdot \cdot 6. $

(8)

其中，f(r₀)为求得的最小冠层阻力分布；r_{0, i}为分段最小冠层阻力（分段初值设置为60 s m⁻¹）；f(LAI)、f(R_n)、f(p_VPD)、f(T)、f(θ)分别为整体胁迫函数，其参数a₁至a₄取值为整体参数拟合结果；f_i(LAI)、f_i (R_n)、f_i (p_VPD)、f_i (T)、f_i (θ)为分段胁迫函数，其参数取值a_{1, i}至a_{4, i}为分段拟合出的参数结果。

参数a₁至a₄初始、限制区间和拟合值见表 1。基于迭代出的最小冠层阻力分布、气象及MODIS叶面积指数数据采用P-M模型计算潜热通量。其中，2012~2015年数据用于最小冠层阻力分布拟合，2016年数据用于独立验证。

（2）高质量通量数据遴选

遴选高质量数据旨在为数据拟合年均最小冠层阻力分布及数据模拟及验证提供可靠性。文章采用潜热和显热阈值筛选、降雨时次剔除和基于湍流发展和平稳性的质量分级方案遴选宁乡通量观测站已有通量数据。根据局地特点，对潜热通量阈值设置为－50~600 W m⁻²，显热通量阈值设置为－100~60 W m⁻²；依照雨量筒数据按时次剔除降雨时段通量数据；基于Foken and Wichura（1996）给出的湍流发展性和湍流平稳性分级方案，选择“0级”（最优）通量数据。

2.2 评估方法

采用多个统计量用以评估模拟性能，包括相关系数（R）、均方根误差（RMSE）和一致性指数（IA）。相关性系数（R）：

$ R = \frac{{\sum\limits_n {\left( {o - \overline o } \right)\left( {s - \overline s } \right)} }}{{\sqrt {\sum\limits_n {{{\left( {o - \overline o } \right)}^2}} \sum\limits_n {{{\left( {s - \overline s } \right)}^2}} } }}, $

(9)

均方根误差（RMSE）：

$ {\rm{RMSE}} = \sqrt {\frac{1}{n}\sum\limits_n {{{\left( {o - s} \right)}^2}} } , $

(10)

一致性指数（IA）：

$ {\rm{IA}} = 1 - \frac{{\sum\limits_n {{{\left( {s - o} \right)}^2}} }}{{\sum\limits_n {{{\left( {\left| {s - \overline s } \right| + \left| {o - \overline o } \right|} \right)}^2}} }}, $

(11)

其中，o为观测值，s为模拟值，o为观测均值，s为模拟均值。

3 模型验证与应用 3.1 站点概况与数据

通量观测站位于湖南省宁乡县（28°20′N，112°34′E）属东亚季风区，海拔高度110 m，年均温17.6℃，年降水1458 mm。观测下垫面相对平坦均质以阔叶和针叶林为主，冠层高度约7.5 m。该通量观测站连续通量观测始于2012年8月，对通量及气象观测达4 a。

通量观测塔高20 m，距地17.5 m安装三维超声风速仪（CSAT3，Campbell公司，USA）和红外气体分析仪（EC150，Campbell公司，USA）监测三维风速及水汽和CO₂浓度，原始数据以10 Hz频率记录于CR3000数据采集器（Campbell公司，USA）。同高度设置四分量辐射计（CNR4，Kipp & Zonen公司，NL）同步观测净辐射。在观测塔9.5 m和13.5 m两个高度放置普通风速计和温湿度传感器同步观测风速和温湿度。塔顶安置雨量筒观测降水状况，塔下测量不同土壤深度（0.05 m、0.1 m、0.2 m、0.4 m、0.6 m、0.8 m、1 m）温湿度，并在地下0.05 m深度放置土壤热通量板记录土壤热通量，上述数据以半小时间隔存储于CR1000（Campbell公司，USA）数据采集器。

选用EddyPro软件将原始涡度相关数据处理为半小时的通量数据。先后经过仪器表征量剔除、5倍方差异常值剔除、坐标旋转、频率损失订正、超声虚温订正和水汽浓度修正等质量控制（Lee et al., 2004）。通量数据R_n－G₅－S（净辐射、5cm土壤热通量和储热项的差值）与H+LE（显热和潜热之和）相关系数为0.93，平均回归斜率1.18；能量闭合度日均值为0.52，略低于一般林地的能量闭合度，能量残差（$\delta = {R_{\rm{n}}} - {G_5} - S - {\rm{LE}} - H$）日均范围为2.44±21.9 W m⁻²。

观测站点并无叶面积实测数据，故以宁乡通量站为中心，1 km范围的MODIS叶面积反演数据代替。该数据时间间隔为8 d，对每年的3~11月采用多项式拟合，1~2月和12月采用数据插值得到叶面积的半小时变化数据（图 3）。本文数据采集过程中由于可规避的人为疏忽，导致部分净辐射、土壤热通量和储热项存在缺失。本文采用该站点CLDAS（中国气象局陆面数据同化系统）数据（Shi et al., 2014）以代替。而关于输入数据不确定性产生的模拟误差将在本文3.4节误差分析中阐述。

图 4为宁乡通量观测站2012年8月至2016年6月通量观测数据质量比例的日均分布图；受仪器不确定性及人为失误的影响，缺测数据占比超过总量的20%，降水期间的数据占比也较大。可用数据约为50%左右，夜间受弱湍流影响占比更低，具有白天高夜间低的日变化特征。

3.2 通量估计和模型验证

经两次迭代终止，最小冠层阻力年均分布见图 5其分布规律存在夏低、冬高的季节变化特征，其中春、秋季分别呈快速下降和上升趋势，变化区间在20.8~90.1 s m⁻¹，5次多元回归拟合曲线如下：

$ \begin{array}{c} y = - 1.2 \times {10^{ - 9}}{x^5} + 1 \times {10^{ - 6}}{x^4} - 3.1 \times {10^{ - 4}}{x^3} + \\ 0.037{x^2} - 1.6x + 85, \end{array} $

(12)

其中，x为日序，y为年均最小冠层阻力。

以2012~2016年改进前后的模拟与实测潜热通量数据的1:1线可知，改进后的模拟效果好于原模型方案（图 6）。基于最小冠层阻力分布的模拟值与实测数据拟合回归系数由原来的1.11降低至0.99。全天数据对比来看，原方案整体呈高估现象，改进后潜热通量高估现状有所缓减，白天数据对比亦是如此，同时改进后模型对夜间模拟效果也有一定改善。从绝对误差来看，改进前后全天误差分别为36.74、35.60 W m⁻²，改进后模拟效果稍好于原始方案的模型。

改进后的各项统计量相比于原始方案也有一定程度提高（表 2），且白天模拟效果好于夜间。改进后整体数据的相关系数（R）为0.86，均方根误差（RMSE）为58.46，一致性指数（IA）0.91，而昼夜R、RMSE和IA分别为0.83、71.8、0.89（昼）和0.1、28.62、0.26（夜）。2016年的独立数据验证也存在白天潜热通量模拟效果好于夜间现象，其中夜间受潜热通量量级小故RMSE较小，其他统计量数值上夜间小于白天表明夜间潜热通量模拟效果不如白天。误差水平来看，总体绝对误差在50 W m⁻²左右，同样的是受量级影响白天绝对误差大于夜间。

对2012~2016年潜热通量日变化按照季节平均后，对比改进前后模拟与实测（图 7）可知，改进后的模拟值能更好刻画潜热通量的夏高冬低的季节变化特征，模拟值相比于实测潜热（春季到冬季潜热实测依次为：53.44、97.02、49.01、23.26 W m⁻²），平均季节潜热模拟误差（春季到冬季平均误差依次为：11.85、13.01、21.67、12.75 W m⁻²）小于原方案（春季到冬季平均误差依次为：12.94、17.27、22.49、12.76 W m⁻²），其中秋季绝对误差较大，而夏季误差偏离程度相对较大。日平均的变化趋势方面两者也非常相似，能表现潜热白天高夜间低的日变化特征，较实测夏天的白天潜热模拟值略高估而其他季节白天有低估现象。此外，对个别时段分析也具有高估现象，主要集中在2013年夏季，可能受当地骤旱影响导致模拟性能降低（Xie et al., 2016），而P-M模型无法及时刻画短时气象要素骤变对植被带来的负影响。年际模拟来看，2014和2015年模拟与实测绝对误差较低，且统计指标优异，表明具有好的模拟效果，2016年独立数据验证平均误差稍大于2013~2015模拟期，但同样也表明具有较好效果。

3.3 缺失数据插补应用

应对站点潜热通量缺测问题，采用本文给出的最小冠层阻力分布方案，基于P-M模型模拟用于通量缺失数据插补。选取2014年7月18日至8月7日连续通量数据，此数据段较完整优质，分别验证单数据缺失和多数据连续缺失的插补性能，各插补方法见表 3。单数据缺失插补性能比较中，从数据段始末依次剔除单个数据，分析误差统计量来评估插补性能；多数据连续插补性能比较中，采用相似方法，从数据段始末连续滑动删除多个数据，并对比连续缺失数据的插补性能。

本文方案插补单点及连续缺失数据方面都表现优异，与实测最为接近。从插补单点缺失数据性能对比来看（图 8），非线性插值、日变化趋势插补法、查表法与本文方案皆能有效捕捉潜热通量日变化特征，以2014年8月4日为例，单点插补都能表现当日潜热的双峰现象；查表法受样本数量限制，插补表现较为离散。四种方法插值误差统计见图 9。对比表明，各方法标准差都在160 W m⁻²附近，本文插补方案RMSE和R好于其他3种方法，相关系数达0.93，RMSE为60.7 W m⁻²，单点缺失插补效果最优；日变化趋势插补法各统计量表现也较好，仅次于本文方法。

从连续缺失数据插补方面来看，本文方法优于其他插补方案。连续缺失数据个数从2升至8，各统计量变化趋势见图 10。结果可知，本文方法和查表法插补性能稳定，各统计量不随连续缺失数据的改变而改变；而非线性插值和日变化随连续缺失数据的增多，RMSE增大，R和IA呈快速下降趋势，插补性能降低，且非线性插值法降低更明显。当连续缺失点超过2个时，本方法插补效果皆保持最优；当连续缺失点超4个时，查表法插补效果好于非线性插值和日变化趋势插补法，但仍次于本文方法。

3.4 数据不确定性分析

本文对潜热通量模拟及插补有赖于其他通量数据（包含净辐射、土壤热通量和土壤储热）及阻力项（冠层阻力和边界层阻力）数据，这些数据的不确定性会给潜热通量的模拟和插补带来不同程度的影响，同时，本文模型具有明确表达式，这里对P-M模型采用微分误差分析方法来预估输入数据的不确定性产生的误差贡献。具体方法如下：

分别对净辐射、土壤热通量、储热项、冠层阻力和边界层阻力求导（刘斌等，2015），并将结果中微分以差分形式代替，然后基于宁乡通量站的相关量级关系（$R \approx 2{E_{\rm{a}}}$，$\Delta \approx 3\gamma $，${r_{\rm{a}}} \approx {r_{\rm{s}}}$，${G_5} \approx 0.005{R_{\rm{n}}}$，$S \approx 0.18{R_{\rm{n}}}$，其中$R = {R_{\rm{n}}} - {G_{\rm{5}}} - S$，${E_a} = \frac{{\rho \cdot {c_p} \cdot ({e_s} - e)}}{{\gamma \cdot {r_a}}}$，各变量含义同上），得出净辐射、土壤热通量、储热项和阻力项引起最大潜热相对误差如下：

$ \begin{array}{c} {\left. {{\delta _{{\rm{input}}}}} \right|_{\max }} = 1.05\left| {\frac{{\Delta {R_{\rm{n}}}}}{{{R_{\rm{n}}}}}} \right| + 0.05\left| {\frac{{\Delta {G_5}}}{{{G_5}}}} \right| + 0.19\left| {\frac{{\Delta S}}{S}} \right| + \\ 0.25\left| {\frac{{\Delta {r_{\rm{s}}}}}{{{r_{\rm{s}}}}}} \right| + 0.06\left| {\frac{{\Delta {r_{\rm{a}}}}}{{{r_{\rm{a}}}}}} \right|. \end{array} $

(13)

潜热模拟误差来源与各输入分量息息相关，随各分项相对误差的增大而增大。若模型各分量相对误差相同时，数据输入项中的净辐射、冠层阻力、储热项、空气动力学阻力和土壤热通量对插补的潜热误差贡献依次降低，其贡献比例依次为0.65、0.16、0.12、0.04和0.03。净辐射的对潜热误差贡献最大，应尽量保持净辐射数据的连续性，数据替代的做法需要更加谨慎。此外，本文净辐射主要通过观测方法获得，在这里我们并不考虑仪器的观测误差。忽略净辐射对潜热通量的误差贡献后，冠层阻力对潜热误差贡献最大，其连续观测较难直接获得，且对模型有高敏感性（刘斌等，2014；Ershadi et al., 2015），故本文基于最小冠层阻力多年平均分布改善了冠层阻力的季节响应的做法是合理的。根据既有资料，土壤热通量影响较小，忽略土壤热通量对潜热通量误差影响不及储热项和边界层阻力项。

4 结论与讨论

本文遴选2012~2015年多年高质量潜热通量和气象数据，采用整体与分段相结合的拟合迭代算法，反向利用P-M模型拟合最小冠层阻力的季节分布，其呈夏低冬高的季节变化特征，变化范围在20.8~90.1 s m⁻¹。

基于Jarvis冠层阻力方案中的最小冠层阻力年均分布曲线以规避其设为定值导致的模拟误差，通过比较2016年模拟与实测的相关系数、均方根误差和一致性指数，表明方法的合理性并能改进了潜热通量估计。并应用于宁乡通量观测站潜热通量缺失数据的插补，比较其与非线性插值法、日变化趋势插补法和查表法等统计方法插补效果，结果表明：就单点和连续缺失数据的插补效果而言，本方法表现皆优于统计方法，且连续缺失插补结果的稳定性不随缺失数据个数的增加而降低。

利用微分误差分析方法量化输入数据对模拟和插补结果的不确定性，误差估计表明：当模型各输入分量相对误差相同时，净辐射、冠层阻力、储热项、空气动力学阻力和土壤热通量对插补的潜热误差贡献依次降低。

本文方法在模型应用和通量估计这两方面都具有扩展性。模型应用方面，考虑了更多的潜热贡献源，物理过程也更完善的双层模型（Shuttleworth and Wallace, 1985）和多源模型（Mu et al., 2011），在采用Javis阻力参数化方案过程中仍涉及到最小冠层阻力为定值的问题，可借助本文方法改进其模拟性能；通量估计方面，基于潜热通量的插补结果及波文比关系本方法可应用于显热通量估计（Irmak et al., 2014），但也还需关注的是，受潜热插补和波文比观测双重不确定性影响，显热通量估计效果会弱于潜热。此外，本研究还面临着一些不足，一方面，本文最小冠层阻力分布曲线是基于单个通量站点拟合得到，受观测时长限制及下垫面代表性的限制，该曲线在其他站点应用还存在不确定性；另一方面，本文改进的阻力方案是基于统计分布而未考虑植被生理过程，未来可在最小冠层阻力分布的过程表述中进行更深入研究。