基于MATLAB的控制网平差程序设计--第四章源代码
基于MATLAB的控制网平差程序设计--第六章源代码
近似坐标计算的函数-calcux0y0函数（126页）function [x0,y0]=calcux0y0(x0,y0,e,d,sid,g,f,dir,s,t,az,pn,xyknow,xyunknow,point,aa,bb,cc)%本函数的作用是计算待定点的近似坐标format short; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% time=0;prelength=length(xyknow);non_orient=0;point_angle=0;while length(xyunknow)>0%考虑的计算方法有：1.极坐标;2.前方交会；3.测边交会；4.后方交会；%5.无定向导线的两种情况：（1）已知两个点；（2）分离的已知点与方位角；基本思路：%采用循环的方法逐一对每一个未知点进行以上各种方法条件的搜索，满足后即解算。
aa0=[];bb0=[];cc0=[];%记录搜索到两条观测边但需用户给顺序的点，注意要放在while 里面。
time=time+1; % 用于统计循环次数。
way=0;for i=xyunknow %依次循环向量中的各元素%============================================================= =====%方法1.极坐标条件搜索与计算-->way=1,基本思路：找到或求出一个方位角，找出一条边。
temp1=[]; temp2=[]; temp3=[]; temp4=[]; temp5=[]; temp6=[]; temp7=[]; temp8=[];temp9=[]; temp10=[];A=[];B=[];P=[];%第一步：寻找观测条件：两种情况：一是有已知方位角；二是由两个已知点及方向观测值推出方位角。
平差实验报告及完整matlab程序
fprintf(fn,'%6.3f\n',v);
fprintf(fn,'%s\n','未知参数估值x');
fprintf(fn,'%6.3f\n',x1);
fprintf(fn,'%s\n','各点高程平差值');
fprintf(fn,'%6.3f\n',x2);
86.806
3
31.225
344.8
14.846
4
-71.952
149.2
25.706
5
-61.084
142.9
31.216
6
-44.178
250.0
42.626
7
10.847
128.2
8
16.350
98.0
9
11.409
196.1
五、计算过程
1、经典自由网平差（以6号点为例）
由已知可列出误差方程
（1）
fnal=strcat(p1,f);fpath=p1;
fm=fopen(fnal,'r');
while (~feof(fm))
s=fscanf(fm,'%f',inf)
end
%---------------从文件获取数据----------------------------------------------
switch dh_index
case 1
b=b(:,2:6);
x0=x0(2:6,:);
基于MATLAB的水准网和测边网平差程序设计
基于MATLAB的水准网和测边网平差程序设计摘要MATLAB是目前在研究机构广泛应用的一种数值计算及图形工具软件，它的特点是语法结构简明、数值计算高效、图形功能完备，特别适合非专业编程员完成数值计算、科学试验处理等任务。
以往的测量数据处理方法需要编制特定的处理矩阵运算程序，而且程度复杂，难度大。
本文介绍一种基于MATLAB的水准网和测边网的程序设计方法，与其它算法语言相比，具有编程简单，运算速度快的特点。
文中分别阐述了水准网和测边网程序的理论基础、实现步骤和运行结果。
通过实例的分析，总结出利用MATLAB对测量数据处理有很大的应用价值，它缩短了编程的时间，提高工作效率。
关键词：MATLAB；水准网；测边网；程序设计ABSTRAC TMATLAB is one species of numerical-values calculation and graphic tools software which is widely used to apply at research institutions at present. The particularities are: concise grammar-structure、highly efficient in numerical values calculating、complete function of graphs、especially it is adapted to evildoing professional programmer to accomplish the tasks that are numerical-values calculating and scientific experiments treating. The ancient methods of measured data-processing need establishing special proceedings of treating matrices operation, moreover, it is complex and greatly difficult.This article introduces one programming method dealing with leveling and measuring edge network based on MATLAB. Compared with other algorithm language, it has particularities which are simply programming and quickly operating. The article separately expatiate the theories basics、realizing steps and running results at leveling and measuring edge network. With the analysis of examples, it has prodigious application value in measured data-processing by use of MATLAB. Moreover, it shortens programming time and improves working effectiveness.Key words：MATLAB；leveling network；measuring edge network；programming目录绪论 (4)1. MATLAB软件简介 (5)2．MATLAB 在测量平差中的应用 (6)2.1测量平差原理的概述 (6)2.2平差程序总体方案 (7)3.1程序的功能 (8)3.2水准模型网的间接平差 (8)3.2.1 “权”值的确定 (8)3.2.2 水准路线的平差计算 (9)3.2.3 精度评定 (11)3.3水准网间接平差程序信息设计 (11)3.4 水准网程序与使用说明 (12)3.4.1 水准网程序流程图 (12)3.4.2 水准网程序的使用 (12)3.5案例 (13)4. 测边网平差程序设计 (15)4.1数学模型 (15)4.1.1 误差方程和法方程的组成 (15)4.1.2 边长观测的权 (15)4.1.3 解算法方程 (16)4.1.4 精度评定 (19)4.2 测边网平差信息设计 (20)4.2.1 主要的技术要求 (21)4．3利用MATLAB的绘图语句绘制网图 (21)4.4测边网程序和使用说明 (22)4.5 程序代码说明： (23)4.6程序的使用算例 (25)结论 (29)致谢 (30)参考文献 (31)附录一 (32)附录二 (36)附录三 (46)绪论作为一名测量技术人员，如果不掌握一门PC机编程语言与便携计算工具，要想提高测量工作的效率几乎寸步难行。
基于MATLAB的高程控制网平差系统的设计与应用
基于MATLAB的高程控制网平差系统的设计与应用高霞;应建福【摘要】文中利用MATLAB开发的计算机程序代码短小高效,能够快速的实现图形绘制和精度分析,与利用C++等高级编程语言开发的平差处理程序相比,可以在短时间内完成对平差的设计和验证.【期刊名称】《矿山测量》【年(卷),期】2016(044)004【总页数】3页(P63-65)【关键词】MATLAB;高程控制网;条件平差【作者】高霞;应建福【作者单位】甘肃工业职业技术学院测绘学院,甘肃天水741025;天水天正设计咨询有限公司,甘肃天水741000【正文语种】中文【中图分类】P224测绘科学是一门以大规模数据甚至是海量数据处理、分析与应用为基础的学科，其各项具体工作如测量计算、测量平差、大地测量数据处理、GPS数据结算、图像处理、坐标换算等都涉及大量的数值计算与分析。
运用MATLAB解算测绘信息数据处理问题具有其它程序设计语言难以比拟的优越性。
在高程控制网平差的计算中涉及大量的矩阵求逆和方程组的求解等问题，而MATLAB正是解决此类问题较好的数学软件。
对于一些常用的平差模型，可以用MATLAB软件通过调用其函数来完成，从而大大提高工作效率。
(1)简单易用的程序语言。
MATLAB是一个高级的矩阵/列阵语言，它包括控制语句、函数、数据结构、输入输出和面向对象编程特点。
用户可以在命令窗口中将输入语句与执行命令同步，也可以先编好一个较大的复杂的应用程序后一起运行。
而且这种语言可移植性较好，可拓展性极强，这也是MATLAB能够深入到科学研究及工程计算的重要原因。
(2)强大的科学计算、数据处理能力。
MATLAB是一个包含大量计算算法的集合，其拥有600多个工程中要用到的数学运算函数，可以方便的实现用户所需的各种计算功能。
在通常情况下，使用编程工作量会大大的减少。
其函数能解决的问题大致包括矩阵运算和线性方程组的求解、特征向量、微分方程及偏微分方程的求解、工程中的优化问题、多维数组操作以及建模动态仿真等。
实验三-利用matlab程序设计语言完成某工程导线网平差计算
实验三-利用m a t l a b程序设计语言完成某工程导线网平差计算(总11页)本页仅作为文档封面，使用时可以删除This document is for reference only-rar21year.March实验三利用matlab程序设计语言完成某工程导线网平差计算实验数据；某工程项目按城市测量规范（CJJ8-99）不设一个二级导线网作为首级平面控制网，主要技术要求为：平均边长200cm，测角中误差±8,导线全长相对闭合差<1/10000,最弱点的点位中误差不得大于5cm，经过测量得到观测数据，设角度为等精度观测值、测角中误差为m=±8秒，鞭长光电测距、测距中误差为m=±√smm，根据所学的‘误差理论与测量平差基础’提出一个最佳的平差方案，利用matlab完成该网的严密平差级精度评定计算;平差程序设计思路:1采用间接平差方法,12个点的坐标的平差值作为参数.利用matlab进行坐标反算，求出已知坐标方位角；根据已知图形各观测方向方位角；2计算各待定点的近似坐标,然后反算出近似方位角，近似边.计算各边坐标方位角改正数系数；3确定角和边的权,角度权Pj=1；边长权Ps=100/S；4计算角度和边长的误差方程系数和常数项，列出误差方程系数矩阵B，算出Nbb=B’PB,W=B’Pl,参数改正数x=inv(Nbb)*W;角度和边长改正数V=Bx-l；6 建立法方程和解算x,计算坐标平差值, 精度计算；程序代码以及说明:s10=;s20=;s30=;s40=;s50=;s60=;s70=;s80=;s90=;s100=;s110=;s120=;s130=;s140=; %已知点间距离Xa=;Ya=;Xb=;Yb=;Xc=;Yc=;Xd=;Yd=;Xe=;Ye=;Xf=;Yf=; %已知点坐标值a0=atand((Yb-Ya)/(Xb-Xa))+180;d0=atand((Yd-Yc)/(Xd-Xc));f0=atand((Yf-Ye)/(Xf-Xe))+360; %坐标反算方位角a1=a0+(163+45/60+4/3600)-180a2=a1+(64+58/60+37/3600)-180;a3=a2+(250+18/60+11/3600)-180;a4=a3+(103+57/60+34/3600)-180;a5=d0+(83+8/60+5/3600)+180;a6=a5+(258+54/60+18/3600)-180-360;a7=a6+(249+13/60+17/3600)-180;a8=a7+(207+32/60+34/3600)-180;a9=a8+(169+10/60+30/3600)-180;a10=a9+(98+22/60+4/3600)-180;a12=f0+(111+14/60+23/3600)-180;a13=a12+(79+20/60+18/3600)-180;a14=a13+(268+6/60+4/3600)-180;a15=a14+(180+41/60+18/3600)-180; %推算个点方位角aa=[a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a12 a13 a14 a15]'X20=Xb+s10*cosd(a1);X30=X20+s20*cosd(a2);X40=X30+s30*cosd(a3);X50a=X40+s40*cosd(a4);X60=Xd+s50*cosd(a5);X70=X60+s60*cosd(a6);X80=X70+s70*cosd(a7);X90=X80+s80*cosd(a8);X100=X90+s90*cosd(a9);X50c=X100+s100*cosd(a10);X130=Xf+s110*cosd(a12);X140=X130+s120*cosd(a13);X150=X140+s130*cosd(a14);X50e=X150+s140*cosd(a15); %各点横坐标近似值X0=[X20 X30 X40 X60 X70 X80 X90 X100 X130 X140 X150 X50a X50c X50e]'Y20=Yb+s10*sind(a1);Y30=Y20+s20*sind(a2);Y40=Y30+s30*sind(a3);Y50a=Y40+s40*sind(a4);Y60=Yd+s50*sind(a5);Y70=Y60+s60*sind(a6);Y80=Y70+s70*sind(a7);Y90=Y80+s80*sind(a8);Y100=Y90+s90*sind(a9);Y50c=Y100+s100*sind(a10);Y130=Yf+s110*sind(a12);Y140=Y130+s120*sind(a13);Y150=Y140+s130*sind(a14);Y50e=Y150+s140*sind(a15); %个点从坐标近似值Y0=[Y20 Y30 Y40 Y60 Y70 Y80 Y90 Y100 Y130 Y140 Y150 Y50a Y50c Y50e]'P=[X0 Y0];X50=(X50a+X50c+X50e)/3Y50=(Y50a+Y50c+Y50e)/3s4=sqrt((Y40-Y50)^2+(X40-X50)^2);s1=sqrt((Y100-Y50)^2+(X100-X50)^2);s14=sqrt((Y150-Y50)^2+(X150-X50)^2);A1=[cosd(a1) cosd(a2) cosd(a3) cosd(a4) cos(a5) cosd(a6) cosd(a7) cosd(a8) cosd(a9) cosd(a10) cosd(a12) cosd(a13) cosd(a14)cosd(a15)]';B11=[sind(a1) sind(a2) sind(a3) sind(a4) sin(a5) sind(a6) sind(a7) sind(a8) sind(a9) sind(a10) sind(a12) sind(a13) sind(a14) sind(a15)]'; s=blkdiag(s10,s20,s30,s4,s50,s60,s70,s80,s90,s10',s110,s120,s130,s 14);a=*inv(s)*B11b=*inv(s)*A1ab4=atand((Y50-Y40)/(X50-X40))+180;ab10=atand((Y50-Y100)/(X50-X100));ab14=atand((Y50-Y150)/(X50-X150))+360;m4=ab4-a3+180;m10=ab10-a9+180;m11=ab4-ab10;m15=ab14-a14+180;m16=ab10-ab14+360;m04=103+57/60+34/3600;m010=98+22/60+4/3600;m011=94+53/60+50/3600;m015=180+41/60+18/3600;m016=ab10-ab14+360;l=[0 0 0 m4-103-57/60-34/3600 0 0 0 0 0 m10-98-22/60-4/3600m11-94-53/60-50/3600 0 0 0 m15-180-41/60-18/3600 m16-103-23/60-8/3600 0 0 0 s40-s4 0 0 0 0 0 s100-s1 0 0 0 s140-s14]';e1=(abs(X20-Xb))/s10;e2=(abs(X30-X20))/s20;e3=(abs(X40-X30))/s30;e4=(abs(X50-X40))/s4;e5=(abs(X60-Xd))/s50;e6=(abs(X70-X60))/s60;e7=(abs(X80-X70))/s70;e8=(abs(X90-X80))/s80;e9=(abs(X100-X90))/s90;e10=(abs(X50-X100))/s1;e11=(abs(X130-Xf))/s110;e12=(abs(X140-X130))/s120;e13=(abs(X150-X140))/s130;e14=(abs(X50-X150))/s14; e=[e1 e2 e3 e4 e5 e6 e7 e8 e9 e10 e11 e12 e13 e14]'m1=(abs(Y20-Yb))/s10;m2=(abs(Y30-Y20))/s20;m3=(abs(Y40-Y30))/s30;m4=(abs(Y50-Y40))/s4;m5=(abs(Y60-Yd))/s50;m6=(abs(Y70-Y60))/s60;m7=(abs(Y80-Y70))/s70;m8=(abs(Y90-Y80))/s80;m9=(abs(Y100-Y90))/s90;m10=(abs(Y50-Y100))/s1;m11=(abs(Y130-Yf))/s110;m12=(abs(Y140-Y130))/s120;m13=(abs(Y150-Y140))/s130;m14=(abs(Y50-Y150))/s14; m=[m1 m2 m3 m4 m5 m6 m7 m8 m9 m10 m11 m12 m13 m14]' %以上为求得误差方程系数B=[ 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 00 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0] %系数矩阵BP=blkdiag(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,100/s10,100/s20,100/s30,1 00/s40,100/s50,100/s60,100/s70,100/s80,100/s90,100/s100,100/s 110,100/s120,100/s130,100/s140); %定义权矩阵Nbb=B'*P*BW=B'*P*l;x=inv(Nbb)*WV=B*x-l;inv(Nbb);Y=V'*P*V;O=sqrt(Y/6)*3600 %精度评定计算结果：平差值坐标X： +003 *Qx1= Qy1= Qx2= Qy2= ……Qx15= Qy15=。
基于MATLAB的控制网平差程序设计--第五章源代码 - 副本
观测数据读入程序，rddat1函数（85页）global net ed dd sd dd1 pn x0 y0 m1 m2 m3 ms pp e d sid md g f dir ni si ma s t az aa bb cc rt rr tt global pathname filenamex0=[];y0=[];e=[];d=[];sid=[];g=[];f=[];dir=[];si=[];ni=[];s=[];t=[];az=[];pn=[];[filename,pathname]=uigetfile('*.txt','请选择原始数据');fit1=fopen(strcat(pathname,filename),'rt');if(fit1==-1)msgbox('Input File or Path is not correct','Warning','warn');return;endnet=fscanf(fit1,'%d',1);[a]=fscanf(fit1,'%d',3);ed=a(1);dd=a(2);dd1=a(3);sd=ed+dd;[pn]=fscanf(fit1,'%d',sd);[a]=fscanf(fit1,'%f',2*ed);for i=1:edx0(i)=a(2*i-1);y0(i)=a(2*i);end[a]=fscanf(fit1,'%d',3);m1=a(1);m2=a(2);m3=a(3);isid=0;[a]=fscanf(fit1,'%f',2);ms=a(1);pp=a(2);[a]=fscanf(fit1,'%d %d %f',3*m1);for i=1:m1e(i)=a(3*i-2);d(i)=a(3*i-1);sid(i)=a(3*i);end[e,i1]=chkdat(sd,pn,e);[d,i2]=chkdat(sd,pn,d);i3=0;isid=i1+i2+i3;idir=0;md=fscanf(fit1,'%f',1);[a]=fscanf(fit1,'%d %d %f',3*m2);for i=1:m2n1(i)=a(3*i-2);n2(i)=a(3*i-1);unk(i)=a(3*i);end[n1,i1]=chkdat(sd,pn,n1);[n2,i2]=chkdat(sd,pn,n2);i3=0;ik=1;si(1)=1;for i=1:sdii=0;for j=1:m2if(n1(j)==j)ii=ii+1;g(ik)=n1(j);f(ik)=n2(j);dir(ik)=unk(j);ik=ik+1;endendni(i)=ii;si(i+1)=si(i)+ni(i);endidir=i1+i2+i3;iaz=0;if(m3>0)ma=fscanf(fit1,'%f',1);[a]=fscanf(fit1,'%d %d %f',3*m3);for i=1:m3s(i)=a(3*i-2);t(i)=a(3*i-1);az(i)=a(3*i);end[s,i1]=chkdat(sd,pn,s);[t,i2]=chkdat(sd,pn,t);i3=0;iaz=i1+i2+i3;endkk=isid+idir+iaz;if(kk>0)msgbox('Error by function rddat1','Warning','warn');return;endfclose('all');open(strcat(pathname,net_name,b_datafile));return误差方程与法方程的组成函数-obnorm函数（90页）function obnormglobal ed dd dd1 ni si e d g f s tglobal m1 m2 m3 ms pp md ma x0 y0 sid dir az c fit1 fit2 global a q1 pa3 qls wlo=2062.648062470964;m=m1+m2+m3;n=2*dd;sum=n*(n+1)/2.0;sd=ed+dd;a(1:m,1:9)=0.0;for i=1:sdii=4*(ni(i)+1);pa3(i,1:ii)=0.0;endc(1:sum)=0.0;w(1:n)=0.0;for i=1:m1 %边长观测误差方程dx=x0(d(i))-x0(e(i));dy=y0(d(i))-y0(e(i));ss=sqrt(dx*dx+dy*dy);cosa=dx/ss;sina=dy/ss;a(i,1)=2*e(i)-1-2*ed+1.0e-9;a(i,2)=-cosa;a(i,3)=a(i,1)+1;a(i,4)=-sina;a(i,5)=2*d(i)-1-2*ed+1.0e-9;a(i,6)=cosa;a(i,7)=a(i,5)+1;a(i,8)=sina;a(i,9)=100.0*(ss-sid(i));q1(i)=(ms^2+(ss*pp*0.0001)^2);endq1(m1+1:m2+m1)=md*md;for i=1:sdif(ni(i)==0)continue;endjj=5;z0=0.0;zal=0;for j=si(i):s(i)+ni(i)-1dx=x0(f(j))-x0(g(j));dy=y0(f(j))-y0(g(j));a0=alfa(dx,dy);z1=a0-dir(j);if(z1<=0.0)z1=z1+2.0*pi;endzal=zal+1./ql(m1+j);z0=z0+z1/q1(m1+j);endz0=z0/zal;for j=si(i):si(i)+ni(i)-1dx=x0(f(j))-x0(g(j));dy=y0(f(j))-y0(g(j));ss=dx*dx+dy*dy;a0=alfa(dx,dy);ai=-dy/ss*lo;bi=dx/ss*lo;ii=m1+j;a(ii,1)=2*g(j)-1-2*ed+1.0e-9;a(ii,2)=-ai;a(ii,3)=a(ii,1)+1;a(ii,4)=-bi;a(ii,5)=2*f(j)-1-2*ed+1.0e-9;a(ii,6)=ai;a(ii,7)=a(ii,5)+1;a(ii,8)=bi;ss=dir(j)+z0;if(ss>=2.0*pi)ss=ss-2.0*pi;enda(ii,9)=(a0-ss)*lo*100.0;pa3(i,jj)=a(ii,5);pa3(i,jj+1)=a(ii,6)/q1(ii);pa3(i,jj+2)=a(ii,7);pa3(i,jj+3)=a(ii,8)/q1(ii);pa3(i,2)=pa3(i,2)+a(ii,2)/q1(ii);pa3(i,4)=pa3(i,4)+a(ii,4)/q1(ii);jj=jj+4;endpa3(i,1)=a(ii,1);pa3(i,3)=a(ii,3);qls(i)=-zal;endfor i=1:m3dx=x0(t(i))-x0(s(i));dy=y0(t(i))-y0(s(i));a0=alfa(dx,dy,a0);ss=dx*dx+dy*dy;ai=-dy/ss*lo;bi=dx/ss*lo;ii=m1+m2+i;a(ii,1)=2*s(i)-1-2*ed+1.0e-9;a(ii,2)=-ai;a(ii,3)=a(ii,1)+1;a(ii,4)=-bi;a(ii,5)=2*t(i)-1-2*ed+1.0e-9;a(ii,6)=ai;a(ii,7)=a(ii,5)+1;a(ii,8)=bi;if((a0-az(i))>pi)a(ii,9)=(a0-az(i)-2.0*p)*lo*100;elsea(ii,9)=(a0-az(i))*lo*100.0;endq1(ii)=ma*ma;endfor i=1:m %形成法方程for j=1:4jj=fix(a(i,2*j-1));if(jj<=0)continue;endw(jj)=w(jj)+a(i,2*j)*a(i,9)/q1(i);di=(jj-1)*(n-jj/2.0);for k=1:4kk=fix(a(i,2*k-1));if(kk<=0|jj>kk)continue;endc(di+kk)=c(di+kk)+a(i,2*k)*a(i,2*j)/q1(i);endendendif(m2>0) %和误差方程形成法方程for i=1:sdif(ni(i)==0)continue;endfor j=1:2*(ni(i)+1)jj=fix(pa3(i,2*j-1));if(jj<=0)continue;enddi=(jj-1)*(n-jj/2);for k=1:2*(ni(i)+1)kk=fix(pa3(i,2*k-1));if(kk<=0|jj>kk)continue;endc(di+kk)=c(di+kk)+pa3(i,2*k)*pa3(i,2*j)/qls(i);endendendendreturn平差值与精度评定（94页）global net ed dd sd dd1 pn x0 y0 m1 m2 m3 ms pp e d sid md g f dir ni si ma s t az global aa bb cc rt rr ttglobal a q1 pa3 qls w c x y uw0global pathname net_name s_datafile a_datafile;fit2=fopen(strcat(pathname,net_name,a_datafile,'wt');if(fit2==-1)msgbox('Input File or Path is not correct','Warning','warn');return;endk=1;while(k)m=m1+m2+m3;obnorm;c=invsqr(c,2*dd);[uw0,k]=adjxy(fit2);endellipse(uw0,fit2);n=2*dd;sum=n*(n+1)/2.0;n1=2*(ed+dd);sum1=n1*(n1+1)/2.0;for i=sum1:-1:sum1-sum+1c(i)=c(i-(sum1-sum));endfor i=1:2*eddi=(i-1)*(n1-i/2.0);for j=i:n1if(j==i)c(di+j)=0.00000001;elsec(di+j)=0.0;endendendif(m1>0)adjs(uw0,fit2);endif(m2>0)adjd(uw0,fit2);endif(m3>0)adja(uw0,fit2);endfclose(fit2);open(strcat(pathname,net_name,a_datafile));坐标改正数计算及单位权中误差计算函数-adjxy函数（95页）function [uw0,k]=adjxy(fit2)global ed dd dd1 ni si e d g f s t pn x yglobal m1 m2 m3 ms md ma x0 y0 sid dir az cglobal a q1 pa3 qls wsd=ed+dd;n=2*dd;k=0;for i=1:ndxy(i)=0.0;di=(i-1)*(n-i/2.0);for j=1:ndj=(j-1)*(n-j/2.0);if(j<i)dxy(i)=dxy(i)-c(dj+i)*w(j);elsedxy(i)=dxy(i)-c(di+j)*w(j);endendif(abs(dxy(i))>1.0)k=1;enddxy(i)=dxy(i)/100.0;endx(1:ed)=x0(1:ed);y(1:ed)=y0(1:ed);for i=1:ddx(ed+i)=x0(ed+i)+dxy(2*i-1);y(ed+i)=y0(ed+i)+dxy(2*i);endx0(1:sd)=x(1:sd);y0(1:sd)=y(1:sd);for i=1:sdif(i<=ed)vx(i)=0.0;vy(i)=0.0;elsevx(i)=dxy(2*(i-ed)-1);vy(i)=dxy(2*(i-ed));endendfprintf(fit2,' adjusted coordinates\n');fprintf(fit2,' pn vx x vy y\n');for i=1:sdfprintf(fit2,' %6d %8.4f %14.4f %8.4f %14.4f\n',pn(i),vx(i),x(i),vy(i),y(i));endpvv=0.0;for i=1:npvv=pvv+w(i)*dxy(i)*100.0;endm=m1+m2+m3;for i=1:mpvv=pvv+a(i,9)*a(i,9)/ql(i);endif(m2>0)ii=0;for i=1:sdif(ni(i)~=0;ii=ii+1;endendenduw0=sqrt(pvv/((m-n-ii)*1.0e0));return计算各点误差椭圆-ellipse函数（98页）function ellipse(uw0,fit2)glosbal ed dd pn c x y ai bi fifprintf(fit2,' parameter of error ellipse\n');fprintf(fit2,' pn(i) mx my mm a b fi\n'); n=2.0*dd;maxmm=0.0;smm=0.0;for i=1:ddii=ed+i;di=(2*i-2)*(n-(2*i-1)/2.0);dj=(2*i-1)*(n-i);q1=c(di+2*i-1);q2=c(dj+2*i);q3=c(di+2*i);d1=sqrt(abs(q1+q2-q3));d=uw0*dl;xx=q1-q2;yy=2*q3;zz=q1+q2;mx1=sqrt(q1);my1=sqrt(q2);mm1=sqrt(zz);if(abs(xx)<1d-10)fi(i)=sign(90.0,q3);elsefi(i)=atan(yy/xx)*57.2958;endif(xx>=0&yy>=0)fi(i)=fi(i)/2.0;elseif(xx>=0&yy<=0)fi(i)=(fi(i)+360)/2.0;elseif(xx<0)fi(i)=(fi(i)+180)/2.0;endww=sqrt(xx*xx+yy*yy);a1=sqrt((zz+ww)/2.0);b1=sqrt((zz-ww)/2.0);ab1=a1-b1;mx=uw0*mx1;my=uw0*my1;mm=uw0*mm1;if(mm>maxmm)maxmm=mm;i1=ii;endsmm=smm+mm;ai(i)=uw0*a1;bi(i)=uw0*b1;ab=ai(i)-bi(i);fprintf(fit2,' %10d %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f\n',pn(ii),mx,my,mm,ai(i),bi(i),fi(i));endsmm=smm/dd;fprintf(fit2,' mse of unit weight= %9.6f\n',uw0);fprintf(fit2,' the maximum station error mm= %8.3f(cm) pn= %4d\n',maxmm,pn(i1));fprintf(fit2,' the average station error mm= %8.3f\n',smm);return边长观测值平差值改正数及精度评定-adjs函数（101页）function adjs(uw,fit2)global ed dd sd pn m1 e d sid x yfprintf(fit2,'adjusted sides\n');fprintf(fit2,'i pn(e) pn(d) side vs(cm) side+vs ms\n');for i=1:m1dx=x(d(i))-x(e(i));dy=y(d(i))-y(e(i));if(e(i)<d(i))[maa,mss]=trel(uw,e(i),d(i));else[maa,mss]=trel(uw,d(i),e(i));endss=dx*dx+dy*dy;ss=sqrt(ss);vs=(ss-sid(i))*100;sid1=sid(i)+vs/100;fprintf(fit2,' %3d %8d %8d %15.4f %10.4f %15.4f %8.2f\n',i,pn(e(i)),pn(d(i)),sid(i),vs,sid1,mss);endreturn方向观测值的平差与精度评定-adjd函数（102页）function adjd(uw,fit2)global ed dd dd1 sd pn ni si e d g f s t netglobal m1 m2 m3 ms pp md ma x0 y0 x y sid dir az cglobal a q1 pa3 qls wfprintf(fit2,'adusted directions and their accuracy\n');fprintf(fit2,' i pn(g) pn(f) dir vd(") dir+vd\n');lo=206264.8062470964;for j=1:m2q1(j+m1)=md*md;endfor i=1:sdzi=0.0;zal=0.0;for j=si(i):si(i)+ni(i)-1dx=x(f(j))-x(g(j));dy=y(f(j))-y(g(j));a0=alfa(dx,dy);z1=a0-dir(j);if(z1<0.0)z1=z1+2.0*pi;endzal=zal+1./ql(m1+j);zi=zi+z1/ql(m1+j);endif(ni(i)~=0)zi=zi/zal;endfor j=si(i):si(i)+ni(i)-1dx=x(f(j))-x(g(j));dy=y(f(j))-y(g(j));a0=alfa(dx,dy);ss=dir(j)+zi;if(ss>=2.0*pi)ss=ss-2*pi;endvd=(a0-ss)*lo;dir1=dir(j)+vd/lo;dir1=rad_dms(dir1);dir(j)=rad_dms(dir(j));fprintf(fit2,' %3d %8d %8d %14.5f %9.2f %16.5f\n', j,pn(g(j)),pn(f(j)),dir(j),vd,dir1);endendfprintf(fit2,'adusted directions\n');fprintf(fit2,' i pn(g) pn(f) dir az ma\n');for i=1:sdfor j=si(i):si(i)+ni(i)-1dx=x(f(j))-x(g(j));dy=y(f(j))-y(g(j));if(f(j)<g(j))[maa,mss]=trel(uw,f(j),g(j));else[maa,mss]=trel(uw,g(j),f(j));enda0=alfa(dx,dy);if(j==si(i))a00=a0;endss=a0-a00;if(ss<0)ss=ss+2.0*pi;endss=rad_dms(ss);a0=rad_dms(a0);fprintf(fit2,' %3d %10d %10d %14.5d %14.5d %10.4f\n',j,pn(g(j)),pn(f(j)),ss,a0,maa);endendreturn方位角观测值改正数与精度评定-adja函数（104页）function adja(uw,fit2)global m3 x y s t azl0=206264.8062470964;fprintf(fit2,'adjusted azimath \n i pn(s) pn(t) az va(") az+va ma');for i=1:m3dx=x(t(i))-x(s(i));dy=y(t(i))-y(s(i));a0=alfa(dx,dy);if(t(i)<s(i))[maa,mss]=trel(uw,t(i),s(i));else[maa,mss]=trel(uw,s(i),t(i));endif((a0-az(i))>pi)va=(a0-az(i)-2*pi)*lo;elseva=(a0-az(i))*lo;endaz1=az(i)+va/lo;az1=wg(az1);az(i)=wg(az(i));fprintf(fit2,' %3d %8d %8d %13.5f %8.2f %12.5f %8.2f',i,pn(s(i)),pn(t(i)),az(i),va,az1,maa);endreturn调用的的trel函数function [maa,mss]=trel(uw,rr,tt)global ed dd x0 y0 pn cn=2*(dd+ed);dx=x0(tt)-x0(rr);dy=y0(tt)-y0(rr);ss=sqrt(dx*dx+dy*dy);a0=alfa(dx,dy);aij=2062.648*sin(a0)/ss;bij=-2062.648*cos(a0)/ss;b=a+1;f=2*tt-1;g=f+1;da=(a-1)*(n-a/2.0);db=(b-1)*(n-b/2.0);df=(f-1)*(n-f/2.0);dg=(g-1)*(n-g/2.0);q1=c(da+a)+c(df+f)-2.0*c(da+f);q2=c(db+b)+c(dg+g)-2.0*c(db+g);q3=c(da+b)+c(df+g)-c(da+g)-c(db+f);x=q1-q2;y=2.0*q3;z=q1+q2;qs=q1*cos(a0)^2+q2*sin(a0)^2+q3*sin(2.0*a0); qa=q1*aij*aij+q2*bij*bij+2.0*q3*aij*bij;ms=sqrt(qs)*uw;maa=sqrt(qa)*uw;return。
matlab在测量控制网优化设计与平差中的应用
Ｎ －ＳＥ Ｗ Ｅ ．应属于 Ｗ向 Ｗ Ｎ 构造休系． 二期构造压性结构而 （ 褶皱轴面和压性断层） Ｎ 呈Ｎ Ｅ 向， 张性结构Ｉ走向Ｎ ｂ ｌ ＷＷ。 组扭性面走向Ｎ Ｅ Ｅ ；另一组
ｘ二Ｂｓｎ ｎｉ ｉ ＋Ａｃ ｓ ｏ ｉ ｉ ７ｓ ｏ Ｔｃ ｓ
Ｙ二 一 ｃｓ ｎ ＋ ｃｓ ｃｓ ｏ ７ ｉ ｉ Ａ ｏ ７ ｏ ｉ Ｂ ｓ 式， １ ，
ｉ － 误差椭圆的参数变觉， 一 般取０．＂３０。，ｉ －６ 角增量大小取决于绘制椭圆的粕度： Ａ Ｂ— 分别为椭圆的长、短半轴； ．
了 ９ 小 ， 中 为 长半轴方位角值 。 ＝ ０。一 ． ３ 粗俘与框图
Ｍ ＴＡ 有两种常用的＿作方式：一种是直接交互的 ＡＬＢ ［ 指令行操作方式；另 一 种是Ｍ文件的编程方式。在前一种 工作方式下，Ｍ Ｔ Ａ Ａ Ｌ Ｂ被当作一种高级的数学演算和计算 可视器来使用。对于简单问题，在ＭＡ Ｌ Ｂ的提示符下直 ＴＡ 接输入命令是快速有效的。然而，当命令数量增加或希望 改变一个或几个变量的值或某个命令需要重复多次时，直 接输入就非常麻烦。这时可把许多可执行的ＭＡ Ｌ Ｂ命令 ＴＡ 放在Ｍ文件中，只要在ＭＡ Ｌ Ｂ提示符下愉入Ｍ文件的 ＴＡ 文件名，即可执行多个命令。以下是针对＿述问题由 Ｌ ＭＡ Ｌ Ｂ软件编辑器而建立的Ｍ文件． ＴＡ
由干系敬阵 Ａ列亏 ．位的．小二乘解就不唯一 为此
ＡＰＸ ＡＰ Ｔ ＝Ｔ Ａ Ｌ 即Ｎ Ａ Ａ 二Ｔ 则 Ｐ
协因数阵为：
（） ２
采 在 小乘 －ｉ 加最 范 气 一ｉ 用最 二 内Ｖ 和 权 小 批ＴＸｍ ｍ ｎ ｎ
的准则下，来求未知参致的最佳估值ｘ． 即在 Ｑ ＝ Ｉ（Ｐ）１ ， Ｎ Ａ Ａ一 ｕ ＝ Ｔ （ ３ ） Ｖ Ｖ ｍｎ Ｐ＝ ｉ 若现增加一组观测值Ｌ 二其权为Ｐ．相应的法方程系数 ｉ 下组成法方程
matlab常见经典平差程序
matlab常见经典平差程序希望本人编写的经典平差可以对matlab初学者以及测绘专业的学生有用By cowboyfunction void()sprintf('请选择平差类型')sprintf('1: 参数平差\n2: 条件平差\n3：具有条件的参数平差\n4: 具有参数的条件平差\n') a=input('请输入对应号码\n');switch (a)case 1% 参数平差V=AX-Lsprintf('1:参数平差V=AX-L')n=input('请输入n=');t=input('请输入t=');A=input('请输入系数矩阵A=');L=input('请输入常数项矩阵L=');P=diag(input('请输入权阵P='))%P=input('请输入权阵P=');N=A'*P*A;U=A'*P*L;sprintf('开始计算')X=inv(N)*(U);V=A*X-L;Qx=inv(N);u=sqrt((V'*P*V)/(n-t));fid=fopen('data.txt','wt'); %写入文件路径，文件输出fprintf(fid,'=====================参数平差: V=AX-L =====================\n');fprintf(fid,'输出n,t\n') ;fprintf(fid,'%12.5f %12.5f\n',n ,t) ;%n,t fprintf(fid,'输出矩阵A\n')[m,n]=size(A); %Afor i=1:1:mfor j=1:1:nif j==nfprintf(fid,'%12.5f\n',A(i,j));elsefprintf(fid,'%12.5f\t',A(i,j));endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵L\n') ;[m,n]=size(L); %Lfor i=1:1:mfor j=1:1:nif j==nfprintf(fid,'%12.5f\n',L(i,j));elsefprintf(fid,'%12.5f\t',L(i,j));endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵P\n') ;[m,n]=size(P); %P for i=1:1:mfor j=1:1:nif j==nfprintf(fid,'%12.5f\n',P(i,j));elsefprintf(fid,'%12.5f\t',P(i,j));endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵X\n'); [m,n]=size(X); %X for i=1:1:m for j=1:1:nif j==nfprintf(fid,'%12.5f\n',X(i,j));elsefprintf(fid,'%12.5f\t',X(i,j));endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵V\n'); [m,n]=size(V); %V for i=1:1:m for j=1:1:nif j==nfprintf(fid,'%12.5f\n',V(i,j));elsefprintf(fid,'%12.5f\t',V(i,j));endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵Qx\n'); [m,n]=size(Qx); %Qx for i=1:1:mfor j=1:1:nif j==nfprintf(fid,'%12.5f\n',Qx(i,j)); elsefprintf(fid,'%12.5f\t',Qx(i,j));endendendfprintf(fid,'\n');fprintf(fid,'输出单位权中误差u\n'); [m,n]=size(u); %u for i=1:1:m for j=1:1:nif j==nfprintf(fid,'%12.5f\n',u(i,j));elsefprintf(fid,'%12.5f\t',u(i,j));endendendfprintf(fid,'\n');fclose(fid);case 2% 条件平差BV+W=0sprintf('2:条件平差BV+W=0')n=input('请输入n=');t=input('请输入t=');B=input('请输入系数矩阵B=');W=input('请输入自由项向量W='); P=input('请输入权阵P='); sprintf('开始计算')K=-inv((B*inv(P)*B'))*W;V=inv(P)*B'*K;u=sqrt((V'*P*V)/(n-t));Ql=inv(P)-inv(P)*B'*inv(B*inv(P)*B')*B*inv(P);fid=fopen('data.txt','wt'); %写入文件路径，文件输出fprintf(fid,'=====================条件平差: BV+W=0 =====================\n'); fprintf(fid,'输出n,t\n') ;fprintf(fid,'%12.5f %12.5f\n',n ,t) ;%n,tfprintf(fid,'输出矩阵B\n');[m,n]=size(B); %Bfor i=1:1:mfor j=1:1:nif j==nfprintf(fid,'%12.5f\n',B(i,j));elsefprintf(fid,'%12.5f\t',B(i,j));endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵W\n');[m,n]=size(W); %Wfor i=1:1:mfor j=1:1:nif j==nfprintf(fid,'%12.5f\n',W(i,j));elsefprintf(fid,'%12.5f\t',W(i,j));endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵P\n'); [m,n]=size(P); %Pfor i=1:1:mfor j=1:1:nif j==nfprintf(fid,'%12.5f\n',P(i,j)); elsefprintf(fid,'%12.5f\t',P(i,j));endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵V\n'); [m,n]=size(V); %Vfor i=1:1:mfor j=1:1:nif j==nfprintf(fid,'%12.5f\n',V(i,j)); elsefprintf(fid,'%12.5f\t',V(i,j)); endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵Qx\n'); [m,n]=size(Ql); %Ql for i=1:1:m for j=1:1:nif j==nfprintf(fid,'%12.5f\n',Ql(i,j));elsefprintf(fid,'%12.5f\t',Ql(i,j));endendendfprintf(fid,'\n');fprintf(fid,'输出单位权中误差u\n');[m,n]=size(u); %u for i=1:1:mfor j=1:1:nif j==nfprintf(fid,'%12.5f\n',u(i,j));elsefprintf(fid,'%12.5f\t',u(i,j));endendendfprintf(fid,'\n');fclose(fid);case 3% 具有条件的参数平差% 参数平差V=AX-L% 条件平差BX+W=0sprintf('3:具有条件的参数平差V=AX-L;BX+W=0') n=input('请输入n=');t=input('请输入t=');r=input('请输入r=');A=input('请输入误差方程系数矩阵A=');L=input('请输入常数项矩阵L=');B=input('请输入条件方程系数矩阵B=');W=input('请输入自由项向量W=');P=input('请输入权阵P=');sprintf('开始计算')K=-(inv(B*inv(A'*P*A)*B'))*(W+B*inv(A'*P*A)*A'*P*L)X=(inv(A'*P*A))*(B'*K+A'*P*L)N=(A'*P*A)M=B*inv(N)*B'Qx=inv(N)-inv(N)*B'*inv(M)*B*inv(N)u=sqrt((V'*P*V)/(n-t+r))fid=fopen('data.txt','wt'); %写入文件路径，文件输出fprintf(fid,'=====================具有条件的参数平差: V=AX-L;BX+W=0 ===========================\n');fprintf(fid,'输出n,t,r\n') ;fprintf(fid,'%12.5f %12.5f %12.5f\n',n ,t,r); %n,t,rfprintf(fid,'输出矩阵A\n');[m,n]=size(A); %Afor i=1:1:mfor j=1:1:nif j==nfprintf(fid,'%12.5f\n',A(i,j));elsefprintf(fid,'%12.5f\t',A(i,j));endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵L\n');[m,n]=size(L); %Lfor i=1:1:mfor j=1:1:nif j==nfprintf(fid,'%12.5f\n',L(i,j)); elsefprintf(fid,'%12.5f\t',L(i,j)); endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵B\n'); [m,n]=size(B); %Bfor i=1:1:mfor j=1:1:nif j==nfprintf(fid,'%12.5f\n',B(i,j)); fprintf(fid,'%12.5f\t',B(i,j)); endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵W\n'); [m,n]=size(W); %W for i=1:1:m for j=1:1:nif j==nfprintf(fid,'%12.5f\n',W(i,j)); elsefprintf(fid,'%12.5f\t',W(i,j)); endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵P\n'); [m,n]=size(P); %P for i=1:1:m for j=1:1:nif j==nfprintf(fid,'%12.5f\n',P(i,j)); elsefprintf(fid,'%12.5f\t',P(i,j)); endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵X\n'); [m,n]=size(X); %X for i=1:1:m for j=1:1:nif j==nfprintf(fid,'%12.5f\n',X(i,j)); elsefprintf(fid,'%12.5f\t',X(i,j)); endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵V\n'); [m,n]=size(V); %Vfor j=1:1:nif j==nfprintf(fid,'%12.5f\n',V(i,j)); elsefprintf(fid,'%12.5f\t',V(i,j));endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵Qx\n'); [m,n]=size(Qx); %Qx for i=1:1:m for j=1:1:nif j==nfprintf(fid,'%12.5f\n',Qx(i,j)); elsefprintf(fid,'%12.5f\t',Qx(i,j)); endendendfprintf(fid,'\n');fprintf(fid,'输出单位权中误差\n'); [m,n]=size(u); %u for i=1:1:m for j=1:1:nif j==nfprintf(fid,'%12.5f\n',u(i,j));elsefprintf(fid,'%12.5f\t',u(i,j));endendendfprintf(fid,'\n');fclose(fid);case 4% 具有参数的条件平差% 具有参数的条件方程式BV+AX+W=0sprintf('4:具有参数的条件平差BV+AX+W=0 ')n=input('请输入n=');t=input('请输入t=');r=input('请输入r=');B=input('请输入条件方程中V的系数矩阵B=');A=input('请输入条件方程中X的系数矩阵A=');W=input('请输入条件方程自由向量W=');P=input('请输入权阵P=');sprintf('开始计算')N=(B*inv(P)*B')M=(A'*inv(N)*A)X=-inv(M)*A'*inv(N)*WK=-inv(N)*(A*X+W)V=inv(P)*B'*KQx=inv(M)u=sqrt((V'*P*V)/(r-t))fid=fopen('data.txt','wt'); %写入文件路径，文件输出fprintf(fid,'=====================具有参数的条件平差: BV+AX+W=0 ================================\n');fprintf(fid,'输出n,t,r\n') ;fprintf(fid,'%12.5f %12.5f %12.5f\n',n ,t,r) ; %n,t,rfprintf(fid,'输出矩阵B\n');[m,n]=size(B); %Bfor i=1:1:mfor j=1:1:nif j==nfprintf(fid,'%12.5f\n',B(i,j));elsefprintf(fid,'%12.5f\t',B(i,j)); endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵A\n'); [m,n]=size(A); %Afor i=1:1:mfor j=1:1:nif j==nfprintf(fid,'%12.5f\n',A(i,j)); elsefprintf(fid,'%12.5f\t',A(i,j)); endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵W\n'); [m,n]=size(W); %Wfor i=1:1:mfor j=1:1:nif j==nfprintf(fid,'%12.5f\n',W(i,j)); elsefprintf(fid,'%12.5f\t',W(i,j)); endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵P\n');[m,n]=size(P); %P for i=1:1:m for j=1:1:nif j==nfprintf(fid,'%12.5f\n',P(i,j)); elsefprintf(fid,'%12.5f\t',P(i,j)); endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵X\n'); [m,n]=size(X); %X for i=1:1:m for j=1:1:nif j==nfprintf(fid,'%12.5f\n',X(i,j)); elsefprintf(fid,'%12.5f\t',X(i,j)); endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵V\n'); [m,n]=size(V); %V for i=1:1:m for j=1:1:nif j==nfprintf(fid,'%12.5f\n',V(i,j)); elsefprintf(fid,'%12.5f\t',V(i,j)); endendendfprintf(fid,'\n');fprintf(fid,'输出矩阵Qx\n'); [m,n]=size(Qx); %Qx for i=1:1:m for j=1:1:nif j==nfprintf(fid,'%12.5f\n',Qx(i,j)); elsefprintf(fid,'%12.5f\t',Qx(i,j)); endendendfprintf(fid,'\n');fprintf(fid,'输出单位权中误差\n'); [m,n]=size(u); %u for i=1:1:m for j=1:1:nif j==nfprintf(fid,'%12.5f\n',u(i,j));elsefprintf(fid,'%12.5f\t',u(i,j));endendendfprintf(fid,'\n');fclose(fid);otherwisedisp('代号输入有误')end。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
chkdat函数（72页）function [n1,k]=chkdat(sd,pn,n1)n=length(n1);k=0;for i=1:ni1=0;for j=1:sdif(n1(i)==pn(j))i1=1;n1(i)=j;break;endendif(i1==0)% fprintf(fit2,'%5d %5d\n',i,n1(i)k=1;endendreturnreadlevelnetdata函数（73页）function [ed,dd,sd,gd,pn,h0,k1,k2,h1,s]=readlevelnetdata global filename filepath;global ed dd sd pn gd h0 k1 k2 h1 s k11 k12;k1=[];k2=[];h=[];s=[];[filename,filepath]=uigetfile('*.txt','选择高程数据文件');fid1=fopen(strcat(filepath,filename),'rt');if(fid1==-1)msgbox('Input File or Path is not correct','Warning','warn');return;ended=fscanf(fid1,'%f',1);dd=fscanf(fid1,'%f',1);sd=ed+dd;gd=fscanf(fid1,'%f',1);pn=fscanf(fid1,'%f',sd);h0=fscanf(fid1,'%f',ed);h0(dd+1:ed+dd)=h0(1:ed);heightdiff=fscanf(fid1,'%f',[4,gd]);heightdiff=heightdiff';k1=heightdiff(:,1);%起点k2=heightdiff(:,2);%终点k11=heightdiff(:,1);%起点k12=heightdiff(:,2);%终点h1=heightdiff(:,3);%高差s=heightdiff(:,4);%距离fclose('all');%点号转换[k1,k01]=chkdat(sd,pn,k1);[k2,k02]=chkdat(sd,pn,k2);h0(1:dd)=20000;ie=0;while(1)%计算近似高程for k=1:gdi=k1(k);j=k2(k);if(h0(i)<1e4&h0(j)>1e4)h0(j)=h0(i)+h1(k);ie=ie+1;endif(h0(i)>1e4&h0(j)<1e4)h0(i)=h0(j)-h1(k);ie=ie+1;endendif(ie==dd)break;endendh0=reshape(h0,length(h0),1); returnbm1函数（75页）function id=bm1(gd,dd,k1,k2)%计算一维压缩存放的数组idid=[];for i=1:ddk=i;for j=1:gdi1=k1(j);i2=k2(j);if(i1==i&i2<k)k=i2;endif(i2==i&i1<k)k=i1;endendid(i)=k;endfor i=2:ddid(i)=id(i-1)+i-id(i)+1;endreturn一维压缩存储法方程平差（76页）global pathname filenameglobal ed dd sd pn gd h0 k1 k2 h1 s dh;p=1./s;id=bm1(gd,dd,k1,k2);mm=id(dd);a(1:mm)=0;b(1:dd)=0;for k=1:gd %形成法方程i=k1(k);j=k2(k);h1=h1(k)+h0(i)-h0(j);if(i<=dd)ii=id(i)-i;a(ii+i)=a(ii+i)+p(k);b(i)=b(i)-h1*p(k);endif(j<=dd)jj=id(j)-j;a(jj+j)=a(jj+j)+p(k);b(j)=b(j)+h1*p(k);if(i<=dd)if(i>=j)a(ii+j)=a(ii+j)-p(k);elsea(jj+i)=a(jj+i)-p(k);endendendenda=gs5(dd,a,id);%变带宽下三角紧缩存储高斯消元法dh=cy6(a,b,id,dd,1);%常数项约化与回代子程序dh(dd+1:ed+dd)=0;hm(dd+1:ed+dd)=0;for i=1:sdh(i)=h0(i)+dh(i);endvv=0;for i=1:gdL(i)=h(k2(i))-h(k1(i));v(i)=h(k2(i))-h(k1(i))-h1(i);vv=vv+v(i)*v(i)/s(i);endu=sqrt(vv/(gd-dd));for i=1:ddb(1:dd)=0;b(i)=1.0;b=cy6(a,b,id,dd,i);hm(i)=sqrt(b(i))*u;endwritelevelnetdata(pn,k1,k2,h1,v',L',h0,dh',h',hm',u);gs5函数（77页）function a=gs5(dd,a,id)%变带宽高斯消去法for i=1:ddii=id(i)-i;if(i-1==0)li=1-ii;elseli=id(ii-1)-ii+1;endfor j=li:ijj=id(j)-j;if(j-1)==0lj=1-jj;elselj=id(j-1)-jj+1;endlk=li;if(li<lj)lk=lj;endfor k=lk:j-1kk=id(k);a(ii+j)=a(ii+j)-a(ii+k)/a(kk)*a(jj+k);endendendreturncy6函数（78页）function b=cy6(a,b,id,dd,k1)%常数项约化与回代子程序for i=k1:ddii=id(i)-i;if(i==1)nd=id(i);elsend=id(i)-id(i-1);ende=0;for k=1:i-1if((i-k)<nd)e=e+a(ii+k)*b(k);endendb(i)=(b(i)-e)/(a(ii+i);endfor i=dd-1:-1:k1ii=id(i);for k=i+1:ddkk=id(k)-k;nk=id(k)-id(k-1);if(k-i<nk)b(i)=b(i)-a(kk+i)/a(ii)*b(k);endendendreturn上三角存储法方程平差程序（79页）mm=(dd+1)*dd/2;a(1:mm)=0;b(1:dd)=0;for k=1:gdi=k1(k);j=k2(k);h1=h1(k)+h0(i)-h0(j);if(i<=dd)ii=(i-1)*(dd-i/2);a(ii+i)=a(ii+i)+1/s(k);b(i)=b(i)+1./s(k)*h1;endif(j<=dd)jj=(j-1)*(dd-j/2);a(jj+j)=a(jj+j)+1/s(k);b(j)=b(j)-1./s(k)*h1;if(i<=dd)if(i<j)a(ii+j)=a(ii+j)-1/s(k);elsea(jj+i)=a(jj+i)-1/s(k);endendendenda=invsqr(a,dd);for i=1:dddh(i)=0;di=(i-1)*(dd-i/2);for j=1:dddj=(j-1)*(dd-j/2);if(j<i)dh(i)=dh(i)-a(dj+i)*b(j);elsedh(i)=dh(i)-a(di+j)*b(i);endendenddh(dd+1:ed+dd)=0;hm(dd+1:ed+dd)=0;for i=1:sdh(i)=h0(i)+dh(i);endvv=0;for i=1:gdL(i)=h(k2(i))-h(k1(i));v(i)=h(k2(i))-h(k1(i))-h1(i);vv=vv+v(i)*v(i)/s(i);enduw0=sqrt(vv/(gd-dd));for i=1:ddii=(i-1)*(dd-i/2);hm(i)=sqrt(a(ii+i))*uw0;endreturn输出数据函数（79页）function writelevelnetdata(pn,k1,k2,h1,v,L,h0,dh,h,hm,uw0)disp('待定点高程平差值及中误差:')disp('---点号----近似高程(m)-高程改正(m)-高程平差值(m)-中误差')[pn,h0,dh,h,hm]disp('高差观测值平差值:')disp('---点号------点号----观测高差(m)---高差改正(m)-平差高差(m)')[pn(k1),pn(k2),h1,v,L][filename1,pathname1]=uigetfile('*.txt','请选择输出文件');fid2=fopen(strcat(pathname1,filename1),'wt');if(fid2==-1)msgbox('Error by Opening Output File','Warning','warn');return;endfprintf(fid2,'待定点高程平差值及中误差：\n 点号--近似高程(m)--高程改正(m)-高程平差值(m)-中误差\n');fprintf(fid2,'%5d %10.4f %10.4f %10.4f %10.4f\n',[pn,h0,dh,h,hm]');fprintf(fid2,'高差观测值平差值:\n -点号---点号--观测高差(m)--高差改正(m)-平差高差(m)\n'); fprintf(fid2,'%5d %5d %10.4f %10.4f %10.4f\n',[pn(k1),pn(k2),h1,v,L]');fprintf(fid2,'单位权中误差：%10.4fm\n',uw0);% open(strcat(pathname1,filename1));fclose(fid2);return利用Matlab矩阵运算的平差程序（81页）function level3ticdisp('平差已经开始---->>>>')global ed dd sd pn gd h0 k1 k2 h1 s dh;[ed,dd,sd,gd,pn,h0,k1,k2,h1,s]=readlevelnetdata;[dh,h,V,L,uw0,uwh,uwl]=calculatelevelnet(ed,dd,sd,gd,pn,h0,k1,k2,h1,s);writelevelnetdata(pn,k1,k2,h1,V,L,h0,dh,h,uwh,uw0); %输出水准网解算结果yunxing=toc;disp(['平差过程的运行时间为',num2str(yunxing),'秒。