考虑地表反照率时域重建的海河流域蒸散发模拟研究

段 浩，韩 昆，赵红莉，徐浩玮，文 铭，靳晓辉，蒋云钟

（1.中国水利水电科学研究院，北京 100048；2.黄河水利科学研究院，河南 郑州 450003）

1 研究背景

地表蒸散发（ET）是联系水、能量和碳循环的中间环节［1］，是陆面生态过程［2］和水资源管理［3］的关键参数。准确获取ET时空分布，对理解陆气间水热交换、加强水资源管理具有重要意义。

基于遥感监测技术，可获取ET在大时空尺度上的变化特征，与传统的站点观测相比具有显著优势。在Penman-Monteith（P-M）公式基础上发展而来的Penman-Monteith-Leuning（PML）模型［4］，是基于遥感数据计算ET的方法之一，该方法通过引入地表导度Gs，可基于植被叶面指数和P-M公式直接计算ET，已在全球不同气候条件下得到了应用。王海波等［5］基于PML模型，模拟了黑河流域ET；李红霞等［6］和Zhang等［7］分别利用此模型对澳大利亚的ET进行了模拟和验证；Li等［8］利用该模型分析了植被变化对ET的影响；Zhang等［9］还基于PML模型分析了全球范围内气象要素对土壤蒸发和植被散发的影响。这些研究使PML模型在ET模拟上逐渐成熟，成为开展大尺度ET监测的重要手段［10］。

地表反照率（ac）是反映地表对太阳辐射反射能力的指标，它的细微变化都会对辐射（Rn）收支分配产生显著影响［11］，最终影响ET的模拟精度。在利用PML模型进行ET模拟时，通常通过两种方式对ac进行处理，一种是直接使用8天合成的ac产品，如Zhang等［12］对MODIS的8天合成ac产品进行了质量控制和平滑处理，在重采样后得到0.05°分辨率的数据；该作者在计算全球ET时［9］同样是使用了Global Land Surface Satellite（GLASS）8天0.05°的ac产品。这种处理方式只能得到每8天的ET结果。另一种方式是对8天合成的ac产品做线性插值，得到逐日数据，如Li等［8］对8天合成的ac进行线性插值，分析了子牙河流域ET对植被变化的响应；Mu等［13］在计算ET时，同样使用了线性插值的方法，得到逐日的ac产品。基于线性插值的方法获取时间连续的ac产品，处理方式简单，但该方法将ac在时间上的演变简化为线性关系，机理性不足。ac的变化受到土壤含水量［14］、下垫面条件［15］、积雪［16］、温度［17］等多种要素的综合影响，具有显著的时空异质性特征，线性插值难以体现其变化过程。张振宇等对SEBAL蒸散发模型的研究表明［18］，通过不同方法得到的ac间的不同，是导致日蒸散量差异的主要原因。

为强化PML模型模拟中ac产品的遥感机理，本文引入充分考虑遥感数据时空信息的谐波分析与泊松方程协同的时域重建算法，在充分考虑ac时空变异特征的基础上对海河流域ac产品进行重建，并基于时域重建后的ac产品进行流域蒸散发反演。同时基于线性插值方法得到的ac，分析不同方法得到的ac对PML模型中ET结果的影响。

2 研究方法

2.1 基于叶面指数的PML蒸散发模型 PML模型计算ET的基本思想是，在潜在蒸发计算方法基础上引入“表面阻抗”的概念［19］，得到非饱和下垫面蒸散发的计算形式，其计算过程可拆分为土壤蒸发和植被蒸腾两部分，具体可表示为［4］：

其中：ET为蒸散发；λ为汽化潜热；ε为温度-饱和水汽压曲线斜率与干湿表常数的比值；ρa为空气的密度；Cp为空气定压比热；Da为参考高度饱和水汽压差；A为可用能量，代表净辐射与土壤热通量的差值；Ga为空气动力学导度；Gs为地表导度；Gc为冠层导度；f为土壤蒸散发系数；Ac和As分别为可用能量中被冠层和土壤吸收的部分。

冠层导度的计算需要利用最大气孔导度gsx和LAI的关系得到：

