杀毒软件360
需积分: 0 67 浏览量
更新于2014-05-06
收藏 15.81MB RAR 举报
unction [Err,Ai,Bj,perO,perN,Enew]=balancing(A,E,Flag,J,bJ)
%计算测网焦点差值的平差值。
%输出:Ai-主测网平差值,应用(主测线-Ai);Bj-联络测网平差值,应用(联络测网+Bj)
% perO-平差前测网精度,perN-平差后测网精度
% Err-平差是否正确完成,‘0’-正确,其他-出错信息
%输入:A-测网焦点构造的矩阵,E-测网焦点的差异值(主-辅)
% Flag-平差算法选择,0-规则测网平差(要求交点都存在),1-不规则测网最小二乘平差
% (J,bJ)当Flag==0时无用,当Flage==1时,如果J=0表示最小方差平差,bJ无用
% 当0<J<=(M+N)时,bJ表示第J条测线(主[0<J<=M]或辅[M<J<=M+N])给定的调差值(一般bJ=0),以此调差其他测线值
Ai=0;Bj=0;perO=0;perN=0;
Err='0';
Enew=0;
Z=0;
[M,N]=size(E);
for i=1:M
Z=Z+A(i,i);
end
%%
%规则平差算法,要求测网必须完备(无虚拟焦点)
if Flag==0
if Z~=M*N
Err='非规则测网,不能使用规则测网算法平差';
fprintf(Err);
return;
end
W=sum(sum(E.^2));
perO=sqrt(W/Z/2);%平差前测网精度计算
b=calE(E);
Ai=zeros(1,M);
Bj=zeros(1,N);
for i=1:M
Ai(i)=b(i)/N;
for j=1:N
E(i,j)=E(i,j)-Ai(i);
end
end
b=calE(E);
for j=1:N
Bj(j)=b(j+M)/M;
for i=1:M
E(i,j)=E(i,j)-Bj(j);
end
end
W=sum(sum(E.^2));
perN=sqrt(W/Z/2);%平差后测网精度计算
%%
%不规则测网,依据最小方差原则平差
else
[M,N]=size(E);
W=sum(sum(E.^2));
perO=sqrt(W/Z/2);%平差前测网精度计算
if J==0
for i=1:M
A(N+M,i)=A(i,i);
end
for j=1:N
A(N+M,j+M)=-A(j+M,j+M);
end
b=calE(E);
b(N+M)=0;
if rank(A)<M+N
Err='测网被分为两个以上,不能计算';
return;
end
AB=A\b;
Ai=AB(1:M)';
Bj=AB(M+1:M+N)';
%%
%不规则测网,依据指定第条J测线调差值bJ平差注意主测线是减去该值,联络测线是加上该值(一般bJ=0)
else
if J<1||J>(M+N)
Err='测线号指定错误,超出测线范围。';
fprintf(Err);
return;
end
b=calE(E);
j=0;
b1=zeros(M+N-1,1);
for i=1:M+N
if i~=J
j=j+1;
b1(j)=b(i)-A(i,J)*bJ;
end
end
A(J,:)=[];
A(:,J)=[];
AB=A\b1;
if J<=M
Ai=zeros(1,M);
j=0;
Ai(J)=bJ;
for i=1:M-1
j=j+1;
if i==J
j=j+1;
end
Ai(j)=AB(i);
end
Bj=AB(M:M+N-1)';
else
Ai=AB(1:M)';
Bj=zeros(1,N);
j=0;
Bj(J-M)=bJ;
for i=M+1:M+N-1
j=j+1;
if i==J
j=j+1;
end
Bj(j)=AB(i);
end
end
end
for i=1:M
for j=1:N
if A(i,M+j)==1
E(i,j)=E(i,j)-Ai(i)-Bj(j);
end
end
end
W=sum(sum(E.^2));
perN=sqrt(W/Z/2);%平差后测网精度计算
end
Enew=E;
end
function b=calE(E)
%计算AM=b的b值,M-主测线条数;N-联络测线条数
[M,N]=size(E);
b=zeros(M+N,1);
for i=1:M
b(i)=sum(E(i,:));
end
for i=1:N
b(i+M)=sum(E(:,i));
end
end