駱麗華，覃 輝

（1.桂林理工大學 廣西礦冶與環(huán)境科學實驗中心，廣西 桂林 541004；2. 廣西空間信息與測繪重點實驗室，廣西 桂林 541004；3.桂林理工大學 測繪地理信息學院，廣西 桂林 541004；4.廣東科學技術職業(yè)學院 , 廣東 珠海 519090）

Matlab程序設計在GPS高程擬合中的應用

駱麗華1,2,3，覃 輝4

（1.桂林理工大學 廣西礦冶與環(huán)境科學實驗中心，廣西 桂林 541004；2. 廣西空間信息與測繪重點實驗室，廣西 桂林 541004；3.桂林理工大學 測繪地理信息學院，廣西 桂林 541004；4.廣東科學技術職業(yè)學院 , 廣東 珠海 519090）

介紹了二次曲面擬合法、多面函數(shù)擬合法以及Matlab自帶的griddata擬合插值函數(shù)法。計算結果表明，Matlab強大的數(shù)據(jù)處理能力運用在GPS高程擬合中，處理方式簡潔，計算量少，精度可靠，同時具有立體可視特征。

Matlab；GPS高程擬合；二次曲面；多面函數(shù)；griddata插值

GPS坐標是以地球質心為坐標原點的地心坐標，測量獲得的高程是以WGS84橢球面為基準面的大地高[1]。而在我國工程中采用的高程通常是以似大地水準面為基準面的正常高，大地高與正常高之差稱為高程異常ζ。因此，在實際測量工作中，大地高和正常高之間的轉換就變得十分有意義，只要求出高程異常就可以方便地求得正常高。通常求取高程異常的方法是通過某種數(shù)學模型根據(jù)已知的GPS三維點位模擬求出未知點的高程異常。由于Matlab突破了高級編程語言要編寫大量循環(huán)語句的思想，它只需通過一個或幾個簡單的命令就可以完成矩陣或數(shù)組的運算，所以文中結合實例利用Matlab進行程序設計，對局部區(qū)域的高程異常進行擬合，實現(xiàn)了數(shù)據(jù)的自動處理。

1 GPS高程擬合原理

1.1 曲面擬合法

曲面擬合法，就是將高程異常近似地看作是一定范圍內各點坐標的曲面函數(shù)，認為高程異常在此范圍內變化平緩，可以用已聯(lián)測水準的GPS點的高程異常來擬合這一函數(shù)。在求得函數(shù)的擬合系數(shù)之后，就可用這一擬合函數(shù)來計算其他GPS點的高程異常和正常高[2]。設某公共點的高程異常ζ與其平面坐標x、y之間有以下關系：

式中，ζ為高程異常；f (x，y)為ζ的趨勢值；ε為誤差值。

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

寫成矩陣形式有：

誤差方程為：

對于每一個已知點，均可以列出上式方程，在∑ε2=min條件下，可求解系數(shù)陣：

再由己知高程異常的權陣情況，上式可寫為：

系數(shù)求出后，再按式（2）求出待求點的正常高。

1.2 多面函數(shù)曲面擬合法

多面函數(shù)法的基本思想是對于任何數(shù)學表面或任何不規(guī)則的圓滑表面，總可以用一系列有規(guī)則的數(shù)學表面的總和，以任意精度逼近[3]。根據(jù)這樣的思想，高程異常函數(shù)可以表示為：

式中，ai為待定系數(shù)；Q（x，y，xi，yi）為核函數(shù)。核函數(shù)是關于x、y的函數(shù)，核函數(shù)的中心在（xi，yi）處。理論上核函數(shù)是可以任意構造的，在實際應用中，通常用以下幾種函數(shù)來充當核函數(shù)。

錐面：

雙曲面：

倒曲面：

在上述各式中，x、y表示內插點坐標；xi、yi表示已知點的坐標；核函數(shù)中的表示內插點到已知點的水平距離；δ為光滑系數(shù)。

以核函數(shù)為雙曲面為例，說明多項式曲面擬合法的具體求解過程。設測區(qū)內已知點個數(shù)為n，求解式（7）中的系數(shù)（a1，a2，a3，…，an），其矩陣形式為：

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

求解未知點的高程異常值，根據(jù)式（8）和式（9）即可得到求解公式：

在選擇多面函數(shù)求解高程異常值的時候，需要注意核函數(shù)的選取問題。由于是自主取值，為了達到最佳擬合效果，就要逐步試驗進而改進，然后選定一個最佳取值。同時，通過改變平滑因子[4]，可以使計算結果達到最理想的狀態(tài)。因此，沒有必要去關心平滑因子值和數(shù)據(jù)的具體聯(lián)系，只需要找到一個合適的值使得計算結果處于最佳狀態(tài)。

下文實驗中確定平滑因子的具體方法是：在程序中將平滑因子設計為變量，初始值設為0，步長為0.000 01（用戶可以自定義步長）。當所求的外符合精度不在要求精度范圍內，程序會繼續(xù)調節(jié)平滑因子。每次調節(jié)平滑因子后，系統(tǒng)進行一次運算，并求得該次運算后的外符合精度，如此循環(huán)運算，直到外符合精度達到要求。實驗中平滑因子確定為6 498.860 11。

首先,對60例患者的病因構成及一般資料統(tǒng)計分析中顯示,60例非炎癥性腸病大腸潰瘍患者中,腸結核患者16例,感染性腸炎患者10例,吻合口發(fā)生潰瘍患者4例,潰瘍性大腸癌患者19例,缺血性腸炎患者5例,放射性腸炎患者2例,內痔術后潰瘍患者2例,急性闌尾炎患者2例;其中,以潰瘍性大腸癌以及腸結核患者數(shù)量最多,其次分別為感染性腸炎、缺血性腸炎以及吻合口潰瘍等情況患者,這也是導致非炎癥性腸病大腸潰瘍的主要病因。

