北航惯性导航综合实验一实验报告

实验一

陀螺仪关键参数测试与分析实验

加速度计关键参数测试与分析实验

二零一三年五月十二日

实验一 陀螺仪关键参数测试与分析实验

一、   实验目的

通过在速率转台上的测试实验,增强动手能力和对惯性测试设备的感性认识;通过对陀螺仪测试数据的分析,对陀螺漂移等参数的物理意义有清晰的认识,同时为在实际工程中应用陀螺仪和对陀螺仪进行误差建模与补偿奠定基础。

二、   实验内容

利用单轴速率转台,进行陀螺仪标度因数测试、零偏测试、零偏重复性测试、零漂测试实验和陀螺仪标度因数与零偏建模、误差补偿实验。

三、   实验系统组成

单轴速率转台、MEMS 陀螺仪(或光纤陀螺仪)、稳压电源、数据采集系统与分析系统。

四、   实验原理

1.   陀螺仪原理 

陀螺仪是角速率传感器,用来测量载体相对惯性空间的角速度,通常输出与角速率对应的电压信号。也有的陀螺输出频率信号(如激光陀螺)和数字信号(把模拟电压数字化)。以电压表示的陀螺输出信号可表示为:

             (1-1)

式中是与比力有关的陀螺输出误差项,反映了陀螺输出受比力的影响,本实验不考虑此项误差。因此,式(1-1)简化为

                   (1-2)

由(1-2)式得陀螺输出值所对应的角速度测量值:

                        (1-3)

对于数字输出的陀螺仪,传感器内部已经利用标度因数对陀螺仪模拟输出进行了量化,直接输出角速度值,即:

                      (1-4)

是是陀螺仪的零偏,物理意义是输入角速度为零时,陀螺仪输出值所对应的角速度。且

                            (1-5)

精度受陀螺仪标度因数、随机漂移、陀螺输出信号的检测精度和的影响。通常表现为有规律性,可通过建模与补偿方法消除,表现为随机特性,可通过信号滤波方法抵制。因此,准确标定是实现角速度准确测量的基础。

五、   陀螺仪测试实验步骤

1)   标度因数和零偏测试实验

a. 接通电源,预热一定时间;

b. 陀螺工作稳定后,测量静止情况下陀螺输出并保存数据;

c. 转台正转,测试陀螺仪输出,停转;转台反转,测试陀螺仪输出,停转。在正转和反转时测试陀螺仪输出量,并分别保存数据;

d. 改变转台输入角速率重复步骤c,正负角速率的速率档分别不少于5 个(按军标要求是11 个);

e. 转速结束后,当转台静止时,采集陀螺仪输出数据,并保存。

f.根据最小二乘法公式

                  (1-6)

                      (1-7)

计算陀螺标度因数和零偏。

2)   零漂测试(零偏稳定性)

在静止下采集陀螺仪数据,并由测试数计算陀螺仪零偏稳定性。军标中通常的测试时间是1 小时,并对所采集的数据进行1 秒、10 秒及100秒等不同时间的平滑。本实验中可采集数据10 分钟左右,并分别进行1 秒、10 秒及100 秒平滑。

按如下公式

                   (1-8)

计算陀螺仪零偏稳定性,并进行比较。

3)   零偏重复性测试

a. 令转台某角速度200下进行正转,转速平稳后,采集陀螺输出数据,并保存。

b. 令转台某角速度-200下进行反转,转速平稳后,采集陀螺输出数据,并保存。

c. 按计算陀螺零偏;

d. 关掉陀螺电源,并重新启动,重复步骤a、b;

e. 重复步骤d 进行3-5 次,共得到陀螺零偏5-7 个;

f.  对5-7个陀螺零偏按下式(1-9)

                      (1-9)

求均方差,得零偏重复性指标。

六、   实验结果

1.  数据处理

将原始数据剔除后绘图如下

2.  计算陀螺标度因数和零偏

根据陀螺在10°/s,20°/s,40°/s,60°/s ,80°/s角速率下正反转的输出,分别求得正转下陀螺的标度因数和零偏,及反转下陀螺的标度因数和零偏,然后求的均值。

                =  0.9901

                = 0.0358

3.  零偏稳定性对所采集的数据进行1 秒、10 秒及100秒等不同时间的平滑,如下图。

