【RS Notes】Moran et al., 1994 利用地表气温和植被指数的关系估算作物水分亏缺

现有研究的问题

Idso et al.和Jackson et al.于1981年提出的作物水分胁迫指数（Crop Water Stress Index, CWSI）利用叶片温度和气温的差值来估算作物水分胁迫，是一个常用的作物水分胁迫监测指数，但在实际应用中，对于部分植被覆盖的地表，土壤的贡献难以去除，传感器很难直接测量得到叶片温度。

解决办法：植被指数/温度（VIT）梯形（Trapezoid）

图1 梯形特征空间

作出 $T_\text{s}-T_\text{a}$ 和植被覆盖度的散点图，散点的外轮廓将会呈现出类似图1中梯形的形状。

基本假设

$T_\text{s}-T_\text{a}=\mathrm{FVC}\cdot(T_\text{v}-T_\text{a})+(1-\mathrm{FVC})\cdot(T_\text{o}-T_\text{a})$

也即温度的线性叠加。实际上根据Stefan-Boltzman公式，应该是温度的四次方（能量）满足线性叠加关系。这里直接认为温度可以线性叠加，是进行了简化。

$\begin{aligned} & T_\text{v}-T_\text{a} = a +b\cdot\Gamma \\ &T_\text{o}-T_\text{a} = a' +b'\cdot E \end{aligned}$

也即在给定 $R_\text{n}$ 和 $\mathrm{VPD}$ 的条件下， $T_\text{v}-T_\text{a}$ 和 $T_\text{o}-T_\text{a}$ 与蒸发 $E$ 和蒸腾 $\Gamma$ 之间是线性关系。

水分亏缺指数（Water Deficit Index, WDI）

$\begin{aligned} \mathrm{WDI} &=1-\lambda E_{\Gamma} / \lambda E_{\Gamma_p} \\ &=\left[\left(T_\text{s}-T_\text{a}\right)_\min-\left(T_\text{s}-T_\text{a}\right)_{r}\right] /\left[\left(T_\text{s}-T_\text{a}\right)_\min-\left(T_\text{s}-T_\text{a}\right)_\max\right] \end{aligned}$

可以看到WDI的定义和CWSI基本是一致的，区别在于，CWSI中考虑的是蒸腾和潜在蒸腾；而WDI中考虑的是蒸散和潜在蒸散，因为WDI针对的是部分覆盖的情形。

考虑到土壤的蒸发和植被的蒸腾并不是独立进行的，而是有着相互的影响，提出了多成分模型。

$\begin{array}{c} H=c_\text{a,v}\left\{\left(T_\text{v}-T_\text{a}\right) / r_\text{v}+\left(T_\text{o}-T_\text{a}\right) / r_\text{o}\right\} / \\ {\left[1+r_\text{a} / r_\text{v}+r_\text{a} / r_\text{o}\right]} \end{array}$

这里，认为部分覆盖情形下植被冠层阻抗和土壤阻抗分别与完全覆盖时的植被冠层阻抗和裸土阻抗成线性关系。

接下来，考虑遥感的可用性，选择一个认为与 $\mathrm{FVC}$ 具有线性关系的植被指数，代替 $\mathrm{FVC}$ 进行特征空间的构建。本文选择的是土壤调节植被指数（Soil-Adjusted Vegetation Index, SAVI）。从而，由梯形的两腰就可以得到最干和最湿两种情况下温度差 $T_\text{s}-T_\text{a}$ 与SAVI之间的线性关系，进而可以计算SDI。

实际使用中的假设

净辐射 $R_\text{n}$ 可视为定值。
粗糙度长度 $z_\text{o}$ 和位移高度 $d_\text{o}$ 与元素/作物高度满足线性的经验关系。
土壤热流 $G$ 可由 $R_\text{n}$ 和 $\mathrm{FVC}$ （或SAVI）确定。
$kB^{-1}$ 与风速和温度差 $T_\text{s}-T_\text{a}$ 的乘积满足线性关系。

输入参数

最大植被高度和最小土壤粗糙度
完全覆盖和裸土情形下的最大和最小SAVI取值
最大LAI
气孔阻抗的最大和最小值

从梯形到三角形

图3 从梯形到三角形

如图所示，利用梯形的四个顶点，构成四个三角形分别考虑。

接下来是本文的核心内容。

图4 梯形中的“沙漏”

