计算随机误差

参考文献

Finkelstein, P. L. and P. F. Sims (2001). "Sampling error in eddy correlation flux measurements." Journal of Geophysical Research: Atmospheres 106(D4): 3503-3509.

论文中提到，在观测数据为10Hz的情况下，先取m=200(20s)，后测试发现m取值在100-400之间时，结果仅变化1-2%。如果数据中没有时间趋势，则时间长于400时结果不会改变。

function[var] = cal_random_func2(sub_w,sub_mass,n,lag)

%%%%  lag:  Peter L. Finkelstein  2001 中的P值

Ew = mean(sub_w);
Em = mean(sub_mass);
    
Xw = sub_w - Ew;
Xm = sub_mass - Em;

sw = sum(Xw.*Xw)/n;
sm = sum(Xm.*Xm)/n;
swm = sum(Xm.*Xw)/n;


result = zeros(lag,4)*nan;

for h = 1:lag
    
num = n-h;
sub_result = zeros(num,4);  

    for i = 1:(num)
          sub_result(i,1) = Xw(i)*Xw(i+h);
          sub_result(i,2) = Xm(i)*Xm(i+h);
          sub_result(i,3) = Xm(i)*Xw(i+h);
          sub_result(i,4) = Xw(i)*Xm(i+h);
     end
r = sum(sub_result,1)/n;
result(h,:)= r;

end

 var = [sw*sm+swm*swm+(sum(result(:,1).*result(:,2))*2+sum(result(:,3).*result(:,4))*2)]/n;

end

load("E:\data.mat")
w_cut   = 9000;  %%% 因为w从17:00开始，而mass从17:30开始。
lag     = 9;  %%% 延迟时间 s
var_wbc = zeros(1036,2)*nan;

n = 9000;
p = 200

for i = 1:100
    disp(i);
    
    %%% 以w为准
    index_st = (i-1)*n+1;
    index_ed = i*n;
    
    st_bc = index_st+lag*5;
    ed_bc = index_ed+lag*5;
    
    
    sub_w    = wind(index_st+w_cut:index_ed+w_cut,3);
    sub_mass = BC(st_bc:ed_bc,2);
    sub_num  = BC(st_bc:ed_bc,14);
    
    if (numel(find(isnan(sub_w)))>0) || (numel(find(isnan(sub_mass)))>0)
         
        continue;
        
    end

    var_wbc(i,1) = sqrt(cal_random_func2(sub_w,sub_mass,n,200));
    var_wbc(i,2) = sqrt(cal_random_func2(sub_w,sub_num,n,200));
    
end