matlab多面函数拟合法

%-----------------多面函数拟合法,至少有五个公共点-------------

clc;clear all;

[P1,F1]=uigetfile('*.txt','打开GPS水准点平面坐标文件');

fnal1=strcat(F1,P1);fpath=P1;

fid1=fopen( fnal1,'r');

lineNum=0;

zbjz=[];

format long e

while (~feof(fid1))

curline = fgetl(fid1);

lineNum = lineNum + 1;

zbjz(lineNum,1)=1;

zbjz(lineNum,2)=str2num(curline(5:16));

zbjz(lineNum,3)=str2num(curline(18:28));

zbjz(lineNum,4)=str2num(curline(30:36));

end

zbjz;

Y=zbjz(:,2);

X=zbjz(:,3);

Z=zbjz(:,4);

o=zbjz(:,1);

n=length(X);

[P2,F2]=uigetfile('*.txt','打开GPS水准点大地高文件');

fnal2=strcat(F2,P2);fpath1=P2;

fid2=fopen( fnal2,'r');

LineNum=0;

ddg=[];

while (~feof(fid2))

LineNum = LineNum + 1;

urline = fgetl(fid2);

ddg(LineNum,1)=str2num(urline(5:10));

end

ddg;

gcyc=ddg-Z;

jl1=[];

for i=1:n

for j=1:n

jl1(j,i)=sqrt((X(i)-X(j))^2+(Y(i)-Y(j))^2);

end

end

jl1;

aa=max(jl1);

aaa=max(aa);

b=aaa^2;

k=0;

for i=1:n

for j=1:n

B(j,i)=1-((X(i)-X(j))^2+(Y(i)-Y(j))^2+k^2)/b;

end

end

B;

% A1=pinv(S'*S)*S'*gcyc; %水准点和已知点相同A=pinv(B)*gcyc;

A;

for i=1:n

for j=1:n

Qp(j,i)=1-((X(i)-X(j))^2+(Y(i)-Y(j))^2+k^2)/b;

end

end

Qp;

yzdgcyc=Qp*A;

yzdnih=ddg-yzdgcyc;

v1=(Z-yzdnih)*1000;

shuizhun='已知水准点拟合高程(m)为';

fprintf('%s\n',shuizhun)

fprintf('%7.4f\n',yzdnih)

cancha1='已知点拟合残差(mm)为';

fprintf('%s\n',cancha1)

fprintf('%6.1f\n',v1)

neifuhe='内符合精度(mm)';

fprintf('%s\n',neifuhe)

miu=sqrt(sum(v1.*v1)/(n-1));

fprintf('%6.1f\n',miu)

[P,F]=uigetfile('*.txt','打开检核点平面坐标文件');

fnal=strcat(F,P);

fpath=P;

fid=fopen( fnal,'r');

lineNum=0;

jhzbjz=[];

format long e

while (~feof(fid))

curline = fgetl(fid);

lineNum = lineNum + 1;

jhzbjz(lineNum,1)=1;

jhzbjz(lineNum,2)=str2num(curline(5:16));

jhzbjz(lineNum,3)=str2num(curline(18:28));

jhzbjz(lineNum,4)=str2num(curline(30:36)); end

jhzbjz;

Y4=jhzbjz(:,2);

X4=jhzbjz(:,3);

Z4=jhzbjz(:,4);

o4=jhzbjz(:,1);

n2=length(X4);

[P4,F4]=uigetfile('*.txt','打开检核点大地高文件'); fnal4=strcat(F4,P4);

fpath4=P4;

fid4=fopen( fnal4,'r');

LineNum=0;

jhddg=[];

while (~feof(fid4))

LineNum = LineNum + 1;

urline = fgetl(fid4);

jhddg(LineNum,1)=str2num(urline(5:10)); end

jhddg;

gcyc1=jhddg-Z4;

jl11=[];

for i=1:n

for j=1:n

jl11(j,i)=sqrt((X(i)-X(j))^2+(Y(i)-Y(j))^2);

end

end

jl11;

bb=max(jl11);