大探实验_涡度相关观测系统计算程序

2019级 505组 源代码

最终的计算结果部分由python程序给出，另一部分由matlab计算得到
https://www.yuque.com/books/share/a0b1ee7d-cb7b-4b75-9bee-f8782e7b9ba8?

python程序如下

1. 实验工具：
   使用Python语言完成，主要应用numpy库和math库完成计算和存储
2. 实验代码：

import numpy as np
import math
import pandas as pd
import os

path  = r'D:\学习\大二\大探\涡度\05\25m'
llist = os.listdir(path)
filelist = []

for text in llist:
    if text.endswith('.dat') == True:
        filelist.append(text)

def wind_data_clean(start, end, n1, n2, n3):
    U, V, W = [], [], []
    u0 = np.array(df[n1][start:end])
    v0 = np.array(df[n2][start:end])
    w0 = np.array(df[n3][start:end])

    l1 = np.nanmean(u0) + 3 * np.nanstd(u0)
    r1 = np.nanmean(u0) - 3 * np.nanstd(u0)
    U = df[n1][start:end]
    
    l2 = np.nanmean(v0) + 3 * np.nanstd(v0)
    r2 = np.nanmean(v0) - 3 * np.nanstd(v0)
    V = df[n2][start:end]
    
    l3 = np.nanmean(w0) + 3 * np.nanstd(w0)
    r3 = np.nanmean(w0) - 3 * np.nanstd(w0)
    W = df[n3][start:end]
    
    for i in range(start , end):
        if U[i] > l1 or U[i] < r1:
            U[i] = np.nan
            V[i] = np.nan
            W[i] = np.nan
    
    for i in range(start , end):
        if V[i] > l2 or V[i] < r2:
            U[i] = np.nan
            V[i] = np.nan
            W[i] = np.nan
    
    for i in range(start , end):
        if W[i] > l3 or W[i] < r3:
            U[i] = np.nan
            V[i] = np.nan
            W[i] = np.nan
    
    U = np.array(U)
    V = np.array(V)
    W = np.array(W)
    
    return U, V, W

def change_NAN(start, end, name):
    for i in range(start, end):
        if df[name][i] == 999:
            df[name][i] = np.nan
    return 0

def data_clean(start, end, name):
    u0 = np.array(df[name][start:end])
    l = np.nanmean(u0) + 3 * np.nanstd(u0)
    r = np.nanmean(u0) - 3 * np.nanstd(u0)
    
    for i in range(start , end):
        if df[name][i] > l or df[name][i] < r:
            df[name][i] = np.nan
            
    U = 0
    return U

def average_data(start, end, data_):
    average_data_ = np.nanmean(data_[start:end])
    return average_data_

def change_xyz(start, end, Ux_, Uy_, Uz_):
    u0 = Ux_[start:end]
    v0 = Uy_[start:end]
    w0 = Uz_[start:end]
    
    cos_alpha = np.nanmean(u0)/math.sqrt(np.nanmean(u0) ** 2 + np.nanmean(v0) ** 2)
    sin_alpha = np.nanmean(v0)/math.sqrt(np.nanmean(u0) ** 2 + np.nanmean(v0) ** 2)
    
    u1 = u0*cos_alpha + v0 * sin_alpha
    v1 = -u0*sin_alpha + v0 * cos_alpha
    w1 = w0
    
    cos_beta = np.nanmean(u1)/math.sqrt(np.nanmean(u1) ** 2 + np.nanmean(w1) ** 2)
    sin_beta = np.nanmean(w1)/math.sqrt(np.nanmean(u1) ** 2 + np.nanmean(w1) ** 2)
    
    u2 = u1*cos_beta + w1 * sin_beta
    v2 = v1
    w2 = -u1*sin_beta + w1 * cos_beta
    
    return u2, v2, w2

def ganre(start, end):
    p0 = df['amb_press'][start:end]*1000
    t0 = df['Ts'][start:end]+273.15
    H2O = df['H2O'][start:end]/1000
    rou0 = p0/(287.059*t0) + H2O
    cp = 1004.67*(1+0.84*H2O/rou0)
    t0 = np.array(t0)
    Uz_ = np.array(Uz)

    WT = np.nanmean(np.multiply((t0 - np.nanmean(t0)),(Uz_ - np.nanmean(Uz_))))
    
    
    H = rou0*cp*WT
    H = np.array(H)
    return np.mean(H)

def qianre(start, end):
    t0 = df['Ts'][start:end]+273.15
    lamda = (2.501-0.00237*t0)*(10**6)
    H2O = df['H2O'][start:end]/1000
    Uz_ = np.array(Uz)
    H2O_ = np.array(H2O)
    WR = np.nanmean((H2O_ - np.nanmean(H2O_))*(Uz_ - np.nanmean(Uz_)))
    
    return np.mean(lamda*WR)

