基于MATLAB的控制网平差程序设计--第六章源代码

合集下载

MATLAB在测量控制网优化设计与平差中的应用

测绘出版社，１９９５［２］国家测绘局．全球定位系统（ＧＰＳ）测量规范．北京：
测绘出版社，１９９２
第一作者简介詹家民，男，测绘高级工程师，现任教于徐州师范大学工学院，曾在《测绘通报》、《矿山测量》《江苏测绘》等刊物上发表论文多篇。
参考文献［１］刘基余，李征航等．全球定位系统原理及应用．北京：
MATLAB 有两种常用的工作方式：一种是直接交互的
观测数据增删和观测权改变的依据。在这一过程中，其误指令行操作方式；另一种是M 文件的编程方式。在前一种
差椭圆的绘制需要将测量坐标系化为数学坐标系，并以各工作方式下，MATLAB 被当作一种高级的数学演算和计算
点为中心的坐标系绘制误差椭圆，其函数模型为：
幕供设计人员修改。该过程有两种：a 调整观测值；b 调
题之一。根据修改后的数据进一步进行优化，使其达到所整观测值的权。反复调整到符合要求为止，打印输出观测
要求的技术指标。
表和设计图。
绘制控制网图与误差椭圆图是该阶段的主要任务，通
（5）MATLAB 部分源程序：
过绘制的网图进行控制网设计的调整，并将修改结果作为
第一作者简介左廷英，女，硕士，现在中南大学信息地理工程学院从事教学工作，发表多篇论文。
（收稿日期：2003年1月6日）
63
况下，误差方程中的系数A 产生列亏。这种没有足够起始数据的平差问题，称为秩亏自由网平差。
可根据最小二乘准则，组成法方程： ATPAX^ = ATPL
即 N = ATPA 则
（2）
协因数阵为： Q xx = N－ 1 = (ATPA )－ 1 （3）
由于系数阵 A 列亏 , 参数的最小二乘解就不唯一。为此

基于MATLAB的控制网平差程序设计--第四章源代码

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),'秒。

实验三-利用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初学者以及测绘专业的学生有用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 编程程序

导线网严密平差完整程序代码及注释已知数据已知点坐标：XA=730.024 ，YA=126.040 ，111.855=XB，232.172=Y B应得结果 = 5.30表七：实际结果点号X(m) Y(m) 中误差X(mm) 中误差Y(mm) 点位中误差1 678.1641 287.3411 3.3868 5.8737 7.60782 564.6919 446.8827 6.0102 7.5137 9.19373 533.8468 629.2749 8.7985 7.3927 10.05964 475.6776 535.8728 7.3949 8.4987 9.96675 499.1408 374.4231 5.8076 7.6840 9.18276 594.6190 176.7976 5.2308 3.4068 7.34757 826.3848 305.7569 2.7439 5.4473 7.15518 767.8778 401.8510 4.8371 6.3810 8.37349 745.7230 531.0410 6.8992 6.6276 9.194710 671.6256 650.8637 9.0695 6.0533 9.722011 794.2833 600.4688 8.4994 6.1800 9.578412 859.0528 511.5349 6.7399 6.6874 9.160813 898.2695 399.0763 4.7580 6.4207 8.358614 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:28det(i)=sqrt(M(i,i));end% 点位中误差for i=1:14ddet(i)=2.5*sqrt(det(i*2-1)+det(i*2));endddet=ddet';% 整理数据，输出最终结果for i=1:14XX(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);endjieguo(:,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:24S3(i,1)=S1(i,1);S3(i,2)=(S1(i,2)+S1(i,3)/60+S1(i,4)/3600);End====================================================% 计算坐标近似值function [X0]=zbjs(S2,S3)% 起始数据X10=[730.024;126.040;855.111;172.232];X0=zeros(1,28);% 预分配内存% 已知边方位角t=180+atand((X10(2)-X10(4))/(X10(1)-X10(3)));% 下方闭合% 起始方位角t1=t;for i=1:1:6% 计算各边方位角if S2(i,5)>0 % i=4t1=t1+S3(S2(i,4),2)+S3(S2(i,5),2)-180;% 边4方位角+角6+角7-180 elset1=t1+S3(S2(i,4),2)-180;endif t1<0t1=t1+360;end% 计算各点近似坐标if i<2X0(1)=X10(1)+S2(1,2)*cosd(t1);X0(2)=X10(2)+S2(1,2)*sind(t1);endif i>=2X0(2*i-1)=X0(2*(i-1)-1)+S2(i,2)*cosd(t1);X0(2*i)=X0(2*(i-1))+S2(i,2)*sind(t1);endend% 上方闭合% 起始方位角t2=t;t2=t2-S3(12,2);if t2<0t2=t2+360;endfor i=8:1:15% 计算个边方位角if i<9 % i=8X0(2*(i-1)-1)=X10(3)+S2(i,2)*cosd(t2);X0(2*(i-1))=X10(4)+S2(i,2)*sind(t2);endif i>8t2=t2+S3(S2(i,4),2)-180;if t2<0t2=t2+360;end% 计算各点近似坐标X0(2*(i-1)-1)=X0(2*(i-2)-1)+S2(i,2)*cosd(t2); X0(2*(i-1))=X0(2*(i-2))+S2(i,2)*sind(t2);endend==================================================== % 计算方位角function [M]=fangweijiao(a,b,X0)X10=[730.024;126.040;855.111;172.232];% 确定边对应两点类型及坐标。

用MATLAB解决-条件平差和间接平差

B L2
L1 A
L3 C
9
clc Disp(‘条件平差示例２’) Disp(‘三角形内角观测值’) L1 = [42 12 20] L2 = [78 9 9] L3 = [59 38 40] L = [L1; L2; L3] Disp(‘将角度单位由度分秒转换为弧度’) LL = dms2rad(mat2dms(L))
设误差Δ和参数X的估计值分别为V 和 Xˆ
15
则有
V AXˆ l
为了便于计算，通常给参数估计一个充分接近的近似值 X 0
Xˆ X 0 xˆ 则误差方程表示为
V Axˆ l
其中常数项为
l L (AX 0 d)
16
由最小二乘准则，所求参数的改正数应该满足
V T PV min
目标函数对x求一阶导数，并令其为零
if(sum(LL) == pi) disp(‘检核正确’)
else
V = A'*Ka
disp(‘检核错误’) end
11 例《误差理论与测量平差基础》P75
在下图中，Ａ、Ｂ为已知水准点，其高程为 HA=12.013m, HB = 10.013m, 可视为无误差。为了确定点Ｃ及Ｄ点的高程，共观测了四个高差，高差观测值及相应的水准路线的距离为：
disp(‘系数矩阵B’)
B = [1 0; 0 1; 1 0; 0 1; -1 1; -1 0]
l = [0; 0; 4; 3; 7; 2]
disp(‘C是单位权观测高差的线路公里数，S是线路长度’)
C = l*ones(1,6)
23
h6
E h3
h7 B
S = [1.1, 1.7, 2.3, 2.7, 2.4, 4.0]

基于Matlab的水准网间接平差程序设计

基于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)求最或然值、精度评定。

MATLAB平差程序

