Matlab程序设计在GPS高程拟合中的应用

骆丽华，覃 辉

（1.桂林理工大学 广西矿冶与环境科学实验中心，广西 桂林 541004；2. 广西空间信息与测绘重点实验室，广西 桂林 541004；3.桂林理工大学 测绘地理信息学院，广西 桂林 541004；4.广东科学技术职业学院 , 广东 珠海 519090）

Matlab程序设计在GPS高程拟合中的应用

骆丽华1,2,3，覃 辉4

（1.桂林理工大学 广西矿冶与环境科学实验中心，广西 桂林 541004；2. 广西空间信息与测绘重点实验室，广西 桂林 541004；3.桂林理工大学 测绘地理信息学院，广西 桂林 541004；4.广东科学技术职业学院 , 广东 珠海 519090）

介绍了二次曲面拟合法、多面函数拟合法以及Matlab自带的griddata拟合插值函数法。计算结果表明，Matlab强大的数据处理能力运用在GPS高程拟合中，处理方式简洁，计算量少，精度可靠，同时具有立体可视特征。

Matlab；GPS高程拟合；二次曲面；多面函数；griddata插值

GPS坐标是以地球质心为坐标原点的地心坐标，测量获得的高程是以WGS84椭球面为基准面的大地高[1]。而在我国工程中采用的高程通常是以似大地水准面为基准面的正常高，大地高与正常高之差称为高程异常ζ。因此，在实际测量工作中，大地高和正常高之间的转换就变得十分有意义，只要求出高程异常就可以方便地求得正常高。通常求取高程异常的方法是通过某种数学模型根据已知的GPS三维点位模拟求出未知点的高程异常。由于Matlab突破了高级编程语言要编写大量循环语句的思想，它只需通过一个或几个简单的命令就可以完成矩阵或数组的运算，所以文中结合实例利用Matlab进行程序设计，对局部区域的高程异常进行拟合，实现了数据的自动处理。

1 GPS高程拟合原理

1.1 曲面拟合法

曲面拟合法，就是将高程异常近似地看作是一定范围内各点坐标的曲面函数，认为高程异常在此范围内变化平缓，可以用已联测水准的GPS点的高程异常来拟合这一函数。在求得函数的拟合系数之后，就可用这一拟合函数来计算其他GPS点的高程异常和正常高[2]。设某公共点的高程异常ζ与其平面坐标x、y之间有以下关系：

式中，ζ为高程异常；f (x，y)为ζ的趋势值；ε为误差值。

二次多项式曲面f (x，y)可表示为：

写成矩阵形式有：

误差方程为：

对于每一个已知点，均可以列出上式方程，在∑ε2=min条件下，可求解系数阵：

再由己知高程异常的权阵情况，上式可写为：

系数求出后，再按式（2）求出待求点的正常高。

1.2 多面函数曲面拟合法

多面函数法的基本思想是对于任何数学表面或任何不规则的圆滑表面，总可以用一系列有规则的数学表面的总和，以任意精度逼近[3]。根据这样的思想，高程异常函数可以表示为：

式中，ai为待定系数；Q（x，y，xi，yi）为核函数。核函数是关于x、y的函数，核函数的中心在（xi，yi）处。理论上核函数是可以任意构造的，在实际应用中，通常用以下几种函数来充当核函数。

锥面：

双曲面：

倒曲面：

在上述各式中，x、y表示内插点坐标；xi、yi表示已知点的坐标；核函数中的表示内插点到已知点的水平距离；δ为光滑系数。

以核函数为双曲面为例，说明多项式曲面拟合法的具体求解过程。设测区内已知点个数为n，求解式（7）中的系数（a1，a2，a3，…，an），其矩阵形式为：

其 中，ζ = （ζ1，ζ2，ζ3， …，ζn）T； B= （a1，a2，a3，…，an）T。由此方程组可解得系数（a1，a2，a3，…，an）的唯一解：

求解未知点的高程异常值，根据式（8）和式（9）即可得到求解公式：

在选择多面函数求解高程异常值的时候，需要注意核函数的选取问题。由于是自主取值，为了达到最佳拟合效果，就要逐步试验进而改进，然后选定一个最佳取值。同时，通过改变平滑因子[4]，可以使计算结果达到最理想的状态。因此，没有必要去关心平滑因子值和数据的具体联系，只需要找到一个合适的值使得计算结果处于最佳状态。

下文实验中确定平滑因子的具体方法是：在程序中将平滑因子设计为变量，初始值设为0，步长为0.000 01（用户可以自定义步长）。当所求的外符合精度不在要求精度范围内，程序会继续调节平滑因子。每次调节平滑因子后，系统进行一次运算，并求得该次运算后的外符合精度，如此循环运算，直到外符合精度达到要求。实验中平滑因子确定为6 498.860 11。