零漂计算结果(1000s平滑):Bs= 0.0144

4.  零偏重复性

以角速度40°/s正反转,共采集5组数据

零偏重复性:0.009309

七,实验小结

由零漂平滑后的结果可知,对采集的数据平滑时间长可以提高零偏的稳定性。

八,源程序


%%%%加载数据%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Gyro_0end=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Gyro_data\1 标度因数和零偏测试\Gyro_0end.txt');

Gyro_0start=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Gyro_data\1 标度因数和零偏测试\Gyro_0start.txt');

Gyro_10n=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Gyro_data\1 标度因数和零偏测试\Gyro_10n.txt');

Gyro_10p=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Gyro_data\1 标度因数和零偏测试\Gyro_10p.txt');

Gyro_20n=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Gyro_data\1 标度因数和零偏测试\Gyro_20n.txt');

Gyro_20p=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Gyro_data\1 标度因数和零偏测试\Gyro_20p.txt');

Gyro_40n=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Gyro_data\1 标度因数和零偏测试\Gyro_40n.txt');

Gyro_40p=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Gyro_data\1 标度因数和零偏测试\Gyro_40p.txt');

Gyro_60n=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Gyro_data\1 标度因数和零偏测试\Gyro_60n.txt');

Gyro_60p=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Gyro_data\1 标度因数和零偏测试\Gyro_60p.txt');

Gyro_80n=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Gyro_data\1 标度因数和零偏测试\Gyro_80n.txt');

Gyro_80p=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Gyro_data\1 标度因数和零偏测试\Gyro_80p.txt');

%%%%%%%%%%%%剔除不合格数据%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Gyro_10p=Gyro_10p(find((Gyro_10p>9)&(Gyro_10p<11)));

Gyro_20p=Gyro_20p(find((Gyro_20p>15)&(Gyro_20p<25)));

Gyro_40n=Gyro_40n(find((Gyro_40n>-50)&(Gyro_40n<0)));

Gyro_40p=Gyro_40p(find((Gyro_40p>35)&(Gyro_40p<45)));

Gyro_60p=Gyro_60p(find((Gyro_60p>50)&(Gyro_60p<70)));

Gyro_80p=Gyro_80p(find((Gyro_80p>70)&(Gyro_80p<90)));

for i=1:11145

    k(i)=i;

end

plot(k,Gyro_0end(1:11145,1),'r',k,Gyro_0start(1:11145,1),'r',k,Gyro_10n(1:11145,1),'r',k,Gyro_10p(1:11145,1),'r',k,Gyro_20n(1:11145,1),'r',k,Gyro_20p(1:11145,1),'r',k,Gyro_40n(1:11145,1),'r',k,Gyro_40p(1:11145,1),'r',k,Gyro_40n(1:11145,1),'r',k,Gyro_60n(1:11145,1),'r',k,Gyro_60p(1:11145,1),'r',k,Gyro_80n(1:11145,1),'r',k,Gyro_80p(1:11145,1),'r');

title('剔除数据后','fontsize',12);

xlabel('时间t(s)','fontsize',12);

ylabel('度/秒','fontsize',12);

%%%%%%%%%%%%%5555计算标度因数%%%%%%%%%%%%%%%%%%%%%%%%5

Gyro_0end1=mean(Gyro_0end);

Gyro_0start1=mean(Gyro_0start);

Gyro_10n1=mean(Gyro_10n);

Gyro_10p1=mean(Gyro_10p);

Gyro_20n1=mean(Gyro_20n);

Gyro_20p1=mean(Gyro_20p);

Gyro_40n1=mean(Gyro_40n);

Gyro_40p1=mean(Gyro_40p);

Gyro_60n1=mean(Gyro_60n);

Gyro_60p1=mean(Gyro_60p);

Gyro_80p1=mean(Gyro_80p);

Gyro_80n1=mean(Gyro_80n);

%%%%%%%%求正转标度因数%%%%%%

F=[Gyro_10p1  Gyro_20p1  Gyro_40p1  Gyro_60p1  Gyro_80p1];

W=[10 20 40 60 80];

J=[Gyro_10p1*10   Gyro_20p1*20  Gyro_40p1*40   Gyro_60p1*60  Gyro_80p1*80];

KG0=(sum(J)-(sum(F)*sum(W))/5)/(sum(W.^2)-(sum(W)*sum(W))/5); %%%%0.9905

