测量数据处理实验报告

院系____________________班级____________________姓名____________________学号____________________目录（页码根据自己需要填写）实验项目须知实验一秩亏自由网上机实验实验二极大验后滤波、推估上机实验实验三最小二乘配置上机实验实验四平差系统的统计假设检验上机实验编写实验报告要求实验项目须知一、课程介绍本课程是测绘工程类专业的专业选修课程之一，是在学生学习了误差理论与测量平差知识之后又开设的一门专业的数据处理课程，该课程主要讲授近现代测量数据处理理论、模型和方法，同时紧密结合计算机编程，从而提高学生的计算机应用能力、测量数据处理能力，培养学生养成严谨的科学态度，提高分析和解决问题的能力，更好的适应今后工作和学习。二、测量数据处理实验教学要求1.掌握稳健估计方法处理测量数据；2.掌握极大验后滤波、推估理论和方法处理测量数据；3.初步掌握最小二乘配置理论和方法处理测量数据；4.初步掌握工程或科研项目中测量数据处理的理论、步骤、方法和技巧；5.掌握测量数据处理中常用的计算机语言和软件；6.服从实验指导教师的规定，认真、按时、独立完成任务每次实验结束，提交书写工整或按照规范打印实验报告。7.在实验过程中，还应遵守纪律，爱护实验室设备、离开后将自己的随身物品带走，并记录仪器使用状况、将凳子整理好并放到电脑桌的下方。8.上机实验后，提交实验报告，成绩评定重点考核算法、处理结果与结果分析。三、实验选择要求可根据课堂讲课具体内容，完成4个实验，并要求写出实验报告。实验一秩亏自由网上机实验实验目的掌握各种秩亏自由网平差的函数模型和随机模型；熟悉各种秩亏自由网平差解算的若干公式；结合算例能利用matlab进行编程计算。实验内容在下图水准网中，观测高差、距离和各待定点高程近似值见教材P43表2-1，分别进行下列自由网平差：以6号点为固定点的经典自由网平差；以重心基准的自由网平差（Pχ=E）；(3)以1，2，5，6四个点为拟稳基准的拟稳平差.参考程序：（1）首先建立误差方程B=[-11000;-10010;-10001;0-1100;0-1010...;0-1000;00-110;00-101;0000-1];l=[-1;4;5;8;16;2;-13;-20;-1];S=[204.1188.7344.8149.2142.9250128.298196.1];P=diag(10./S);N=B'*P*B;x=inv(N)*B'*P*lV=B*x-lQx=inv(N)(2)首先建立误差方程B=[-110000;-100100;-100010;0-11000;0-10100...;0-10001;00-1100;00-1010;0000-11];l=[-1;4;5;8;16;2;-13;-20;-1]S=[204.1188.7344.8149.2142.9250128.298196.1];P=diag(10./S);S1=[1;1;1;1;1;1];W=B'*P*l;Xr=inv(B'*P*B+S1*S1')*WV=B*Xz-lQxr=inv(B'*P*B+S1*S1')-inv(B'*P*B+S1*S1')*S1*S1’*inv(B'*P*B+S1*S1')(3)首先建立误差方程B=[-110000;-100100;-100010;0-11000;0-10100...;0-10001;00-1100;00-1010;0000-11];l=[-1;4;5;8;16;2;-13;-20;-1]S=[1;1;1;1;1;1];S1=[204.1188.7344.8149.2142.9250128.298196.1];P=diag(10./S1);N=B'*P*BW=B’*P*l;px=[110011]Px=diag(px)Ss=Px*SXs=inv(B'*P*B+Ss*Ss')*WV=B*Xs-lQxs=inv(B'*P*B+Ss*Ss')*N*inv(B'*P*B+Ss*Ss')实验心得和体会实验二极大验后滤波、推估上机实验实验目的和要求掌握极大验后滤波与推估的基本原理和计算公式；能根据观测值和给定的数学模型计算滤波信号的协方差矩阵的估计矩阵；3.能使用matlab进行有关的计算；实验内容:1.（基本要求题目）沿A、B连线在A、1、2、3、4等五个点测定了大气温度，获得18组观测值,其结果见下表，已知各点的距离分别为4.511,10.747,16.753,22.220（km）.假定AB连线上的大气温度，是一个以各点至A的距离S为因素的平稳随机过程，试估计此随机过程的协方差函数。观测数据温度测站序号A1234116.415.916.315.815.4216.316.116.215.915.5316.816.415.916.316.3416.616.216.316.216.3516.816.716.116.216.3616.616.616.416.216.3716.716.617.016.316.2816.816.516.916.516.9916.916.416.116.616.21016.116.616.217.117.01116.516.716.716.816.51217.116.916.616.816.71316.816.916.817.117.01417.016.916.716.917.01517.116.916.716.917.01616.216.016.316.516.41716.115.916.216.416.41816.115.916.016.116.51916.316.016.616.517.02016.116.015.916.516.3参考程序：X=[16.415.916.315.815.4;16.316.116.215.915.5;16.816.415.916.316.3;...16.616.216.316.216.3;16.816.716.116.216.3;...16.616.616.416.216.3;16.716.617.016.316.2;...16.816.516.916.515.9;16.916.416.116.616.2;...16.116.616.217.117.0;16.516.716.716.816.5;...17.116.916.616.816.7;16.816.916.817.117.0;...17.016.916.716.917.0;17.116.916.716.917.0;...16.216.016.316.516.4;16.115.916.216.416.4;...16.115.916.016.116.5;16.316.016.516.617.0;16.116.015.916.516.3]Dx1=cov(X)s1=[0;4.510;10.746;16.754;22.230]A=[ones(5,1),s1,s1.^2,s1.^3,s1.^4]s2=[0;6.236;12.242;17.709]B=[ones(4,1),s2,s2.^2,s2.^3,s2.^4]s3=[0;6.006;11.473]C=[ones(3,1),s3,s3.^2,s3.^3,s3.^4]s4=[0;5.467]D=[ones(2,1),s4,s4.^2,s4.^3,s4.^4]E=[10000]Y=[A;B;C;D;E]L=[0.1168;0.0958;0.0506;0.0466;0.0394;0.1362;0.0705;0.1023;0.1047;0.1073;0.0567;0.0346;...0.1529;0.1561;0.2203]b=inv(Y’*Y)*Y’*Lb=[b(1),b(2)*10,b(3)*100,b(4)*1000,b(5)*10000]%将km单位化为10km单位S=[04.51110.74716.75322.220;4.51106.23612.24217.709;10.7476.23606.00611.473;...16.75312.2426.00605.467;22.22017.70911.4735.4670]b=fliplr(b)fori=1:25S(i)=polyval(b,S(i)/10);endS%即为滤波信号协方差阵的估计矩阵Warning:UnabletocreatepersonalMATLABworkfolder:C:\yonghu\????\MATLABWarning:Userpathmustbeanabsolutepathandmustexistondisk.>>X=[16.415.916.315.815.4;16.316.116.215.915.5;16.816.415.916.316.3;...16.616.216.316.216.3;16.816.716.116.216.3;...16.616.616.416.216.3;16.716.617.016.316.2;...16.816.516.916.515.9;16.916.416.116.616.2;...16.116.616.217.117.0;16.516.716.716.816.5;...17.116.916.616.816.7;16.816.916.817.117.0;...17.016.916.716.917.0;17.116.916.716.917.0;...16.216.016.316.516.4;16.115.916.216.416.4;...16.115.916.016.116.5;16.316.016.516.617.0;16.116.015.916.516.3]X=16.400015.900016.300015.800015.400016.300016.100016.200015.900015.500016.800016.400015.900016.300016.300016.600016.200016.300016.200016.300016.800016.700016.100016.200016.300016.600016.600016.400016.200016.300016.700016.600017.000016.300016.200016.800016.500016.900016.500015.900016.900016.400016.100016.600016.200016.100016.600016.200017.100017.000016.500016.700016.700016.800016.500017.100016.900016.600016.800016.700016.800016.900016.800017.100017.000017.000016.900016.700016.900017.000017.100016.900016.700016.900017.000016.200016.000016.300016.500016.400016.100015.900016.200016.400016.400016.100015.900016.000016.100016.500016.300016.000016.500016.600017.000016.100016.000015.900016.500016.3000>>Dx1=cov(X)Dx1=0.12030.10230.05590.03890.03040.10230.14050.07160.08800.08420.05590.07160.10940.05040.03690.03890.08800.05040.13710.14280.03040.08420.03690.14280.2157>>s1=[0;4.510;10.746;16.754;22.230]s1=04.510010.746016.754022.2300>>A=[ones(5,1),s1,s1.^2,s1.^3,s1.^4]A=1.0e+05*0.000000000.00000.00000.00020.00090.00410.00000.00010.00120.01240.13330.00000.00020.00280.04700.78790.00000.00020.00490.10992.4421>>s2=[0;6.236;12.242;17.709]s2=06.236012.242017.7090>>B=[ones(4,1),s2,s2.^2,s2.^3,s2.^4]B=1.0e+04*0.000100000.00010.00060.00390.02430.15120.00010.00120.01500.18352.24600.00010.00180.03140.55549.8350>>s3=[0;6.006;11.473]s3=06.006011.4730>>C=[ones(3,1),s3,s3.^2,s3.^3,s3.^4]C=1.0e+04*0.000100000.00010.00060.00360.02170.13010.00010.00110.01320.15101.7326>>s4=[0;5.467]s4=05.4670>>D=[ones(2,1),s4,s4.^2,s4.^3,s4.^4]D=1.000000001.00005.467029.8881163.3982893.2979>>E=[10000]Y=[A;B;C;D;E]E=10000Y=1.0e+05*0.000000000.00000.00000.00020.00090.00410.00000.00010.00120.01240.13330.00000.00020.00280.04700.78790.00000.00020.00490.10992.44210.000000000.00000.00010.00040.00240.01510.00000.00010.00150.01830.22460.00000.00020.00310.05550.98350.000000000.00000.00010.00040.00220.01300.00000.00010.00130.01510.17330.000000000.00000.00010.00030.00160.00890.00000000>>L=[0.1168;0.0958;0.0506;0.0466;0.0394;0.1362;0.0705;0.1023;0.1047;0.1073;0.0567;0.0346;...0.1529;0.1561;0.2203]b=inv(Y’*Y)*Y’*LL=0.11680.09580.05060.04660.03940.13620.07050.10230.10470.10730.05670.03460.15290.15610.2203b=inv(Y’*Y)*Y’*L|Error:TheinputcharacterisnotvalidinMATLABstatementsorexpressions.>>b=inv(Y'*Y)*Y'*Lb=0.1468-0.0008-0.00290.0003-0.0000>>b=[b(1),b(2)*10,b(3)*100,b(4)*1000,b(5)*10000]S=[04.51110.74716.75322.220;4.51106.23612.24217.709;10.7476.23606.00611.473;...16.75312.2426.00605.467;22.22017.70911.4735.4670]b=0.1468-0.0080-0.28600.2774-0.0706S=04.511010.747016.753022.22004.511006.236012.242017.709010.74706.236006.006011.473016.753012.24206.006005.467022.220017.709011.47305.46700>>b=fliplr(b)fori=1:25S(i)=polyval(b,S(i)/10);endSb=-0.07060.2774-0.2860-0.00800.1468S=0.14680.10760.05810.07920.04000.10760.14680.08720.05890.08230.05810.08720.14680.08980.05780.07920.05890.08980.14680.09600.04000.08230.05780.09600.1468>>2.同第1题，已知测得A、1、2、3、4五个测站上的大气温度为利用推估公式试推估B点大气温度。参考程序：L=[19.2;20.0;18.6;19.3;19.9]Mutg=Mul=mean(L)S1=[2318.48912.2536.2470.78];fori=1:5S1(i)=polyval(b,S1(i)/10);endS1%DtgLTg=Mutg+S1*inv(S)*(L-Mul)三、实验心得与体会实验三最小二乘配置上机实验一、目的与要求1、掌握最小二乘配置的数学模型；2、掌握最小二乘配置的估值公式及其精度评定，能利用matlab软件进行有关计算。二、实验内容对于观测值是重力异常的情况,要求搞清观测方程的列立,其参数和信号分别是什么?根据讲课内容利用matlab进行解算.L=[-0.55-0.230.58-1.8]';Dx=[10.4640.0990.306;0.46410.2020.331;0.0990.20210.147;0.3060.3310.1471];Do=eye(4)*0.09;%观测噪声协方差阵E=eye(4);Dx2=[10.962;0.9621];%推估信号的先验协方差阵Dx2x=[0.4350.7450.2050.581;0.3810.7940.2380.500];%推估信号关于滤波信号的互协方差阵Dx2=1.00000.96200.96201.0000Dx2x=0.43500.74500.20500.58100.38100.79400.23800.5000X=[640440140620]';Y=[480400140180';X0=mean(X);Y0=mean(Y);X1=X-X0;Y1=Y-Y0;G=[ones(4,1),X1./100,Y1./100];%为何要除以100?G=[11.81.8;1-0.21.0;1-3.2-1.6;11.6-1.2];X2=[500460]’;Y2=[300300]’;X2=X2-X0;Y2=Y2-Y0;Gp=[ones(2,1),X2./100,Y2./100];Gp=[10.40;100];inv(Dx+Do);G'*inv(Dx+Do)*G;Ycan=inv(G'*inv(Dx+Do)*G)*G'*inv(Dx+Do)*L;%滤波信号的先验期望为0Dycan=inv(G'*inv(Dx+Do)*G);Ycan=-0.4767-0.49090.3815Dycan=0.46180.00920.01890.00920.0783-0.04840.0189-0.04840.1440Xcan=Dx*inv(Dx+Do)*(L-G*Ycan);X2can=Dx2x*inv(Dx+Do)*(L-G*Ycan);?Dxcan=Dx-Dx*inv(Dx+Do)*(E-G*Dycan*G'*inv(Dx+Do))*Dx;Dx2can=Dx2-Dx2x*inv(Dx+Do)*(E-G*Dycan*G'*inv(Dx+Do))*Dx2x';Dxcanycan=-Dx*inv(Dx+Do)*G*Dycan;Dx2canycan=-Dx2x*inv(Dx+Do)*G*Dycan;Dxcan=0.92080.62460.03260.36470.62460.67450.33660.21200.03260.33660.94430.19620.36470.21200.19620.9565Dx2can=0.84010.77910.77910.7907Dxcanycan=-0.4857-0.0576-0.1738-0.46190.0449-0.1561-0.37740.14910.0510-0.4323-0.17310.2031Dx2canycan=-0.4653-0.0350-0.0252-0.4483-0.0035-0.0473Ogcan=G*Ycan+Xcan;Ogpcan=Gp*Ycan+X2can;Dogcan=[Geye(4)]*[DycanDxcanycan';DxcanycanDxcan]*[G';eye(4)];Dogpcan=[Gpeye(2)]*[DycanDxoycan';DxoycanDx2can]*[Gp';eye(2)];Ogcan=-0.5734-0.20030.5679-1.7942Ogpcan=-0.8155-0.6397Dogcan=?0.08570.0055-0.00220.00110.00550.08310.0028-0.0014-0.00220.00280.08880.00060.0011-0.00140.00060.0897Dogpcan=?0.36310.32950.32950.35602、已知，观测方程为求信号和倾向参数的估值及其估计误差方差。L=[013]';ux=[00]';Dx=[20;02];Do=eye(3)*2;Dxo=[-110;000];B=[-1-1;-10;01];G=[01-1]';Ycan=inv(G'*inv(B*Dx*B'+Do+B*Dxo+Dxo'*B')*G)*G'*inv(B*Dx*B'+Do+B*Dxo+Dxo'*B')*(L-B*ux);Xcan=ux+(Dx*B'+Dxo)*inv(B*Dx*B'+Do+B*Dxo+Dxo'*B')*(L-G*Y-B*ux);Ycan=-0.3333Xcan=-0.66671.3333Dycan=inv(G'*inv(B*Dx*B'+Do+B*Dxo+Dxo'*B')*G);Dxcan=Dx-(Dx*B’+Dxo)*inv(B*Dx*B'+Do+B*Dxo+Dxo'*B')*(eye(3)-G*Dycan*G’*inv(B*Dx*B'+Do+B*Dxo+Dxo'*B'))*(Dxo’+B*Dx);注意:如果直接拷贝过去有错误,可手工输入命令行.在相距均为1.445的5个点处测得F(u)的函数值Li(i=1,2,…5),观测数据列于下表。B=eye(5);G=[10;11.445;12.890;14.335;15.780];L=[0.61081.08632.90344.59256.2714]';u=G(:,2);D=0.1260*exp(-0.36*u.^2)';Do=0.01*eye(5);Lx=zeros(5,1);Dx=[D;D(2)D(1:4);D(3)D(2)D(1:3);D(4:-1:2)D(1:2);D(5:-1:1)];Ycan=inv(G'*inv(B*Dx*B'+Do)*G)*G'*inv(B*Dx*B'+Do)*(L-B*Lx);Xcan=Lx+Dx*B'*inv(B*Dx*B'+Do)*(L-G*Ycan-B*Lx);R=[0.72250.72252.16753.61255.0575;2.16750.72250.72252.16753.6125;2.16750.72250.72252.1675;5.05753.61252.16750.72250.7225];Dx2x=0.1260*exp(-0.36*(R.^2));Lx2=zeros(4,1);X2can=Lx2+Dx2x*B'*inv(B*Dx*B'+Do)*(L-G*Ycan-B*Lx);F=[ones(4,1),R(:,1)]*Ycan+X2canYcan=0.325218914405570.98910879733786Xcan=0.23026995338496-0.58781954440652-0.29810342411034-0.004825597687310.20616191925934Dx2x=0.104413570950870.104413570950870.023219249557560.001148234628170.000012627120680.023219249557560.104413570950870.104413570950870.023219249557560.001148234628170.001148234628170.023219249557560.104413570950870.104413570950870.023219249557560.000012627120680.001148234628170.023219249557560.104413570950870.10441357095087X2can=-0.22209452936566-0.55898642725186-0.105187790070440.10820133195839三、实验心得和体会实验四平差系统的统计假设检验上机实验实验目的要求:通过此实验要求掌握简单的粗差剔出方法，并应用Matlab编程实现其方法，获得较为满意的平差系统前提。在测量中，按照误差的特点与性质，误差可分为：系统误差，粗差和随机误差。在假定不含有系统误差的情况下，可借助MATLAB对测量数据进行处理，使处理过程快速、结果可靠。实验内容：处理测量数据的过程如下：(1)按测量的先后顺序记录下每个测量值；(2)计算算术平均值；(3)计算残余误差；(4)校核算术平均值及残余误差；(5)判断是否有粗大误差，若有，剔除；(6)计算单次测量的标准差；(7)计算算术平均值的标准差：(8)计算算术平均值的极限误差；(9)列出测量结果。误差处理时常用的MATLAB函数函数名调用格式作用1absB=abs(a)求绝对值2sqrtB=sqrt(a)对向量中的值依次开平方3meanb=mean(a)求平均值4stdb=std(a)求标准差5cova=cov(x,y)求协方差6normrndW=normrnd(,,,)生成正态分布的向量7normstat[E,D]=(mu,sigma)计算正态分布的期望与方差8normfit[muhat,sigmut,muci,sigmaci]=normfit(X,Alpha)已知数据符合正态分布，对参数进行点估计和区间估计其算法流程图如下：开始开始输入数据计算平均值计算残余误差计算单次测量标准差判断是否含粗大误差计算算术平均值标准差S写出计算结果否剔除含粗大误差的数据具体内容：现对某被测量进行20次测量，得到测量序列x，其中第1个数为粗大误差，需运用莱以特准则将其剔除，再对数据进行分析计算，具体程序如下：closeallclearclcx=[28.005724.997424.996224.997024.985224.997725.001225.003125.014424.996525.006225.008025.009424.990125.002125.002424.989924.992625.010824.9987];%含有粗大误差的测量值序列aver=mean(x)%求该序