2024年5月10日发(作者：)

chkdat函数（72页）

function [n1,k]=chkdat(sd,pn,n1)

n=length(n1);

k=0;

for i=1:n

i1=0;

for j=1:sd

if(n1(i)==pn(j))

i1=1;

n1(i)=j;

break;

end

end

if(i1==0)

% fprintf(fit2,'%5d %5dn',i,n1(i)

k=1;

end

end

return

readlevelnetdata函数（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;

end

ed=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:gd

i=k1(k);

j=k2(k);

if(h0(i)<1e4&h0(j)>1e4)

h0(j)=h0(i)+h1(k);

ie=ie+1;

end

if(h0(i)>1e4&h0(j)<1e4)

h0(i)=h0(j)-h1(k);

ie=ie+1;

end

end

if(ie==dd)

break;

end

end

h0=reshape(h0,length(h0),1);

return

bm1函数（75页）

function id=bm1(gd,dd,k1,k2)

%计算一维压缩存放的数组id

id=[];

for i=1:dd

k=i;

for j=1:gd

i1=k1(j);

i2=k2(j);

if(i1==i&i2

k=i2;

end

if(i2==i&i1

k=i1;

end

end

id(i)=k;

end

for i=2:dd

id(i)=id(i-1)+i-id(i)+1;

end

return

一维压缩存储法方程平差（76页）

global pathname filename

global 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);

end

if(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);

else

a(jj+i)=a(jj+i)-p(k);

end

end

end

end

a=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:sd

h(i)=h0(i)+dh(i);

end

vv=0;

for i=1:gd

L(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);

end

u=sqrt(vv/(gd-dd));

for i=1:dd

b(1:dd)=0;

b(i)=1.0;

b=cy6(a,b,id,dd,i);

hm(i)=sqrt(b(i))*u;

end

writelevelnetdata(pn,k1,k2,h1,v',L',h0,dh',h',hm',u);

gs5函数（77页）

function a=gs5(dd,a,id)

%变带宽高斯消去法

for i=1:dd

ii=id(i)-i;

if(i-1==0)

li=1-ii;

else

li=id(ii-1)-ii+1;

end

for j=li:i

jj=id(j)-j;

if(j-1)==0

lj=1-jj;

else

lj=id(j-1)-jj+1;

end

lk=li;

if(li

lk=lj;

end

for k=lk:j-1

kk=id(k);

a(ii+j)=a(ii+j)-a(ii+k)/a(kk)*a(jj+k);

end

end

end

return

cy6函数（78页）

function b=cy6(a,b,id,dd,k1)

%常数项约化与回代子程序

for i=k1:dd

ii=id(i)-i;

if(i==1)

nd=id(i);

else

nd=id(i)-id(i-1);

end

e=0;

for k=1:i-1

if((i-k)

e=e+a(ii+k)*b(k);

end

end

b(i)=(b(i)-e)/(a(ii+i);

end

for i=dd-1:-1:k1

ii=id(i);

for k=i+1:dd

kk=id(k)-k;

nk=id(k)-id(k-1);

if(k-i

b(i)=b(i)-a(kk+i)/a(ii)*b(k);

end

end

end

return

上三角存储法方程平差程序（79页）

mm=(dd+1)*dd/2;

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=(i-1)*(dd-i/2);

a(ii+i)=a(ii+i)+1/s(k);

b(i)=b(i)+1./s(k)*h1;

end

if(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

a(ii+j)=a(ii+j)-1/s(k);

else

a(jj+i)=a(jj+i)-1/s(k);

end

end

end

end

a=invsqr(a,dd);

for i=1:dd

dh(i)=0;

di=(i-1)*(dd-i/2);

for j=1:dd

dj=(j-1)*(dd-j/2);

if(j

dh(i)=dh(i)-a(dj+i)*b(j);

else

dh(i)=dh(i)-a(di+j)*b(i);

end

end

end

dh(dd+1:ed+dd)=0;

hm(dd+1:ed+dd)=0;

for i=1:sd

h(i)=h0(i)+dh(i);

end

vv=0;

for i=1:gd

L(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);

end

uw0=sqrt(vv/(gd-dd));

for i=1:dd

ii=(i-1)*(dd-i/2);

hm(i)=sqrt(a(ii+i))*uw0;

end

return

输出数据函数（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;

end

fprintf(fid2,'待定点高程平差值及中误差：n 点号--近似高程(m)--高程改正(m)-高程平差值

(m)-中误差n');

fprintf(fid2,'%5d %10.4f %10.4f %10.4f %10.4fn',[pn,h0,dh,h,hm]');

fprintf(fid2,'高差观测值平差值:n -点号---点号--观测高差(m)--高差改正(m)-平差高差(m)n');

fprintf(fid2,'%5d %5d %10.4f %10.4f %10.4fn',[pn(k1),pn(k2),h1,v,L]');

fprintf(fid2,'单位权中误差：%10.4fmn',uw0);

% open(strcat(pathname1,filename1));

fclose(fid2);

return

利用Matlab矩阵运算的平差程序（81页）

function level3

tic

disp('平差已经开始---->>>>')

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

return

水准网平差过程（81页）

function [dh,h,V,L,uw0,uwh,uwl]=calculatelevelnet(ed,dd,sd,gd,pn,h0,k1,k2,h1,s)

%水准网平差

A=sparse(zeros(sd,gd));%求解系数阵

b=(0:(gd-1))*sd;

A(k1'+b)=-1;

A(k2'+b)=1;

A=A';

A=A(:,1:dd);

disp('-------水准网间接平差误差方程的系数阵:------')

A=full(A)

l=zeros(gd,1); %求解常数项

disp('-------误差方程的常数项:-------')

l=h0(k1)-h0(k2)+h1

disp('-------权阵:-------')

p=diag(1./s)%权阵

dh=inv(A'*p*A)*A'*p*l ;%高差改正数

h00=h0(dd+1:sd);

h0=h0(1:dd); %待定点高差观测值

h=h0+dh;%待定点高程平差值

V=A*dh-l; %高差观测值改正数

L=h1+V; %高差观测值平差值

%精度评定

uw0=sqrt(V'*p*V/(gd-dd));%单位权中误差

Qxx=inv(A'*p*A);

uwh=uw0*sqrt(diag(Qxx));%待定点高程平差值中误差

uwh(dd+1:ed+dd)=0.0;

Qff=A*Qxx*A';

uwl=uw0*sqrt(diag(Qff));%高差平差值中误差

h=[h;h00];%所有点高程

h0=[h0;h00];

dh=[dh;zeros(ed,1)];

return