def wind_direction(start, end):
    u0 = Ux
    v0 = Uy
    cos_alpha = np.nanmean(u0)/math.sqrt(np.nanmean(u0) ** 2 + np.nanmean(v0) ** 2)
    direction = math.acos(cos_alpha)*180/math.pi
    return direction

def CO(start, end):
    CO2 = df['CO2'][start:end]*0.000001
    Uz_ = np.array(Uz)
    CO2_ = np.array(CO2)
    WC = np.nanmean((CO2_ - np.nanmean(CO2_))*(Uz_ - np.nanmean(Uz_)))
    
    return np.mean(WC)

data = []
H_data_item = []
E_data_item = []
direction_item = []
C_data_item = []
for x in range(0,len(filelist)):

    df = pd.read_csv(filelist[x],
                    header=None,
                    names=['TIMESTAMP','RECORD','Ux','Uy','Uz','Ts','diag_sonic','CO2','H2O','diag_irga','amb_tmpr','amb_press','CO2_sig_strgth','H2O_sig_strgth'],
                    skiprows=lambda x: (x%36000)//600<30,
                    engine='python')
    
    aUx, aUy, aUz = [], [], []
    H_data = []
    E_data = []
    wind_d = []
    C_data = []
    for i in range(0, len(df["Ux"]), 18000):
        a = change_NAN(i, i+18000, 'Ux')
        a = change_NAN(i, i+18000, 'Uy')
        a = change_NAN(i, i+18000, 'Uz')
        a = change_NAN(i, i+18000, 'Ts')
        a = change_NAN(i, i+18000, 'H2O')
        a = change_NAN(i, i+18000, 'amb_press')
        a = data_clean(i, i+18000, 'CO2')
        Ux, Uy, Uz = wind_data_clean(i, i+18000, 'Ux', 'Uy', 'Uz')
        wind_d.append(wind_direction(i, i+18000))
        Ux, Uy, Uz = change_xyz(0, 18000, Ux, Uy, Uz)
        aUx.append(average_data(0, 18000, Ux))  
        aUy.append(average_data(0, 18000, Uy)) 
        aUz.append(average_data(0, 18000, Uz)) 
        H_data.append(ganre(i, i+18000))
        E_data.append(qianre(i, i+18000))
        C_data.append(CO(i, i+18000))
    
    H_data_item.append(H_data)
    E_data_item.append(E_data)
    direction_item.append(wind_d)
    C_data_item.append(C_data)
    
data.append(H_data_item)
data.append(E_data_item)
data.append(direction_item)
data.append(C_data_item)
np.save('25m_10_13-15.npy', data)

matlab相关程序

1.DealData.m

clc;clear;close all;

Len=18000;   % 每30min数据有18000行
     data0temp=importdata("D:\exp\exp3\05\data\25m1013.dat") ;   % 加载数据  importdata
     data0=data0temp.data;      % data0为存放原始数据的矩阵
   

for ii=1:48    %1day - 24hours - 48 halfhours
    
      Rawdata0(:,1:4)=data0(1+(ii-1)*Len:Len*ii,2:5);   
      %  2，3，4列为风速在3个方向上的分量，第5列为气温数据   在Rawdata数据中，1-3列存风速数据
      
      Rawdata0(:,5:6)=data0(1+(ii-1)*Len:Len*ii,7:8);
      %  第7列数据为CO2密度，第8列为水汽密度D
      
       
%进行对数据的质量控制,依旧3倍标准差原则筛选*******************************************************
      clear ave; clear std_data;
      for j=1:6
          ave(j) = mean(Rawdata0(:,j),'omitnan');
          std_data(j) =  std(Rawdata0(:,j),'omitnan');
      end
      ave = ave';  std_data = std_data';   
      
      for j=1:6
          for i= 1:18000
             if Rawdata0(i,j)>ave(j)+3*std_data(j)|Rawdata0(i,j)<ave(j)-3*std_data(j)
                    Rawdata0(i,j) = nan;
            end
          end
      end
      clear Rawdata_1
      Logi2=isnan(Rawdata0(:,1))|isnan(Rawdata0(:,2))|isnan(Rawdata0(:,3))|isnan(Rawdata0(:,4))|isnan(Rawdata0(:,5))|isnan(Rawdata0(:,6));  
       
    for j=1:6
           Rawdata_1(:,j)=Rawdata0(~Logi2,j);   % 剔除数据中含有NaN的行（该行对应的布尔值为0）
    end 
    
    density_of_air = 1.293;
%进行坐标旋转*******************************************************************************************
       data_1=rotate1(Rawdata_1);
       datamean0_1(ii,:)=mean(Rawdata_1);   
       % % mean函数作用于矩阵求的就是各列的平均值
   
       
       datamean_1(ii,:)=mean(data_1);
       % 计算旋转后平均值，以datamean_1的结果为准！！！
       