首先,对60例患者的病因构成及一般资料统计分析中显示,60例非炎症性肠病大肠溃疡患者中,肠结核患者16例,感染性肠炎患者10例,吻合口发生溃疡患者4例,溃疡性大肠癌患者19例,缺血性肠炎患者5例,放射性肠炎患者2例,内痔术后溃疡患者2例,急性阑尾炎患者2例;其中,以溃疡性大肠癌以及肠结核患者数量最多,其次分别为感染性肠炎、缺血性肠炎以及吻合口溃疡等情况患者,这也是导致非炎症性肠病大肠溃疡的主要病因。

1.3 Matlab的griddata插值函数法

Matlab中的griddata函数可以将位于同一空间坐标系下的散点插值为规则格网，提供了以下4种方法：线性内插法、三次多项式内插法、最邻近点内插法和Matlab自带的一种圆滑插值法（V4）,可方便地实现结合邻近离散点分布特征的光滑曲面拟合[5]。

2 Matlab程序实现

以下程序分别是二次曲面拟合法、多面函数拟合法以及griddata插值法的主要程序部分[6-8]。其中的x1,y1,p1,q1,x2,y2,p2,q2已通过Excel输入并调入Matlab中备用。

2.1 二次曲面拟合法

%二次曲面拟合法（文件名：ecqmnh.m）

function F=myfunn1(beta,x)

x=x1-mean(x1);y=y1-mean(y1);%数据中心化

myfunn1=inline('beta(1)+beta(2)*x(:,1)+beta(3)*x(:,2 )+beta(4)*x(:,1).*x(:,1)+beta(5)*x(:,1).*x(:,2)+beta(6)*x(:, 2).*x(:,2)','beta','x');%二次曲面高程模型

z=p-q; %已知点高程异常矩阵

x=[x y];

beta=nlinf i t(x,z,myfunn1,[1 1 1 1 1 1]);%模型参数

zi=beta(1)+beta(2)*xi+beta(3)*yi+beta(4)*xi.*xi+bet a(5)*xi.*yi+beta(6)*yi.*yi);%待拟合点高程异常

z2=[z-zi]；%残差

2.2 多面函数拟合法

%多面函数拟合合法（文件名：dmhsnh.m）

z=p-q;%已知点高程异常

k=6498.86011;%平滑因子

W12=sqrt((x(1)-x(2))^2+(y(1)-y(2))^2+k^2);

W13=sqrt((x(1)-x(3))^2+(y(1)-y(3))^2+k^2);……

W110=sqrt((x(1)-x(10))^2+(y(1)-y(10))^2+k^2);

W23=sqrt((x(2)-x(3))^2+(y(2)-y(3))^2+k^2);

W24=sqrt((x(2)-x(4))^2+(y(2)-y(4))^2+k^2);…….

W210=sqrt((x(2)-x(10))^2+(y(2)-y(10))^2+k^2);……

W89=sqrt((x(8)-x(9))^2+(y(8)-y(9))^2+k^2);

W810=sqrt((x(8)-x(10))^2+(y(8)-y(10))^2+k^2);

W910=sqrt((x(9)-x(10))^2+(y(9)-y(10))^2+k^2);

W11=k;

W21=W12;W22=k;

W31=W13;W32=W23;W33=k;

……

W 1 0 1=W 1 1 0;W 1 0 2=W 2 1 0;……W108=W810;W109=W910;W1010=k;

B=[W11 W12 W13 W14 W15 W16 W17 W18 W19 W110;

W21 W22 W23 W24 W25 W26 W27 W28 W29 W210;

……

W101 W102 W103 W104 W105 W106 W107 W108 W109 W1010];

A=inv(B'*B)*B'*z;%模型参数

f11=sqrt((xi-x(1))^2+(yi)-y(1))^2+k^2);

f12=sqrt((xi-x(2))^2+(yi-y(2))^2+k^2);

……

f110=sqrt((xi)-x(10))^2+(yi-y(10))^2+k^2);

Q=[f11 f12 f13 f14 f15 f16 f17 f18 f19 f110];

I=Q*A %任一点高程异常

2.3 Griddata插值法

%griddata插值法（文件名：griddata.m）

[xi,yi]=meshgrid(min(x): dx:max(x),min(y): dy:max(y)); %生成规格化格网

[xi,yi,zi]=griddata(x,y,z, xi,yi,’method’); % 高 程异常曲面拟合

zi=griddata(x,y,z,xi,yi,’method’); %求待拟合点的高程异常