function BJPC2ZB=[];ZJ=[];BC=[];%%format long[filename,filepath]=uigetfile('*.txt','选择平差文件:'); name=[filepath filename];fid1=fopen(name,'rt');if(fid1==-1)disp('文件未打开，请重试！');return;endn=fscanf(fid1,'%f',1);%输入边的个数b=fscanf(fid1,'%f',1);%输入多余观测个数E=fscanf(fid1,'%f',1);%输入测角中误差t=fscanf(fid1,'%f',1);%输入坐标个数ZB=fscanf(fid1,'%f',[2 t]);dms1=fscanf(fid1,'%f',[3 n+1]);BC=fscanf(fid1,'%f',[1 n]);T0=fscanf(fid1,'%f',1);T1=fscanf(fid1,'%f',1);fclose(fid1);dms=dms1';ZJ=dms2degrees(dms);%%if t==2T0=T0;%若告诉起始方位角就直接输入endif t==4 %若没有告诉起始方位角由坐标反算Xa=ZB(1,1);Xb=ZB(1,2);Ya=ZB(2,1);Yb=ZB(2,2);Tx=Xb-Xa; Ty=Yb-Ya;T0=atan(Ty/Tx)*180/pi;%对方位角的讨论if ((Tx>0)&&(Ty>0))T0;endif (((Tx<0)&&(Ty>0))||((Tx<0)&&(Ty<0)))T0=T0+180;endif ((Tx>0)&&(Ty<0))T0=T0+360;endif ((Tx==0)&&(Ty>0))T0=90;endif ((Tx==0)&&(Ty<0))T0=270;endif((Ty==0)&&(Tx>0))T0=0;endif ((Ty==0)&&(Tx<0))T0=180;endenda0=T0;x0=ZB(1,t/2);y0=ZB(2,t/2);A=[0,0];FWJ=[];%开始计算近似方位角J=1;while(J<=n+1)belta=ZJ(J);%讨论方位角if J==1a=a0+belta;elsea=a+belta;endif a>=180a=a-180;elsea=a+180;endif a>=360a=a-360;endFWJ(J)=a;J=J+1;endFWJ;%%J=1;%开始计算近似坐标Awhile(J<=n)if J==1A(1,1)=x0+BC(1)*cos(FWJ(1)*pi/180);A(2,1)=y0+BC(1)*sin(FWJ(1)*pi/180);elseA(1,J)=A(1,J-1)+BC(J)*cos(FWJ(J)*pi/180); A(2,J)=A(2,J-1)+BC(J)*sin(FWJ(J)*pi/180);endJ=J+1;endA;%%W=[];%开始计算闭合差if t==2T1=T1;endif t==4Xa=ZB(1,3);Xb=ZB(1,4);Ya=ZB(2,3);Yb=ZB(2,4);Tx=Xb-Xa; Ty=Yb-Ya;T1=atan(Ty/Tx)*180/pi;if ((Tx>0)&&(Ty>0))T1;endif (((Tx<0)&&(Ty>0))||((Tx<0)&&(Ty<0)))T1=T1+180;endif ((Tx>0)&&(Ty<0))T1=T1+360;endif ((Tx==0)&&(Ty>0))T1=90;endif ((Tx==0)&&(Ty<0))T1=270;endif((Ty==0)&&(Tx>0))T1=0;endif ((Ty==0)&&(Tx<0))T1=180;endendW(1,1)=-(FWJ(n+1)-T1)*3600;%以秒为单位W(2,1)=-(A(1,n)-ZB(1,t/2+1))*100;%以厘米为单位W(3,1)=-(A(2,n)-ZB(2,t/2+1))*100;W;%%XS=zeros(3,2*n+1);%开始计算系数阵，XS是系数阵e=n+1;while(e<=2*n+1)XS(1,e)=1;e=e+1;endf=1;while(f<=n)XS(2,f)=cos(FWJ(f)*pi/180);%纵坐标边长改正数的系数 XS(3,f)=sin(FWJ(f)*pi/180);%横坐标边长改正数的系数 f=f+1;endu=1;y0=A(2,n);x0=A(1,n);Y=[];X=[];y1=[];x1=[];while(u<=n+1)if u==1Y(1)=ZB(2,t/2);X(1)=ZB(1,t/2);elseY(u)=A(2,u-1);X(u)=A(1,u-1);endu=u+1;endy1=-1/2062.65*(y0-Y);%纵坐标转角改正数的系数x1=1/2062.65*(x0-X);%横坐标转角改正数的系数y1;x1;d=n+1;z=1;while(d<=2*n+1)XS(2,d)=y1(1,z);XS(3,d)=x1(1,z);d=d+1;z=z+1;endXS;%%r=[];h=1;while (h<=n)s=(5+5*BC(h)/1000)/10;%求中误差厘米为单位r(h)=E^2/s^2;%计算距离的权h=h+1;endaa=[];y=1;while(y<=n)aa(y)=r(y);%将距离的权赋给矩阵aay=y+1;endc=n+1;while(c<=2*n+1)aa(c)=1;%转角的权都是1c=c+1;endP=diag(aa);%形成权阵P;%%Q=inv(P);N=XS*Q*XS';N1=inv(N);K=inv(N)*W;V=Q*XS'*K;V;BC1=BC'+V(1:n)/100;%将改正数加到观测值以求平差边长（除一百统一单位）ZJ1=ZJ+V(n+1:2*n+1)/3600;%同上BC1;i=1;while (i<=n+1)ZJ2(i,1)=fix(ZJ(i));ZJ2(i,2)=fix((ZJ(i)-fix(ZJ(i)))*60);ZJ2(i,3)=((ZJ(i)-fix(ZJ(i)))*60-fix((ZJ(i)-fix(ZJ(i)))*60))*60;i=i+1;endZJ2;%%A1=[0,0];FWJ1=[];J=1;while(J<=n+1)belta=ZJ1(J);if J==1a=a0+belta;elsea=a+belta;endif a>=180a=a-180;elsea=a+180;endif a>=360a=a-360;endFWJ1(J)=a;J=J+1;endi=1;while (i<=n+1)FWJ3(i,1)=fix(FWJ(i));FWJ3(i,2)=fix((FWJ(i)-fix(FWJ(i)))*60);FWJ3(i,3)=((FWJ(i)-fix(FWJ(i)))*60-fix((FWJ(i)-fix(FWJ(i)))*60))*60;FWJ2(i,1)=fix(FWJ1(i));FWJ2(i,2)=fix((FWJ1(i)-fix(FWJ1(i)))*60);FWJ2(i,3)=((FWJ1(i)-fix(FWJ1(i)))*60-fix((FWJ1(i)-fix(FWJ1(i)))*60))*60; i=i+1;endFWJ3;FWJ2;%%J=1;%开始计算平差坐标A1while(J<=n)if J==1A1(1,1)=ZB(1,t/2)+BC1(1)*cos(FWJ1(1)*pi/180);A1(2,1)=ZB(2,t/2)+BC1(1)*sin(FWJ1(1)*pi/180);elseA1(1,J)=A1(1,J-1)+BC1(J)*cos(FWJ1(J)*pi/180);A1(2,J)=A1(2,J-1)+BC1(J)*sin(FWJ1(J)*pi/180);endJ=J+1;endA1;m=sqrt(V'*P*V/b);%测量精度单位权中误差。

导线网平差程序源代码

程序源代码C语言程序：#include <windows.h>#include <stdarg.h>#include <stdio.h>#include <stdlib.h>#include<math.h>#define PI 3.1415926535898#define p 206264.806247#define MAX 50//矩阵的乘法运算MatrixMutiply(double Matrix1[MAX][MAX],double Matrix2[MAX][MAX],double MatrixResult[MAX][MAX],int m1,int m2,int m3){int i,j,k;double Sum;for(i=0;i<m1;i++){for(j=0;j<m3;j++){/*按照矩阵乘法的规则计算结果矩阵的i*j元素*/Sum=0;for(k=0;k<m2;k++)Sum+=Matrix1[i][k]*Matrix2[k][j];MatrixResult[i][j]=Sum;}}//return MatrixResult;}//矩阵的转置运算MatrixT(double Matrix1[MAX][MAX],double T[MAX][MAX],int m1,int m2) {//m1,m2分别是矩阵的行列数int i,j;//double T[50][50];for(i=0;i<m1;i++){for(j=0;j<m2;j++){T[j][i]=Matrix1[i][j];}}//return T;}//矩阵的逆运算void swap(double a,double b){double c=a;a=b;b=c;}double DinV(double A[50][50],int n) //n代表阶数{int i,j,k;double d;int JS[50],IS[50];for (k=0;k<n;k++){d=0;for (i=k;i<n;i++)for (j=k;j<n;j++){if (fabs(A[i][j])>d){d=fabs(A[i][j]);IS[k]=i;JS[k]=j;}}if (d+1.0==1.0)return 0;if (IS[k]!=k)for (j=0;j<n;j++)swap(A[k][j],A[IS[k]][j]);if (JS[k]!=k)for (i=0;i<n;i++)swap(A[i][k],A[i][JS[k]]);A[k][k]=1/A[k][k];for (j=0;j<n;j++)if (j!=k) A[k][j]=A[k][j]*A[k][k];for (i=0;i<n;i++)if (i!=k)for (j=0;j<n;j++)if (j!=k) A[i][j]=A[i][j]-A[i][k]*A[k][j];for (i=0;i<n;i++)if (i!=k) A[i][k]=-A[i][k]*A[k][k]; }for (k=n-1;k>=0;k--){for (j=0;j<n;j++)if (JS[k]!=k) swap(A[k][j],A[JS[k]][j]);for (i=0;i<n;i++)if (IS[k]!=k)swap(A[i][k],A[i][IS[k]]);}}//将度分秒连写的角度化为弧度值double jd_hd(double D,double F,double M){//int dd=(int)((int)B/10000); //提取度值//int ff=(int)(((int)B-dd*10000)/100); //提取分值//double mm=B-dd*10000-ff*100;//提取秒值double B;B=(D*3600.0+F*60.0+M)/p;//角度化弧度return B;}main(){double D,F,M,sigma_beta=2.0,sigma_s;//scanf("%f%f%f",&D,&F,&M);//printf("%f",jd_hd(D,F,M));doublebeta[50],alf[50],alfo[50],s[50],so[50],Xo[50],Yo[50],B[50][50]={0.0},L[50][50],P[50][50 ]={0.0},c[50][50]={0.0},wx=0.0;//此处为Xo Yo，B矩阵赋初值为零beta代表夹角，alf方位角，alfo方位角近似值，s距离观测值，so距离近似值//Xo Yo坐标近似值，B[50][50]误差矩阵，L[50]为L矩阵，P[50][50]为权阵,c[1][30]代表限制条件的系数阵,w代表限制条件常数项doubleNbb[50][50],Ncc[50][50],W[50][50],Ks,xgu[50][50],Xgu[30][1],Ygu[30][1],V[50][50],sigma_gu,Q[50][50],sigma_xy[50][50];alf[0]=PI;alfo[0]=PI;int i,j,m1,m2,m3;//将测量数据导入，并存入相应数组FILE *fp3;char strline[100]; //读取文件每行的bufferint du[100],fen[100],miao[100];double bian[100];i=0,j=1;if((fp3=fopen("D:\\111\\测量数据.txt","r"))==NULL) //文件位置和文件名{printf("文件不存在！");return 0;}while(!feof(fp3)) //判断文件是否已到末尾{fgets(strline,100,fp3); //读取一行sscanf(strline,"%d %d %d %lf",&du[i],&fen[i],&miao[i],&bian[i]); //从文件读取到的一行数据分别存放在两个数组中i++;}fclose(fp3);while(1){//printf("\n%d\t%d\t%d\t%lf",du[j],fen[j],miao[j],bian[j]);beta[j-1]=jd_hd(du[j],fen[j],miao[j]);s[j-1]=bian[j];j++;if(j>=i){break;}}Xo[0]=5000.0;Yo[0]=5000.0;Xo[15]=5000.0;Yo[15]=5000.0;so[0]=s[0];//用来求未知点坐标近似值for(i=1;i<15;i++){alf[i]=alf[i-1]+beta[i]-PI;if(alf[i]>=(2*PI)){alf[i]=alf[i]-2*PI;}Xo[i]=Xo[i-1]+s[i-1]*cos(alf[i-1]);Yo[i]=Yo[i-1]+s[i-1]*sin(alf[i-1]);//printf("X=%f\t",Xo[i]);}for(i=1;i<15;i++){//求近似距离so[i]=sqrt((Yo[i+1]-Yo[i])*(Yo[i+1]-Yo[i])+(Xo[i+1]-Xo[i])*(Xo[i+1]-Xo[i]));//求近似方位角，分象限进行讨论if((Yo[i+1]-Yo[i])>0&&(Xo[i+1]-Xo[i])>0)//第一象限{alfo[i]=atan((Yo[i+1]-Yo[i])/(Xo[i+1]-Xo[i]));}else if((Yo[i+1]-Yo[i])>0&&(Xo[i+1]-Xo[i])<0)//第二象限{alfo[i]=atan((Yo[i+1]-Yo[i])/(Xo[i+1]-Xo[i]))+PI;}else if((Yo[i+1]-Yo[i])<0&&(Xo[i+1]-Xo[i])>0)//第三象限{alfo[i]=atan((Yo[i+1]-Yo[i])/(Xo[i+1]-Xo[i]))+2*PI;}else //((Yo[i]-Yo[i-1])<0&&(Xo[i]-Xo[i-1])<0)//第四象限{alfo[i]=atan((Yo[i+1]-Yo[i])/(Xo[i+1]-Xo[i]))+PI;}//printf("alf=%f\t",alfo[i]);//printf("so=%f\t",so[i]);}//求B矩阵//将第一个角度的系数单独算出j=0;//B[0][j]=0.0;((Yo[14]-Yo[0])/(so[14]*so[14]))*p/1000.0;//B[0][j+1]=0.0;-1*((Xo[14]-Xo[0])/(so[14]*so[14]))*p/1000.0;B[0][j]=-1*(Yo[1]-Yo[0])/(so[0]*so[0])*p/1000.0;//B[0][j+1]=(Xo[1]-Xo[0])/(so[0]*so[0])*p/1000.0;B[0][26]=((Yo[14]-Yo[0])/(so[14]*so[14]))*p/1000.0;B[0][27]=-1*(Xo[14]-Xo[0])/(so[14]*so[14])*p/1000.0;//将第二个角度的系数单独算出B[1][j]=((Yo[2]-Yo[1])/(so[1]*so[1])-(Yo[0]-Yo[1])/(so[0]*so[0]))*p/1000.0;//B[1][j+1]=-1*((Xo[2]-Xo[1])/(so[1]*so[1])-(Xo[0]-Xo[1])/(so[0]*so[0]))*p/1000.0 ;B[1][j+2]=-1*((Yo[2]-Yo[1])/(so[1]*so[1]))*p/1000.0;B[1][j+3]=((Xo[2]-Xo[1])/(so[1]*so[1]))*p/1000.0;//求其他角度改正的系数for(i=2;i<15;i++){if(i<14){B[i][j]=((Yo[i-1]-Yo[i])/(so[i-1]*so[i-1]))*p/1000.0;B[i][j+1]=-1*((Xo[i-1]-Xo[i])/(so[i-1]*so[i-1]))*p/1000.0;B[i][j+2]=((Yo[i+1]-Yo[i])/(so[i]*so[i])-(Yo[i-1]-Yo[i])/(so[i-1]*so[i-1]))*p/1000.0;B[i][j+3]=-1*((Xo[i+1]-Xo[i])/(so[i]*so[i])-(Xo[i-1]-Xo[i])/(so[i-1]*so[i-1]))*p/1000. 0;B[i][j+4]=-1*((Yo[i+1]-Yo[i])/(so[i]*so[i]))*p/1000.0;B[i][j+5]=(Xo[i+1]-Xo[i])/(so[i]*so[i])*p/1000.0;}else{B[i][j]=(Yo[i-1]-Yo[i])/(so[i-1]*so[i-1])*p/1000.0;B[i][j+1]=-1*((Xo[i-1]-Xo[i])/(so[i-1]*so[i-1]))*p/1000.0;B[i][j+2]=((Yo[i+1]-Yo[i])/(so[i]*so[i])-(Yo[i-1]-Yo[i])/(so[i-1]*so[i-1]))*p/1000.0;B[i][j+3]=-1*((Xo[i+1]-Xo[i])/(so[i]*so[i])-(Xo[i-1]-Xo[i])/(so[i-1]*so[i-1]))*p/1000. 0;}j=j+2;}B[2][1]=0.0;//求边长改正的系数j=0;//将第一个边长的系数单独算出B[i][j]=(Xo[1]-Xo[0])/so[0];//B[i][j+1]=(Yo[1]-Yo[0])/so[0]; //(Yo[1]-Yo[0])/(so[0]*so[0]); i=i+1;for(i;i<30;i++){if(i<29){B[i][j]=-1*(Xo[i-14]-Xo[i-15])/so[i-15];B[i][j+1]=-1*(Yo[i-14]-Yo[i-15])/so[i-15];B[i][j+2]=-1*B[i][j];B[i][j+3]=-1*B[i][j+1];}else{B[i][j]=-1*(Xo[i-14]-Xo[i-15])/so[i-15];B[i][j+1]=-1*(Yo[i-14]-Yo[i-15])/so[i-15];}j=j+2;}B[16][1]=0.0;for(j=1;j<27;j++){for(i=0;i<30;i++){B[i][j]=B[i][j+1];}}//求L矩阵L[0][0]=(beta[0]-(alfo[0]-alfo[14]+PI))*p;for(i=1;i<30;i++){if(i<15){L[i][0]=(beta[i]-(alfo[i]-alfo[i-1]+PI))*p;if(L[i][0]>PI*p){L[i][0]=(L[i][0]-2*PI*p);}L[i][0]=L[i][0];}else{L[i][0]=(s[i-15]-so[i-15])*1000.0;}}//求权阵Pfor(i=0;i<30;i++){if(i<15){P[i][i]=1;}else{sigma_s=5+10*0.000001*s[i-15]*1000; //单位为（''/mm）的平方P[i][i]=sigma_beta*sigma_beta/(sigma_s*sigma_s);}}double Tb[50][50],Tc[50][50],MatrixResult[50][50],TV[50][50];doubletemp1[MAX][MAX],temp2[MAX][MAX],temp3[MAX][MAX],temp4[MAX][MAX],temp5[MAX][MAX];//计算Nbb矩阵,W矩阵m1=30;m2=28;MatrixT(B,Tb,m1,m2);m3=30;MatrixMutiply(Tb,P,temp1,m2,m1,m3);m3=m2;m1=m2;m2=30;MatrixMutiply(temp1,B,Nbb,m1,m2,m3);m1=27;m2=30;m3=1;MatrixMutiply(temp1,L,W,m1,m2,m3);//矩阵输出D:\\111\\导线网输出数据3027.txt文本FILE *fp;fp=fopen("D:\\111\\导线网输出数据3027.txt","w");if(fp!=NULL){fprintf(fp,"距离近似值so(单位：m):\t");fprintf(fp,"方位角近似值alfo(单位：弧度):\n");for(i=0;i<15;i++){fprintf(fp,"%4.12lf\t",so[i]);fprintf(fp,"%.12lf\n",alfo[i]);}fprintf(fp,"Xo(单位：mm)：\t");fprintf(fp,"Yo(单位：mm)：");fprintf(fp,"\n");for(i=0;i<15;i++){fprintf(fp,"%.6lf\t",Xo[i]);fprintf(fp,"%.6lf",Yo[i]);fprintf(fp,"\n");}fprintf(fp,"B矩阵：");fprintf(fp,"\n");for(i=0;i<30;i++){for(j=0;j<27;j++){fprintf(fp,"%.6f ",B[i][j]);}fprintf(fp,"\n");}fprintf(fp,"L矩阵(单位：秒和mm)："); fprintf(fp,"\n");for(i=0;i<30;i++){fprintf(fp,"%.6lf ",L[i][0]);fprintf(fp,"\n");}fprintf(fp,"P矩阵：");fprintf(fp,"\n");for(i=0;i<30;i++){for(j=0;j<30;j++){fprintf(fp,"%.6f ",P[i][j]);}fprintf(fp,"\n");}fprintf(fp,"Nbb矩阵：");fprintf(fp,"\n");for(i=0;i<27;i++){for(j=0;j<27;j++){fprintf(fp,"%.12f ",Nbb[i][j]);}fprintf(fp,"\n");}fclose(fp); //写入完毕，关闭文件}DinV(Nbb,27); //MatrixResult=c * Nbb的逆,此时Nbb已经变成Nbb的逆//计算x^m1=27;m2=27;m3=1;MatrixMutiply(Nbb,W,xgu,m1,m2,m3);double xgu28[50][50],sigma_xy28[50][50];xgu28[0][0]=xgu[0][0];for(i=1;i<27;i++){xgu28[i+1][0]=xgu[i][0];}xgu28[1][0]=0.0;//计算X^(即Xgu估值)Xgu[0][0]=5000.0;Ygu[0][0]=5000.0;for(i=0;i<14;i++){Xgu[i+1][0]=Xo[i+1]+xgu28[2*i][0]/1000.0;Ygu[i+1][0]=Yo[i+1]+xgu28[2*i+1][0]/1000.0; }//精度评定m1=30;m2=27;m3=1;MatrixMutiply(B,xgu,temp5,m1,m2,m3);for(i=0;i<30;i++){if(i<15){V[i][0]=(temp5[i][0]-L[i][0]);}else{V[i][0]=(temp5[i][0]-L[i][0]);}}m1=30;m2=1;MatrixT(V,TV,m1,m2);m1=1;m2=30;m3=30;MatrixMutiply(TV,P,temp4,m1,m2,m3);m1=1;m2=30;m3=1;MatrixMutiply(temp4,V,temp4,m1,m2,m3);sigma_gu=sqrt(temp4[0][0]/3); //单位权中误差double vv=0.0;for (i=0;i<15;i++){for (j=0;j<1;j++){vv=vv+V[i][j];}//puts("");}printf("%lf\t",vv);for(i=0;i<27;i++){for(j=0;j<27;j++){Q[i][j]=Nbb[i][j];}}for(i=0;i<27;i++){sigma_xy[i][0]=sqrt(Q[i][i])*sigma_gu; //坐标平差值中误差//printf("%lf\n",sigma_xy[i][0]);}sigma_xy28[0][0]=sigma_xy[0][0];for(i=1;i<27;i++){sigma_xy28[i+1][0]=sigma_xy[i][0];}sigma_xy28[1][0]=0.0;//printf("%.10lf\n",Ncc[0][0]);FILE *fp1;fp1=fopen("D:\\111\\导线网输出数据3027.txt","a"); fprintf(fp1,"Nbb的逆：");fprintf(fp1,"\n");for(i=0;i<27;i++){for(j=0;j<27;j++){fprintf(fp1,"%lf ",Nbb[i][j]);}fprintf(fp1,"\n");}fprintf(fp1,"W：");fprintf(fp1,"\n");for(i=0;i<27;i++){for(j=0;j<1;j++){fprintf(fp1,"%.12lf ",W[i][j]);}fprintf(fp1,"\n");}fprintf(fp1,"x^(单位：mm)：");fprintf(fp1,"\n");for(i=0;i<28;i++){for(j=0;j<1;j++){fprintf(fp1,"%.10lf ",xgu28[i][j]);}fprintf(fp1,"\n");}fprintf(fp1,"V(单位：秒和mm)：");fprintf(fp1,"\n");for(i=0;i<30;i++){for(j=0;j<1;j++){fprintf(fp1,"%.10lf ",V[i][j]);}fprintf(fp1,"\n");}fprintf(fp1,"X^(单位：m)：\t");fprintf(fp1,"Y^(单位：m)：");fprintf(fp1,"\n");for(i=0;i<15;i++){for(j=0;j<1;j++){fprintf(fp1,"%.6lf\t",Xgu[i][j]);fprintf(fp1,"%.6lf",Ygu[i][j]);}fprintf(fp1,"\n");}fprintf(fp1,"单位权中误差(单位：mm)：\n"); fprintf(fp1,"%.10lf\n",sigma_gu);fprintf(fp1,"X坐标平差值中误差(单位：mm):\t");fprintf(fp1,"Y坐标平差值中误差(单位：mm)：");fprintf(fp1,"\n");for(i=0;i<13;i++){fprintf(fp1,"%lf\t\t",sigma_xy28[2*i][0]);fprintf(fp1,"%lf",sigma_xy28[2*i+1][0]);fprintf(fp1,"\n");}}C语言画图程序//导线网概略图#include <stdio.h>#include <windows.h>#include <math.h>#define NUM 30LRESULT CALLBACK WndProc (HWND, UINT, WPARAM, LPARAM) ;int WINAPI WinMain (HINSTANCE hInstance, HINSTANCE hPrevInstance, PSTR szCmdLine, int iCmdShow){static char szAppName[] = "SineWave" ;HWND hwnd ;MSG msg ;WNDCLASSEX wndclass ;wndclass.cbSize = sizeof (wndclass) ;wndclass.style = CS_HREDRAW | CS_VREDRAW ;wndclass.lpfnWndProc = WndProc ;wndclass.cbClsExtra = 0 ;wndclass.cbWndExtra = 0 ;wndclass.hInstance = hInstance ;wndclass.hIcon = LoadIcon (NULL, IDI_APPLICATION) ;wndclass.hCursor = LoadCursor (NULL, IDC_ARROW) ;wndclass.hbrBackground = (HBRUSH) GetStockObject (WHITE_BRUSH) ; wndclass.lpszMenuName = NULL ;wndclass.lpszClassName = szAppName ;wndclass.hIconSm = LoadIcon (NULL, IDI_APPLICATION) ;RegisterClassEx (&wndclass) ;hwnd = CreateWindow (szAppName, "控制网图",WS_OVERLAPPEDWINDOW,CW_USEDEFAULT, CW_USEDEFAULT,CW_USEDEFAULT, CW_USEDEFAULT,NULL, NULL, hInstance, NULL) ;ShowWindow (hwnd, iCmdShow) ;UpdateWindow (hwnd) ;while (GetMessage (&msg, NULL, 0, 0)){TranslateMessage (&msg) ;DispatchMessage (&msg) ;}return msg.wParam ;}LRESULT CALLBACK WndProc (HWND hwnd, UINT iMsg, WPARAM wParam, LPARAM lParam){static int cxClient, cyClient ;HDC hdc ;int i ;PAINTSTRUCT ps ;POINT pt [NUM] ;switch (iMsg){case WM_SIZE:cxClient = LOWORD (lParam) ;cyClient = HIWORD (lParam) ;return 0 ;case WM_PAINT:hdc = BeginPaint (hwnd, &ps) ;doubley[16]={500.0,363.2211,250.8730,245.2811,343.9816,399.8183,506.0272,596.9804,6 90.1846,686.1801,688.9648,550.4837,552.8629,499.4146,497.8342,500.0}; doublex[16]={500.0,500.0,506.038,269.917,264.999,268.990,346.426,388.125,405.382,466. 024,553.062,555.201,651.049,652.317,591.707,500.0};for (i = 0 ; i < 16 ; i++){pt[i].x = x[i];pt[i].y = y[i];}Polyline (hdc, pt, 16);return 0 ;case WM_DESTROY:PostQuitMessage (0) ;return 0 ;}return DefWindowProc (hwnd, iMsg, wParam, lParam) ;}C#程序：Form1:using System;using System.Collections.Generic;using ponentModel;using System.Data;using System.Drawing;using System.Linq;using System.Text;using System.Threading.Tasks;using System.Windows.Forms;using System.IO;using System.Collections.Generic;namespace导线控制网{public partial class Form1 : Form{public Form1(){InitializeComponent();}public string jisuanjieguo;public const double p = 180*3600/Math.PI;//定义角度转弧度的函数public double jiao_hu(double du,double fen,double miao){return (du*3600.0+fen*60.0+miao)/p;}//定义矩阵运算的类（包括矩阵的加、减、乘、转置和求逆）public class matrix_yusuan{//矩阵相加public static double[,] matrix_jia(double[,] Arry, double[,] Arry1) {int m = Arry.GetLength(0);int n = Arry.GetLength(1);int s = Arry1.GetLength(0);int t = Arry1.GetLength(1);double[,] temp = new double[m, n];double[,] tem = {{0}};if (m == s && n == t){for (int i = 0; i < m; i++){for (int j = 0; j < n; j++){temp[i, j] = Arry[i, j] + Arry1[i, j];}}return temp;}else{Console.WriteLine("两个矩阵大小不同");return tem ;}}//矩阵相减public static double[,] matrix_jian(double[,] Arry, double[,] Arry1) {int m = Arry.GetLength(0);int n = Arry.GetLength(1);int s = Arry1.GetLength(0);int t = Arry1.GetLength(1);double[,] temp = new double[m, n];double[,] tem = { { 0 } };if (m == s && n == t){for (int i = 0; i < m; i++){for (int j = 0; j < n; j++){temp[i, j] = Arry[i, j] - Arry1[i, j];}}return temp;}else{Console.WriteLine("两个矩阵大小不同");return tem;}}//矩阵转置public static double[,] matrix_t(double[,] Arry){int m = Arry.GetLength(0);int n = Arry.GetLength(1);double[,] temp = new double[n, m];for (int i = 0; i < n; i++){for (int j = 0; j < m; j++){temp[i,j] = Arry[j,i];}}return temp;}//矩阵相乘public static double[,] matrix_cheng(double[,] Arry, double[,] Arry1) {int m = Arry.GetLength(0);//矩阵Arry的行数int r = Arry.GetLength(1);//矩阵Arry的列数int k = Arry1.GetLength(0);//矩阵Arry的行数int n = Arry1.GetLength(1);//矩阵Arry1的列数double[,] temp = new double[m, n];double[,] tem = { { 0 } };if (r == k){for (int i = 0; i < m; i++){for (int j = 0; j < n; j++){for (int t = 0; t < r; t++){temp[i, j] += Arry[i, t] * Arry1[t, j];}}}return temp;}else{Console.WriteLine("两个矩阵无法相乘");return tem;}}//求矩阵Arry的逆矩阵public static double[,] matrix_ni(double[,] Arryni){int Level = Arryni.GetLength(1);double[,] NArry = new double[Level, Level];double HLS = matrix_yusuan.matrix_hls(Arryni);double[,] BArry = matrix_bansui(Arryni);for (int i = 0; i < Level; i++){for (int j = 0; j < Level; j++){NArry[i, j] = BArry[i, j] / HLS;}}return NArry;}//求矩阵Arry的伴随矩阵public static double[,] matrix_bansui(double[,] Arryni){int Level = Arryni.GetLength(1);double[,] BArry = new double[Level,Level];for (int m = 0; m < Level; m++){for (int n = 0; n < Level; n++){BArry[m, n] = matrix_yusuan.matrix_yuzi(Arryni, n, m);}}return BArry;}//求矩阵Arry的元素Arry[i,j]的余子式public static double matrix_yuzi(double[,] Arryni, int i, int j)//第i行，第j列，起始值为0{int Level = Arryni.GetLength(1);double[,] Arry1 = new double[Level-1,Level-1];for (int m = 0; m < Level - 1; m++){for (int n = 0; n < Level - 1; n++){if (m < i && n < j){Arry1[m, n] = Arryni[m, n];}else if(m < i && n >= j){Arry1[m, n] = Arryni[m, n + 1];}else if (m >= i && n < j){Arry1[m, n] = Arryni[m + 1, n];}else if (m >= i && n >= j){Arry1[m, n] = Arryni[m + 1, n + 1];}}}//根据矩阵元素的不同位置返回不同的值if ((i + j) % 2 != 0){return (-1)*matrix_yusuan.matrix_hls( Arry1);}else{return matrix_yusuan.matrix_hls(Arry1);}}//求行列式public static double matrix_hls(double[,] Arryni){int Level = Arryni.GetLength(1);//简单来说，对于常用的二维数组，取0或者1，分别获取列和行的长度(行数和列数)double[,] dArry = new double[Level, Level];for (int i = 0; i < Level; i++){for (int j = 0; j < Level; j++){dArry[i, j] = Arryni[i, j];}}int sign = 1;for (int i = 0, j = 0; i < Level && j < Level; i++, j++){//判断该行的正对角元素dArry[i,j]是否为0，若为0，则寻找i行以下的行m（m>i,且dArry[m,j]!=0）进行交换if (dArry[i, j] == 0){//判断是否达到了矩阵的最大行if (i == Level - 1){return 0;}int m = i + 1;//获取一个dArry[m,j]不为零的行for (; dArry[m, j] == 0; m++){if (m == Level - 1){return 0;}}//把i行和m行互换double temp;for (int n = j; n < Level; n++){temp = dArry[i, n];dArry[i, n] = dArry[m, n];dArry[m, n] = temp;}sign *= (-1);}//把当前行以下的行所对应的列变成0double tmp;//for (int s = Level - 1; s > i; s--)for (int s = i + 1; s < Level; s++){tmp = dArry[s, j];for (int t = j; t < Level; t++){dArry[s, t] -= dArry[i, t] * (tmp / dArry[i, j]); }}}double result = 1;for (int i = 0; i < Level; i++){if (dArry[i, i] != 0){result *= dArry[i, i];}else{return 0;}}return sign * result;}}public class aa//用于判断文件导入和计算是否完成{public static int panduan1 = 0;//panduan1用来判断是否进行了计算public static int panduan2 = 0;//panduan2用来判断是否导入文件}int i;int j;//定义数组double[,] hudu = new double[30, 1];double[,] bian = new double[30, 1];double[] alf = new double[30];double[] alfo = new double[30];double[] so = new double[30];double[] Xo = new double[30];double[] Yo = new double[30];double[,] B = new double[30, 27];double[,] L = new double[30, 1];double[,] P = new double[30, 30];double[,] W = new double[27, 1];double[,] temp = new double[30, 30];//此处为Xo Yo，B矩阵赋初值为零 hudu代表夹角，alf方位角，alfo方位角近似值，s 距离观测值，so距离近似值//doubleNbb[50][50],xgu[50][50],Xgu[30][1],Ygu[30][1],V[50][50],sigma_gu,Q[50][50],sigma_xy[50] [50];double[,] Nbb = new double[27, 27];double[,] xgu = new double[27, 1];double[,] xgu28 = new double[28, 1];double[,] Xgu = new double[30, 1];double[,] Ygu = new double[30, 1];double[,] V = new double[30, 1];double[,] Q = new double[27, 27];double[,] sigma_xy = new double[30, 30];double[,] sigma_xy28 = new double[30, 30];double sigma_gu;double D, F, M, sigma_hudu = 2.0, sigma_s;//用来求未知点坐标近似值private void toolStripMenuItem7_Click(object sender, EventArgs e) {if (aa.panduan1 == 1){textBox2.Text = "";textBox2.Text = "近似方位角alfo" + "\t\t" + "\r\n";for (i = 0; i < 15; i++){textBox2.Text += Math.Round(alfo[i], 6);textBox2.Text += "\r\n\r\n";}}else{MessageBox.Show("请先进行计算！", "系统提示");}}private void toolStripMenuItem10_Click(object sender, EventArgs e) {if (aa.panduan1 == 1){textBox2.Text = "";textBox2.Text = "P矩阵" + "\t\t" + "\r\n";for (i = 0; i < P.GetLength(0); i++){for (j = 0; j < P.GetLength(1); j++){textBox2.Text += Math.Round(P[i, j], 6) + "\t";}textBox2.Text += "\r\n\r\n";}}else{MessageBox.Show("请先进行计算！", "系统提示");}}private void toolStripMenuItem9_Click(object sender, EventArgs e) {if (aa.panduan1 == 1){textBox2.Text = "B矩阵" + "\t\t" + "\r\n";for (i = 0; i < B.GetLength(0); i++){for (j = 0; j < B.GetLength(1); j++){textBox2.Text += Math.Round(B[i, j], 6) + "\t";}textBox2.Text += "\r\n\r\n";}}else{MessageBox.Show("请先进行计算！", "系统提示");}}private void button15_Click(object sender, EventArgs e){if (aa.panduan2 == 1){alf[0] = Math.PI;alfo[0] = Math.PI;Xo[0] = 5000.0;Yo[0] = 5000.0;Xo[15] = 5000.0;Yo[15] = 5000.0;so[0] = bian[0, 0];for (i = 1; i < 15; i++){alf[i] = alf[i - 1] + hudu[i, 0] - Math.PI;if (alf[i] >= (2 * Math.PI)){alf[i] = alf[i] - 2 * Math.PI;}Xo[i] = Xo[i - 1] + bian[i - 1, 0] * Math.Cos(alf[i - 1]);Yo[i] = Yo[i - 1] + bian[i - 1, 0] * Math.Sin(alf[i - 1]);}jisuanjieguo += "近似距离(单位：m)：" + "\r" + "近似方位角(单位：弧度)：" + "\r\n";jisuanjieguo += Math.Round(so[0], 3) + "\t\t" + Math.Round(alfo[0], 6) + "\r\n";for (i = 1; i < 15; i++){//求近似距离so[i] = Math.Sqrt((Yo[i + 1] - Yo[i]) * (Yo[i + 1] - Yo[i]) + (Xo[i + 1] - Xo[i]) * (Xo[i + 1] - Xo[i]));//求近似方位角，分象限进行讨论if ((Yo[i + 1] - Yo[i]) > 0 && (Xo[i + 1] - Xo[i]) > 0)//第一象限{alfo[i] = Math.Atan((Yo[i + 1] - Yo[i]) / (Xo[i + 1] - Xo[i])); }else if ((Yo[i + 1] - Yo[i]) > 0 && (Xo[i + 1] - Xo[i]) < 0)//第二象限{alfo[i] = Math.Atan((Yo[i + 1] - Yo[i]) / (Xo[i + 1] - Xo[i])) + Math.PI;}else if ((Yo[i + 1] - Yo[i]) < 0 && (Xo[i + 1] - Xo[i]) > 0)//第三象限{alfo[i] = Math.Atan((Yo[i + 1] - Yo[i]) / (Xo[i + 1] - Xo[i])) + 2 * Math.PI;}else//((Yo[i]-Yo[i-1])<0&&(Xo[i]-Xo[i-1])<0)//第四象限{alfo[i] = Math.Atan((Yo[i + 1] - Yo[i]) / (Xo[i + 1] - Xo[i])) + Math.PI;}jisuanjieguo += Math.Round(so[i], 3) + "\t\t"+ Math.Round(alfo[i], 6) + "\r\n";}//将第一个角度的系数单独算出j = 0;//B[0][j]=0.0;((Yo[14]-Yo[0])/(so[14]*so[14]))*p/1000.0;//B[0][j+1]=0.0;-1*((Xo[14]-Xo[0])/(so[14]*so[14]))*p/1000.0;B[0, j] = -1 * (Yo[1] - Yo[0]) / (so[0] * so[0]) * p / 1000.0;//B[0][j+1]=(Xo[1]-Xo[0])/(so[0]*so[0])*p/1000.0;B[0, 25] = ((Yo[14] - Yo[0]) / (so[14] * so[14])) * p / 1000.0;B[0, 26] = -1 * (Xo[14] - Xo[0]) / (so[14] * so[14]) * p / 1000.0;//将第二个角度的系数单独算出B[1, j] = ((Yo[2] - Yo[1]) / (so[1] * so[1]) - (Yo[0] - Yo[1]) / (so[0] * so[0])) * p / 1000.0;//B[1][j+1]=-1*((Xo[2]-Xo[1])/(so[1]*so[1])-(Xo[0]-Xo[1])/(so[0]*so[0]))*p/1000.0;B[1, j + 1] = -1 * ((Yo[2] - Yo[1]) / (so[1] * so[1])) * p / 1000.0;B[1, j + 2] = ((Xo[2] - Xo[1]) / (so[1] * so[1])) * p / 1000.0;//i = i + 1;//将第三个角度的系数单独算出i = 2;B[i, j] = ((Yo[i - 1] - Yo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;//B[i, j + 1] = -1 * ((Xo[i - 1] - Xo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;B[i, j + 1] = ((Yo[i + 1] - Yo[i]) / (so[i] * so[i]) - (Yo[i - 1] - Yo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;B[i, j + 2] = -1 * ((Xo[i + 1] - Xo[i]) / (so[i] * so[i]) - (Xo[i - 1] - Xo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;B[i, j + 3] = -1 * ((Yo[i + 1] - Yo[i]) / (so[i] * so[i])) * p / 1000.0; B[i, j + 4] = (Xo[i + 1] - Xo[i]) / (so[i] * so[i]) * p / 1000.0;//求其他角度改正的系数for (i = 3; i < 15; i++){if (i < 14){B[i, j + 1] = ((Yo[i - 1] - Yo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;B[i, j + 2] = -1 * ((Xo[i - 1] - Xo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;B[i, j + 3] = ((Yo[i + 1] - Yo[i]) / (so[i] * so[i]) - (Yo[i - 1] - Yo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;B[i, j + 4] = -1 * ((Xo[i + 1] - Xo[i]) / (so[i] * so[i]) - (Xo[i - 1] - Xo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;B[i, j + 5] = -1 * ((Yo[i + 1] - Yo[i]) / (so[i] * so[i])) * p / 1000.0;B[i, j + 6] = (Xo[i + 1] - Xo[i]) / (so[i] * so[i]) * p / 1000.0; }else{B[i, j + 1] = (Yo[i - 1] - Yo[i]) / (so[i - 1] * so[i - 1]) * p / 1000.0;B[i, j + 2] = -1 * ((Xo[i - 1] - Xo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;B[i, j + 3] = ((Yo[i + 1] - Yo[i]) / (so[i] * so[i]) - (Yo[i - 1] - Yo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;B[i, j + 4] = -1 * ((Xo[i + 1] - Xo[i]) / (so[i] * so[i]) - (Xo[i - 1] - Xo[i]) / (so[i - 1] * so[i - 1])) * p / 1000.0;}j = j + 2;}//求边长改正的系数j = 0;//将第一个边长的系数单独算出B[i, 0] = (Xo[1] - Xo[0]) / so[0]; //(Xo[1]-Xo[0])/(so[0]*so[0]);i = i + 1;//将第二个边长的系数单独算出B[i, j] = -1 * (Xo[i - 14] - Xo[i - 15]) / so[i - 15];//B[i, j + 1] = -1 * (Yo[i - 14] - Yo[i - 15]) / so[i - 15];B[i, j + 1] = -1 * B[i, j];B[i, j + 2] = (Yo[i - 14] - Yo[i - 15]) / so[i - 15];for (i = i + 1; i < 30; i++){if (i < 29){B[i, j + 1] = -1 * (Xo[i - 14] - Xo[i - 15]) / so[i - 15];B[i, j + 2] = -1 * (Yo[i - 14] - Yo[i - 15]) / so[i - 15];B[i, j + 3] = -1 * B[i, j + 1];B[i, j + 4] = -1 * B[i, j + 2];}else{B[i, j + 1] = -1 * (Xo[i - 14] - Xo[i - 15]) / so[i - 15];B[i, j + 2] = -1 * (Yo[i - 14] - Yo[i - 15]) / so[i - 15];}j = j + 2;}jisuanjieguo += "B矩阵：" + "\r\n";for (i = 0; i < 30; i++){for (j = 0; j < 27; j++){jisuanjieguo += Math.Round(B[i, j], 6) + "\t";}jisuanjieguo += "\r\n";}//求L矩阵，角度和边长分别求解L[0, 0] = (hudu[0, 0] - (alfo[0] - alfo[14] + Math.PI)) * p;for (i = 1; i < 30; i++){if (i < 15){L[i, 0] = (hudu[i, 0] - (alfo[i] - alfo[i - 1] + Math.PI)) * p;if (L[i, 0] > Math.PI * p){L[i, 0] = (L[i, 0] - 2 * Math.PI * p);}else{L[i, 0] = L[i, 0];}}else{L[i, 0] = (bian[i - 15, 0] - so[i - 15]) * 1000.0;}}jisuanjieguo += "L矩阵(单位：秒和mm)：" + "\r\n";for (i = 0; i < 30; i++){for (j = 0; j < 1; j++){jisuanjieguo += Math.Round(L[i, j], 6) + "\t\t";}jisuanjieguo += "\r\n";}//求权阵Pfor (i = 0; i < 30; i++){if (i < 15){P[i, i] = 1;}else{sigma_s = 5 + 10 * 0.000001 * bian[i - 15, 0] * 1000; //单位为（''/mm）的平方P[i, i] = sigma_hudu * sigma_hudu / (sigma_s * sigma_s);}}jisuanjieguo += "P矩阵：" + "\r\n";for (i = 0; i < 30; i++){for (j = 0; j < 30; j++){jisuanjieguo += Math.Round(P[i, j], 6) + "\t";}jisuanjieguo += "\r\n";}Nbb =matrix_yusuan.matrix_cheng(matrix_yusuan.matrix_cheng(matrix_yusuan.matrix_t(B), P), B); //计算Nbb矩阵double[,] temp2 = new double[27, 30];。

基于MATLAB的单三角锁严密平差程序设计

一
ｉｓ＝；０＝ｌ％当待定点为第奇数个时，待ｌ＝ｅｏ（；ＯｚｒｓＯ）ｒｏ＝：ｆｒｉｌＨ定点近似坐标的求法；３ｘ（（１ｃｔ１２＂ｘ一）ｃｔ１２＂ｙｆｒｊ１；ｘ￣ｉ）ｏ（－，＋Ｏ２ｏ（一， — ｏ＝：ｉ一ｄ（１ｉｄ（３－ｉａ１ｙ— ）一）＇＋；ｉｌｊ０ｆＯ）＝＞ｆ）ｃｔ（２）ｃｔａ２肛￣＝ｏ０－，— ｏ０一 ’ ｉｄｉ１－ｄ３１，１（）０ｊ０ｉ；６）０，＝ｊｅｓｌｅｙ＝０１ｃｔ１２＂ｙ —）ｏ０－３））一）ｏ（一，＋６２ｃｔ６２）ｘ－ｄ（ｌｙｉｄ＋
⑤计算法方程式的系数及权的编程如下：
ｌ＝２６６．０６４０６；３０；ｏ０２４８２７９４Ａ＝６０，
“ ２＋１＿）一ｆｉｃｔ（－，）ｃｔＯｉ２）ｆ￣ｏｄ１２）ｏｄ（，）（＝（３＋ｉ－１：ｙ（ｙ一）ｏ（一’ ａ１ｃｔｏ２ ” ｙ＝（２ｃｔ１２ｙ— ）ｏ０一，＋ｉｉｄ（３）ｉｄ１（２－（１ｉ）－ｉ）－ｘ一；列出误差方程＝ —ｆ；ｘ－ｘ）（（ｘ（ｆ）ｉｉｆ）／ｉ（）由误差方程系数Ｂ和自由项ｌ３组成法ｙ）ｙ（ｆｉ（ｙ）ｆｉｉ（－／）ｅｄｎ方程Ｎ — ：Ｗ０，法方程个数等于参数的个数ｔ；ｅｄｎ
Ｂ＝ｅｏ０ｎｎ２；Ｙｚｒ＊，（）ｓ２＋）始化矩阵Ｂ；Ｙ

Matlab源程序代码

正弦波的源程序：(一)，用到的函数1,f2t函数function x=f2t(X)global dt df t f T N%x=f2t(X)%x为时域的取样值矢量%X为x的傅氏变换%X与x长度相同并为2的整幂%本函数需要一个全局变量dt(时域取样间隔) X=[X(N/2+1:N),X(1:N/2)];x=ifft(X)/dt;end2,t2f函数。

function X=t2f(x)global dt df N t f T%X=t2f(x)%x为时域的取样值矢量%X为x的傅氏变换%X与x长度相同,并为2的整幂。

%本函数需要一个全局变量dt(时域取样间隔) H=fft(x);X=[H(N/2+1:N),H(1:N/2)]*dt;end(二),主程序。

1，%（1）绘出正弦信号波形及频谱global dt df t f Nclose allk=input('取样点数＝2^k, k取10摆布');if isempty(k), k=10; endf0=input('f0=取1(kz)摆布');if isempty(f0), f0=1; endN=2^k;dt=0.01; %msdf=1/(N*dt); %KHzT=N*dt; %截短期Bs=N*df/2; %系统带宽f=[-Bs+df/2:df:Bs]; %频域横坐标t=[-T/2+dt/2:dt:T/2]; %时域横坐标s=sin(2*pi*f0*t); %输入的正弦信号S=t2f(s); %S是s的傅氏变换a=f2t(S); %a是S的傅氏反变换a=real(a);as=abs(S);subplot(2,1,1) %输出的频谱plot(f,as,'b');gridaxis([-2*f0,+2*f0,min(as),max(as)])xlabel('f (KHz)')ylabel('|S(f)| (V/KHz)') %figure(2)subplot(2,1,2)plot(t,a,'black') %输出信号波形画图gridaxis([-2/f0,+2/f0,-1.5,1.5])xlabel('t(ms)')ylabel('a(t)(V)')gtext('频谱图')最佳基带系统的源程序：(一)，用到的函数f2t函数和t2f函数。

学位论文—基于matlab测量平差程序设计-创新实践报告

创新实践报告实践名称：基于MATLAB测量平差程序设计系部名称：测绘工程学院专业班级：测绘工程11-6班学生姓名：学号：指导教师：xxx工程学院教务处制注：1、此报告为参考格式，各栏项目可根据实际情况进行调整；实验成绩以优（90~100）、良（80~89）、中（70~79）、及格（60~69）、不及格（60以下）五个等级评定。

附录1 外文参考文献原文参考文献一Improving oil fields recovery through real-time water flooding optimizationPamela Alessia Chiara MarescalcoPolitecnico di Torino, Land, Environment and GeoTechnologies Engineering Department,24, C.so Duca degli Abruzzi, 10129 Torino, ItalySUMMARYIncreasing oil recovery from reservoirs is a strong urge. One of the most effective ways to get the result iswater flooding and that’s why its application is nowadays widely used in the petroleum industry. Obviously,water flooding efficiency strongly depends on reservoir properties; this makes simulating a water injectionprocess a priori an extremely important step of the reservoir production strategy. Simulation is commonlydone adopting a finite difference (FD) simulation approach.This paper explores a different and complementary approach, represented by streamline-based simulation,coupled with a tool to optimize water flooding campaigns and to help quick decision making. Inthe present study, water flooding simulation is performed via two commercial software: an FD and astreamline-based simulator, to highlight advantages and disadvantages of both simulation techniques indescribing a water injection campaign and to exploit the two appr oaches’ uniqueness in parallel.The final goal of iteratively converging to the optimal water flooding scheme, which is the core of thepresent work, is achieved through a customized Matlab script. The generated automatic procedure showsits effectiveness in improving oil recovery, expediting decision making and saving time and FD simulationruns. A three steps workflow is outlined to get the best water flooding scheme for the examples shownbelow. Copyright q 2008 John Wiley & Sons, Ltd.Received 13 March 2008; Revised 1 September 2008; Accepted 8 November 2008KEY WORDS: water displacement optimization; FD simulation; streamline simulation; Matlab automatedroutine; real-time decision making; quantitative and iterative adjustment of water rates tobe injectedINTRODUCTIONIn the last decades water flooding has been widely applied in the petroleum field, both in mature and in newly developed fields, and its attractiveness lies in supporting the entire field pressureduring depletion and in improving the final oil recovery [1–3]. The techniqueconsists in injecting water with the purpose of displacing and therefore producing oil, especially if the reservoir lacks an underlying aquifer able to counterbalance the depletion and to drive oil to the producer wells.To really understand a water flooding process and to predict its efficiency, it is useful to simulate it a priori. This is usually done via finite difference (FD) simulation, which is able to describe any kind of existing reservoir very accurately, but unfortunately is not able to give enough details about the way the flow occurs throughout the field. Recent works [4, 5] have proposed a newly developed approach, based on streamline simulation, whose main attraction lies in providing information not obtainable from FD simulation and useful for the purpose of improving reservoir performances.In the present study, water flooding simulation is performed via two commercial software: an FDand a streamline-based simulator, to highlight advantages and disadvantages of both simulationtechniques in describing a water injection campaign and to exploit the two approaches’ uniquenessin parallel. If FD simulation is essential to checking the streamline simulation results, the streamlinebasedsimulator is, on the other hand, an ideal tool to perform a procedure able to optimizewater injection. Then, the core of the work and its innovative approach lies in exploiting the twosoftware features in conjunction with the application of a customized Matlab code developed inorder to elaborate all the streamline simulation outputs, to calculate the changes to be made to theproduction/injection constraints for the subsequent simulation run, so as to iteratively converge tothe optimal injection scheme for the reservoir under study. A three steps workflow is outlined toget the best water flooding scheme for the examples shown below.FD APPROACH VERSUS STREAMLINE SIMULATIONECLIPSE is a complete and complex simulator, whose attractiveness resides in being able todescribe any kind of reservoir, including geological complexity, formation features, and fluidproperties. The simulator is based on a time and spatial discretization and solves a three-dimensionalequation by assigning to all the parameters involved in the simulation a unique value, associated withthe entire cell, for every grid block. For models with a large number of cells, using FD relaxationcan be computationally heavy; therefore, a good balance must be kept between having sufficientaccuracy in describing the reservoir and keeping simulation time within reasonable limits [6].3DSL, on the other side, solves two different equations on two different grids and this is usuallyfaster than FD simulation, especially for large models: the pressure equation is solved implicitlyon the background grid (or pressure grid), whereas the saturation equation is solved explicitly onthe streamline grid. This involves a minor effect of grid refinement on the resultsof the simulation,time-step limitations not as severe, thanks to a better stability of the geometrical grid, numericaldiffusion easier to control, faster solutions with respect to FD approach [7].Streamline simulation solves mono-dimensional equations along streamlines, which means itsolves multiple streamlines in parallel, and the fluid transport, which for FD approach occursbetween grid cells [8], occurs along streamlines: this gives an immediate answer in terms of howthe streamlines (connecting injector and producer wells—to say, well pairs) are distributed, sothat the fluid trajectories and their rates at the wells are known at every time step [9, 10]. Thanksto the available information, the distribution of injected water volumes can be modified, and amore effective production strategy can be planned to maximize oil recovery [4, 5]. Same as for FDsimulation, when using streamline-based simulation a good balance must be maintained betweenpressure updates and computational speed.In conclusion, 3DSL is simpler from a reservoir description (includes both flow physics andpetro-physical/geological information) point of view, but it cannot take into account importantparameters such as capillary pressure or cross flow effects: moreover, fluid compressibility cannot beeasily taken into account and mass conservation errors may occur while mapping between pressuregrid and streamline grid therefore being incomplete or inadequate for most of real reservoirs. Everyfurther detail can be found in previous works [11]. THE IDEA OF THE STUDYIn this paper the use of the two software as complementary tools to optimize water injectionis shown: ECLIPSE appeared to be the most reliable software to simulate reservoir models andto be essential to verify 3DSL’s results accuracy upstream and downstream of the methodologydeveloped for this research, while 3DSL has been shown to be simpler and usually faster, providinginformation not obtainable by means of FD simulation, and which could be exploited to save atrial-and-error process trying to find the best injection schemeThe idea of the study resulted from previous works [12].The research presented here follows a three step workflow. In the first part of the study, thetwo software were used to simulate simple synthetic models. Once the results obtained from thetwo simulators were checked for consistency, a method was developed to exploit 3DSL’s features.The most interesting pieces of information available from 3DSL are the connections between wellpairs (injector-producer) and the rates at the wells and the connections. With this informatioavailable, the streamline simulator, coupled with Matlab, was involved into an automatic procedurethat reallocated water injected volumes in order to optimize water injection campaigns, for all the examined cases.This was gained by a customized Matlab script, writtenfor this specific purpose, which automaticallyinteracted with the streamline simulator, processed the input data supplied by the softwaregiving back new input rates to be used for the following simulation run. This was done for a fixednumber of runs to iteratively converge to the optimal injection scheme. From preliminary analyses[13] the procedure was shown to be effectiveEventually, in the third experimental phase, the last rates calculated from Matlab were inputinto ECLIPSE to check the procedure’s effectiveness and the added val ue in terms of oil recoveryimprovement.THE AUTOMATED PROCEDURE AND THE EXPERIMENTAL METHOD Let’s now focus on the second and main part of the experimental study which, as we said, consistedin writing a customized Matlab code. Streamline answers in terms of streamlines distribution(connecting injector and producer wells—or well pairs) and fluid rates at the wells are writtenat every time step in a 3DSL output file named ∗.WAF, which is crucial for the entire codingdescribed below and for its application.The endeavor of the Matlab customized code is the optimization of the displacement processthrough a gradual reallocation of injected water rates, carried out by increasing the water volumeinjected at the highly efficient connections and decreasing it at the poorly efficient connections.The term injection efficiency (for a well pair connection) stands for the ratio between offset oilproduced thanks to water displacement and the amount of water injected at a certain well. In ananalogue manner the injection efficiency of an injector well is the ratio between offset oil at allproducers connected to it and the total water injected at the same well. The approach used forthe current application was aimed at increasing the oil production through a better use of a givenwater volume available for injection. The Matlab code written for this purpose mainly works asfollows: the data needed for the procedure are read from 3DSL’s output ∗.WAF file and loadedinto Matlab environment, then they are processed throughout the code; eventually new rates forthe next simulation run are output to 3DSL. This is all done automatically with an approach aimedat maximizing oil production through a more efficient use of a given water volume available forthe injection [13]. The code steps are here summarized [12]:1st step consists in: establishing the number of iterations to be run for the procedure of ratereallocation and their duration; fixing the volume of water to be injected and the target liquid rate.2nd step consists in: running a ‘do nothing’ 3DSL simulation in order to choose a properstarting time for the entire procedure of rate reallocation (usually the reallocation startswhereas aproduction plateau starts).3rd step consists in: reading (within the Matlab environment, by means of functions coded onpurpose) from the ∗.WAF file: rates (total rate for each injector well, and partial rate for eachproduction well connected to the injector); time; name and number of injector or producer wellsand connections between them; number of existing connections in the model at any time step ofthe simulation.4th step consists in: calculating the injection efficiency for each well pair, for each injector well,and for the field (average injection efficiency).5th step consists in: calculating for each iteration and for each injector well a new rate:Qnew i=(1+wi )·qoldi (1)where i stands for the well, wi for the weight it has been assigned, and qoldi for the rate injectedat the well at the previous step.The average reference field injection efficiency (referred to as signed e) is a mean value:depending on its value, positive or negative weights are assigned to the efficient or inefficient connections.Once the maximum weight, minimum weight, minimum injection efficiency allowed, maximum injection efficiency awaited, and _ (which is the grade of the polynomial that interpolates the relation weight-efficiency) are fixed, the weights are evaluated as follows:where wmin is the minimum weight at the least efficiency, wmax is the maximum weight at thehighest efficiency, emin is the least acceptable injection efficiency, emax is the highest awaitedinjection efficiency.6th step consists in: calculating new rates for each well pair connection and, consequently, foreach injector; calculating for each injector the differences between the last known rates and thenew ones; determining the new total volume of water injected as a summation of the new rates atthe injectors; checking that the constraints on the total water volume available (x) are satisfied byFigure 1. Weighting functions (as from Equations (2) and (3)) for different values of exponent_.introducing a correction factor c:c=_x qnew i(4)7th step consists in: calculating new liquid rates (for each producing well the novel rate iscalculated by adding/subtracting the same amounts of water rate, _q±, calculated for the injectorwells connected to it), imposing that the constraints on the total liquid volume (y) are satisfiedthrough the correction factor c1:c1= _y qnew i(5)8th step consists in: making Matlab write a text file to be included in the 3DSL dataset,summarizing all the new rates and time information.9th step consists in: running the reallocation procedure, following the eight steps listed above,for as many iterations as initially fixed (well rates and time are updated at every _t).Both _t and the set of parameters chosen at the 5th step affect the final results.For all the examples shown here, to signed e was allotted the average field efficiency value, inorder to have a case sensitive value. At the same time, _ was put equal to 2 for all the exam inedcases, so that the weighting function would be nonlinear and the most significant changes in injection rates would occur far-off from the average efficiency, while only small changes would take place around the mean value (Figure 1).For the examples shown below, parameter sensitivities were performed in order to get the bestset for the application of the Matlab code.RESERVOIR MODELS DESCRIPTIONIn this paper two synthetic models and a real one are presented.The numbers of cells for the synthetic models are 20×20×1 for the first case, and 123×54×1for the second case. A uniform discretization of the volume was chosen to accurately describe the phases present in both models are water and oil, which is supposed to be dead oil. Both water and oil have the same viscosity, equal to 1 cp. Oil and water density are, respectively, 780 and1029kg/m3. Datum pressure is equal to 30 bars. For both models, water and liquid constraints are set equal, to let the displacement occur by voidagereplacement, and correspond to 80.000rm3/day and to 150.000rm3/day, respectively. For both models the reallocation procedure starts at day 730;the number of steps, whose duration is 730 days each, is equal to 10.Similarly a real field was tested. The model was available from ECLIPSE; a translation work was done to obtain a model written for 3DSL. The two models were then checked to verify the coherence of the main parameters. Once the match was gained, 3DSL model was used to apply the procedure. The real model, as shown in Figure 4(a, b), is divided into three regions by the presence of faultsFigure 4. (a) Faults divide model three into three regions and (b) model three has three different regions.Figure 5. Permeability (a) and porosity map (b). Model three.The number of cells is 68×71×23. The grid cells are not uniform. Permeability and porosity arehighly heterogeneous and range from 0 to 1315mD and from 4 to 27%, respectively (Figure 5(a, b)).The other model properties are summarized in Tables I and II.Table II. Model three: main parameters.Swi (three regions) 0.13; 0.29; 0.06RESULTSThe three models examined were all used for the three steps workflow related before.First of all, the match and coherence between the two models were always obtained. At thispoint no relevant discrepancies were observed between the two sets of results for thethree models.Figure 6(a, b) shows, respectively, oil and water cumulative production from the third model only,as an example.At a second stage, Matlab code was applied to each 3DSL model, starting from a certain moment in time, fixed at day 730 for all the three examples. Ten iterations were run, each of them lasting730 days, and the total time for the procedure to run was fixed so as to be compared with the base case ‘do nothing’ simulation run.As for the frequency of the injection update, the longer the time step is, the higher the rates injected/produced and the major the changes concerning the values of well pairs injection efficiencies during the time steps. For the present work the well pairs injection efficiencies were assumed to remain constant throughout the single time step. Shorter time steps would maybe more suitable to this hypothesis but, on the other hand, they would imply a much higher computational time. In a real field application it might be useful to try with medium time steps (i.e. 6 months),eventually shortening them if, by analyzing the real-time field data, an intervention would appear to be needed.Each model was then compared with its respective optimized case, as shown in Figures 7–9(a, b).The figures again show, respectively, oil and water cumulative production from the field for the three models in sequence. For all the models the generated methodology was shown to give good results and to carry advantages with respect to the mere base case.The two synthetic cases showed good results both in terms of a relevant improvement in oil recovery and of a significant reduction in water production (the results are summarized in Tables III–V). Figure 10(a, b) shows the displacement of model two at the end of the simulation,for the base case and for the reallocated rates case, respectively.The real case, which at the time of this research was a pure forecast example, since the reservoir was not producing, so no production history was available, showed to gain a certain amount of oil production, but in parallel a higher water production.At the third stage the final rates observed from the 3DSL ‘optimized’ case were input in ECLIPSE, to double check the effectiveness of the optimization for the specific case. Figure 11(a, b)P. A. C. MARESCALCOFigure 6. Match ECLIPSE versus 3DSL: oil (a) and water cumulative production (b). Model three.Figure 7. Base case versus reallocated rates case: oil (a) and watercumulative production (b). Model one.Figure 8. Base case versus reallocated rates case: oil (a) and watercumulative production (b). Model two.RECOVERY THROUGH REAL-TIME WATER FLOODING OPTIMIZATIONFigure 9. Base case versus reallocated rates case: oil (a) and water cumulative production (b). Model three shows oil and water cumulative production for the second synthetic model, obtained with ECLIPSE before and after the reallocation procedure.For the synthetic models, whose ECLIPSE results downstream of Matlab procedure are as good as expected, the methodology proved to be valid.Figure 10. Oil displaced at the end of the simulation run: base case(a) and reallocated rates case (b). Model two.Figure 11. Reallocated rates case versus base case, ECLIPSE check simulation run, oil(a) and water cumulative production (b). Model two As for the real case, the results at this stage are not as good as expected. In the next paragrapha few comments on this will be remarked.DISCUSSION AND CONCLUSIONSLooking at the results, it can be generally said that the Matlab code written for this study has added value to the reallocation procedure, making it reliable and attractive for itsapplication in the petroleum industry. Previous works [12] conducted for fields under production and the two simpler examples reported here undoubtedly proved that reallocating water to optimize waterinjection is effective both for models with simpler reservoir features (incompressible fluids, simple geology) and for more complex cases (compressible fluids, heterogeneous formation): an in creased oil production was observed together with a reduced water production.As stated in the previous paragraph, for the real model, the optimized case with respect tothe base case shows a good improvement in oil recovery, but a parallel higher water production:the field in fact had not yet been deployed, so no real-time production neither geological no rpetrophysical information was provided. The reservoir features and the specific values assigned to the main reservoir parameters were then simply estimated: the values assigned to model parameters were a choice among other equally valid assessments. Last but not the least, ECLIPSE’s results certainly depend on the large amount of modifications done with respect to the original ECLIPSE model, modification necessary in order to comply with 3DSL’s binding simplifying hypotheses at the time the model ‘translation’ was doneThe aim of the research work of developing a completely automatic procedure to find the best scheme for the volumes to be injected and of using 3DSL coupled with Matlab and ECLIPSE as complementary tools proved to be effective and valid. The entire work represents a real help when a quick decision on a water flooding arrangement has to be made, both for injection campaigns yet to plan or for injection operations to be adjusted.If a reallocation procedure was to be performed only by means of ECLIPSE, a trial-and-error procedure would be needed, with a considerable waste of computational time and simulation runs.On the other hand, 3DSL allows for more simulations in parallel, saving computational time and redundant simulations.The present research in the future can be extended to a real case with a known production history to calibrate the best injection scheme.附录2 外文参考文献翻译提高油田采收率通过实时水驱优化帕米拉·阿卢嘉勒都灵理工大学摘要：提高油藏的采收率是强烈的欲望。

通信原理基于matlab的计算机仿真源代码

例错误!文档中没有指定样式的文字。

-1%周期信号（方波）的展开,fb_jinshi.mclose all;clear all;N=100; %取展开式的项数为2N＋1项T=1;fs=1/T;N_sample=128; %为了画出波形，设置每个周期的采样点数dt = T/N_sample;t=0:dt:10*T-dt;n=-N:N;Fn = sinc(n/2).*exp(-j*n*pi/2);Fn(N+1)=0;ft = zeros(1,length(t));for m=-N:Nft = ft + Fn(m+N+1)*exp(j*2*pi*m*fs*t);endplot(t,ft)例错误!文档中没有指定样式的文字。

-4利用FFT计算信号的频谱并与信号的真实频谱的抽样比较。

脚本文件T2F.m定义了函数T2F，计算信号的傅立叶变换。

function [f,sf]= T2F(t,st)%This is a function using the FFT function to calculate a signal's Fourier %Translation%Input is the time and the signal vectors,the length of time must greater %than 2%Output is the frequency and the signal spectrumdt = t(2)-t(1);T=t(end);df = 1/T;N = length(st);f=-N/2*df:df:N/2*df-df;sf = fft(st);sf = T/N*fftshift(sf);脚本文件F2T.m定义了函数F2T，计算信号的反傅立叶变换。

function [t st]=F2T(f,sf)%This function calculate the time signal using ifft function for the input %signal's spectrumdf = f(2)-f(1);Fmx = ( f(end)-f(1) +df);dt = 1/Fmx;N = length(sf);T = dt*N;%t=-T/2:dt:T/2-dt;t = 0:dt:T-dt;sff = fftshift(sf);st = Fmx*ifft(sff);另写脚本文件fb_spec.m如下：%方波的傅氏变换, fb_spec.mclear all;close all;T=1;N_sample = 128;dt=T/N_sample;t=0:dt:T-dt;st=[ones(1,N_sample/2), -ones(1,N_sample/2)]; %方波一个周期subplot(211);plot(t,st);axis([0 1 -2 2]);xlabel('t'); ylabel('s(t)');subplot(212);[f sf]=T2F(t,st); %方波频谱plot(f,abs(sf)); hold on;axis([-10 10 0 1]);xlabel('f');ylabel('|S(f)|');%根据傅氏变换计算得到的信号频谱相应位置的抽样值sff= T^2*j*pi*f*0.5.*exp(-j*2*pi*f*T).*sinc(f*T*0.5).*sinc(f*T*0.5);plot(f,abs(sff),'r-')例错误!文档中没有指定样式的文字。

基于Matlab的导线网平差程序设计

基于Matlab的导线网平差程序设计
李建章
【期刊名称】《兰州交通大学学报》
【年(卷),期】2010(029)004
【摘要】导线网数据量大,网形复杂多变,其数据处理过程大多涉及到矩阵的计算.利用VC、VB等编程语言进行导线网程序的开发,算法比较复杂.基于Matlab平台,利用其强大的矩阵处理能力,设计出导线网数据结构,此基础上进行导线网平差程序的设计与开发,减小了代码编写的工作量.
【总页数】3页(P88-90)
【作者】李建章
【作者单位】兰州交通大学,土木工程学院,甘肃,兰州,730070
【正文语种】中文
【中图分类】P209
【相关文献】
1.基于Excel和MATLAB的导线网平差 [J], 孙庆华;熊伟;张宪柱;王磊
2.基于 MATLAB 的水准网平差程序设计 [J], 王鹏磊;刘长星
3.基于MATLAB的改进导线网平差程序设计 [J], 王鹏磊
4.基于 MATLAB 的导线网平差软件设计及误差椭圆的绘制 [J], 翟敏;陶秋香
5.一种基于MATLAB的改进的水准网平差程序设计与实现 [J], 李亮亮;郭恒林;王利华
因版权原因，仅展示原文概要，查看原文内容请购买。

唐昌建（MATLAB编程基础及应用）第六章

唐昌建（MATLAB编程基础及应⽤）第六章第六章数值计算⽅法与数据第⼀节⾮线性数值计算⼀、⾮线性微分⽅程的求解ordinary differential equation （ODE ） 1.基本问题（1）⾮线性问题。

（2）基本命令：⾮刚性常微分⽅程（组）ode23,ode45,ode113 2.基本解法：（1）建⽴标准微分⽅程（组）(,)d yF y t dx =⽬的：将所有不同变量归为同⼀变量。

⾼阶微分⽅程归为⼀阶微分⽅程。

例如：y''+2y'-y=0 令 y 1=y;y 2=y'原式变为： y 2'=y 1-2y 2;y 1'=y 2（2）建⽴ode 相应的函数m ⽂件，格式如function ff=fun(t,y ) ff=[(,)F y t ]（3）调⽤fun ⽤ode 函数求解可在命令窗进⾏也可以在m ⽂件中进⾏。

ode 的调⽤格式：[t,y ]=ode45('fun',[⾃变量的范围],[初值列阵]) （4）例⼦练习 6-1（1）2222dyy x xdx =-++ y(0)-1(fun1 exno19)（2）捕⾷者与被捕⾷者函数关系（x 是被捕⾷者）0.01dx x xy dt =+ 0.02dyx xy dt =-+ x(0)=30 y(0)=20(fun2exno20)（3）222(1)0d y dy y y dx dx µ=-+=(fun3 exno21)(0)0dydx = y(0)=1变换为：12dy y dx = 22121(1)dy y y y dx µ=-+y 1(0)=1 y 2(0)=0fun1function dy=fun1(x,y) dy=[-2*y+2*x.^2+2*x]e19[x,y]=ode23('fun1',[0,0.5],1)%⼀步⼀步的解，相当占⽤资源 plot(x,y)fun2function f=fun2(t,y)f=[y(1)-0.01*y(1)*y(2);-y(2)+0.02*y(1)*y(2)]e20ts=0:0.1:20 y0=[30,20][t,y]=ode45('fun2',ts,y0) [t,y]subplot(221)plot(t,y(:,1),'r',t,y(:,2),'b'),gtext('x'),gtext('y') subplot(222)plot(y(:,1),y(:,2)),gtext('x'),gtext('y')fun3function f=fun3(t,y)f=[y(2);2*(1-y(1)^2)*y(2)-y(1)]e21（⾃编）ts=0:0.1:20y0=[1,0][t,y]=ode113('fun3',ts,y0)[t,y]subplot(221)plot(t,y(:,2),'r',t,y(:,1),'b')⼆、⾮线性函数的最⼩值1.单变量（可代参数）a)[x,y,po1,po2]=fminbnd(fun,x1,x2,options,p1,p2,...)b)fun:可带参数的⽬标函数c)x2,x2:⾃变量的求值范围d)options：输⼊的优化参数，包括四个域display,maxfunevals,maxlter,tolx.这四个域由函数optimset()来传递。

