DDKC.rar_导线平差_间接平差模型代码Matlab

clc clear %n=19,t=14;r=5; 10个角度，9个边长 p_miao=206265; S=xlsread('shuju','A1:A9'); %读入边长数据 S=S'; L_d=xlsread('shuju','B1:B10');%读入度，分，秒 L_d=L_d'; L_f=xlsread('shuju','C1:C10'); L_f=L_f'; L_m=xlsread('shuju.xlsx','D1:D10'); L_m=L_m'; L_hd=L_d/180*pi+L_f/60/180*pi+L_m/3600/180*pi; %角度化为弧度 %起始方位角 T0=175/180*pi+16/60/180*pi+35/3600/180*pi; %近似方位角 L_A=zeros(10,1); L_A(1)=T0+L_hd(1)-pi; for i=2:10 L_A(i)=L_A(i-1)+L_hd(i)-pi; if L_A(i)>2*pi L_A(i)=L_A(i)-2*pi; elseif L_A(i)<0 L_A(i)=L_A(i)+2*pi; end end %近似坐标值 x_j=zeros(10,1); y_j=zeros(10,1); x_j(1)=164366.5283;y_j(1)=122993.9747; x_j(10)=164366.5283;y_j(10)=122993.9747; test=S(1)*sin(L_A(1)); for i=2:9 x_j(i)=x_j(i-1)+S(i-1)*cos(L_A(i-1)); y_j(i)=y_j(i-1)+S(i-1)*sin(L_A(i-1)); end x_j(7)=164548.3383;y_j(7)=122978.9521; %近似边长 S0=zeros(9,1); for i=1:9 S0(i)=sqrt((x_j(i+1)-x_j(i))^2+(y_j(i+1)-y_j(i))^2); end %权阵 P=zeros(19,1); for i=1:19 if i<11 P(i)=1; else P(i)=4/S(i-10); end end P=diag(P); %边长改正数常数ls ls=zeros(9,1); for i=1:9 ls(i)=(S(i)-S0(i))*1000; end %角度观测近似值 L_0=zeros(10,1); L_0(1)=L_A(1)-T0-pi; if L_0(1)<0 L_0(1)=L_0(1)+2*pi; end for i=2:10 L_0(i)=L_A(i)-L_A(i-1)+pi; if(L_0(i)<0) L_0(i)=L_0(i)+2*pi; elseif(L_0(i)>2*pi) L_0(i)=L_0(i)-2*pi; end end %角度改正常数ll ll=zeros(10,1); for i=1:10 ll(i)=(L_hd(i)-L_0(i))*3600*180/pi; end %坐标方位角改正系数 B_Aa=zeros(10,1); B_Ab=zeros(10,1); for i=1:9 B_Aa(i)=p_miao*sin(L_A(i))/S0(i)/1000; B_Ab(i)=-p_miao*cos(L_A(i))/S0(i)/1000; end B_Aa(10)=p_miao*sin(L_A(10))/S0(1)/1000; B_Ab(10)=-p_miao*cos(L_A(10))/S0(1)/1000; %边长改正数系数 B_Sa=zeros(9,1); B_Sb=zeros(9,1); for i=1:9 B_Sa(i)=(x_j(i+1)-x_j(i))/S0(i); B_Sb(i)=(y_j(i+1)-y_j(i))/S0(i); end %误差方程系数 B=zeros(19,14); %角度 B(1,1)=-B_Aa(1);B(1,2)=-B_Ab(1); B(2,1)=-B_Aa(1)+B_Aa(2);B(2,2)=B_Ab(2)+B_Ab(2); B(2,3)=-B_Aa(2);B(2,4)=-B_Ab(2); for i=2:4 B(i+1,2*i-3)=-B_Aa(i); B(i+1,2*i-2)=-B_Ab(i); B(i+1,2*i-1)=B_Aa(i)+B_Aa(i+1);B(i+1,2*i)=B_Ab(i)+B_Ab(i+1); B(i+1,2*i+1)=-B_Aa(i+1);B(i+1,2*i+2)=-B_Ab(i+1); end B(6,7)=-B_Aa(5);B(6,8)=-B_Ab(5); B(6,9)=B_Aa(5)+B_Aa(6);B(6,10)=B_Ab(5)+B_Ab(6); B(7,9)=-B_Aa(6);B(7,10)=-B_Ab(6); B(7,11)=-B_Aa(7);B(7,12)=-B_Ab(7); B(8,11)=B_Aa(7)+B_Aa(8);B(8,12)=B_Ab(7)+B_Ab(8); B(8,13)=-B_Aa(8);B(8,14)=-B_Ab(8); B(9,11)=-B_Aa(8);B(9,12)=-B_Ab(8); B(9,13)=B_Aa(8)+B_Aa(9);B(9,14)=B_Ab(8)+B_Ab(9); B(10,13)=-B_Aa(9);B(10,14)=-B_Ab(9); B(10,1)=-B_Aa(10);B(10,2)=-B_Ab(10); %边长 B(11,1)=B_Sa(1);B(11,2)=B_Sb(2); for i=2:5 B(i+10,2*i-3)=-B_Sa(i);B(i+10,2*i-2)=-B_Sb(i); B(i+10,2*i-1)=B_Sa(i);B(i+10,2*i)=B_Sb(i); end B(16,9)=-B_Sa(6);B(16,10)=-B_Sb(6); B(17,11)=B_Sa(7);B(17,12)=B_Sb(7); B(18,11)=-B_Sa(8);B(18,12)=-B_Sb(8); B(18,13)=B_Sa(8);B(18,14)=B_Sb(8); B(19,13)=B_Sa(1);B(19,14)=B_Sb(2); %求l阵 l=[ll;ls]; %求NBB的逆矩阵 NBB_1=inv(B'*P*B); W=B'*P*l; %参数平差值 xy=NBB_1*W; XY=zeros(14,1); for i=1:5 XY(2*i-1)=x_j(i+1); XY(2*i)=y_j(i+1); end XY(11)=x_j(8);XY(12)=y_j(8); XY(13)=x_j(9);XY(14)=y_j(9); XY=XY+xy/1000; %待定点的协因数阵 QXX=zeros(7,1); QYY=zeros(7,1); QXY=zeros(7,1); V=B*(xy)-l; cgema0=sqrt((V'*P*V)/5); %单位权中误差 cgema=zeros(7,1); %各待定点的中误差 for i=1:7 QXX(i)=NBB_1(2*i-1,2*i-1); QYY(i)=NBB_1(2*i,2*i); QXY(i)=NBB_1(2*i-1,2*i); cgema(i)=cgema0*sqrt(QXX(i)+QYY(i)); end % % %10个待定点误差椭圆的绘制 % K=zeros(10,1); % QEE=zeros(10,1); % QFF=zeros(10,1); % E=zeros(10,1); % F=zeros(10,1); % faiE=zeros(10,1); % % XY_line=[x_j(1);y_j(1);XY; % x_j(12);y_j(12);x_j(1);y_j(1);]; % x_line=zeros(13,1); % y_line=zeros(13,1); % for i=1:13 % x_line(i)=XY_line(2*i-1); % y_line(i)=XY_line(2*i); % end % plot(y_line,x_line); % hold on; % % for i=1:10 % K(i)=sqrt((QXX(i)-QYY(i))^2+4*(QXY(i)^2)); % QEE(i)=0.5*(QXX(i)+QYY(i)+K(i)); % QFF(i)=0.5*(QXX(i)+QYY(i)-K(i)); % E(i)=cgema0*sqrt(QEE(i))/1000; % F(i)=cgema0*sqrt(QFF(i))/1000; % faiE(i)=atand((QEE(i)-QXX(i))/QXY(i)); % if(faiE(i)<0) % faiE(i)=faiE(i)+180; % end % t=0:0.01:2*pi; % x1=170*E(i)*sin(t); %为了显示，乘以比例系数170 % y1=170*F(i)*cos(t); % dei=pi/2-faiE(i)/180*pi; % x=XY(2*i)+cos(dei)*x1-sin(dei)*y1; % y=XY(2*i-1)+sin(dei)*x1+cos(dei)*y1; % plot(x,y); % hold on; % xlabel('Y'); % ylabel('X') % title('误差椭圆') % end % % %数据储存 % xlswrite('V',V); % xlswrite('L_A',L_A); % xlswrite('x_j',x_j); % xlswrite('y_j',y_j); % xlswrite('S0',S0); % xlswrite('P',P); % xlswrite('l',l); % B_A=[B_Aa,B_Ab]; % xlswrite('B_A',B_A); % B_S=[B_Sa,B_Sb]; % xlswrite('B_S',B_S); % xlswrite('B',B); % xlswrite('NBB_1',NBB_1); % xlswrite('XY',XY); % xlswrite('cgema',cgema); % xlswrite('xy_',xy); % xlswrite('E',E); % xlswrite('F',F); % xlswrite('faiE',faiE); % %