涡度

大气探测实验报告涡度相关观测系统[作者姓名]2016-5涡度相关观测系统2/7实验报告：涡度相关观测系统一、实验目的利用观测数据进行涡度相关计算，包括平均速度、标准差、平均风向、湍流强度、雷诺应力和摩擦速度等。二、实验原理1.涡度以及涡度相关：气流可以描述为很多湍涡的水平流动。这些湍涡具有3维结构，携带着不同的下垫面或者环境属性，比如大气成分浓度，温度和湿度。涡度相关：通过快速测定大气的物理量(如温度，湿度，CO2浓度等)与垂直风速的协方差来计算湍流通量的一种方法，它是基于大气湍流理论和数据统计分析相结合的一种技术。目前已经被广泛应用于地气交换的测量。通量：单位时间单位面积内通过了多少物质或能量垂直风速和某种属性的协方差。通量取决于通过面积的物质或能量的量、物质或能量通过的面积大小、通过关注面积的时间是多长湍流：如果湍流可以看做一系列关于平均值的波动，即雷诺平均，那么任何瞬时物理量的值都等于给定时间(比如30min)内的均值加瞬时分量（脉动值）。两个物理量的乘积等于物理量均值的乘积加一个协方差，''yxyxxy，任何脉动量的均值等于0，0'x，两个物理量之和的均值等于物理量平均值的和yxyx，'uuu，'，'TTT。涡动相关通量测量的主要假设：在A点测量应该代表从A点的整个上风向区域；测量在感兴趣的边界层内进行（常通量层）；对于某一测量，源区面积要足够大，这样测量通量才仅代表感兴趣的区域；净垂直通量是通过湍涡产生的；仪器可以以很高的频率探测非常小的变化；在一个适当的时间间隔内，没有干空气的垂直质量流（即平均垂直风速为零）。显热通量''TwcHpρ=空气密度cp=比气体热容潜热通量''vwLLEρv=水汽密度L=汽化潜热CO2通量''ccwFρc=CO2密度4/12''2''*wvwuu2.观测仪器涡度相关观测系统3/7HMP45C空气温湿度探头由Vaisala公司生产的空气温湿度探头，探头使用电容聚合体作为湿度传感器，PRT作为温度传感器。相对湿度：范围：0.8～100%，精度：±2%；温度：范围：-40～+60℃，精度：±0.1℃。CNR4净辐射传感器CNR4分开测量太阳光辐射和远红外辐射。太阳光辐射：太阳光辐、反射的太阳光辐射、反照率。远红外辐射：来自天空的远红外辐射、来自土壤表面的辐射HFP01土壤特通量板量程：±100Wm-2；精度：读数的3%CS616土壤水分传感器TCAV土壤平均温度传感器CS150三维超声及CO2/H2O分析仪3.平均时间和采样频率平均时间的确定：涡度相关系统所采集的数据是连续不断的，为了计算不同时间段的通量（通量的日变化），必须要对数据进行分时间段分别计算。研究表明，不同的时间划分方法会得到不完全同的结果，一般情况下，平均时间取10－60min是合适的，30min是目前最常用的。采样频率的确定：采样频率就是指要求采集器按照什么速度记录数据，理论上，采样频率应该是越高越好，但毕竟仪器和采集器的响应时间是有限的。确定采样频率要考虑采集器和仪器的反应时间，数据存储设备容量，精度要求等影响，同时要考虑仪器设备的成本。5－20Hz是可以接受的，现在一般是10Hz（每秒钟记录10组数据），研究它的目的就是要明白降低采样频率可能造成的影响有多大，采样频率最低数是多少。4.数据校正的必要性虽然涡度相关方法是目前最好的测定地气交换的微气象方法之一，但涡度相关理想的观测条件是难以达到的，如地形平坦，地势开阔，下垫面均匀，风向不能从仪器支臂后过来等，即使都满足也难以保证结果完全准确。另外，大气层结条件，大气湍流强弱，仪器精度和电子噪声等都会影响数据的采集和计算结果。所以，原始数据和通量数据的校正是必须的。5.原始数据的校正仪器的标定（确定仪器系数，仪器性能漂移校正等）要按照仪器所要求的时间和方法进行标定，不同的要求方法不同。经常检查原始数据，发现问题及时解决。异常数据的剔除和插补：方法：极限值法（绝对值法）；标准差法（相对值法）3倍。6.密度校正（WPL校正）WPL校正就是空气密度变化的校正。一般情况下空气密度作为常数来处理，对于感热通量是可以的，但对于象CO2,H2O等微量痕量气体的交换，密度变化对结果很敏感。因为红外仪器现在测量的是空气的CO2的绝对浓度(mg/m3)而不是相对浓度(ppm，百万比)，同样的绝对浓涡度相关观测系统4/7度在不同的温度下其相对浓度是不一样的，空气温度的微弱变化引起密度的变化对于象CO2气体这样量级的物质来说是非常敏感的，而计算公式要求的是相对浓度s和垂直速度w的协方差。7.坐标旋转改变坐标系统的方向，目的是通过旋转改变原来的坐标系统，使得我们更容易计算或者表达更简洁直观。在新坐标系统，垂直方向上的风速要使得它尽可能的小，甚至为零，通过旋转坐标就可以达到目的。第一次：目的是使得X轴平行于平均气流（风向）,即水平横向风速为0。第一次变换后的3个方向的风速u2,v2和w2分别为：sincos112vuu，sincos112uvv，12ww，)/arctan(11uv。第二次：第2次旋转的目的是使得垂直风速为零)/arctan(sincossincos2223223223uwvvu第三次：经过两次旋转后，其风速仪在新平面的方向与最初的方向有关.虽然此时经向的平均风速和垂直平均风速等于零,但他们的协方差不一定为零.为了避免这种模糊现象可能对动量计算产生的影响,McMillen建议再进行一次旋转,即第3次旋转,使得,其计算公式如下:)2arctan(21sincossincos23233333433434wvwvv8.能量闭合LE为潜热通量,H为显热通量,Rn为净辐射,G为土壤热通量,S为冠层热储量,Q为附加能量源汇的总和,因Q项值很小常常被忽略.QSGRnHLE三、实验仪器有编程功能的计算机。涡度相关观测系统5/7四、数据剔除计算数据的标准差和平均值，利用3sigma原则，剔除所有距平值大于3sigma的数值，注意筛选数据时有一个方向上不符合即全部筛去。数据列表略去，前30分钟舍去1112个数据，后30分钟舍去461个。五、均值、均方差、平均风向计算niiunu11，niivnv11，niiwnw11；niiuuun121，niivvvn121，niiuvw前30分钟平均值1.3729-0.34850.0532标准差0.41140.27400.1853后30分钟平均值1.7167-0.22160.0692标准差0.52720.57180.2464Theta1=-0.2486，Theta2=-0.1284六、坐标旋转后的均值、均方差计算uvw前30分钟平均值1.41740.00000.0000标准差0.40400.28340.1872后30分钟平均值1.73230.00000.0000标准差0.51670.51140.3703七、湍流强度UIuu，UIvv，UIwwIxIyIz前30分钟0.28500.20000.1321后30分钟0.29830.29520.2138八、雷诺应力niiivunvu1''''1，niiiwunwu1''''1，niiivwnvw1''''1x方向y方向z方向前30分钟0.0200-0.0157-0.0018后30分钟0.0403-0.00320.1221涡度相关观测系统6/7九、摩擦速度4/12''2''*wvwuu1*u=0.2745，2*u=0.3460十、附MATLAB程序数据剔除：functiondata=sigma(data)datmean=mean(data);sigma3=3*std(data);sgm=zeros(18000,3);fori=1:3sgm(:,i)=data(:,i)-datmean(i);data(find(sgm(:,i)sigma3(i)|sgm(:,i)-sigma3(i)))=NaN;enddata=data(find(~isnan(data(:,1))&~isnan(data(:,2))&~isnan(data(:,3))),:);end第一次旋转：function[dat,theta]=peach1(data)dat=data;mn=mean(data);theta=atan(mn(2)/mn(1));dat(:,1)=data(:,1)*cos(theta)+data(:,2)*sin(theta);dat(:,2)=data(:,2)*cos(theta)-data(:,1)*sin(theta);end第二次旋转：functiondat=peach2(data)dat=data;mn=mean(data);theta=atan(mn(3)/mn(1));dat(:,1)=data(:,1)*cos(theta)+data(:,3)*sin(theta);dat(:,3)=data(:,3)*cos(theta)-data(:,1)*sin(theta);end第三次旋转：functiondat=peach3(data)dat=data;mn=mean(data);theta=0.5*atan(2*mn(2)*mn(3)/(mn(2)^2-mn(3)^2));dat(:,2)=data(:,2)*cos(theta)+data(:,3)*sin(theta);涡度相关观测系统7/7dat(:,3)=data(:,3)*cos(theta)-data(:,2)*sin(theta);end数据处理：clearimportdata3data=cell2mat(ans);data=[str2num(data(:,29:35)),str2num(data(:,39:45)),str2num(data(:,49:55))];dat1=data(1:18000,:);dat2=data(18001:36000,:);dat1=sigma(dat1);dat1mean1=mean(dat1);dat1std1=std(dat1);v1=(dat1mean1(1)^2*dat1mean1(3)^2+dat1mean1(2)^2*dat1mean1(3)^2)^0.25;[dat1,theta1]=peach1(dat1);dat1=peach2(dat1);dat1=peach3(dat1);dat1mean2=mean(dat1);dat1std2=std(dat1);I1=dat1std2/dat1mean2(1);ln1=[mean(dat1(:,1).*dat1(:,2)),mean(dat1(:,1).*dat1(:,3)),mean(dat1(:,3).*dat1(:,2))];dat2=sigma(dat2);dat2mean1=mean(dat2);dat2std1=std(dat2);v2=(dat2mean1(1)^2*dat2mean1(3)^2+dat2mean1(2)^2*dat2mean1(3)^2)^0.25;[dat2,theta2]=peach1(dat2);dat2=peach2(dat2);dat2=peach3(dat2);dat2mean2=mean(dat2);dat2std2=std(dat2);I2=dat2std2/dat2mean2(1);ln2=[mean(dat2(:,1).*dat2(:,2)),mean(dat2(:,1).*dat2(:,3)),mean(dat2(:,3).*dat2(:,2))];