%%%%%%求反转标度因数%%%%%%%%%%%

F1=[Gyro_10n1  Gyro_20n1  Gyro_40n1  Gyro_60n1  Gyro_80n1];

W1=[10 20 40 60 80];

J1=[Gyro_10n1*(10)   Gyro_20n1*(20)  Gyro_40n1*(40)   Gyro_60n1*(60)  Gyro_80n1*(80)];

KG1=-(sum(J1)-(sum(F1)*sum(W1))/5)/(sum(W1.^2)-(sum(W1)*sum(W1))/5);   %%%%%0.9895

KG=(KG0+KG1)/2;   %%%%%0.9901

%%%%%%%%%%%%%%%%%%%%%%%%%%求零偏%%%%%%%%%%%%%%%%%%%%55

F0=-(sum(F1)/5+KG*sum(W1)/5);

F01=sum(F)/5-KG*sum(W)/5;

F0=(F0+F01)/2;          %%%%%%%%%%%%%%%%%% F0=0.3580

%%%%%%%%%%%加载静止时的数据%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Gyro_result=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Gyro_data\2 零偏稳定性测试\Gyro_result.txt');

Gyro_result=Gyro_result(find((Gyro_result>-0.8)&(Gyro_result<0.8))); %%%%%%%%%%剔除数据

Gyro_result=smooth(Gyro_result,128000);%利用移动平均法做平滑处理1000s

Gyro_result1=smooth(Gyro_result,1280);  %利用移动平均法做平滑处理10s

Gyro_result2=smooth(Gyro_result,12800);%利用移动平均法做平滑处理100s

for i=1:206224

    u(i)=i;

end

figure;                                                                                       %新建一个图形窗口

plot(u,Gyro_result,'g');                                                                             %绘制加噪波形图

hold on;

plot(u,Gyro_result1,'r');    %绘制平滑后波形图

hold on;

plot(u,Gyro_result2,'k');

xlabel('时间t(s)','fontsize',12);

ylabel('零漂平滑后结果','fontsize',12);

legend('1000s平滑','100s平滑','10s平滑');

b0=mean(Gyro_result); %%%

b1=mean(Gyro_result1);  %%%  -0.0479

b2=mean(Gyro_result2); %%%% -0.0482

c0=sum((Gyro_result-b0).^2)/206223;

B0=sum(c0.^0.5)/KG;   %%%%B0= 0.0144

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%零偏重复性测试

Gyro_result11=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Gyro_data\3 零偏重复性测试\1_40n.txt');

Gyro_result22=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Gyro_data\3 零偏重复性测试\2_40p.txt');

Gyro_result33=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Gyro_data\3 零偏重复性测试\3_40p.txt');

Gyro_result44=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Gyro_data\3 零偏重复性测试\4_40n.txt');

Gyro_result55=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Gyro_data\3 零偏重复性测试\5_40n.txt');

Gyro1=Gyro_result11(1:5000,:);

Gyro2=Gyro_result22(1:5000,:);

Gyro3=Gyro_result33(1:5000,:);

Gyro4=Gyro_result44(1:5000,:);

Gyro5=Gyro_result55(1:5000,:);

Gyro1=Gyro1(find((Gyro1>-50)&(Gyro1<0))); %%%%%%%%%%剔除数据

Gyro3=Gyro3(find((Gyro3>35)&(Gyro3<50)));

b1=mean(Gyro1); %%%

b2=mean(Gyro2);  %%%  -0.0479

b3=mean(Gyro2);%%%% -0.0482

b4=mean(Gyro4);  %%%  -0.0479

b5=mean(Gyro5);%%%% -0.0482

c1=sum((Gyro1-b1).^2)/4959;

BS1=sum(c1.^0.5)/KG;   %%%%B0= 0.0144

c2=sum((Gyro2-b2).^2)/4959;

BS2=sum(c2.^0.5)/KG;

c3=sum((Gyro3-b3).^2)/4959;

BS3=sum(c3.^0.5)/KG;

c4=sum((Gyro4-b4).^2)/4959;

BS4=sum(c4.^0.5)/KG;

c5=sum((Gyro5-b5).^2)/4959;

BS5=sum(c5.^0.5)/KG;

BSS=[BS1 BS2 BS3 BS4 BS5];

