测绘 matlab 编程 程序

导线网严密平差 完整程序代码及注释

已知数据

已知点坐标：

X

A

=730.024 ，

Y

A

=126.040 ，

111.855=X

B

，232.172=Y B

应得结果 = 5.30

表七：

实际结果

点号X(m) Y(m) 中误差X(mm) 中误差Y(mm) 点位中误差

1 678.1641 287.3411 3.3868 5.8737 7.6078

2 564.6919 446.8827 6.0102 7.5137 9.1937

3 533.8468 629.2749 8.7985 7.3927 10.0596

4 475.6776 535.8728 7.3949 8.4987 9.9667

5 499.1408 374.4231 5.807

6 7.6840 9.1827

6 594.6190 176.7976 5.2308 3.4068 7.3475

7 826.3848 305.7569 2.7439 5.4473 7.1551

8 767.8778 401.8510 4.8371 6.3810 8.3734

9 745.7230 531.0410 6.8992 6.6276 9.1947

10 671.6256 650.8637 9.0695 6.0533 9.7220

11 794.2833 600.4688 8.4994 6.1800 9.5784

12 859.0528 511.5349 6.7399 6.6874 9.1608

13 898.2695 399.0763 4.7580 6.4207 8.3586

14 919.3459 276.9409 3.1903 4.8006 7.0670

单位权中误差： 5.3393

详细代码

% 贡XX导线严密平差程序

clear;clc;

% 从文件中输入数据

% S1 观测角度数据

% 【角号度分秒左点角点右点】-1为已知点B，0为已知点A fid1=fopen('S1.txt','rt');

[S1,~]=fscanf(fid1,'%f ',[7 24]);

fclose(fid1);

S1=S1';

% S2 边长观测值

% 【边号距离测边中误差对应角号】

fid2=fopen('S2.txt','rt');

[S2,~]=fscanf(fid2,'%f ',[5 17]);

fclose(fid2);

S2=S2';

% S3：度分秒到角度

[S3]=dfm2jd(S1);

% 计算坐标近似值

% X0=[x1 y1 x2 y2 ......]

[X0]=zbjs(S2,S3);

% 循环，未知参数改正量最大值mx<0.01mm停止

mx=100;

while mx>0.01

% 计算B，l，P 矩阵

[B,l,P]=BBll(X0,S2,S1,S3);

% P68页

x=(B'*P*B)\B'*P*l;% 改正量

mx=max(x);

v=B*x+l;

[rr,cc]=size(B);

m0=sqrt((v'*P*v)/(rr-cc));% 单位权中误差

Q=inv(B'*P*B);% 协因数阵

M=m0^2*Q;% 未知数方差矩阵

% 平差后坐标

X=X0'+x/1000;

X0=X';

end

% 坐标中误差

for i=1:28

det(i)=sqrt(M(i,i));

end

% 点位中误差

for i=1:14

ddet(i)=2.5*sqrt(det(i*2-1)+det(i*2));

end

ddet=ddet';

% 整理数据，输出最终结果

for i=1:14

XX(i,1)=X(2*i-1);

XX(i,2)=X(2*i);

DD(i,1)=det(2*i-1);

DD(i,2)=det(2*i);

end

jieguo(:,1)=1:14;

jieguo(:,2)=XX(:,1);

jieguo(:,3)=XX(:,2);

jieguo(:,4)=DD(:,1);

jieguo(:,5)=DD(:,2);

jieguo(:,6)=ddet;

disp('计算完成！')

fid=fopen('输出结果.txt','wt');

fprintf(fid,'点号 X(m) Y(m) 中误差X(mm) 中误差Y(mm) 点位中误差 \n');

fprintf(fid,' %2d %11.4f %11.4f %5.4f %5.4f %5.4f\ n',(jieguo)');

fprintf(fid,'单位权中误差：%8.4f',(m0)');

fclose(fid);

====================================================

% 单位转换

function [S3]=dfm2jd(S1)

for i=1:24

S3(i,1)=S1(i,1);

S3(i,2)=(S1(i,2)+S1(i,3)/60+S1(i,4)/3600);

End