如图4所示，任意进行一次测量，可以根据这一测量点落在梯形中的哪一个子区域，决定将其与梯形的哪两个顶点相连，从而确定一个“沙漏”范围，继而确定该测量点所对应的胁迫范围。图4d说明在进行一次测量后，结合 $T_\text{o}-T_\text{a}$ 就可以得到 $T_\text{v}-T_\text{a}$ 。

误差源分析

云状况的快速变化将导致严重的误差
阴影植被和光照植被，以及阴影土壤和光照土壤之间都存在差异
气象数据的误差也是重要的误差来源
未考虑土壤和植被冠层之间的耦合效应

验证

模拟

使用一个两组分的模型进行了模拟，发现 $H$ 和 $\lambda E_\Gamma$ 与 $\mathrm{FVC}$ 之间呈现出明显的非线性关系。同样，WDI与 $\lambda E_\Gamma / \lambda E_{\Gamma_p}$ 之间也存在明显的非线性关系。

图6 模拟结果（一）

图7 模拟结果（二）

实验

图8 实验测点散点图，实现为理论梯形

图9 四种灌溉条件的WDI

图10 WDI与CWSI和基线CWSI的比较

回顾：CWSI原理

在附录中，作者回顾了CWSI的原理和应用，写得也很好，值得学习。

基本原理

$R_\mathrm{n}=G+H+\lambda E_\Gamma$

对于稠密作物，有如下的简易表达式：

$H=c_\text{a,v}(T_\text{v}-T_\text{a})/r_\text{a}$

其中 $c_\text{a,v}$ 为空气的体积热容（ $\mathrm{J\cdot^\circ\hspace{-0.3em}C^{-1}\cdot m^{-3}}$ ）， $r_\mathrm{a}$ 为空气动力学阻抗（ $\mathrm{s\cdot m^{-1}}$ ）。

$\lambda E_\Gamma=c_\text{a,v}\cdot\mathrm{VPD}/[\gamma(r_\mathrm{a}+r_\mathrm{v})]$

其中 $\mathrm{VPD}$ 为蒸气压差（ $\mathrm{kPa}$ ）， $\gamma$ 为湿度线常数（Psychrometric Constant， $\mathrm{kPa\cdot^\circ\hspace{-0.3em}C^{-1}}$ ）， $r_\mathrm{v}$ 为冠层（对水汽输运的）阻抗（ $\mathrm{s\cdot m^{-1}}$ ）。

对于植被完全覆盖的情形， $G$ 可以忽略不计。从而可以得到：

$\begin{aligned} \left(T_\text{v}-T_\text{a}\right)=&\left[r_\text{a} R_\text{n} / c_\text{a,v}\right]\left[\gamma\left(1+r_\text{v} / r_\text{a}\right) /\right.\\ &\left.\left\{\Delta+\gamma\left(1+r_\text{v} / r_\text{a}\right)\right\}\right] \\ &-\left[\operatorname{VPD} /\left\{\Delta+\gamma\left(1+r_\text{v} / r_\text{a}\right)\right\}\right] \end{aligned}$

其中 $\Delta$ 为饱和蒸气压对温度的斜率（ $\mathrm{kPa\cdot^\circ\hspace{-0.3em}C^{-1}}$ ）。

CWSI的建立

CWSI定义为：

$\begin{aligned} \mathrm{CWSI} &=1-\lambda E_{\Gamma} / \lambda E_{\Gamma\mathrm{p}} \\ &=\left[\gamma\left(1+r_\text{v} / r_\text{a}\right)-\gamma^{*}\right] /\left[\Delta+\gamma\left(1+r_\text{v} / r_\text{a}\right)\right] \end{aligned}$

其中 $\gamma^*=\gamma(1+r_\text{vp}/r_\text{a})$ ， $r_\text{vp}$ 为潜在冠层阻抗。

CWSI的求解

第一种办法是求解 $r_\text{v} / r_\text{a}$ ：

$\begin{aligned} r_\text{v} / r_\text{a}=&\left\{\left[\gamma r_\text{a} R_\text{n} / c_\text{a,v}\right]-\left[\left(T_\text{v}-T_\text{a}\right)(\Delta+\gamma)\right]\right.\\ &-\mathrm{VPD}\} /\left\{\gamma\left[\left(T_\text{v}-T_\text{a}\right)-r_\text{a} R_\text{n} / c_\text{a,v}\right]\right\} \end{aligned}$

注意这里并没有消去 $r_\text{a}$ 。

第二种办法是：