1.3 Matlab的griddata插值函數(shù)法

Matlab中的griddata函數(shù)可以將位于同一空間坐標系下的散點插值為規(guī)則格網(wǎng)，提供了以下4種方法：線性內插法、三次多項式內插法、最鄰近點內插法和Matlab自帶的一種圓滑插值法（V4）,可方便地實現(xiàn)結合鄰近離散點分布特征的光滑曲面擬合[5]。

2 Matlab程序實現(xiàn)

以下程序分別是二次曲面擬合法、多面函數(shù)擬合法以及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);%數(shù)據(jù)中心化

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]);%模型參數(shù)

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 多面函數(shù)擬合法

%多面函數(shù)擬合合法（文件名：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;%模型參數(shù)

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)); %生成規(guī)格化格網(wǎng)

[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 實驗數(shù)據(jù)來源

本文數(shù)據(jù)源于文獻[9]某個沿江且地勢平坦地區(qū)，面積約為8 km2。該GPS控制網(wǎng)的17個水準高程點都是等精度且不含粗差，點位分布如圖1所示。本實驗將這17個點分為2組，前10個點（藍色）作為實驗的模型擬合點，后7個點（紅色）數(shù)據(jù)作為檢核數(shù)據(jù)點。利用Matlab分別用二次曲面擬合法、多面函數(shù)擬合法和griddata插值法編程實現(xiàn)，各種方法計算的結果如表1～4所示。

3.2 數(shù)據(jù)精度分析

圖1 點位分布圖

表1 二次曲面擬合結果

表2 多面函數(shù)擬合結果

表3 Matlab三次多項式內插（cubic）擬合結果

表4 Matlab’V4’法擬合結果

由于高程擬合大多采用數(shù)學模型進行插值計算，對于模型誤差很難確定，所以GPS水準計算的精度一般利用內、外符合精度進行評定。Griddata函數(shù)進行高程擬合，由于擬合曲面通過各擬合點，故擬合點的內符合精度為0。本文為了便于與二次曲面及多面函數(shù)擬合法比較，取插值點的單點最大偏差及外部符合精度為評價指標。

根據(jù)核檢點ζi與擬合值ζi'之差，按下式可計算出外符合精度m[10]：

式中，v為擬合殘差；n 為檢核點數(shù)。經(jīng)過計算，各方法的精度指標如表5所示。

表5 幾種擬合方法的精度對比

由表5可以看出，二次曲面擬合檢核點最大殘差為3.12 mm，多面函數(shù)擬合檢核點最大殘差為9.09 mm，Matlab三次多項式內插法最大殘差為6.13 mm，Matlab自帶的V4插值法最大殘差為3.83 mm。二次曲面擬合效果最好，能達到三等水準測量的精度要求；多面函數(shù)法擬合效果最差，其原因可能是核函數(shù)選取和平滑因子取值不是最理想狀態(tài)，但能達到四等水準測量精度；同樣Griddata插值函數(shù)擬合法也能達到四等水準測量的精度要求。由于似大地水準面是一個不規(guī)則的連續(xù)曲面，二次曲面擬合的高程異常成果較好。從表5也反映出Matlab的griddata函數(shù)中V4法相對三次多項式內插擬合法精度也有一定提高。根據(jù)外符合精度值可以看出，用Matlab的擬合插值函數(shù)進行GPS高程擬合也具有較好的外推能力，完全可以滿足大多數(shù)工程的要求。 Matlab的圖像處理功能[11]可以清楚地顯示擬合區(qū)域的高程異常分布。圖2是以Matlab的3次多項式擬合的高程異常空間分布圖。

圖 2 高程異常擬合曲面示意圖

4 結 語

運用Matlab進行高程異常數(shù)據(jù)處理，可以免除列誤差方程和法方程再求解的過程，大大提高了計算效率。同時Matlab自帶的griddata插值函數(shù)運用在GPS高程擬合中完全能代替?zhèn)鹘y(tǒng)的模型擬合方法，具有程序編制簡便、轉換精度可靠的特點，Matlab的圖像處理功能為我們提供清晰可靠的立體視覺效果。

[1] 李征航，黃勁松. GPS 測量原理與數(shù)據(jù)處理[M]. 武漢:武漢大學出版社，2005

[2] 王旭，劉文生. GPS高程擬合方法的研究[J]. 測繪科學，2010(35):28-31

[3] 李秀海，韓冰.基于多面函數(shù)模型的GPS高程擬合精度分析[J]. 測繪與空間地理信息，2010(1):12-14

[4] 邱斌，朱建軍. 基于模糊集合理論的GPS高程多面函數(shù)擬合模型平滑因子優(yōu)選[J]. 測繪工程，2004,13(3):35-39

[5] 求是科技. Matlab7.0從入門到精通[M]. 北京：人民郵電出版社，2006

[6] 白征東.Matlab在測量平差教學中的應用[J].測繪通報，2009(11):73-76

[7] 石博強，趙金. Matlab數(shù)學計算與工程范例教程[M]. 北京:中國鐵道出版社， 2005

[8] 陳本富，王貴武. 基于Matlab的數(shù)據(jù)處理方法在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]. 北京：電子工業(yè)出版社，2002

P208

B

1672-4623（2015）01-0099-03

10.3969/j.issn.1672-4623.2015.01.033

駱麗華，碩士，研究方向為測繪信息采集與數(shù)據(jù)處理。

2014-01-15。

項目來源：國家自然科學基金資助項目（41071294）；廣西自然科學基金資助項目（0991023）。