基于MATLAB的高程控制网平差系统的设计与应用

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

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

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

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

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

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

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

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

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

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

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

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

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

近似坐标计算的函数-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=[];%第一步：寻找观测条件：两种情况：一是有已知方位角；二是由两个已知点及方向观测值推出方位角。

temp7=find(t==i);if length(temp7)>0temp8=find(xyknow==s(temp7(1)));if length(temp8)>0temp9=find(e==s(temp7(1))&d==t(temp7(1))|e==t(temp7(1))&d==s(temp7(1)));if length(temp9)>0 %第一种情况：有已知方位角（一般适用于闭合导线）S=mean(sid(temp9));Alfa0=az(temp7(1));x0(i)=x0(s(temp7(1)))+S*cos(Alfa0);y0(i)=y0(s(temp7(1)))+S*sin(Alfa0);way=1;method(i)=11;%----->由已知方位角算出endendelsetemp1=find(f==i);%第二种情况：由两个已知点及方向观测值推出方位角（适用于附合、支导线）if length(temp1)>0for j=xyknowtemp2=find(g(temp1)==j);if temp2>0temp3(end+1)=temp1(temp2);endend %找到所有以已知点为起点，以i为终点的方向值位置，并存放于temp3中if length(temp3)>0for j=temp3temp4=find(e==g(j)&d==f(j)|e==f(j)&d==g(j));if temp4>0 %找到一条观测边temp5=g(j);for k=xyknowtemp6=find(g==temp5&f==k);if temp6>0 %找到另一个方向值位置A=k;B=temp5;S=mean(sid(temp4));dir1=dir(temp6);dir2=dir(j);way=1;method(i)=12;%由两个已知点及方向观测值算得break;endendendendendend%第二步：按极坐标公式计算坐标if way==1deltax=x0(B)-x0(A);deltay=y0(B)-y0(A);a1=alfa(deltax,deltay);a2=a1+dir2-dir1;if a2<0a2=a2+2*pi;endif a2>2*pia2=a2-2*pi;endx0(i)=x0(B)+S*cos(a2);y0(i)=y0(B)+S*sin(a2);endend%%================================================================ ======%方法2：前方搜索与计算way=2%思路：找出合适的A,B点，然后对A,B,P进行逆时针排序，照公式计算坐标if way==0%第一步：搜索合适的A,B点temp1=[];temp2=[];temp3=[];temp4=[];temp5=[];temp6=[];temp7=[];temp8=[];A=[];B=[];P=[];temp1=find(f==i); %找出终点为i的方向值位置for j=xyknowtemp2=find(g(temp1)==j);if temp2>0temp3(end+1)=temp1(temp2);endend %找出起点为已知点，终点为目标点未知点的方向值位置，记录在temp3中if length(temp3)>1for j=g(temp3)for k=g(temp3)temp4=find(g==j&f==k);if temp4>0temp5(end+1)=temp4;%temp5中记录了所有与目标未知点有观测的已知点之间方向观测值的位置endendif length(temp5)>0for j=temp5temp6=find(g(temp5)==f(j)&f(temp5)==g(j));if temp6>0way=2;method(i)=2;P=i;A=g(i);B=f(j);break;endendend %找出一组合适的A,B即可，退出循环endif way==2 %第二步：对A,B,P进行逆时针排序[A,B,P]=reorderunknow(A,B,P,1,g,f,dir);bet1=beta(P,A,B,g,f,dir);bet2=beta(P,B,A,g,f,dir);%%第三步：代入公式计算坐标if(bet1~=pi/2)&(bet2~=pi/2)x0(P)=(x0(A)*tan(bet1)+x0(B)*tan(bet2)+(y0(B)-y0(A))*tan(bet1)*tan(bet2))/(tan(bet1)+tan(bet2 ));y0(P)=(y0(A)*tan(bet1)+y0(B)*tan(bet2)+(x0(A)-x0(B))*tan(bet1)*tan(bet2))/(tan(bet1)+tan(bet2 ));elsesidab=sqrt((x0(A)-x0(B))^2+(y0(A)-y0(B))^2);alfaab=alfa(x0(B)-x0(A),y0(B)-y0(A));alfaap=alfaab-bet1;if alfaap<0alfaap=alfaap+2*pi;endif bet1==pi/2sidap=tan(bet2)*sidab;elsesidap=sidab/cos(bet1);endx0(P)=x0(A)+sidap*cos(alfaap);y0(P)=y0(A)+sidap*sin(alfaap);endend %for if wey==2end % for if way==0%%============================================================ ======%方法3 ：测边交会与计算way=3,思路：三种算法：一是找到一个端点为所求目标未知点的三条边%二是寻求方向观测信息以对A,B,P进行逆时针排序，三是给出A,B,P的排列顺序。

