弹性系数系数在水文气象中的应用及其MATLAB实现

       前言，在写这个博客之前，本人经过了很长时间的纠结，看过很多关于弹性系数的文献，也请教过很多人，但一直没有彻底搞明白。所以，只要看到文章中提高“弹性系数”，“因子贡献率”，“全微分”或者“敏感性分析”，心里都发憷，而且因为一直没有彻悟而揪心，简直到了寝食难安的地步。今天经过再次思考，突然茅塞顿开，抓紧整理出来以备同样困惑的朋友学习。废话少数，开始讲解！

1、弹性系数概念及推导

       弹性系数，或者是敏感性，指的是因变量变化率与自变量变化率的比值。是不是不好理解？那我就举个通俗的例子，比如人的幸福感（因变量）由工资收入，职位高低和孩子多少3个因素决定（自变量），如果保持职位，孩子多少不变，仅仅把你的工资提高1%看看你的幸福感提高多少，如果变化很明显，那么说明工资收入对幸福感很敏感；同理，如果说你再多要一个孩子你想死的心都有，那么说明孩子多少也很敏感，只是负影响而已。如果敏感性为0.1，那么就表明自变量增加10%，因变量则会随之增加1%。公式推导如下

其中Sx是指的敏感性指数，y是因变量（如上例子中的幸福指数），x是自变量（如例子中的工资收入）。这个公式中只显示出来了对x自变量的敏感性，当然对其他因子也是类似，各个自变量之间要为独立关系，那么∂y/∂x就是偏微分。

2、气象中干燥指数的弹性系数及MATLAB实现

   在气象中常用到一个指数是干燥指数，指的是潜在蒸散发ET₀（是辐射，温度，水汽压，风速的函数，本文以Penmen-Montieth公式为例）与降水P的得比值（如下）。干燥指数是用来衡量一个地区气候干湿状况的一个重要指标，在气候、水文、植被等领域得到了广泛的应用。

   问题就来了，这么重要的一个指数，对一个特定区域它会受哪些因素影响显著，或者对哪些因素敏感呢？这就涉及到了AI对气象因子的敏感性分析了。

比如分析AI对2米高风速的敏感性，根据敏感系数定义就有：

利用MATLAB的diff求偏微分的函数很容易得到ET0对u2的偏微分公式,以及ET₀对VP（水汽压），Rs（辐射）,Ta（温度）的偏导。如下所示：

 

3、水文中径流的弹性系数以及MATLAB实现

水文研究中常提高一个问题就是分析径流的变化，到底是哪些因素影响的，并且各个因素的贡献率是多少？一般而言，径流变化了包括了气候变化（降水P，潜在蒸散发ET0），植被变化等。以杨大文老师一篇文章为例“基于水热耦合平衡方程的黄河流域径流变化归因分析”。文中数据如下：

首先，对一个流域水量平衡方程为：

但对于一个闭合流域，长时间尺度（比如10年），土壤湿度△ S一般不会变化，可以看做0。故而，多年尺度一个流域的实际蒸散发可以用

     

其次，根据布迪克水热平衡理论Bydyko,一定的气候和植被条件下, 流域长期的水文气候特征服从水分和能量平衡原理。由于Budyko有多公式，本文采用杨大文老师的

其中P为年降雨，ET0为潜在蒸散发，n 是反映流域下垫面特征的参数, 流域下垫面特征包括地形、土壤和植被等。综合以上两个公式，可以得出

这样，就可以根据流域多年降水P，潜在蒸散发ET0，径流R求出反应流域特征的参数n值。利用MATLAB的solve函数解方程，求出参数n。

如上代码求得结果n=1.51，与文中的n值一致，其他也尝试无误。

最后，由于径流R是由降水P，潜在蒸散发ET₀和植被变化n引起，

即有径流R：

所以年径流R对P，ET₀，n的敏感性分别为：

        

利用MATLAB较容易获取偏微分：

（我原先多次尝试用手推出偏微分，想想太滑稽了，有了工具就好了！）

剩下的就是利用MATLAB根据多年的径流R，降水P，潜在蒸散发ET0和植被参数n来计算敏感性了：

注意：弹性系数的精髓就是时间一定要对应，如果P,ET0，R是相对转折点前的，那么对应的弹性系数也是转折点前的。如果是全时段，那么弹性系数也应该是全时段。

综上，可以看出计算结果，与文中给出的弹性系数值一致。据此关于弹性系数或者敏感性分析介绍完了。