%计算协方差矩阵*******************************************************************************
       datacov_1=zeros(2,0);  
       
       datacov_1=cov(data_1); % data_1为坐标旋转后的最终矩阵
       %datacov 为协方差方针6*6
       
       flux_1(ii,:)=datacov_1(3,:);  % 垂直方向上的通量
       
       sigma_1(ii,1)=datacov_1(1,1)^0.5;    sigma_1(ii,2)=datacov_1(2,2)^0.5;    
       sigma_1(ii,3)=datacov_1(3,3)^0.5;    sigma_1(ii,4)=datacov_1(4,4)^0.5;
       
       % sigma_1存储标准差数据
     
       
%计算湍流强度，摩擦速度，动量通量**********************************************************
       U(ii)=sqrt( datamean_1(ii,1)^2+datamean_1(ii,2)^2+datamean_1(ii,3)^2 ); 
       %calculate average wind speed
       
       for j=1:3
           tur_inten(ii,j) = sigma_1(ii,j)/U(ii);
       end
       
       U = U';
       
       %计算磨擦速度 用fric_speed表示
       fric_speed(ii) =  (datacov_1(1,3)^2+datacov_1(2,3)^2)^0.25; 
       
       % 计算摩擦速度，也可直接利用垂直方向通量数据：fric_speed(ii) = (flux_1(ii,1)^2+flux_1(ii,2)^2)^0.25;
       
       
       %计算切应力（动量通量），用mom_flux表示
       mom_flux = -density_of_air*fric_speed;
       
      
end


fric_speed = fric_speed';
mom_flux = mom_flux';
mean_wind = datamean_1(2:2:48,1:3);

for i = 1:24
    frictionspeed(i) = fric_speed(2*i); 
    tur(i,:) = tur_inten(2*i,:);
    windstd(i,1:3) = sigma_1(2*i,1:3);
    momen(i) = mom_flux (2*i);
    
    
end
frictionspeed = frictionspeed';

  momen = momen';
  flux_1=flux_1(2:2:48,:);

2.rotate1.m

function data=rotate1(Rawdata)
%该程序的目的是做坐标变换，TR
%input:Rawdata
%output:data
% 

% 只通过旋转改变传入参数矩阵前3列的风速数据，不影响第4-6列的其他数据

datamean=mean(Rawdata);
U=datamean(1);V=datamean(2);W=datamean(3);


%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%旋转，使得U沿着平均速度方向，而V＝0，
YawMatrix=[U/(sqrt(U^2+V^2)),-V/(sqrt(U^2+V^2)),0;V/(sqrt(U^2+V^2)),U/(sqrt(U^2+V^2)),0;0,0,1];
YawData=Rawdata(:,1:3)*YawMatrix;
datamean=mean(YawData);
U=datamean(1);V=datamean(2);W=datamean(3);
%通过乘变换方阵实现旋转

%旋转，使得U沿着平均速度方向，而W＝0，
PitchMatrix=[U/(sqrt(U^2+W^2)),0,-W/(sqrt(U^2+W^2)); 0,1,0;W/(sqrt(U^2+W^2)),0,U/(sqrt(U^2+W^2))];
PitchData=YawData*PitchMatrix;
datamean=zeros(2,0);datamean=mean(PitchData);
datacov=zeros(2,0);datacov=cov(PitchData);

data=Rawdata;
data(:,1:3)=PitchData;    

3.程序A（基于向量的数据处理程序）

DealDataFormer.m

    i=1;
     for(j=18001:36000:846001)
        Ux_raw_ave(i) = mean(Ux_raw(j:j+17999),'omitnan');
        Ux_raw_std(i) = std(Ux_raw(j:j+17999),'omitnan');
        
        Uy_raw_ave(i) = mean(Uy_raw(j:j+17999),'omitnan');
        Uy_raw_std(i) = std(Uy_raw(j:j+17999),'omitnan');
        
        Uz_raw_ave(i) = mean(Uz_raw(j:j+17999),'omitnan');
        Uz_raw_std(i) = std(Uz_raw(j:j+17999),'omitnan');
        
        i=i+1;
     end
     
        Ux = Ux_raw;
        Uy = Uy_raw;
        Uz = Uz_raw;