if way==0temp1=[];temp2=[];temp3=[];temp4=[];temp5=[];temp6=[];temp7=[];temp8=[];i1=0;A=[];B=[];C=[];P=[];%第一步：搜索合适的A,B,C点或带有方向观测值的A，B点temp1=find(e==i|d==i);if length(temp1)>0for j=xyknowtemp2=find(e(temp1)==j|d(temp1)==j);if temp2>0temp3(end+1)=temp1(temp2);endendif length(temp3)>2 %第一种: 找到至少带有已知点的3条边，不需排序即可进行计算way=3;i1=1;P=i;if e(temp3(1))==iA=d(temp3(1)); %AelseA=e(temp3(1));endif e(temp3(2))==iB=d(temp3(2)); %BelseB=e(temp3(2));endif e(temp3(3))==iC=d(temp3(3)); %CelseC=e(temp3(3));ends1=sid(temp3(1));s2=sid(temp3(2));s3=sid(temp3(3));elseif length(temp3)==2 %只能找到两条带有已知点的边P=i;if e(temp3(1))==iA=d(temp3(1)); %AelseA=e(temp3(1));endif e(temp3(2))==iB=d(temp3(2)); %BelseB=e(temp3(2));end%第二种：需进一步查找方向观测值信息，以对A,B,P进行逆时针排序n1=find(g==A&f==P);n2=find(g==A&f==B);n3=find(g==B&f==P);n4=find(g==B&f==A);n5=find(g==P&f==A);n6=find(g==P&f==B);if n1>0&n2>0way=3;i1=2;[A,B,P]=reorderunknow(A,B,P,1,g,f,dir);elseif n3>0&n4>0way=3;i1=2;[B,A,P]=reorderunknow(B,A,P,1,g,f,dir);elseif n5>0&n6>0way=3;i1=2;[A,B,P]=reorderunknow(A,B,P,2,g,f,dir);else %第三种：需给出A,B,P的排列顺序temp4=find(cc==P);if length(temp4)>0A=aa(temp4(1));B=bb(temp4(1));way=3;i1=2;elseaa0(end+1)=A;bb0(end+1)=B,cc0(end+1)=P;endendendend%第二步：坐标计算if way==3&i1==1 %三个已知点，三条边长method(i)=31;deltx=x0(B)-x0(A);delty=y0(B)-y0(A);alfaAB=alfa(deltx,delty);ss=sqrt(deltx*deltx+delty*delty);ee=(s1*s1+ss*ss-s2*s2)/(2*ss);ff=sqrt(s1*s1-ee*ee);xp1=x0(A)+ee*cos(alfaAB)+ff*sin(alfaAB);yp1=y0(A)+ee*sin(alfaAB)-ff*cos(alfaAB);xp2=x0(A)+ee*cos(alfaAB)-ff*sin(alfaAB);yp2=y0(A)+ee*sin(alfaAB)+ff*cos(alfaAB);s4=sqrt((xp1-x0(C))^2+(yp1-y0(C))^2);s5=sqrt((xp2-x0(C))^2+(yp2-y0(C))^2);if abs(s4-s3)<abs(s5-s3)x0(P)=xp1;y0(P)=yp1;elsex0(P)=xp2;y0(P)=yp2;endelseif way==3&i1==2 %由方向观测值确定逆时针顺序后坐标计算q1=find(e==A&d==P|e==P&d==A);s1=sid(q1);q2=find(e==B&d==P|e==P&d==B);s2=sid(q2);method(i)=32;deltx=x0(B)-x0(A);delty=y0(B)-y0(A);alfaAB=alfa(deltx,delty);ss=sqrt(deltx*deltx+delty*delty);ee=(s1*s1+ss*ss-s2*s2)/(2*ss);ff=sqrt(s1*s1-ee*ee);x0(P)=x0(A)+ee*cos(alfaAB)+ff*sin(alfaAB);y0(P)=y0(A)+ee*sin(alfaAB)-ff*cos(alfaAB);endend %if way==0%=============================================================== =====%方法4：后方交会搜索与计算way=4%思路：1.找到三个已知点；2.判断是否是危险圆；3.坐标计算%第一步：找点if way==0temp1=[];temp2=[];temp3=[];temp4=[];temp5=[];temp6=[];A=[];B=[];C=[];P=[];temp1=find(g==i);%找出起点为P的方向值位置if length(temp1)>2for j=xyknowtemp2=find(f(temp1)==j);if temp2>0temp3(end+1)=temp1(temp2);%记录起点为P终点为已知点的方向值位置endendif length(temp3)>2 %找到三个或者三个以上的点A=f(temp3(1));B=f(temp3(2));P=i;for j=3:length(temp3)C=f(temp3(j));[A,B,C]=reorderknow(A,B,C,x0,y0);%将A,B,C排成顺时针方向sa=sqrt((x0(C)-x0(B))^2+(y0(C)-y0(B))^2);sb=sqrt((x0(C)-x0(A))^2+(y0(C)-y0(A))^2);sc=sqrt((x0(A)-x0(B))^2+(y0(A)-y0(B))^2);angleA=acos((sb*sb+sc*sc-sa*sa)/(2*sb*sc));angleB=acos((sa*sa+sc*sc-sb*sb)/(2*sa*sc));angleC=acos((sa*sa+sb*sb-sc*sc)/(2*sa*sb));temp4=find(g==i&f==A);RA=dir(temp4);temp5=find(g==i&f==B);RB=dir(temp5);temp6=find(g==i&f==C);RC=dir(temp6);alfa1=RC-RB;beta1=RA-RC;gama1=RB-RA;alfa2=beta(B,P,C,g,f,dir);beta2=beta(A,P,C,g,f,dir);gama2=beta(A,P,B,g,f,dir);%%第二步：判断四点是否共圆if(angleA==alfa2&angleB==beta2)|(angleA==alfa2&angleC==gama2)|(angleC==gama2&angleB ==beta2)continue %危险圆！！！！！！！！！！！！！！！！else%第三步：非危险圆，坐标计算way=4;method(i)=4;PA=tan(alfa1)*tan(angleA)/(tan(alfa1)-tan(angleA));PB=tan(beta1)*tan(angleB)/(tan(beta1)-tan(angleB));PC=tan(gama1)*tan(angleC)/(tan(gama1)-tan(angleC));x0(P)=(PA*x0(A)+PB*x0(B)+PC*x0(C))/(PA+PB+PC);y0(P)=(PA*y0(A)+PB*y0(B)+PC*y0(C))/(PA+PB+PC);break;endend % for j=3:length(temp3)end %if length(temp3)>2end %if length(temp1)>2end %if way==0%=============================================================== =====%求出一个未知点的近似坐标后终止循环，更新已知点点号和未知点点号数组if way>0temp=[];xyknow(end+1)=i;temp=find(xyunknow==i);xyunknow(temp:end-1)=xyunknow(temp+1:end);xyunknow=xyunknow(1:end-1);break; %%%%%%% for i=xyunknowend % if way>0end %for i=xyunknow%================================================================= =====%方法5.倘若对所有未知点循环一次的过程中一个也未求出，则考虑用无定向导线计算%way==5%第一种情况：已知条件为两个或者两个以上点的坐标if(length(xyknow)==prelength)&(way==0)&(length(xyknow)>=2)siddir=0;seeksid=0;seekdir=0;for j=xyknowfor k=xyunknowtemp1=[];temp1=find(e==j&d==k|e==k&d==j);if(length(temp1)>0) %找到与已知点联测的一条边的另一端点seeksid=1;siddir=1;break;endendif seeksid==1break;endendif seeksid==0 %未找到与已知点联测的一条边的另一端点，则找方向观测的另一个端点for j=xyknowfor k=xyunknowtemp1=[];temp1=find(g==j&f==k|g==k&f==j);if(length(temp1)>0) %找到方向观测的另一端点seekdir=1;seekdir=1;break;endendif seekdir==1break;endendendif siddir==1 % (length(temp1)>0)trueX0=x0;trueY0=y0;x0=[];y0=[];xyknow=[];xyunknow=[];x0(j)=0.0;y0(j)=0.0; % 能用此方法解决，建立新的坐标系if seeksid==1x0(k)=mean(sid(temp1));y0(k)=0;elseif seekdir==1x0(k)=1000;y0(k)=0;endnon_orient=1; %更新已知点数组和未知点数组xyknow=[j k];point(find(point==k|point==j))=[];xyunknow=point';end%============================================================== =====%第二种情况：已知条件为一个已知点和一个方位角，但分散开。