其中：kQ为短波辐射衰减系数，常取0.6；Qh为冠层上方的可见光辐射通量；Q50和D50为当气孔导度gs=gsx/2（gsx是gs的最大值）时的可见光辐射通量和水汽压差，通常取值分别为2.6 MJ/（m2/d）和0.8 kPa。

可用能量A为净辐射Rn与土壤热通量G之差，而Rn则由ac及太阳短波辐射和净长波辐射得到，计算方法为：

其中Rs为太阳短波辐射（MJ·m-2·d-1），Rnl为净长波辐射（MJ·m-2·d-1）。

在对f和gsx两个参数的处理上，gsx的参数化参考Zhang等的研究，根据土地利用类型确定［9］；在f的参数化方面，本文采用Zhang［12］提出的基于累积降雨与平衡蒸发速率比值的方案来确定f的取值。

2.2 地表反照率时域重建方法 ac数据对ET的模拟有重要影响，为充分考虑ac的时空变异性，增强其处理过程的物理机理，本文引入谐波分析与泊松方程协同的ac时域重建算法［20］，计算得到海河流域逐日的ac产品。时域重建的主要过程如下：

（1）初始值确定。指定目标年（模拟年份）和参考年的ac数据序列，本文选取MODIS的ac产品进行分析，对每一个目标年数据，以其前三年ac数据为参考年数据。首先根据ac的数据质量信息对全部原始数据做质量控制，区分有效区域和数据缺失区域；然后利用矩匹配的方法对参考年数据做归一化，并分析与目标数据的相关性，得到参考年的权值；最后通过多年参考数据的加权处理初步计算缺失区域的ac。

关于参考年数据的归一化处理，又可分为以下三种情况：

任意年份参考数据与目标数据的有效区域均有重叠，此时以目标年有效区域重叠部分的均值和标准差为基准，进行归一化；任意年份参考数据与目标数据的有效区域均无重叠网格，此时无法计算权值，基于多年均值进行数据缺失区域的填补；存在参考年数据与目标数据有效区域都不存在重叠的部分情况，根据缺失区域的质量标识和在上述两种情况下得到的临时像元值来确定。

在上述各情况下初始值的具体计算过程，可参见文献［20］的相关描述。

（2）时域插值。当同一网格在参考年的不同年份均为空值，则通过时域插值来填补步骤（1）未能填补的区域，采用真值谐波分析（HANTS）法进行，其基本算法如下：

3 研究数据

3.1 遥感数据 研究所用到的遥感影像数据包括MODIS的ac数据MCD43A3、LAI数据MCD15A2H，覆盖范围为海河流域，数据时间为2005至2018年（其中2005—2007的ac数据作为2008年ac时域重建的参考年数据），以上数据均可从EARTHDATA（https：//search.earthdata.nasa.gov/）下载，所有影像均通过MRT工具转换为0.01°的分辨率。MCD43A3数据在时间上虽然为逐日序列，但在空间上覆盖不完整，可通过时域重建方式得到时空连续的逐日数据产品。此外，本研究中土地利用数据采用LUCC（Land-Use and Land-Cover Change）地物分类标准。

为检验本研究得到的ET结果，本文还使用了MODIS的蒸散发产品（MOD16）。MOD16数据由美国国家航空航天局（NASA）发布，是目前应用最广泛的遥感ET产品之一［21-22］。

3.2 实测数据 由于ET的站点监测成本较高，实测资料缺乏，本研究采集了大兴站2008年、馆陶站和密云站2010年和望都站2018年的涡度相关数据，作为对PML模型模拟结果精度的验证数据。研究所用气象数据来自国家气象站的地面监测数据（www.data.cma.gov.cn），包括用来计算蒸散发的气温、湿度、日照时数、风速和气压等数据，数据序列为2008、2010和2018年三年的逐日观测资料。将国家气象站点的观测数据基于IDW（Inverse Distance Weighted）方法插值成栅格数据，空间分辨率统一为0.01°。海河流域的主要土地利用分类及气象站、涡度相关站点的位置如图1所示，各涡度相关站点的数据情况如表1所示。