$\begin{aligned} \mathrm{CWSI} &=1-\lambda E_{\Gamma} / \lambda E_{\Gamma p} \\ &=\left[\left(T_\text{v}-T_\text{a}\right)-\left(T_\text{v}-T_\text{a}\right)_\min\right] /\left[\left(T_\text{v}-T_\text{a}\right)_\max-\left(T_\text{v}-T_\text{a}\right)_\min\right] \end{aligned}$

第二种方法看起来只包含温度了，但是这里的 $\left(T_\text{v}-T_\text{a}\right)_\min$ 和 $\left(T_\text{v}-T_\text{a}\right)_\max$ 并不是可以直接得到的数值，还需要通过以下方法来计算：

$\begin{aligned} \left(T_\text{v}-T_\text{a}\right)_\min=&\left[r_\text{a} R_\text{n} / c_{\mathrm{a,v}}\right]\left[\gamma\left(1+r_\text{vmin} / r_\text{a}\right) /\right.\\ &\left.\left\{\Delta+\gamma\left(1+r_\text{vmin} / r_\text{a}\right)\right\}\right] \\ &-\left[\mathrm{VPD} /\left\{\Delta+\gamma\left(1+r_\text{vmin} / r_\text{a}\right)\right\}\right] \end{aligned}$

$\begin{aligned} \left(T_\text{v}-T_\text{a}\right)_\max=&\left[r_\text{a} R_\text{n} / c_\text{a,v}\right]\left[\gamma\left(1+r_\text{vmax} / r_\text{a}\right) /\right.\\ &\left.\left\{\Delta+\gamma\left(1+r_\text{vmax} / r_\text{a}\right)\right\}\right] \\ &-\left[\mathrm{VPD} /\left\{\Delta+\gamma\left(1+r_\text{vmax} / r_\mathrm{a}\right)\right\}\right] \end{aligned}$

其中 $r_\text{vmin}=r_\text{vp}$ ， $r_\text{vmax}\rightarrow\infty$ 。此外，也可以利用LAI来计算：

$\begin{aligned} r_\text{vmin}=r_\text{smin}/\mathrm{LAI} \\ r_\text{vmax}=r_\text{smax}/\mathrm{LAI} \end{aligned}$

其中 $r_\text{s}$ 为气孔阻抗， $r_\text{smin}$ 的经验值为25--100（ $s\cdot m^{-1}$ ）， $r_\text{smax}$ 的经验值为1000--1500（ $s\cdot m^{-1}$ ）。

CWSI应用于裸土

上面讨论的都是完全植被覆盖的情形。CWSI也可以被推广到裸土，此时 $\lambda E_{\Gamma} / \lambda E_{\Gamma_p}$ 表示土壤实际蒸发与潜在蒸发之间的比值。对应的，要将叶片温度改为土壤表面温度。

对于完全湿润的裸土，阻抗可以视为0；对于完全干燥的裸土，阻抗视为 $\infty$ 。

空气动力学阻抗的估算

$\begin{array}{c} r_\text{a}=\left\{\left[\ln \left(\left(z-d_{o}\right) / z_{o m}\right)+\ln \left(z_{o m} / z_{o h}\right)-\Psi_{h}\right]\right. \\ \left.\quad \times\left[\ln \left(\left(z-d_{o}\right) / z_{o m}\right)-\Psi_{m}\right]\right\} / k^{2} U \end{array}$

此处可以参考FAO关于Penman-Monteith方程的介绍。

CWSI的基线（Baseline）方法

基线指的是土壤水分不是限制因素时， $T_\text{v}-T_\text{a}$ 与 $\mathrm{VPD}$ 的关系。此时植被按照潜在蒸腾速率进行蒸腾作用。

利用基线的斜率 $f_1$ 和截距 $f_0$ ，可以计算CWSI：

$\begin{aligned} \mathrm{CWSI}=&\left\{\left(f_{0}+f_{1} \mathrm{VPD}\right)-\left(T_\text{v}-T_\text{a}\right)_{r}\right\} / \\ &\left\{\left(f_{0}-f_{1} \mathrm{VPD}\right)-\left(f_{0}+f_{1} \mathrm{VPD}_{x}\right)\right\} \end{aligned}$

文章中并未提及植被部分覆盖时如何应用CWSI——这是正文的梯形方法着力解决的问题。

【RS Notes】Moran et al., 1994 利用地表气温和植被指数的关系估算作物水分亏缺