h0=mean(BSS);

Br=sum((BSS-h0).^2)/4

Br=sum(Br.^0.5);        %%%Br = 8.6662e-005

实验一加速度计关键参数测试与分析实验

一、   实验目的

通过在位置转台上的测试实验,增强学生的动手能力和对惯性测试设备的感性认识;通过对加速度计测试数据的分析,让学生对比力、加速度计标度因数和偏置等参数的物理意义有清晰的认识。为在实际工程中应用加速度计和对加速度计进行误差建模与补偿奠定基础。

二、   实验内容

利用双轴位置转台,进行加速度计标度因数测试、偏置测试、偏置重复性测试和加速度计标度因数与偏置建模、误差补偿实验。

三、   实验原理

加计的误差模型与陀螺仪类似,不在赘述。不同之处在于,加计的输出包含了运动加速度和引力加速度两项。当加计静止或匀速运动时,运动加速度为零,加计输出的是引力加速度。可根据其输出信号与重力矢量的投影间的关系对加速度计的指标参数进行标定。依据公式为

实际操作时,调整的是加计的输入轴与垂直方向的角度,故

加计的单位为(g)的话,

实验目的是标定k和f0.

1)   标度因数和偏置测试实验

a. 把俯仰现调到90 度位置;

b. 调节方位轴,使加速度计敏感轴近似指向地心;

c. 设当然角位置为0°,加速度计通电,记录数据并求出平均值;

d. 在0 到360 度范围内改变输入角度值jθ,记录输出电压数据,并求出平均值jpF 。

通常采取等间隔采样方法,即每次变化角度60 度、30 度或更小间隔。

e.按最小二乘公式

                 (1-10)

                    (1-11)

计算加速度计标度因数和偏置

2)   偏置测试实验

a. 把俯仰现调到近水平位置(任意角度α);

b. 采集加速度计输出数据,保存数据,并计算均值

c. 将位置转台方位轴旋转180 度;

d. 采集加速度计输出数据,保存数据,并计算均值;

e.按(1-12)式

                  (1-12)

计算加速度计偏置。

3)   偏置重复性测试

a. 可按以上两种方法这一计算加速度计偏置

b. 重复步骤a4 到6 次,并保存数据;

c. 按(1-13)式

                   (1-13)

计算加速度计偏置重复性。

四、   实验结果

1.  零偏和标度因数

将原始数据剔除异常后绘图如下:以指向地心为0度;

数据记录:

与平均值在0到360 度范围内变化数据

           

计算结果:

标度因数: -1.0159

零偏:  0.0366

2、偏置测试实验

把俯仰现调到近水平位置;采集加速度计输出数据,保存数据,并计算均值将位置转台方位轴旋转180 度;采集加速度计输出数据,保存数据,并计算均值; 按式Bias计算加速度计偏置。

有实验数据可求Bias=   -0.0112

3零偏重复性计算

在同一条件下计算得4组数据的零偏分别为:

相应的重复性为:Br= 0.0448

五,实验结果分析

   加速度计采用旋转180°计算均值求零偏,方法简单,但是结果变化大,不精确,应当多组求平均,或采用最小二乘法。

六,实验源程序:

%%%%%%%%%加载数据%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

A0=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\1 标度因数和偏置测试\0.txt');

A30=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\1 标度因数和偏置测试\30.txt');

A60=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\1 标度因数和偏置测试\60.txt');

A90=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\1 标度因数和偏置测试\90.txt');

A120=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\1 标度因数和偏置测试\120.txt');

A150=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\1 标度因数和偏置测试\150.txt');

A180=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\1 标度因数和偏置测试\180.txt');

A210=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\1 标度因数和偏置测试\210.txt');

A240=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\1 标度因数和偏置测试\240.txt');

A270=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\1 标度因数和偏置测试\270.txt');

A300=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\1 标度因数和偏置测试\300.txt');

A330=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\1 标度因数和偏置测试\330.txt');

%%%%%%%%%%%%%%%%%%%剔除不合格数据%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

A0=A0(find((A0>-2)&(A0<0)));

A60=A60(find((A60>-1)&(A60<0)));

A90=A90(find((A90>0)&(A90<0.5)));

A120=A120(find((A120>0)&(A120<1)));

A180=A180(find((A180>0)&(A180<2)));

