资源描述:
《通用水准网形秩亏自由网和拟稳平差程序设计》由会员上传分享,免费在线阅读,更多相关内容在教育资源-天天文库。
1、MATLAB设计任何网形的秩亏自由网平差和拟稳平差程序设计长安大学王省超2015.5程序介绍:程序适合于任何网形的水准网平差,原始数据输入到连个excel表格程序界面:原始数据录入表格:(1)DH表(2)GXLB表程序代码:functionxsz=XSZ(num1,num2)%函数功能提取误差系数A[m1,n1]=size(num1);[m2,n2]=size(num2);n=0;fori=1:m1%用来判断参数个数ifnum1(i,2)==1n=n+0;elsen=n+1;endendxsz=zeros(m2,n);%建立系数阵,全为零
2、fori=1:m2%提取系数阵q=num2(i,1);z=num2(i,2);xsz(i,z)=1;xsz(i,q)=-1;endend%---------------常数项L--------------------------------------------------------------functionl=L(num1,num2)[m1,n1]=size(num1);[m2,n2]=size(num2);l=zeros(m2,1);fori=1:m2%计算lq=num2(i,1);z=num2(i,2);l(i,1)=num2
3、(i,4)-num1(z,4)+num1(q,4);endl=l*1000;%把l从米换算为毫米end%-------------求平差权阵P------------------------------------------------------------functionp=P1(num1,num2)[m1,n1]=size(num1);[m2,n2]=size(num2);p=zeros(m2,m2);fori=1:m2p(i,i)=1/num2(i,5);endend%------------秩亏自由网平差-----------
4、----------------------------------------------------------function[v,Qxx]=PTPCclearallclcglobalH1;globalK1;[num1]=xlsread('DH');K=num1(:,4);K1=K';[m1,n1]=size(num1);%用来判断num1得行列[num2]=xlsread('GXLB');[m2,n2]=size(num2);A=XSZ(num1,num2)%误差方程系数r=rank(A);%矩阵的秩d=1;%秩亏数l=L(num
5、1,num2)%误差方程常数项单位mmP=P1(num1,num2);%权N=A'*P*A%法方程系数W=A'*P*lNN=N*NN0=NN(1:r,1:r);NN_=blkdiag(inv(N0),zeros(d,d));%广义逆求逆Nm=N*NN_x=Nm*W%-----------高程平差值------H=x/1000+num1(:,4)%单位统一为MH1=H';%----精度评定-------v=A*x-l%单位mmQxx=N*NN_*N*NN_*Nmsgbox('秩亏自由网普通平差完成')%---------------拟稳平差
6、---------------------------------------------------------------------------------------function[v,Qxx]=NWPCclearallglobalH2;globalK1;[num1]=xlsread('DH');K=num1(:,4);K1=K';;[m1,n1]=size(num1);%用来判断num1得行列[num2]=xlsread('GXLB');[m2,n2]=size(num2);A=XSZ(num1,num2)%误差方程系数r=r
7、ank(A);%矩阵的秩d=m1-r;%秩亏数l=L(num1,num2)%误差方程常数项单位mmP=P1(num1,num2);%权fori=1:m1%找出稳定点ifnum1(i,7)==1f(i)=1;elsef(i)=0;endendc1=find(f==1);%找出稳定点的下标c0=find(f==0);%非稳定点下标A1=A(:,c0)A2=A(:,c1)N11=A1'*P*A1;N12=A1'*P*A2;N21=N12';N22=A2'*P*A2;M=N22-N21*inv(N11)*N12;r=rank(M);d=1;MM=
8、M*M;M0=MM(1:r,1:r);MM_=blkdiag(inv(M0),zeros(d,d));%广义逆求逆Mm_=M*MM_;a=(A2'-N21*inv(N11)*A1')';a_=