end
（4）水准网平差
%间接平差后高程、高差及中误差 %确定系数矩阵A
A=zeros(num_obs,num_unknown); for i=1:num_obs
m=strmatch(start_point(i),id_bm); %从id_bm中寻找和起始一样点号的下标
n=strmatch(end_point(i),id_bm); %从id_bm中寻找和终点一样的点号下标
L(i)=height_dif(i)+height(m)-height(n); if m<=num_unknown %确定系数矩阵A A(i,m)=-1; end
if n<=num_unknown A(i,n)=1; end end
%权系数的确定、条件平差和中误差 p=diag(1./distance); %权矩阵 dh=inv(A'*p*A)*A'*p*L'; v=A*dh-L';
height_dif_v=height_dif'+v; %平差后的高差
height_v=height(1:num_unknown)'+dh; %平差后的未知点高程 xlgema0=sqrt(v'*p*v/(num_obs-num_unknown)); %单位权中误差 Qxx=inv(A'*p*A); %协因数
xlgema_H=xlgema0*sqrt(diag(Qxx)); %未知点高程平差值中误差 Qff=A*Qxx*A';
xlgema_h_dif=xlgema0*sqrt(diag(Qff));%高差平差值中误差
（5）输出结果
%计算成果和精度评定输出成.txt格式的文本文件
[filename,pathname]=uiputfile('.txt','请选择需要输出的文件'); data2=[pathname,filename]; fidout=fopen(data2,'wt');
fprintf(fidout,'一、所有点近似高程(height) \\n');
fprintf(fidout,' %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f\\n',height); %输出所有点近似高程，一行输出 fprintf(fidout,'\\n') %换行 fprintf(fidout,'二、权矩阵(p)\\n');
fprintf(fidout,' %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f\\n',p); %输出权矩阵
fprintf(fidout,'\\n')
fprintf(fidout,'三、平差后单位权中误差（xlgema0）\\n'); %输出平差后单位权中误差，一个值
fprintf(fidout,' %5.4f \\n',xlgema0); fprintf(fidout,'\\n')
fprintf(fidout,'四、平差后未知点高程（height_v）\\n'); %输出平差后未知点高程
fprintf(fidout,' %5.4f %5.4f %5.4f \\n',height_v); fprintf(fidout,'\\n')
fprintf(fidout,'五、平差后未知点精度（xlgema_H）\\n'); %输出平差后未知点精度
fprintf(fidout,' %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f \\n',xlgema_H);
fprintf(fidout,'\\n')
fprintf(fidout,'六、改正后的高差值（height_dif_v）\\n'); %输出改正后的高差值
fprintf(fidout,' %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f \\n',height_dif_v); fclose(fidout);
（6）最终代码
[filename, pathname] = uiputfile('*.txt', '保存数据文件到'); %通过对话窗口找到或创建要存放观测数据的文本文件，获取他们的文件名路径名 file = [pathname,filename]; fid = fopen(file,'wt');
num_known = input('请输入已知水准点个数：'); num_unknown = input('请输入未知水准点个数：'); num_obs = input('请输入观测值个数：');
fprintf(fid,'%5.0f %5.0f %5.0f\\n',num_known,num_unknown,num_obs); for i = 1:num_unknown
id_bm1 = input(['请输入第',num2str(i),'个未知水准点编号'],'s'); id_bm1 = upper(id_bm1); fprintf(fid,' %s',id_bm1); end
for i = 1:num_known
id_bm1 = input(['请输入第',num2str(i),'个已知水准点编号'],'s'); id_bm1 = upper(id_bm1); fprintf(fid,' %s',id_bm1);
height(i) = input(['请输入该已知水准点的高程']); end
fprintf(fid,'\\n'); %换行 for i = 1:num_known
fprintf(fid,' %5.4f',height(i)); end
fprintf(fid,'\\n'); %换行
for i = 1:num_obs
start_point = input(['请输入第',num2str(i),'个测段起始点编号'],'s');%输入num2str（i）个测段起始点编号
start_point = upper(start_point); %大小写字母的转化 end_point = input('终点编号','s');
end_point = upper(end_point); %大小写字母的转化 height_dif = input('测段高差观测值（m）'); distance = input('测段水准路线长度（km）');
fprintf(fid,'%s %s %4.3f %3.2f\\n',start_point,end_point,height_dif,distance); end
fclose('all'); clear;
[filename, pathname]= uigetfile('*.*'); %文件查找窗口 file=fullfile(pathname, filename); %合并路径文件名 fid=fopen(file,'rt'); %rt表示打开的只读文件
num_known = fscanf(fid,'%f',1); %fscanf 函数当能正常读入数据时，返回读入数据的个数，否则返回 %f'以数字格式读取，fid句柄，1.读的数目 num_unknown = fscanf(fid,'%f',1); %每次读一个数 num_obs = fscanf(fid,'%f',1);
num_all = num_known+ num_unknown;
for i = 1: num_all %按字符串方式逐个读取，并存储在字符串数组id_bm中，{}以字符串数组整个BM1存储
id_bm(i) = {fscanf(fid,'%s',1)}; %读出的是点号，字符型 end
height(1: num_unknown)= 0;
height(1+num_unknown: num_all)= fscanf(fid,'%f', num_known);%一一对应高程 for i=1: num_obs
start_point(i) = {fscanf(fid,'%s',1)};%所有起始点在字符串数组中 end_point(i) = {fscanf(fid,'%s',1)}; height_dif(i) = fscanf(fid,'%f',1); distance(i) = fscanf(fid,'%f',1); end
flag(1: num_unknown) = 0; %flag是一个变量
flag(1+num_unknown: num_all) = 1;%0表示该点近似高程未计算,1表示该点近似高程已计算。 while(1)
for i = 1: num_obs
m = strmatch(start_point(i) ,id_bm); %从id_bm中寻找和起始点一样点号的下标
n = strmatch(end_point(i) ,id_bm); %从id_bm中寻找和终点一样点号的下标
if (flag(m)==0 && flag(n)==1)
height(m) = height(n)-height_dif(i); flag(m) = 1;
end
if (flag(n)==0 && flag(m)==1)
height(n) = height(m)+height_dif(i); flag(n) = 1; end end
id_un = find(flag ==0); % 寻找未计算近似高程点的下标 if length(id_un) == 0 %如果无为计算近似高程点，终止循环 break; end end
r=num_obs-num_unknown; %r表示自由度 A=zeros(num_obs,num_unknown); %求解系数阵,初始值化为0 for i=1:num_obs
m=strmatch(start_point(i),id_bm); %从id_bm中寻找和起始一样点号的下标 n=strmatch(end_point(i),id_bm); %从id_bm中寻找和终点一样的点号下标 L(i)=height_dif(i)+height(m)-height(n); %从start_point(i)到end_point(i)高程值 if m<=num_unknown
A(i,m)=-1; %判断未知数的存在与否，开始赋值 end
if n<=num_unknown A(i,n)=1; end end
p=diag(1./distance); %权矩阵，函数diag生成对角矩阵，或提取矩阵对角元素 dh=inv(A'*p*A)*A'*p*L'; v=A*dh-L';
height_dif_v=height_dif'+v; %平差后的高差
height_v=height(1:num_unknown)'+dh; %平差后的未知点高程 xlgema0=sqrt(v'*p*v/r); %单位权中误差 Qxx=inv(A'*p*A);
xlgema_H=xlgema0*sqrt(diag(Qxx)); %未知点高程平差值中误差 Qff=A*Qxx*A';
xlgema_h_dif=xlgema0*sqrt(diag(Qff));%高差平差值中误差 [filename,pathname]=uiputfile('.txt','请选择需要输出的文件'); data2=[pathname,filename]; fidout=fopen(data2,'wt');
fprintf(fidout,'一、所有点近似高程(height) \\n');
fprintf(fidout,' %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f \\n',height); %输出所有点近似高程，一行输出 fprintf(fidout,'\\n') %换行 fprintf(fidout,'二、权矩阵(p)\\n');
fprintf(fidout,' %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f \\n',p); %输出权矩阵
fprintf(fidout,'\\n')
fprintf(fidout,'三、平差后单位权中误差（xlgema0）\\n'); %输出平差后单位权中误差，一个值
fprintf(fidout,' %5.4f \\n',xlgema0); fprintf(fidout,'\\n')
fprintf(fidout,'四、平差后未知点高程（height_v）\\n'); %输出平差后未知点高程 fprintf(fidout,' %5.4f %5.4f %5.4f \\n',height_v); fprintf(fidout,'\\n')
fprintf(fidout,'五、平差后未知点精度（xlgema_H）\\n'); %输出平差后未知点精度
fprintf(fidout,' %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f \\n',xlgema_H); fprintf(fidout,'\\n')
fprintf(fidout,'六、改正后的高差值（height_dif_v）\\n'); %输出改正后的高差值
fprintf(fidout,' %5.4f %5.4f %5.4f %5.4f %5.4f %5.4f \\n',height_dif_v); fclose(fidout);
3、最终成果
《误差理论与测量平差》书上108页例题