基于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的水准网和测边网平差程序设计摘要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开发的计算机程序代码短小高效,能够快速的实现图形绘制和精度分析,与利用C++等高级编程语言开发的平差处理程序相比,可以在短时间内完成对平差的设计和验证.【期刊名称】《矿山测量》【年(卷),期】2016(044)004【总页数】3页(P63-65)【关键词】MATLAB;高程控制网;条件平差【作者】高霞;应建福【作者单位】甘肃工业职业技术学院测绘学院,甘肃天水741025;天水天正设计咨询有限公司,甘肃天水741000【正文语种】中文【中图分类】P224测绘科学是一门以大规模数据甚至是海量数据处理、分析与应用为基础的学科，其各项具体工作如测量计算、测量平差、大地测量数据处理、GPS数据结算、图像处理、坐标换算等都涉及大量的数值计算与分析。

运用MATLAB解算测绘信息数据处理问题具有其它程序设计语言难以比拟的优越性。

在高程控制网平差的计算中涉及大量的矩阵求逆和方程组的求解等问题，而MATLAB正是解决此类问题较好的数学软件。

对于一些常用的平差模型，可以用MATLAB软件通过调用其函数来完成，从而大大提高工作效率。

(1)简单易用的程序语言。

MATLAB是一个高级的矩阵/列阵语言，它包括控制语句、函数、数据结构、输入输出和面向对象编程特点。

用户可以在命令窗口中将输入语句与执行命令同步，也可以先编好一个较大的复杂的应用程序后一起运行。

而且这种语言可移植性较好，可拓展性极强，这也是MATLAB能够深入到科学研究及工程计算的重要原因。

(2)强大的科学计算、数据处理能力。

MATLAB是一个包含大量计算算法的集合，其拥有600多个工程中要用到的数学运算函数，可以方便的实现用户所需的各种计算功能。

在通常情况下，使用编程工作量会大大的减少。

其函数能解决的问题大致包括矩阵运算和线性方程组的求解、特征向量、微分方程及偏微分方程的求解、工程中的优化问题、多维数组操作以及建模动态仿真等。

实验三-利用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=。

观测数据读入程序，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初学者以及测绘专业的学生有用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。

基于Matlab的水准网间接平差程序设计赵亚红;周文国【摘要】设计水准网数据结构，存储在文本中，按照水准网的起点、终点、观测数据相对应关系建立矩阵，利用Matlab强大的矩阵运算功能，通过间接平差方法，按照最小二乘原理，求得任意水准网的未知点的最或然高程值，对平差结果输出存储，程序直观、简便。

并用实例验证了其正确性及通用性。

%On the basis of data structure designed, the relation of point s and lines of level net, the surveying data and the known data are st ored in the text, the matrixs were set up through the relation of the jupping -off points, end -points and the surveying data. And a program is designed in MATLAB to get the value of most probable by its strong abilty of calculating matrix, the result was output and stored. At last,the example proved the programme was right.【期刊名称】《华北科技学院学报》【年(卷),期】2011(008)003【总页数】3页(P58-60)【关键词】水准网;间接平差;Matlab【作者】赵亚红;周文国【作者单位】华北科技学院土木工程系,北京东燕郊101601;华北科技学院土木工程系,北京东燕郊101601【正文语种】中文【中图分类】P207.2水准网间接平差的的具体过程是：(1)根据水准网形进行分析，列误差方程;(2)根据误差方程系数列法方程;(3)解算法方程，求参数X及V;(4)求最或然值、精度评定。

