发布时间 : 星期三 文章通用水准网形秩亏自由网和拟稳平差程序设计更新完毕开始阅读
MATLAB设计任何网形的秩亏自由网平差和拟稳平差程序设计 长安大学王省超 2015.5
程序介绍：程序适合于任何网形的水准网平差，原始数据输入到连个excel表格 程序界面：
原始数据录入表格： （1）DH表
（2）GXLB表
程序代码：
function xsz=XSZ(num1,num2)%函数功能提取误差系数A [m1,n1]=size(num1); [m2,n2]=size(num2); n=0;
for i=1:m1 %用来判断参数个数 if num1(i,2)==1 n=n+0; else
n=n+1; end end
xsz=zeros(m2,n); %建立系数阵，全为零 for i=1:m2 %提取系数阵 q=num2(i,1); z=num2(i,2); xsz(i,z)=1; xsz(i,q)=-1; end end
%---------------常数项L-------------------------------------------------------------- function l=L(num1,num2) [m1,n1]=size(num1); [m2,n2]=size(num2); l=zeros(m2,1);
for i=1:m2 %计算l q=num2(i,1); z=num2(i,2);
l(i,1)=num2(i,4)-num1(z,4)+num1(q,4); end
l=l*1000; %把l从米换算为毫米 end
%-------------求平差权阵P------------------------------------------------------------ function p=P1(num1,num2) [m1,n1]=size(num1); [m2,n2]=size(num2); p=zeros(m2,m2); for i=1:m2
p(i,i)=1/num2(i,5); end
end
%------------秩亏自由网平差--------------------------------------------------------------------- function[v,Qxx]=PTPC clear all clc
global H1; global K1;
[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(num1,num2) %误差方程常数项单位mm P=P1(num1,num2); %权
N=A'*P*A %法方程系数 W=A'*P*l NN=N*N
N0=NN(1:r,1:r);
NN_=blkdiag(inv(N0),zeros(d,d));%广义逆求逆 Nm=N*NN_ x=Nm*W
%-----------高程平差值------
H=x/1000+num1(:,4)%单位统一为M H1=H';
%----精度评定-------
v=A*x-l %单位mm
Qxx=N*NN_*N*NN_*N
msgbox('秩亏自由网普通平差完成')
%---------------拟稳平差--------------------------------------------------------------------------------------- function[v,Qxx]= NWPC clear all global H2; global K1;
[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=m1-r; %秩亏数
l=L(num1,num2) %误差方程常数项单位mm P=P1(num1,num2); %权
for i=1:m1 %找出稳定点 if num1(i,7)==1 f(i)=1; else
f(i)=0; end end
c1=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=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_=Mm_*a';
b_=inv(N11)*(A1'-N12*a_);
%-------结算未知数------------------------------------------------------------- x2=a_*P*l; x1=b_*P*l; x=[x1;x2]
H=x/1000+num1(:,4)%单位统一为M H2=H';
%------计算改正数---------------------------------------------------------------------------