% 后续计算只需要每小时后30分钟的数据，故只对这30分钟数据进行筛选和修正，因此循环不遍历全部数据（如下）       
for i=1:24
    for j=36000*i-17999:36000*i
        if (Ux_raw(j)>Ux_raw_ave(i)+3*Ux_raw_std(i) || Ux_raw(j)<Ux_raw_ave(i)-3*Ux_raw_std(i))
            Ux(j)=NaN;
            Uy(j)=nan; Uz(j)=nan;
        end 
        
        if (Uy_raw(j)>Uy_raw_ave(i)+3*Uy_raw_std(i) || Uy_raw(j)<Uy_raw_ave(i)-3*Uy_raw_std(i))
            Uy(j)=NaN;
            Ux(j) = nan;  Uz(j) = nan;
        end   
        
        if (Uz_raw(j)>Uz_raw_ave(i)+3*Uz_raw_std(i) || Uz_raw(j)<Uz_raw_ave(i)-3*Uz_raw_std(i))
               Uz(j)=NaN;
               Ux(j)=nan; Uy(j)=nan;
        end 
    end
end

for i=1:24
    for j=36000*i-17999:36000*i
        if isnan(Ux(j))  | isnan(Uy(j))  |isnan(Uz(j)) 
            Ux(j) = nan;
            Uy(j) = nan; 
            Uz(j) = nan;
        end
    end
end

%重新计算平均值，生成向量
i = 1;
for(j=18001:36000:846001)
        Ux_ave(i) = mean(Ux(j:j+17999),'omitnan');
        Ux_std(i) = std(Ux(j:j+17999),'omitnan');
        
        Uy_ave(i) = mean(Uy(j:j+17999),'omitnan');
        Uy_std(i) = std(Uy(j:j+17999),'omitnan');
        
        Uz_ave(i)= mean(Uz(j:j+17999),'omitnan');
        Uz_std(i) = std(Uz(j:j+17999),'omitnan');
        
        i=i+1;
end

% 进行第1次坐标旋转
% 
Ux1 = Ux; Uy1 = Uy;
i = 1;
while(i<=24)
        alpha(i)= atan(Uy_ave(i)/Ux_ave(i));
        Ux1(i*36000-17999:i*36000) = Ux(i*36000-17999:i*36000) * cos(alpha(i))+ Uy(i*36000-17999:i*36000)*sin(alpha(i));
        Uy1(i*36000-17999:i*36000) =-1*Ux(i*36000-17999:i*36000) * sin(alpha(i))+ Uy(i*36000-17999:i*36000)*cos(alpha(i));

        i = i+1;
end



%
% 进行第2次坐标旋转
%
% 先算第一次坐标旋转后Ux的平均值（24个时间段，每段半小时）

i=1;
while(i<=24)
        Ux1_ave(i) = mean(Ux1(i*36000-17999:i*36000),'omitnan');
        i=i+1;
end


Ux2 = Ux1;   Uz2 = Uz;
i =1 ;
while(i<=24)

    beta(i) = atan(Uz_ave(i)/Ux1_ave(i));
    
    Ux2(i*36000-17999:i*36000) = Ux1(i*36000-17999:i*36000)*cos(beta(i))+Uz(i*36000-17999:i*36000)*cos(beta(i));
    Uz2(i*36000-17999:i*36000) = -Ux1(i*36000-17999:i*36000)*sin(beta(i))+Uz(i*36000-17999:i*36000)*cos(beta(i));
    
    i=i+1;
    
end


% 计算坐标旋转后风速的平均值，作为最终结果
i = 1;
while(i<=24)
    Ux_ave_final(i) = mean(Ux2(i*36000-17999:i*36000),'omitnan');
    Uy_ave_final(i) = mean(Uy1(i*36000-17999:i*36000),'omitnan');
    Uz_ave_final(i) = mean(Uz2(i*36000-17999:i*36000),'omitnan');
    
    Ux_std_final(i) = std(Ux2(i*36000-17999:i*36000),'omitnan');
    Uy_std_final(i) = std(Uy1(i*36000-17999:i*36000),'omitnan');
    Uz_std_final(i) = std(Uz2(i*36000-17999:i*36000),'omitnan');

    i=i+1;
end

i = 1;
while(i<=24)
    uv(i) = mean((Ux2(i*36000-17999:i*36000)-Ux_ave_final(i)).*(Uz2(i*36000-17999:i*36000)-Uz_ave_final(i)),'omitnan');
    
    i=i+1;
end




U_average = [Ux_ave_final;Uy_ave_final;Uz_ave_final]';
U_std = [Ux_std_final;Uy_std_final;Uz_std_final]';

%计算风向
i = 1;
while(i<=24)
    judge = Uy_ave_final(i)/Ux_ave_final(i);
    if (Uy_ave_final(i)>0)
        D(i) = 270 - atand(judge);
    end
    if (Uy_ave_final(i)<0)
        D(i) = 90 - atand(judge);
    end
    if ((Ux_ave_final(i)>0) && (Uy_ave_final(i)==0))
        D(i) = 180;
    end
    if ((Ux_ave_final(i)<0) && (Uy_ave_final(i)==0))
        D(i) = 0;
    end
    i=i+1;
end