第29卷　第4期2010年8月兰州交通大学学报J ou rnal of Lanzh ou Jiaotong UniversityV ol.29N o.4A ug.2010 文章编号:1001-4373(2010)04-0088-03基于M atlab的导线网平差程序设计＊李建章(兰州交通大学土木工程学院,甘肃兰州　730070)摘　要:导线网数据量大,网形复杂多变,其数据处理过程大多涉及到矩阵的计算.利用VC、VB等编程语言进行导线网程序的开发,算法比较复杂.基于M atlab平台,利用其强大的矩阵处理能力,设计出导线网数据结构,此基础上进行导线网平差程序的设计与开发,减小了代码编写的工作量.关键词:导线网;数据结构;平差;程序设计中图分类号:P209 文献标志码:A 导线网网形灵活多变,在城市测量中应用非常广泛.通常情况下,其外业观测数据量大,数据处理过程中大多涉及到矩阵的计算,且由于导线网网形的不确定性,因此其程序设计非常复杂[1].本文总结导线网的规律,设计出通用数据结构,并基于Matlab 强大的矩阵计算能力,编制了导线网数据处理程序. 1　导线网数据结构设计导线网由导线点、导线边和角度3类要素构成,其中导线边包括起点和终点,角度包括左边和右边.要使程序能对于任意形状的导线网进行处理,首先需要设计数据结构,以存储相关数据.这些数据包括起算数据、观测数据和网形各要素连接关系等,它们都是导线网各要素的属性值.本文用3个表来存储各要素,如表1-3所示.表1　点表数据结构Tab.1　Data structure of points序号字段类型备注1点名称整型2初始纵坐标浮点3初始横坐标浮点4已知点标志整型1为已知点,0为未知点5平差纵坐标浮点6平差横坐标浮点表2　角度数据结构Tab.2　Data structure of angles序号字段类型备注1角度编号整型2左边整型3右边整型4角度浮点度分秒表3　边数据结构Tab.3　Data structure of lines序号字段类型备注1导线边编号整型2起点整型3终点整型4边长浮点5方位浮点弧度6纵坐标增量浮点7横坐标增量浮点 表1-3中分别存储在ptTab,lineTab,angleTab 矩阵中,保存为.mat文件,程序运行时加载[2].2　近似坐标计算近似坐标的计算是导线网平差中关键的一个环节.其精度直接影响到后续平差计算的点位精度和迭代平差工作量大小.近似坐标计算包括近似方位角的计算和近似坐标的计算两个步骤.2.1　近似方位角计算近似方位角的计算以角度为单位,将已知方位传递到网中每一条边.设某角度一边方位已知,而另一边方位未知,由于两边夹角已知,可计算出未知边的方位.图1所示为4种可能情况.图1　方位角计算4种情况Fig.1　Four situation of azimuth angle calculation＊收稿日期:2010-04-06作者简介:李建章(1974-),男,甘肃会宁人,讲师.第4期李建章:基于M atlab 的导线网平差程序设计假定已知左边方位角为fw 1,夹角为α,则以上4种情况下右边的方位角f w r 讨论如下.情况一:fw r =fw 1+α(1)情况二:fw r =fw 1+α(2)情况三:fw r =fw 1±π+α(3)情况四:fw r =fw 1±π+α(4)同理可得,已知右边方位角计算左边方位角的情况,也有4种可能性.程序根据角度两边的端点点名的关系判断以上8种情况,采用相应的计算公式计算未知边方位角[3].程序在获得未知边方位角后,将计算结果保存到边表相应记录中.然后在角度表中搜索相邻角度,搜索条件是:该角度的一边必须是上一角度的一边,而另一边不是上一角度的一边.查询到满足条件的角度后,判断其是否为截至角(两边方位已知),如否则计算出该角度未知边方位,重复前面的步骤直至某一截至角停止.然后在边表中查询有无近似方位未知的边,如有,再次执行以上步骤,直至边表中所有边近似方位计算完毕,这个过程可以通过一个函数自身迭代来实现.程序流程如图2所示.图2　计算近似方位流程图Fig .2　Flow chart of calculating azimuth angle计算导线边方位的子程序如下所示:functio n [ptT ab ,angleT ab ,lineT ab ,ok ]=g etfw0(pt -T ab ,ang le T ab ,lineT ab )[baindex ]=getanglebeg in (angleT ab ,line Tab );%查找起算角[ok ,ptT ab ,line T ab ,ang le T ab ]=caculateFW0(baindex ,ptT ab ,line Tab ,ang le Tab );%由起算角往前传递方位.if no t (fwisok (lineTab ))%判断有没有方位未知的导线边　[ptT ab ,ang leT ab ,lineT ab ,ok ]=getfw 0(ptT ab ,a n -gleT ab ,lineT ab );　%函数迭代计算.end end2.2　近似坐标计算通过2.1计算,所有导线边的近似方位计算完毕,此时可以利用每条边的边长和近似方位计算其坐标增量,这个过程只需要在边表中循环计算即可.然后以导线边为单位,从起始边出发,将已知坐标传递到各未知点.所谓起始边即该边一端点坐标已知,另一端点坐标未知.利用导线边坐标增量计算未知点坐标,然后查找相邻边,判断其是否为截至边(两端点坐标皆知),如否,计算未知点坐标.重复查找计算直至截至边,然后程序在点表中判断有无近似坐标未知的点,如有,则重复以上步骤,否则程序退出.这个过程和2.1中计算方位角的过程是类似的,因此不再列出程序流程图和代码.3　误差方程系数矩阵计算导线网观测值有角度和边长两种类型,一个观测值可列出一个误差方程.因此程序需分别读取角度表和导线边表中每个记录来计算误差方程系数矩阵B 、常数项矩阵l 和权阵P .3.1　利用观测角度计算误差方程系数矩阵设某角度观测值为ang le ,其左边近似方位、左边长、左边起点点名,左边起点点序号,左边终点点名,左边终点点序号,左边坐标增量分别为:lfw ,ls ,lbname ,lbindex ,lename ,leindex ,ldetx ,ldety ,同理对应右边各项为:rfw ,rs ,rbname ,rbindex ,rename ,reindex ,rdetx ,rdety .由图1可知,各观测角度两边的方向有4种情况,为简化程序,首先将其全部转化为情况1.设转换后左导线边起、终点序号分别为tem lbindex ,tem -leindex .右导线边起、终序点号分别为tem rbindex ,temreindex .如下是图1第三种情况的转换代码.if lename ==rbname temlbindex =leinde x ; temleindex =lbinde x ; tem rbindex =rbindex ; tem reindex =reindex ; lde tx =(-1)＊lde tx ; lde ty =(-1)＊lde ty ; rdetx =rdetx ;89兰州交通大学学报第29卷 rdety =rdety ; lfw =lfw +pi (); rfw =rfw ;end设该角度序号为index .又左边起点、终点和右边起点、终点的未知数序号分别为:lbcontrindex ,le -contrindex ,rbco ntrindex ,recontrindex .则B (index ,(2×leco ntrindex -1))=ρΔy L 1000s 2L B (index ,(2×leco ntrindex ))=-ρΔx L1000s 2L B (index ,(2×reco ntrindex -1))=-ρΔy R1000s 2RB (index ,(2×reco ntrindex ))=-ρΔx R1000s 2RB (index ,(2×lbcontrindex -1))=ρΔy R 1000s 2R -ρΔy L1000s 2L B (index ,(2×lbcontrindex ))=-ρΔx R 1000s 2R +ρΔx L1000s 2L angle 0=r f w -l f wl (index ,1)=(angle 0-ang le )×180×3600πP (index ,index )=1[4]上述各点如为已知点,则对应系数为0.未知点序号为所有待求点的排列序号,用于控制误差方程系数在B 矩阵中的列位置.3.2　利用观测边长计算误差方程系数矩阵设某观测边坐标增量为Δx ,Δy ,观测边长为S ,S 0=Δx 2+Δy 2,则B (index ,(2×lbco ntrindex -1))=-ΔxS 0(5)B (index ,(2×lbco ntrindex ))=-ΔyS 0(6)B (index ,(2×leco ntrindex -1))=ΔxS 0(7)B (index ,(2×lecontrindex ))=ΔyS 0(8)l (index ,1)=1000×(S 0-S )(9)P (index ,index )=S100[5](10)其中已知点对应项系数为0,其中index 为观测边序号,观测边序号是该边在边表中序号加上角度观测值个数.lbcontrindex 和leco ntrindex 为该边起点、终点的未知点序号.以上过程需要一个循环语句即可完成.4　结束语本文基于Matlab 开发了导线网平差程序,算法简单,编程工作量小[6].对于测绘专业类似计算问题如三角网、水准网等的程序设计有一定的借鉴意义.由于篇幅的限制,本文没有论述其他相关处理过程的算法,包括利用误差方程迭代计算获取最优解、度分秒和弧度互化、边点要素属性的查询等,这些算法在M atlab 中实现相对比较容易.参考文献:[1]　汪自军,陈圣波,臧立娟,等.导线网数据处理系统关键技术及其实践[J ].微计算机信息,2008(2):216-218.[2]　李建章.基于M a tlab 的水准网平差程序设计[J ].兰州交通大学学报,2009,28(3):29-31.[3]　武艳强,黄立人,江在森.导线网平差中近似坐标的无限定推算方法[J ].测绘通报,2006(12):12-15.[4]　武汉大学测绘学院测量平差学科组.误差理论与测量平差基础[M ].武汉:武汉大学出版社,2003:102-125.[5]　姚德新.土木工程测量学教程[M ].北京:中国铁道出版社,2003:67-68.[6]　薄志义,曹福生.程序寻找支导线网计算路径的研究[J ].测绘科学,2007(9):68-69.Adjustment Programming of Traverse Network on the Basis of MatlabLI Jian -zhang(Sch ool of Civil En gineering ,Lanzhou Jiaotong University ,Lanzh ou 730070,China )A bstract :H aving large data and complex shape ,the data pro cessing of traverse netw o rk involves matrix cal -culation .It 's mo re difficult to prog ram on o ther prog ramming lang uages like VC and VB .In the e ssay ,the data structure o f traverse netw ork is desig ned by utilizing the matrix dispo sal capability of M atlab ,and thenthe prog ram that has the adv antag es of being simple in algo rithm is finished .Key words :trave rse ne tw ork ;data structure ;adjustment ;pro gramming90。