表1 海河流域涡度相关测站情况

图1 海河流域土地利用及气象站、涡度相关站点分布示意图

4 结果及讨论

4.1 海河流域ac重建结果对比 根据时域重建方法和线性插值法可得到的海河流域时空连续的ac产品。图2展示了基于时域重建方法得到的2018年各季节海河流域ac产品。从图中可以看出，2018年ac具有明显的时空变异特征。在空间分布上，2018年ac整体上呈现出北部高，南部低的趋势，同时城市区域的ac相对较低。在季节变化上，夏季和冬季ac的高值区范围相对较小，主要集中在海河流域北部地区；而春季和秋季ac的空间分布相对均衡。在ac的最值方面，冬季的最大值为0.49，同时也是全年的最大值，秋季的最大值较小，为0.31。

图2 时域重建法获取的海河流域2018年各季节ac均值分布图

与时域重建方法相比，基于线性插值方法得到的ac结果整体偏高（见图3）。在2018年各季节中，时域重建的ac结果比线性插值结果偏低0.2以内的地区分别占到88.46%、92.60%、88.31%和82.23%，覆盖了海河流域的大部分地区。在空间分布上，时域重建结果偏高的地区主要集中在海河流域北部及西部山区，其余地区时域重建的ac结果总体偏低。形成这种差异的原因在于两种方法对数据缺失日期插补方法的不同，当某个网格的ac数据出现缺失时，线性插值法基于缺失日期前后的数据建立线性关系，插值结果不会出现明显的低值，当有效数据在时间上的分布较为集中时，还可能会插值出异常高的ac值，樊辉［23］使用不同插值方法对NDVI（Normalized Difference Vegetation Index）数据的重建结果也表明，线性插值的重建结果偏差最大。而基于HANTS的时域重建方法在插补数据时是通过选取均值和显著谐波，利用最小二乘拟合来迭代求解的过程，对缺失数据序列的插补会出现小于缺失日期前后数据的情况。

图3 时域重建法与线性插值法获取的海河流域各季节ac差值分布及统计图

4.2 模型精度评价 基于PML模型对海河流域蒸散发进行反演，将采用时域重建方法得到ac作为输入记为方案“PML_HANTS”，并利用大兴站、馆陶站和密云站的涡度相关数据对模拟结果进行验证（表2）。对模型精度的评价指标选用相关系数（r）、均方根误差（RMSE）和偏离度（Bias）［24］。从图中可以看到，PML_HANTS方案的模拟结果与观测值在三个站点具有较好的相关性，相关系数r分别为0.73、0.77和0.76；在RMSE和Bias指标方面，RMSE指标在2.3 mm/d以下，Bias指标为负值，表明本文的模拟值与观测值相比偏低。与其他基于分布式的f和gsx参数化的模拟结果相比，三个站点的模拟结果的相关系数与Morillas等［25］的研究结果较为一致，表明本文基于PML_HANTS方案在海河流域的ET模拟结果与观测值有较好的相关性。

表2 PML_HANTS方案下大兴、馆陶和密云三个站点ET模拟值与观测值对比

为进一步分析PML_HANTS方案在流域尺度上的精度，选用MOD16的全球蒸散发产品与PML_HANTS方案的模拟结果进行对比。由于MOD16产品是8天合成的ET，在对比时将实测数据与模拟结果统一到与MOD16产品相同的时间尺度上。图4给出了海河流域各涡度相关站点（见表1）8天合成的观测数据与PML_HANTS方案模拟结果及MOD16产品的对比情况，从图中可以看出，PML_HANTS方案模拟结果与涡度相关数据的r是0.68，RMSE是14.23，而MOD16产品与涡度相关数据的r和RMSE分别为0.45和15.09，本文提出的PML_HANTS方案模拟结果要优于MOD16产品。总体来看，PML_HANTS方法在8天尺度上与观测值的RMSE仍有一定偏差，其原因可能有两方面：一是使用基于降雨数据的f参数化方法可能带来不确定性，Morillas等［25］的研究也表明f的参数化方法对ET的模拟精度有一定影响；二是本文采用的气象数据输入是通过对气象站点数据插值后得到，气象数据的插值过程中也会对ET的模拟结果带来一定误差。