A270=A270(find((A270>0)&(A270<0.5)));

A300=A300(find((A300>-1)&(A300<1.5)));

 plot (A120, 'DisplayName', 'A120', 'YDataSource', 'A120'); hold all; plot (A150, 'DisplayName', 'A150', 'YDataSource', 'A150'); plot (A180, 'DisplayName', 'A180', 'YDataSource', 'A180'); plot (A210, 'DisplayName', 'A210', 'YDataSource', 'A210'); plot (A240, 'DisplayName', 'A240', 'YDataSource', 'A240'); plot (A270, 'DisplayName', 'A270', 'YDataSource', 'A270'); plot (A30, 'DisplayName', 'A30', 'YDataSource', 'A30'); plot (A300, 'DisplayName', 'A300', 'YDataSource', 'A300'); plot (A330, 'DisplayName', 'A330', 'YDataSource', 'A330'); plot (A60, 'DisplayName', 'A60', 'YDataSource', 'A60'); plot (A90, 'DisplayName', 'A90', 'YDataSource', 'A90');  hold off; figure(gcf)

title('剔除数据后','fontsize',12);

%%%%%%%%%%%%%%%%%%%%%%%%%%%5555计算标度因数%%%%%%%%%%%%%%%%%%%%%%%%5

f0=mean(A0);

f30=mean(A30);

f60=mean(A60);

f90=mean(A90);

f120=mean(A120);

f150=mean(A150);

f180=mean(A180);

f210=mean(A210);

f240=mean(A240);

f270=mean(A270);

f300=mean(A300);

f330=mean(A330);

%%%%%%%%求正转标度因数%%%%%%

F=[ f0 f30  f60  f90  f120  f150  f180  f210   f240  f270  f300  f330];  %%%%%%%%以指向地心为初始零度

A=[ 1 1*cos(pi/6) 1*cos(pi/3) 1*cos(pi/2) 1*cos(pi*2/3) 1*cos(pi*5/6) 1*cos(pi) 1*cos(pi*7/6) 1*cos(pi*4/3) 1*cos(pi*3/2) 1*cos(5/3*pi) 1*cos(22/12*pi)];

J=[1*f0  1*cos(pi/6)*f30 1*cos(pi/3)*f60  1*cos(pi/2)*f90  1*cos(pi*2/3)*f120  1*cos(pi*5/6)*f150  1*cos(pi)*f180  1*cos(pi*7/6)*f210  1*cos(pi*4/3)*f240 1*cos(pi*3/2)*f270 1*cos(5/3*pi)*f300 1*cos(22/12*pi)*f330 ];

KG0=(sum(J)-(sum(F)*sum(A))/12)/(sum(A.^2)-(sum(A)*sum(F))/12); %%%%KG0=-1.0159

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%求零偏%%%%%%%%%%%%%%%%%%%%55

F01=sum(F)/12-KG0*sum(A)/12;%%%%%%%%%F01=-0.036595

%%%%%%%%%%%%%%加载静止时的数据%%%%%%%%%%%%%%%%%%%%%%%%%%%%

S1=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\2 偏置稳定性测试\偏置测试0度.txt');

S2=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\2 偏置稳定性测试\偏置测试180度.txt');

ff0=(mean(S1)+mean(S2))/2;  %%ff0=-0.011212g

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%零偏重复性测试

D10=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\3 加计偏置重复测试\1_0.txt');

D11=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\3 加计偏置重复测试\1_180.txt');

D20=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\3 加计偏置重复测试\2_0.txt');

D21=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\3 加计偏置重复测试\2_180.txt');

D30=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\3 加计偏置重复测试\3_0.txt');

D31=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\3 加计偏置重复测试\3_180.txt');

D50=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\3 加计偏置重复测试\5_0.txt');

D51=load('E:\惯性器件综合实验\惯性导航试验数据\1\惯导实验1实验数据\Acc_data\3 加计偏置重复测试\5_180.txt');

d1=(mean(D10)+mean(D11))/2;

d2=(mean(D20)+mean(D21))/2;

d3=(mean(D30)+mean(D31))/2;

d5=(mean(D50)+mean(D51))/2;

d=[d1 d2 d3 d5];

h0=mean(d);

Br=sum((d-h0).^2)/3 

Br1=sum(Br.^0.5);      

相关推荐