surf(xi,yi,zi); %高程异常曲面图

3 实例分析

3.1 实验数据来源

本文数据源于文献[9]某个沿江且地势平坦地区，面积约为8 km2。该GPS控制网的17个水准高程点都是等精度且不含粗差，点位分布如图1所示。本实验将这17个点分为2组，前10个点（蓝色）作为实验的模型拟合点，后7个点（红色）数据作为检核数据点。利用Matlab分别用二次曲面拟合法、多面函数拟合法和griddata插值法编程实现，各种方法计算的结果如表1～4所示。

3.2 数据精度分析

图1 点位分布图

表1 二次曲面拟合结果

表2 多面函数拟合结果

表3 Matlab三次多项式内插（cubic）拟合结果

表4 Matlab’V4’法拟合结果

由于高程拟合大多采用数学模型进行插值计算，对于模型误差很难确定，所以GPS水准计算的精度一般利用内、外符合精度进行评定。Griddata函数进行高程拟合，由于拟合曲面通过各拟合点，故拟合点的内符合精度为0。本文为了便于与二次曲面及多面函数拟合法比较，取插值点的单点最大偏差及外部符合精度为评价指标。

根据核检点ζi与拟合值ζi'之差，按下式可计算出外符合精度m[10]：

式中，v为拟合残差；n 为检核点数。经过计算，各方法的精度指标如表5所示。

表5 几种拟合方法的精度对比

由表5可以看出，二次曲面拟合检核点最大残差为3.12 mm，多面函数拟合检核点最大残差为9.09 mm，Matlab三次多项式内插法最大残差为6.13 mm，Matlab自带的V4插值法最大残差为3.83 mm。二次曲面拟合效果最好，能达到三等水准测量的精度要求；多面函数法拟合效果最差，其原因可能是核函数选取和平滑因子取值不是最理想状态，但能达到四等水准测量精度；同样Griddata插值函数拟合法也能达到四等水准测量的精度要求。由于似大地水准面是一个不规则的连续曲面，二次曲面拟合的高程异常成果较好。从表5也反映出Matlab的griddata函数中V4法相对三次多项式内插拟合法精度也有一定提高。根据外符合精度值可以看出，用Matlab的拟合插值函数进行GPS高程拟合也具有较好的外推能力，完全可以满足大多数工程的要求。 Matlab的图像处理功能[11]可以清楚地显示拟合区域的高程异常分布。图2是以Matlab的3次多项式拟合的高程异常空间分布图。

图 2 高程异常拟合曲面示意图

4 结 语

运用Matlab进行高程异常数据处理，可以免除列误差方程和法方程再求解的过程，大大提高了计算效率。同时Matlab自带的griddata插值函数运用在GPS高程拟合中完全能代替传统的模型拟合方法，具有程序编制简便、转换精度可靠的特点，Matlab的图像处理功能为我们提供清晰可靠的立体视觉效果。

[1] 李征航，黄劲松. GPS 测量原理与数据处理[M]. 武汉:武汉大学出版社，2005

[2] 王旭，刘文生. GPS高程拟合方法的研究[J]. 测绘科学，2010(35):28-31

[3] 李秀海，韩冰.基于多面函数模型的GPS高程拟合精度分析[J]. 测绘与空间地理信息，2010(1):12-14

[4] 邱斌，朱建军. 基于模糊集合理论的GPS高程多面函数拟合模型平滑因子优选[J]. 测绘工程，2004,13(3):35-39

[5] 求是科技. Matlab7.0从入门到精通[M]. 北京：人民邮电出版社，2006

[6] 白征东.Matlab在测量平差教学中的应用[J].测绘通报，2009(11):73-76

[7] 石博强，赵金. Matlab数学计算与工程范例教程[M]. 北京:中国铁道出版社， 2005

[8] 陈本富，王贵武. 基于Matlab的数据处理方法在GPS高程拟合中的应用[J]. 昆明理工大学学报，2009(34):1-4

[9] 谢劭峰,王新桥. Matlab在GPS高程转换中的应用[J]. 测绘与空间地理信息，2005,28(6): 75-76

[10] 杨庆振，郭春喜.多种拟合方法在似大地水准面精化中的比较[J].测绘标准化，2009,25 (1):30-32

[11] Magrab E B，高会生. Matlab原理与工程应用[M]. 北京：电子工业出版社，2002

P208

B

1672-4623（2015）01-0099-03

10.3969/j.issn.1672-4623.2015.01.033

骆丽华，硕士，研究方向为测绘信息采集与数据处理。

2014-01-15。

项目来源：国家自然科学基金资助项目（41071294）；广西自然科学基金资助项目（0991023）。