GPS高程拟合计算原理

GPS高程拟合就是利用在范围不大的区域中，应用GPS测得的坐标（大地坐标或者是平面坐标）和应用水准测量技术测得的水准高程，计算测点的高程异常，而高程异常具有一定的几何相关性这一原理，用数值拟合法，拟合出测区似大地水准面，由此根据测区内的其他点坐标内插出其高程异常，从而求出待定点的正常高。

这里观测量为高程异常，未知数为拟合参数。程序设计思路：将这里观测量为高程异常，未知数为拟合参数。坐标构成系数矩阵；将同名点高程异常存入和系数矩阵分别传入特定的数组，按照最小二乘法求出拟合参数，同时按照间接平差法思路求出拟合点高程精度。并将拟合结果按照文本文件输出输出。原始数据可以直接输入创建矩阵，也可以按文本文件读入，这里按文本文件读入矩阵。

把观测数据编辑成文本文件

点名    B(°′″)          L(°′″)      大地高（m） 水准高程（m）

-----------------------------------------------------------------------------------

TJ01 3418.2458    11110.1922    122.638     111.278

TJ02 3418.2723    11107.0242    119.060     107.676

TJ0     33418.3663    11107.5653    204.956     193.595

TJ04 3418.4284    11108.5949    187.392     176.056

TJ05 3418.4796    11109.1467    181.213     169.888

TJ06 3418.3322    11109.5378    172.604     161.282

TJ07 3418.194    111110.3462    148.991       137.678

-----------------------------------------------------------------------------------

程序清单:

%------------------------------------------

%  read_2 GPS 水准高程拟合    

%-----------------------------------------------

clear

clc

[filename,filepath]=uigetfile('*.txt','TEXT file');

fn=num2str(filename);

fid=fopen(fn);

nk=0;

format long

while 1

    line=fgetl(fid)

       if ~isstr(line)

        break

    end

      nk=nk+1

end

nk;

np=nk-3;

fclose(fid)

fid=fopen(fn);

[ff]=fscanf(fid,'%s',5) %%%read  head of file    

    for i=1:np

id(i)= fscanf(fid,'%d',1);

        bg(i)=fscanf(fid,'%f',1);

lg(i)=fscanf(fid,'%f',1)

        Bg(i)=dms_rad(bg(i))  

        Lg(i)=dms_rad(lg(i))

        BLg(i)= Bg(i)* Lg(i)

        Hg(i)=fscanf(fid,'%g',1);

        hg(i)=fscanf(fid,'%g',1);

        yh(i)=Hg(i)-hg(i)      

end%for i=np

%******************************************************************

B0=sum(Bg)/np

L0=sum(Lg)/np

%转换时 回加 xx yy

for i=1:np

Bg(i)=Bg(i)-B0;

Lg(i)=Lg(i)-L0;

end

%--------------------原始数据

X=[ones(np,1)  Bg'  Lg' ]

Yt=yh';

N=X'*X;

bet=inv(N)*X'*Yt;

V=(X*bet-Yt)*100;

[n t]=size(X);

r=n-t;

sgm2=V'*V/r;

ip=1:np

obs=[id' bg' lg' Hg'  hg'   yh'  V]

%----------------------------------

[filename,filepath]=uiputfile('redata.txt','TEXT file');

ftn=num2str(filename)

fit1=fopen(ftn,'w')

%=========================================================保存近似坐标

fprintf(fit1,'                 高程异常计算结果\n');

fprintf(fit1,' 点号 纬度(°′″)  经度(°′″)  大地高 水准高(m) 高程异常(m) 拟合误差 \n');

fprintf(fit1,'----------------------------------------------------------------------\n');

fprintf(fit1,'%5d%12.3f%12.3f%10.3f%10.3f%10.3f%10.3f\n',obs');

fprintf(fit1,'-----------------------------------------------------------------------\n');

fprintf(fit1,'      ----拟合系数----\n');

fprintf(fit1,'%10.3f%10.3f%10.3f\n',bet);

fprintf(fit1,'  ----------------------------\n');

fclose(fit1);

open('redata.txt');