图4 海河流域8天尺度涡度相关数据与PML_HANTS方案模拟结果及MOD16产品对比图

4.3 ac对ET的影响分析 为分析本文两种ac结果对ET的影响，还利用线性插值方法得到的ac对海河流域ET进行了模拟，记为方案“PML_LINEAR”，并对比在两种ac输入条件下，PML模型对望都站模拟结果的差异（表3）。从表2可以看出，两种方案的模拟结果与观测值均有较好的相关性，PML_LINEAR方案的r值为0.52，与PML_HANTS方案差异较小。在其他指标方面，PML_HANTS方案的RMSE比PML_LINEAR方案降低了0.18 mm/d，偏离度也减少了7.61%。因此，对望都站而言，PML_HANTS方案的模拟值与涡度相关观测值的准确性更好。

表3 望都站PML_HANTS与PML_LINEAR方案模拟精度对比

在基于PML模型的ET模拟中，ac主要参与净辐射的计算，进而影响到被冠层和土壤吸收的能量Ac和As，并与净辐射的大小成反比。为更加清楚地反映PML_HANTS和PML_LINEAR方案中ac的差异对ET模拟结果的影响，以望都站2018年5月份（有观测值日期）两种方案下ac变化过程与ET变化过程为例进行对比分析（图5）。从图中可以看出，5月份两种方案得到的ac值变化均较为平缓，但PML_HANTS方案下的ac值比PML_LINEAR方案整体偏低，平均值偏低约0.14。在ET的差异方面，两种方案的模拟结果均较好地反映了ET的变化过程，但由于PML_HANTS方案的ac值相对较小，计算得到的净辐射偏大，即冠层和土壤吸收到的辐射更多，PML_HANTS方案得到的ET结果平均偏高0.32 mm/d。由此可以看出，ac的计算方案的不同，影响了望都站ET的模拟精度。

图5 2018年5月不同方案下ac及ET模拟值与观测值对比图

为进一步分析两种ac重建方案对ET结果在海河流域范围内的影响，将PML_HANTS方案的年蒸散发量与PML_LINEAR方案的结果进行了对比，两种方案条件下2018年海河流域ET总量的差值及频率统计如图6所示。从图中可以看出，PML_HANTS方案得到的年ET总量整体高于PML_LINEAR方案，年ET总量与PML_LINEAR方案的差值在0～200 mm的地区占到92.79%，其中海河流域中部和南部地区偏高较为明显，望都站PML_HANTS方案的结果比PML_LINEAR方案的模拟结果偏高，是这种空间差异特征的具体体现。此外，还有7.15%的地区PML_HANTS方法偏低0～200 mm，主要分布在海河流域北部等地区。

图6 海河流域PML_HANTS与PML_LINEAR模拟结果差值分布及统计图

5 结论

本文基于PML模型对海河流域逐日蒸散发进行了反演研究，通过引入时域重建方法，获取了时空连续的地表反照率数据，并与通过线性插值方法得到ac来计算ET的方法进行对比，探讨了ac对ET模拟结果的影响，主要得到以下结论：

（1）通过与不同年份、4个涡度相关站点的监测数据进行对比，本文提出的PML_HANTS方案模拟结果与涡度相关数据的相关性系数在0.53至0.77之间，相关性较好，但比观测值偏低，需在后续工作中持续改进；同时在8天的时间尺度上，PML_HANTS方案模拟结果与涡度相关数据的相关性也优于MOD16产品。

（2）对海河流域2018年ac的时域重建结果显示，ac整体上呈现出北部高，南部低的分布特征，同时夏季和冬季ac的高值区范围相对较小，而春季和秋季ac的空间分布相对均衡。

（3）对望都站2018年5月的分析表明，基于时域重建的ac结果整体上低于基于线性插值得到的ac数据，蒸散发模拟结果的RMSE降低了0.18 mm/d，偏离度降低了7%以上，可为进一步改善蒸散发精度提供思路。