测量平差实习报告

课程编号：课程性质：选修平差实习学院：测绘学院专业：测绘工程班级：2011级03班学号：2011301610***姓名：******指导教师：*******目录HYPERLINK一、题目内容 2HYPERLINK二、函数模型 3HYPERLINK三、计算框图、程序说明及源程序HYPERLINK1、计算框图： 3HYPERLINK2、程序说明 3HYPERLINK3、源程序 8HYPERLINK三、计算过程 29HYPERLINK四、计算结果 29HYPERLINK五、结果检核HYPERLINK第一题 31HYPERLINK第二题 32HYPERLINK六、实习心得 34一、题目内容第一题：图1为一水准网，A、B、C为已知高程点，P1~P9为9个待定高程点，第一次观测了图中1~15条水准路线的高差，各水准路线的观测高差、距离及已知点高程均列于表4，试求（1）1公里高差平差值中误差；（2）各待定点高程平差值及中误差；如果准备加测5段高差16~20（图中用虚线表示），（3）试估算平差后各待定点的精度；（4）精度最弱的点发生改变没有；（5）哪一点精度提高得最多。图1表1水准网观测数据线路号观测高差(m)距离(km)线路号观测高差(m)距离(km)11.49570.4112.04160.821.78920.5121.53550.735.20620.7131.74500.842.08530.9142.74340.453.41750.4151.37850.561.20801.0160.676.03730.6170.684.33070.6181.893.65100.6190.8102.51200.3201.2已知点高程：、、、第二题：见平差辅助软件函数模型第一题：本题采用间接平差模型即，误差方程为。第二题：本题采用间接平差模型即，误差方程为。计算框图、程序说明及源程序1、计算框图：计算B计算B计算得出计算P(P(i))=1/s(i)计算得出计算P(P(i))=1/s(i)计算l(l=L-L0)计算l(l=L-L0)计算得出计算得出2、程序说明这两个程序由matlab写成。第一题可以适用于任何闭合水准平差，先要把数据写进一个txt文件，格式为：已知点个数未知点数观测值个数增加观测值个数已知点点号已知点高程起点点号终点点号高程距离……注意已知点点号变为数字，按字母顺序接在未知点点号之后。然后运行程序，程序会弹出一个窗口，选择刚刚写的txt文件，程序就可算出1公里高程平差值中误差（q0），各待定点高程平差值及中误差，以及加测5段高差后各待定点精度，精度最弱点的点号，各待定点精度的变化。1水准网数据结构的设计水准网根据其网形可分为附合、闭合和支水准。其中应用最广的为附合和闭合水准网。不管何种网形，一条水准路线包括四个数据:起始点、高差、距离、终点，且H起+h高差=H终，因此在进行数据结构的设计中，必须体现出终点、起点从而判断其与高差关系。在本程序设计中，利用qd［］，zd［］，h［］，jl［］，分别表示起点点号、终点点号、高差、距离所组成的矩阵。2.1水准网所有点高程确定在水准网间接平差计算过程中，按照已知点的个数确定参数个数，有已知点时，参数个数等于未知点个数，没有已知点时参数个数等于未知点个数减1，但是在进行方程列立时，为了计算方便首先要求得参数X的近似值，但是对于用程序实现其自动求取近似值比较麻烦，为了Matlab实现运算简单，因此在本程序中，所有未知点的高程起始值均为0，通过以下代码实现。HHH=HH;%HH为读取文件中已知点高程矩阵，HHH为未知点高程%zd总点数，wz未知点数，fori=1:bHHH(zd－wz+i)=0;%给未知高程点添加高程零生成HHH，b为未知点个数2.2水准网误差方程系数矩阵生成对于水准网间接平差而言，最重要的计算步骤不外乎误差方程的列立，而对于程序设计重点和难点就是自动生成误差方程系数矩阵B和常数项L。在本程序的设计过程中，利用起始点、终点、高差。距离之间的对应关系，利用循环计算常数项L，由于水准网比较特殊，系数阵B也特殊基本上是0、1、－1组成的矩阵，根据这一特点利用对应关系确定其系数阵B。%qd，zd存放的是与各边相对应的点的点号，通过自定义函数实现，方便确定B与L。fori=1:se=qh(i);f=zh(i);L(i)=h(i)-(HHH(f)-HHH(e));endL=L';%常数项Lfori=1:se=qh(i);f=zh(i);if(e<=wd)B(i,e)=-1;endif(f<=wd)B(i,f)=1;endend2.3水准网平差计算与结果输出1)对于水准网平差计算，相对Matlab而言。仅需几行代码可以实现P=diag(1./s);%解权阵x=inv(B’*P*B)*(B’*P*l);%高程改正数x，B系数阵，l为常数项阵v=B*x－l;%观测值改正数q0=sqrt(v’*P*v/(s-wd));Qxx=inv(B’*p*B);q=q0*sqrt(diag(Qxx));%待定点高程平差值中误差第二题也适用于任何已知数据为已知点坐标（其中有两个已知点相邻）的闭合导线网平差，把题目软件输出的观测值表用word处理一下，删掉汉字，前面加上几个数据，分别是边观测值个数，角观测值个数，已知点个数，点总个数，去掉角度编号前面的角标志和后面角度值的度分秒标志，在文本的末尾加上已知点点号和已知点坐标。另外，要把所有0点号改成最大点号再加上1。程序先计算各个待定点坐标的估计值。由观测角出发，按观测顺序检索，若第i个角的末点的坐标未知，则判断起点和中点坐标是否已知，若已知，则根据这两个点的坐标算出后视方向方位角，再结合前视方向边长观测值（根据点号给每个角的两条边按边长观测矩阵顺序编个号，方便引用）及第i个观测角的大小，推算出末点的坐标。再用同样的方法算出未知起点的坐标。然后跳到下一个角进行计算，如此完成一遍计算，制作一个循环，重复刚才的运算直到所有未知点的坐标都估算出来。然后程序根据刚算出的未知点坐标估计值和已知点坐标算出所有边和角的估算值（边和角的各个点可由点号矩阵检索到，然后再结合标矩阵得到各个坐标），以此为基础，计算系数矩阵B。先默认B为0，再针对每个角，结合导线网间接平差的误差方程，求出每个观测值改正数用参数式时参数前面的系数，然后组成Ｂ矩阵。因为对于某些观测边来说，它们的端点为已知点，已知点的坐标不是参数，因此还要把刚生成的矩阵中去掉几列，已知点点号是多少，就去掉第几列，（注意要按照点号大小，从大往小删除。）在求B矩阵时要注意要把估计值中的角化成角度值，边化成以mm为单位。然后算出矩阵ｌ．对于权阵，角度的权全部定为１，边的权定为１４４／边观测值。算出Ｂ,ｌ，Ｐ矩阵后，就可以算出坐标改正数ｘ，观测值改正数ｖ，以及单位权中误差了，代码如下：x=inv(B'*P*B)*B'*P*l;v=B*x-l;%v是改正数，角以秒为单位，边以mm位单位u0=sqrt(v'*P*v/((zg-yg)*2));%单位权中误差注意坐标改正数和边长改正数以mm为单位，而角度改正数以秒为单位。求出改正数后，即可求出坐标平差值和观测值平差值，继而还可求出坐标平差值中误差及各观测值平差值中误差。对平差结果进行整合时，把坐标平差值带入，反算出每个边长和角的大小，与观测值平差值对照，若能对应上，则说明结果符合要求。3、源程序第一题：[f_name1,filepath]=uigetfile('*.txt','pickafileforread');fid=fopen(strcat(filepath,f_name1),'rt');yd=fscanf(fid,'%f',1);%已知点个数wd=fscanf(fid,'%f',1);%未知点数s=fscanf(fid,'%f',1);%观测值个数s2=fscanf(fid,'%f',1);%增加观测值个数yh=fscanf(fid,'%f',yd);%已知点点号HH=fscanf(fid,'%f',yd);%已知点高程jz=fscanf(fid,'%f',[4,s+s2]);%第四行到末尾的观测数矩阵jz=jz';qh=jz(:,1);%起点点号zh=jz(:,2);%终点点号h=jz(:,3);%高差jl=jz(:,4);%距离fclose('all');fori=(wd+1):(wd+yd)HHH(i)=HH(i-wd);endfori=1:wdHHH(i)=0;%给未知高程点添加高程零生成HHHendfori=1:se=qh(i);f=zh(i);l(i)=h(i)-(HHH(f)-HHH(e));endl=l';%常数项lfori=1:se=qh(i);f=zh(i);if(e<=wd)B(i,e)=-1;endif(f<=wd)B(i,f)=1;endendjl1=jl(1:s,:);%未算上第二题加上的几个测段P=diag(1./jl1);xg=inv(B'*P*B)*B'*P*l;%xg是参数改正数;V=B*xg-l;q0=sqrt(V'*P*V/(s-wd));disp('****************各待测点高程平差值为***************');XP=xgdisp('**************************************************');%q0为单位权中误差，xg为参数改正值，由于参数的估计值为0，所以各参数的平差值，即各待定点高程的平差值等于xg%第二小题Qxx=inv(B'*P*B);%协因数阵disp('*****************各待定点高程平差值中误差为：*********');q=q0*sqrt(diag(Qxx))disp('**************************************************');%各待定点的中误差为：k=q(1);fori=1:wdifq(i)>kk=q(i);m=i;%m为精度最弱点的点号endend%寻找精度最弱的点disp('****************精度最弱的点点号为：********************');disp(m);disp('******************************************************');%求出加测段后的B矩阵fori=s+1:s+s2e=qh(i);f=zh(i);if(e<=wd)B(i,e)=-1;endif(f<=wd)B(i,f)=1;endendP2=diag(1./jl);%加测段后的权阵PQxx2=inv(B'*P2*B);%加测段后的协方差disp('*****************平差后各待定点的中误差为：*****************');q2=q0*sqrt(diag(Qxx2))disp('**************************************************');disp('****************精度最弱的点点号为：********************');disp(m);disp('**************************************************');tisheng=q-q2;iftisheng(m)>0disp('最弱点的精度提升了');elsedisp('最弱点的精度没有提升');endk=tisheng(1);fori=1:wdiftisheng(i)>tisheng(1)k=tisheng(i);n=i;endenddisp('*******************精度提升最高的点点号为*********************');disp(n);disp('**************************************************');%结果检核%求出根据各点平差高程算出的各水准路线的高差，与观测高差的平差值比较，若符合，则说明结果符合要求fori=wd+1:wd+ydXP(i)=HH(i-wd);end%加入已知点高程fori=1:se=qh(i);f=zh(i);hy(i)=XP(f)-XP(e);endhy=hy';fori=1:shp(i)=h(i)+V(i);endhp=hp';chabie=hy-hp;第二题：[f_name1,filepath]=uigetfile('*.txt','pickafileforread');fid=fopen(strcat(filepath,f_name1),'rt');bg=fscanf(fid,'%f',1);%边长观测值个数jg=fscanf(fid,'%f',1);%角度观测值个数yg=fscanf(fid,'%f',1);%已知点个数zg=fscanf(fid,'%f',1);%点总数wg=zg-yg;%未知点个数BIAN=fscanf(fid,'%f',[3,bg]);%读入所有边长BIAN=BIAN';%读入后发生了转置，得重新转过来q=BIAN(:,1);%边起点点号m=BIAN(:,2);%边末点点号sg=BIAN(:,3);%边长观测值JIAO=fscanf(fid,'%f',[6,jg]);%读入所有角度JIAO=JIAO';Q=JIAO(:,1);%角度起点点号Z=JIAO(:,2);%角度中点点号M=JIAO(:,3);%角度末点点号DU=JIAO(:,4);%角度观测值度FEN=JIAO(:,5);%角度观测值分MIAO=JIAO(:,6);%角度观测值秒HD=(DU+FEN/60+MIAO/3600)/180*pi;%角度值转换为弧度值在matlab里以弧度计算W=zeros(zg,1);%生成矩阵X=W(:,1);Y=X;%定义两个向量X,Y作为点坐标估计值yh=fscanf(fid,'%f',yg);%已知点点号fori=1:ygX(yh(i))=fscanf(fid,'%f',1);Y(yh(i))=fscanf(fid,'%f',1);end%if条件句或循环语句要以end结尾%求未知点坐标估计值,fori=1:jgforj=1:bgif((q(j)==Z(i)&m(j)==M(i))|(q(j)==M(i)&m(j)==Z(i)))qbh(i)=j;break;endendforj=1:bgif((q(j)==Q(i)&m(j)==Z(i))|(q(j)==Z(i)&m(j)==Q(i)))hbh(i)=j;break;endendend%找到所有角对应的前边编号k=0;%注意，只算角末点坐标不够，因为对有一些图来说，已知点（两相邻的）点号靠后，很多点没法计算，得再加上计算起点坐标while(k<wg)%开始时判断语句为k<zg,死循环！fori=1:jgif(X(M(i))==0&Y(M(i))==0)if((X(Z(i))~=0|Y(Z(i))~=0)&(X(Q(i))~=0|Y(Q(i))~=0))%下标不能为0所以将所有0点号变成zg+1fj=atan((Y(Z(i))-Y(Q(i)))/(X(Z(i))-X(Q(i))));if(X(Z(i))>X(Q(i)))fj=fj+pi;end%fj为方位角k=k+1;X(M(i))=X(Z(i))+sg(qbh(i))*cos(fj+HD(i));Y(M(i))=Y(Z(i))+sg(qbh(i))*sin(fj+HD(i));endendif(X(Q(i))==0&Y(Q(i))==0)if((X(Z(i))~=0|Y(Z(i))~=0)&(X(M(i))~=0|Y(M(i))~=0))%下标不能为0所以将所有0点号变成zg+1fj=atan((Y(Z(i))-Y(M(i)))/(X(Z(i))-X(M(i))));if(X(Z(i))>X(M(i)))fj=fj+pi;end%fj为方位角k=k+1;X(Q(i))=X(Z(i))+sg(hbh(i))*cos(fj+2*pi-HD(i));Y(Q(i))=Y(Z(i))+sg(hbh(i))*sin(fj+2*pi-HD(i));endendendend%求边估计值fori=1:bgsj(i)=sqrt((X(m(i))-X(q(i)))^2+(Y(m(i))-Y(q(i)))^2);%sj为边估计值（由坐标估计值推出）endsj=sj';%求角估计值fori=1:jgfj1=atan((Y(Z(i))-Y(Q(i)))/(X(Z(i))-X(Q(i))));if(X(Z(i))>X(Q(i)))fj1=fj1+pi;end%后视方位角fj2=atan((Y(Z(i))-Y(M(i)))/(X(Z(i))-X(M(i))));if(X(Z(i))>X(M(i)))fj2=fj2+pi;end%前视方位角jj(i)=fj2-fj1;if(jj(i)<0)jj(i)=jj(i)+2*pi;%保证所求夹角为我们需要的角，fj1和fj2的范围为-90度到270度，若fj2>fj1,则fj2-fj1就为所求，否则要加360度endendjj=jj';%计算系数矩阵B%计算角部分fori=1:jgB1(i,Q(i))=(Y(Q(i))-Y(Z(i)))/(sj(hbh(i)))^2;B1(i,M(i))=(Y(M(i))-Y(Z(i)))/sj(qbh(i))^2;B1(i,Z(i))=B1(i,M(i))-B1(i,Q(i));B1(i,Q(i)+zg)=-(X(Q(i))-X(Z(i)))/sj(hbh(i))^2;B1(i,M(i)+zg)=(X(M(i))-X(Z(i)))/sj(qbh(i))^2;B1(i,Z(i)+zg)=-(B1(i,Q(i)+zg)+B1(i,M(i)+zg));end%计算边部分fori=1:bgB2(i,m(i))=(X(m(i))-X(q(i)))/sj(i);B2(i,q(i))=-B2(i,m(i));B2(i,m(i)+zg)=(Y(m(i))-Y(q(i)))/sj(i);B2(i,q(i)+zg)=-B2(i,m(i)+zg);endB=[B1*180*3600/pi/1000;B2];%注意此处B1要化为角度值，以秒为单位，另外要注意统一边长单位，后面都是毫米级，所以要除以1000（*1000/1000^2）fori=1:ygB(:,yh(yg+1-i)+zg)=[];%去掉已知点坐标对应的系数，但注意删除时要从后往前endfori=1:ygB(:,yh(yg+1-i))=[];end%计算矩阵ll1=(HD-jj)*180*3600/pi;%，角度，以秒为单位l2=(sg-sj)*1000;%边，注意此处要以mm为单位l=[l1;l2];%计算权阵Pfori=1:jgP(i,i)=1;endfori=(jg+1):(jg+bg)P(i,i)=144/sg(i-jg);%定权方法见课本141页endx=inv(B'*P*B)*B'*P*l;v=B*x-l;%v是改正数，角以秒为单位，边以mm位单位u0=sqrt(v'*P*v/((zg-yg)*2));%单位权中误差%求角平差值，此处只需把秒部分改正DU1=DU;FEN1=FEN;fori=1:jgMIAO1(i)=MIAO(i)+v(i);MIAO1=MIAO1';if(MIAO1(i)>=60)MIAO1(i)=MIAO1(i)-60;FEN1(i)=FEN(i)+1;ifFEN1(i)==60FEN1(i)=0;DU1(i)=DU(i)+1;endendif(MIAO1(i)<0)MIAO1(i)=MIAO1(i)+60;FEN1(i)=FEN(i)-1;ifFEN1(i)<0FEN1(i)=FEN1(i)+60;DU1(i)=DU(i)-1;endendenddisp('角度平差值为');jp=[DU1FEN1MIAO1]%求边平差值fori=1:bgsp(i)=sg(i)+v(i+jg)/1000;%注意单位enddisp('边平差值为：');sp=sp'%sp为边平差值%求坐标平差值%由于坐标阵X,Y中包括已知点坐标，因此坐标改正数x必须在插上几个0元素%要从前往后插fori=1:wg*2x1(i)=x(i);endfori=1:yg%先把要插入元素位置后面的元素全部后移一位k=wg*2+i-1;%末尾元素位置forj=yh(i):ke=yh(i)+k;x1(e-j+1)=x1(e-j);%倒着移endx1(yh(i))=0;endfori=1:yg%先把要插入元素位置后面的元素全部后移一位k=wg*2+yg+i-1;%末尾元素位置forj=yh(i)+zg:ke=yh(i)+zg+k;x1(e-j+1)=x1(e-j);%倒着移endx1(yh(i)+zg)=0;endfori=1:zgXP(i)=X(i)+x1(i)/1000;%注意单位YP(i)=Y(i)+x1(i+zg)/1000;%注意单位enddisp('待定点坐标平差值为：（此处把已知点包含了进来，直接忽略掉）');XP=XP'YP=YP'disp('各待定点坐标平差值中误差为:(此处还是按点号顺序，但去掉了已知点)');disp('(前面一半是X,后面一半是Y)');%求待定点坐标平差值中误差Qxx=inv(B'*P*B);Qxd=diag(Qxx);uz=u0*sqrt(Qxd)%还是按点号顺序，只是除去了已知点%求各待测点点点位中误差k=u0*sqrt(Qxd(1)+Qxd(1+wg));fori=1:wgup(i)=u0*sqrt(Qxd(i)+Qxd(i+wg));ifup(i)>k;k=up(i);w2=i;%精度最弱点点号endendup=up';%求各观测值平差值中误差QLP=B*inv(B'*P*B)*B';%观测值协因数阵u=u0*sqrt(diag(QLP));%各观测值中误差mu=u(jg+1)/sg(1);%第一个边边长相对中误差fori=jg+1:jg+bgxu(i-jg)=u(i)/sg(i-jg);%各边长相对中误差ifxu(i-jg)>mumu=xu(i-jg);%精度最弱边相对中误差w=i-jg;%精度最弱边边号endenddisp('精度最弱边边号为：');disp(w);disp('精度最弱边相对中误差为：');disp(mu);%检验结果是否闭合%求边检验值fori=1:bgsy(i)=sqrt((XP(m(i))-XP(q(i)))^2+(YP(m(i))-YP(q(i)))^2);%sj为边估计值（由坐标估计值推出）endsy=sy';%求角检验值fori=1:jgfj1=atan((YP(Z(i))-YP(Q(i)))/(XP(Z(i))-XP(Q(i))));if(XP(Z(i))>XP(Q(i)))fj1=fj1+pi;end%后视方位角fj2=atan((YP(Z(i))-YP(M(i)))/(XP(Z(i))-XP(M(i))));if(XP(Z(i))>XP(M(i)))fj2=fj2+pi;end%前视方位角jy(i)=fj2-fj1;if(jy(i)<0)jy(i)=jy(i)+2*pi;%保证所求夹角为我们需要的角，fj1和fj2的范围为-90度到270度，若fj2>fj1,则fj2-fj1就为所求，否则要加360度endendjy=jy';bchabie=sy-sp;HDY=(DU1+FEN1/60+MIAO1/3600)/180*pi;%角度值转换为弧度值在matlab里以弧度计算jchabie=HDY-jy;%求极值Qxx=inv(B'*P*B);fori=1:wgK(i)=sqrt((Qxx(i,i)-Qxx(i+wg,i+wg))^2+4*Qxx(i,i+wg)^2);E(i)=sqrt(u0^2/2*(Qxx(i,i)+Qxx(i+wg,i+wg)+K(i)));F(i)=sqrt(u0^2/2*(Qxx(i,i)+Qxx(i+wg,i+wg)-K(i)));JF(i)=atan(2*Qxx(i,i+wg)/(Qxx(i,i)-Qxx(i+wg,i+wg)))/2;endE=E'/1000%极大值F=F'/1000%极小值JF=JF'/pi*180%极值方向disp('精度最弱点的误差椭圆的E,F,极值方向（弧度制）分别为');disp(E(w2));disp(F(w2));disp(JF(w2));x=x/1000%以m为单位fclose('all');三、计算过程直接运行程序即可得到结果。计算结果第一题计算结果如下：改正后X的值(高程平差值)m:X=高程平差值中误差δx20=（mm）加测后各点Qxx最弱点n=7146.35870.67080.599148.1490.67540.499最弱点精度Dmax=1.2294145.52860.76060.640146.0640.55280.465148.65010.81390.892

精度最弱点没有改变144.99911.07741.776142.48691.10731.779145.90660.95141.391精度提高最多的点m=2144.52830.89491.3741km高差平差值中误差=δ0／√3.6=1.2422mm第二题计算结果如下：点号坐标平差值坐标中误差边编号平差值x（m）y（m）x(mm)y(mm)0--142.3350308.0558191.44180.03530.05351--259.7501299.1832232.83750.02370.02732--356.3402265.6599282.29810.0260.03173--422.2703209.5504287.3980.01760.0554--551.9985147.2956274.2350.01420.06795--650.9106120.6149230.87550.01770.05746--729.554798.3522211.4370.00680.05557--872.4368122.76143.2370.0210.0218--910.2099129.4086135.48850.01890.03389--1056.21810168.668695.24980.01170.048410--1141.49213294.0068164.51780.01320.047411--1278.05014179.8645129.17810.0150.04612--1335.60915177.2636190.65770.01850.049813--030.36816189.6373234.37250.02790.061710--1435.72717141.9249211.78050.03310.070314--1561.534点位精度最低点误差椭圆参数

E=0.557F=0.006极大值方向正切值为3.487915--1645.43216--356.64115--1741.17017--628.613观测角角度观测平差值观测角角度观测平差值。′″。′″∠0122020556.6∠89101833830.1∠10131401940.3∠910112220757.1∠123230414.5∠910142972459∠2341473834.2∠1011122150351.6∠23162543737∠11101475172∠3452482524.1∠1112132091955.7∠4316106592.7∠121301813325.7∠456207257∠1014152004013.8∠5671624046.2∠141516161456.6∠56172594934.7∠1415172363954.1∠6782483234.97617970543∠161517745331.3∠7892011243.79五、结果检核第一题源程序中以下语句为对结果的检核：%结果检核%求出根据各点平差高程算出的各水准路线的高差，与观测高差的平差值比较，若符合，则说明结果符合要求fori=wd+1:wd+ydXP(i)=HH(i-wd);end%加入已知点高程fori=1:se=qh(i);f=zh(i);hy(i)=XP(f)-XP(e);endhy=hy';fori=1:shp(i)=h(i)+V(i);endhp=hp';chabie=hy-hp;若chabie为0或很小，则说明结果符合要求，经检验得结果符合要求。第二题源程序中以下语句为结果的检核：%求边检验值fori=1:bgsy(i)=sqrt((XP(m(i))-XP(q(i)))^2+(YP(m(i)