平均、方差与相关，积分权重导致的计算陷阱：I.空间平均

斥鷃

上一次发帖都已经过去一年了，最近人这个懒啊~另外已向@风情犹在夸下海口，于是趁着月黑风高，我又偷偷摸出了键盘，开始敲新帖啦，这次这帖含金量很高哦，乃是本人经九九八十一难，七七四十九炼，五五二十五夜……哎呀，谁丢的鸡蛋！好了，不背乘法表了，开始正文。

正题：

话说大气科学中常用计算存在几个陷阱，常用计算包括哪些？例如求平均，求差分，求积分以及由这些运算衍生出的运算，比方说求相关系数、相对复杂的EOF。其实这里的陷阱主要指的是权重问题，后续再作进一步阐述，这里先给出这两个常见的（PS：所有以上陷阱问题只针对于用Fortran，MATLAB等软件底层编程，如用NCL、GRADS等内置求差分及平均函数的软件，内置函数则已经很大程度上考虑下述问题）：

1.      水平空间平均——带面积权重的平均，应用：平均温度、NINO指数、EOF等

2.      垂直积分——带质量权重的积分，应用：整层XXX，上层热含量，可降水量等

这里提供一个自编的计算空间相关的函数，大家懒得看内容的可以研究一下程序，相信也可以知道很多东西。

1. >> help scorr
   
2. 计算两序列空间相关的函数,函数形式
   
3.       r=scorr(var1,var2,lat(option))
   
4. var1,var2为两个维数相同的变量,按纬度X经度的方式排列，
   
5. lat为纬度序列，与var1的纬度维相同长度的向量，可选，不填时计算按格点的空间相关

复制代码

scorr.m (831 Bytes, 下载次数: 57)

2016-2-17 21:35 上传

点击文件名下载附件

这一帖先来简要说明水平平均的陷阱：

一、水平空间平均

先给个极端的例子：老师出题，要求80°N以北的温度区域平均，数据给出的是经纬网数据，

	0°	90°E	180°	90°E
90°N	-30	-30	-30	-30
85°N	-20	-20	-20	-20
80°N	-10	-10	-10	-10

好了，这时候小明出现了，直接拿12个数加起来再除12，然后说算得的-20是80°N以北的平均温度，高高兴兴的把结果拿给老师，紧接着老师就给他看了一下这几个点的在地球上的分布图：

显而易见，小明的计算是非常有问题的，因为90°N是一个点，但是由于经纬网取点的原因，被加了4次，而80°N纬圈所包含的点则是前两者的好几倍，但是也只加了4次。问题的实质在于，我们认为地球上每一块相同大小的面积是等权重的，因此面积是求区域平均的单元，而例子中给我们的经纬网并不是等面积的，所以各个数字所代表的区域实际上是不等权的，不能直接求平均，换句话说，如果我们真在这片区域行走，更大可能性出现的情况是在80°N区域温度所影响的区域，所以80°N处温度在求平均时候应当给以更大的权重。

现在来给出一个相对合理的计算结果，因为给出的例子比较特殊，这里按纬度采用线性插值，然后对区域求平均，公式如下：

其中T为温度，只跟与90N中心点的距离r有关，s是面积元，上式可以用极坐标求解，得到-16.67，与小明的-20是有差异的。

因为例子比较特殊，具体计算就不再深究了，如果没弄清楚也没关系，大家只需理解这种由于经纬网存在使得区域平均产生差异的原因就行了。同样的原理，实际上，尤其在处理大尺度大气海洋问题中，区域平均都是不宜用点求和取平均的方式来进行的，原则上都应先从物理量的空间分布情况着手，再通过加权求平均，下面再深入介绍。

二、求平均的一般原理及在水平平均上的应用

在我们处理的大尺度问题中，我们涉及的物理量一般都是连续的物理量，只不过我们都是用离散的方式处理，离散量的求和与连续量的积分相联系，所以本质上我们所做的平均及方差等运算都逃不过积分，对于球面的区域平均问题，平均实质是带面积权重的平均，所以要在这些问题中，如果用的坐标系是经纬网，面积权重就是一个我们必须要注意的问题。

由两种方式可以得到球面在经纬网坐标下的面积权重：

第一个是图形方法，在吕美仲版动力气象P25页开始有具体介绍，用图形解释x,y,z三个坐标元与经度、纬度和高度关系，这里就放一张美图（不美也不要说出来噻……），不赘述了。

第二个是数学方法，大家应该还记得高数中提到过球坐标系和直角坐标系的面积元有对应关系

，

这里简单点说就是把球坐标系的积分元中r换成地球半径a即可，而后积分变成二维积分。两种方式得到的面积权重结果都是a^2*cosψ，最后我们可以得到（前方高能蓝色预警）经纬网下的积分表达式：

这里用<·>表示对某个东西求积分，其中f是求积分的物理量，a是地球半径，λ和ψ是经度和纬度，上式便是通用的空间求积分的式子。于是乎，我们可以得到空间平均和空间方差的表达式，（<1>表示直接对面积元素积分,<>沿用上面的定义，代入即可）：

注意第二式中的f的平均值是第一式中算出来的哦~ 可以证明，用同样定义的积分方式做的协方差和相关系数满足自身故有的性质，比如相关系数的绝对值小于1等。

PS：这一部分内容可以从内积、范数等泛函的角度更进一步认识，这里不再啰嗦，有兴趣的同学可以进一步思考。

三、积分问题的延伸，EOF

我关注这个问题，是来源于NCL官网，example中EOFs里的eof1.ncl。其中有一段代码提到了权重，当时我就很纳闷这个权重是干什么用的，后来想到EOF计算过程中要用到平均值和方差我才恍然大悟，这个权重是为了在计算过程中补足方差计算中的面积权重才放进去的。因为直接的奇异值分解是没有带权重的，所以这里必须要补充这个权重，防止面积权重对奇异值分解的影响。不过程序后半段EOF还原部分除掉的面积权重我暂时没有领会到它的含义，大家可以进一步交流。

由于涉及的内容实在太多，而且打公式实在太烦（求体力！），所以没有更进一步深入，有兴趣的同学有时间我们可以再进一步交流。总结上面的内容，大家以后要是做大尺度的，尤其涉及中高纬问题时，平均或者方差计算要注意一下权重问题哦。

曦飛

不错，mark一下

river

这绝对是大神啊，支持一个！

斥鷃

river 发表于 2016-2-17 21:56
这绝对是大神啊，支持一个！

伦家会害羞的啦

kongfeng0824

太厉害了   

aprden

这是处理数据中的一个关键问题。

风情犹在

嗯 插值的时候同样有这个问题

wxl_2015

专家级别，精英人才

开♂o心

求平均、插值都存在这个问题

KIMO23

  太給力了~~~感謝分享