常宜峰 柴洪洲 王 敏

(解放軍信息工程大學(xué)測(cè)繪學(xué)院，鄭州 450052)

基于正交最小二乘多面函數(shù)構(gòu)建局域海洋地磁場(chǎng)模型研究*

常宜峰 柴洪洲 王 敏

(解放軍信息工程大學(xué)測(cè)繪學(xué)院，鄭州 450052)

針對(duì)多面函數(shù)擬合局域海洋地磁場(chǎng)模型中質(zhì)量控制和節(jié)點(diǎn)選擇問題，提出使用抗差估計(jì)與正交最小二乘原理相結(jié)合的方法。在質(zhì)量控制過程中，將局部有顯著起伏的采樣點(diǎn)視為異常觀測(cè)值，采用抗差估計(jì)減弱異常值點(diǎn)的不利影響。在節(jié)點(diǎn)選擇中，使用基于正交最小二乘原則判斷節(jié)點(diǎn)的貢獻(xiàn)效力大小，按照此原則自適應(yīng)地選擇節(jié)點(diǎn)。結(jié)果表明，基于正交最小二乘原理的多面函數(shù)法的精度和穩(wěn)定性更好，更適于構(gòu)建局域海洋地磁場(chǎng)模型。

局域海洋地磁場(chǎng)模型;抗差估計(jì);正交最小二乘原理;多面函數(shù);節(jié)點(diǎn)

1 引言

地磁場(chǎng)在了解地球內(nèi)部的物質(zhì)分布和輔助導(dǎo)航中有重要用途，其中，地磁場(chǎng)模型起著非常重要的作用。常用的是泰勒多項(xiàng)式模型，優(yōu)點(diǎn)是計(jì)算簡(jiǎn)單，但地磁數(shù)據(jù)存在較大跳變性，難以取得較高的精度。1971年Hardy給出了多面函數(shù)擬合法，1976年用于重力異常和大地水準(zhǔn)面差距計(jì)算，1978年用于地殼形變計(jì)算，2009年趙建虎將此方法引入局域海洋地磁場(chǎng)建模中，并提出分區(qū)建模來提高精度［1］。多面函數(shù)是利用一些規(guī)則的數(shù)學(xué)曲面來逼近實(shí)際曲面，從而獲得較高的精度。本文提出將抗差估計(jì)與正交最小二乘相結(jié)合來擬合局域海洋地磁場(chǎng)模型。

對(duì)于多面函數(shù)模型，核函數(shù)的選擇和節(jié)點(diǎn)的確定至關(guān)重要。對(duì)于核函數(shù)，已有學(xué)者進(jìn)行過深入研究。正雙曲面函數(shù)構(gòu)建的模型在檢查點(diǎn)上的誤差分布和精度統(tǒng)計(jì)參數(shù)均優(yōu)于倒雙曲面函數(shù)，并且正雙曲面核函數(shù)構(gòu)造的模型精度隨平滑因子的變化平緩而相對(duì)穩(wěn)定，倒雙曲面函數(shù)構(gòu)造的模型精度變化劇烈且最優(yōu)平滑因子相對(duì)較大，而過大的平滑因子可能造成局部模型失真，尤其對(duì)于地磁變化相對(duì)較大的數(shù)據(jù)的擬合，在此取正雙曲面函數(shù)，平滑因子取0.1［1］。對(duì)于節(jié)點(diǎn)的選擇，通常是直接取采樣點(diǎn)，可是選擇節(jié)點(diǎn)不同產(chǎn)生的結(jié)果也不同，具有一定的隨意性。由于每個(gè)節(jié)點(diǎn)對(duì)擬合結(jié)果貢獻(xiàn)是不同的，有學(xué)者采用t檢驗(yàn)或逐步回歸的方法來選擇節(jié)點(diǎn)，得到的效果較好［2］。本文利用正交最小二乘原理的基本思想，首先將誤差方程系數(shù)正交化，并計(jì)算觀測(cè)向量與正交化向量之間的夾角，節(jié)點(diǎn)對(duì)曲面擬合的貢獻(xiàn)就用夾角的大小來表示。按照這種原則方法依次選擇合適的節(jié)點(diǎn)，完成自主選擇節(jié)點(diǎn)的目的［3］。此方法的優(yōu)點(diǎn)是既保證了擬合的精度又兼顧了擬合的穩(wěn)定性。

2 基于抗差估計(jì)的質(zhì)量控制

在測(cè)量過程中受到各種因素的影響，可能出現(xiàn)異常觀測(cè)數(shù)據(jù)，一般的質(zhì)量控制是采用判斷測(cè)量值與對(duì)應(yīng)趨勢(shì)面擬合值之間差值大小的趨勢(shì)面法來進(jìn)行的。趨勢(shì)面分析方法是將短尺度的、局部的變化看作隨機(jī)的和非結(jié)構(gòu)的噪聲，因此當(dāng)局部區(qū)域變化比較劇烈時(shí)，可將其視作異常觀測(cè)。但是單純的趨勢(shì)面擬合并不具備抗差的能力，故本文引入抗差估計(jì)理論［4，5］，采用抗差趨勢(shì)面法進(jìn)行質(zhì)量控制。

趨勢(shì)面函數(shù)如下式所示

其中，αij(i=0，…，n，j=0，…，n)為待求系數(shù)，(x，y)為觀測(cè)點(diǎn)坐標(biāo)為對(duì)應(yīng)的磁力觀測(cè)值。如果有m個(gè)觀測(cè)點(diǎn)，應(yīng)用最小二乘準(zhǔn)則即可解得多項(xiàng)式系數(shù)。

設(shè)有誤差方程

根據(jù)抗差估計(jì)理論，以及等價(jià)權(quán)理論可得

3 基于正交最小二乘法節(jié)點(diǎn)選擇

經(jīng)過質(zhì)量控制后，再對(duì)剩下的殘差進(jìn)行擬合。根據(jù)多面函數(shù)原理有，

在已知有限測(cè)量數(shù)據(jù)的情況下求解方程，重要的是節(jié)點(diǎn)的選擇，一般情況下直接取均勻分布在研究區(qū)域內(nèi)的觀測(cè)點(diǎn)，如果選擇所有的采樣點(diǎn)作為節(jié)點(diǎn)，可能會(huì)出現(xiàn)內(nèi)符合精度虛高，推估和內(nèi)插精度無法保證的情況。特別在采樣點(diǎn)受異常誤差污染時(shí)，擬合精度會(huì)更低。為解決以上問題，本文在此引入基于正交最小二乘原理的思想來建立一種自主選擇節(jié)點(diǎn)的方法。

在m個(gè)磁力數(shù)據(jù)中選擇n個(gè)數(shù)據(jù)建模，即

其中，X為模型系數(shù)，A為核函數(shù)矩陣，由最小二乘原理得，

模型精度評(píng)價(jià)指標(biāo)為［6］

其中，n為節(jié)點(diǎn)個(gè)數(shù)，(xi，yi)為節(jié)點(diǎn)坐標(biāo)，βi為模型系數(shù)，Q(x，y，xi，yi)為二次核函數(shù)，一般取正雙曲函數(shù)(如下式)，δ2為模型的平滑因子(本文取0.1)。

先假設(shè)選取所有的采樣點(diǎn)為節(jié)點(diǎn)，即m=n，此時(shí)構(gòu)建的誤差方程系數(shù)陣為

其中，

從式(11)可以看出，每個(gè)ai對(duì)求系數(shù)β的決定效力是有差別的，即每一節(jié)點(diǎn)對(duì)函數(shù)擬合的影響是不同的，故提供了一種評(píng)價(jià)準(zhǔn)則，便于依照節(jié)點(diǎn)貢獻(xiàn)效力的大小來選擇節(jié)點(diǎn)。

下面就引入一種表示節(jié)點(diǎn)貢獻(xiàn)大小的原則，將式(5)表示成向量形式，

假設(shè)ak之間是相互正交的，則左乘βk有，

兩邊求和并移項(xiàng)得

定義η為選擇t(t≤m)個(gè)節(jié)點(diǎn)的列向量時(shí)的總貢獻(xiàn)

η的值在0到1之間。并且當(dāng)t越大，η就越大，表明逼近程度越高;當(dāng)t=m時(shí)，η=1，此時(shí)選擇了所有的采樣點(diǎn)，擬合誤差為零，逼近精度最高。但由于沒有多余觀測(cè)，可靠性較差，為此應(yīng)根據(jù)每個(gè)基矢量的貢獻(xiàn)大小來選擇若干節(jié)點(diǎn)，使η足夠大且有多余觀測(cè)。

要想使η最大，即需使βk取最大值，在式(12)兩邊同時(shí)左乘，可得

顯然，要使βk取得最大值，需向量ak與觀測(cè)向量Z之間的夾角θakZ最大，即

其具體計(jì)算步驟為:

1)將所有的采樣點(diǎn)視為剩余節(jié)點(diǎn)，計(jì)算系數(shù)陣A，m個(gè)列向量可記為a1，a2，…，am;

2)計(jì)算Z與ak之間的夾角，求得最大值

式(18)表明aj對(duì)Z的貢獻(xiàn)最大，于是可將與aj相對(duì)應(yīng)的坐標(biāo)(xj，yj)選為節(jié)點(diǎn)。

3)按照最小二乘原理解算系數(shù)β

4 實(shí)驗(yàn)分析

為驗(yàn)證以上所述方法在局域海洋地磁場(chǎng)模型建立中的有效性，采用某測(cè)區(qū)海洋磁力測(cè)量的數(shù)據(jù)進(jìn)行計(jì)算和分析。該測(cè)區(qū)定位系統(tǒng)采用NGD-60型差分GPS定位系統(tǒng)，船載磁力儀采用美國(guó)Geometrics公司的G-881型海洋磁力儀，共有1 681個(gè)測(cè)量數(shù)據(jù)。經(jīng)過了各項(xiàng)的改正，測(cè)區(qū)的等值線如圖1所示。均勻選取其中的441個(gè)數(shù)據(jù)作為已知的測(cè)量點(diǎn)，將剩余的測(cè)量數(shù)據(jù)作為外部的檢核點(diǎn)，分別使用多面函數(shù)法、正交最小二乘多面函數(shù)法、基于抗差趨勢(shì)面的多面函數(shù)法和基于抗差趨勢(shì)面的正交最小二乘多面函數(shù)法，4種方法進(jìn)行計(jì)算(表1，其中最大值和最小值是指模型計(jì)算結(jié)果與實(shí)測(cè)值的偏離的最大和最小值，平均值是所有偏離量的算數(shù)平均值，均方誤差是對(duì)外部檢核點(diǎn)的綜合評(píng)價(jià)指標(biāo))。

圖1 某測(cè)區(qū)磁力異常分布Fig.1 Distribution of magnetic abnormity in the studied area

表1 不同擬合方法的結(jié)果比較(單位：nT)Tab.1 Comparison among results with different methods (unit：nT)

從表1可以看出，直接使用多面函數(shù)和正交多面函數(shù)進(jìn)行擬合，均方誤差分別可達(dá)17nT和15nT，平均值都在0.6 nT左右，說明模型的集中程度較好，兩者所得到的精度基本相當(dāng)，但正交多面函數(shù)法略好于多面函數(shù)法，并且都有所偏大。加上抗差估計(jì)趨勢(shì)面進(jìn)行質(zhì)量控制之后，同樣比較兩種方法的計(jì)算結(jié)果，均方誤差分別為12 nT和7 nT，精度提高了41%，平均值分別為-8.4 nT和0.2 nT，基于抗差趨勢(shì)面的多面函數(shù)出現(xiàn)了明顯的系統(tǒng)偏差，分析原因是多面函數(shù)法對(duì)節(jié)點(diǎn)沒有選擇性，不是選擇一些特征節(jié)點(diǎn)進(jìn)行建模，而是將所有的節(jié)點(diǎn)作為建模節(jié)點(diǎn)，從而引入了一些系統(tǒng)誤差，進(jìn)而產(chǎn)生了系統(tǒng)性偏差。比較多面函數(shù)和基于抗差趨勢(shì)面的多面函數(shù)，增加了抗差以后，精度由17 nT變?yōu)?2 nT;同樣，對(duì)正交多面函數(shù)增加了抗差以后，精度由15 nT變?yōu)? nT，說明基于抗差趨勢(shì)面的質(zhì)量控制可以減小一部分不良影響。綜合表明基于抗差趨勢(shì)面的正交多面函數(shù)法總體上要優(yōu)于以上其他3種方法。

圖2和圖3分別是使用基于抗差趨勢(shì)面的正交多面函數(shù)法和基于抗差趨勢(shì)面的多面函數(shù)法擬合的檢驗(yàn)點(diǎn)的殘差分布圖，從圖中可以看出，圖2的殘差分布相對(duì)更為集中，且普遍都靠近零點(diǎn)。

為進(jìn)一步驗(yàn)證按照正交多面函數(shù)法選擇節(jié)點(diǎn)的多少與所得模型擬合精度之間的關(guān)系，依次計(jì)算了選擇個(gè)檢驗(yàn)點(diǎn)時(shí)的均方誤差，并將用正交多面函數(shù)選擇的節(jié)點(diǎn)個(gè)數(shù)與擬合的均方誤差的關(guān)系繪制成圖4。

圖2 基于抗差趨勢(shì)面的正交多面函數(shù)法擬合殘差Fig.2 Fitting residual with orthogonal multi-function method based on robust trend-surface

從圖4可以看出，隨著選擇的節(jié)點(diǎn)個(gè)數(shù)的增加，均方誤差逐漸減少。在采用250個(gè)節(jié)點(diǎn)時(shí)，精度已經(jīng)達(dá)到了17 nT，基本上相當(dāng)于選擇所有節(jié)點(diǎn)進(jìn)行多面函數(shù)法建模時(shí)的精度。隨著選擇節(jié)點(diǎn)個(gè)數(shù)的繼續(xù)增加，精度基本保持穩(wěn)定并有所提高。同時(shí)應(yīng)該指出，針對(duì)地磁變化比較復(fù)雜和比較平緩的兩種區(qū)域，數(shù)據(jù)的代表性對(duì)于模型的精度有著重要的影響。因此，在變化比較復(fù)雜的區(qū)域可以適當(dāng)?shù)姆艑挼Ｖ沟呐袛鄺l件，盡可能多的自主選擇特征點(diǎn)進(jìn)行建模。

圖3 基于抗差趨勢(shì)面的多面函數(shù)法擬合殘差Fig.3 Fitting residual with multi-function method based on robust trend-surface

圖4 節(jié)點(diǎn)選擇與擬合誤差關(guān)系Fig.4 Relation between node collecting and error fitting

圖5 兩種方法擬合效果比較Fig.5 Comparion between fitting results with two methods

為了使計(jì)算結(jié)果更直觀，用正交多面函數(shù)法和一般多面函數(shù)法分別計(jì)算了整個(gè)區(qū)域的網(wǎng)格點(diǎn)地磁場(chǎng)總強(qiáng)度值，并利用計(jì)算所得數(shù)據(jù)分別繪制正交多面函數(shù)法和普通多面函數(shù)法磁力異常圖，與該區(qū)域?qū)嶋H磁力異常場(chǎng)圖進(jìn)行比較。圖5(a)是原始的磁力異常分布圖，(b)是使用基于抗差趨勢(shì)面的正交多面函數(shù)法得到的擬合圖，(c)是使用基于抗差趨勢(shì)面的多面函數(shù)法得到的結(jié)果。從圖中可以看出，圖5(b)與(a)等高線的分布較為接近，顯現(xiàn)出了較為細(xì)節(jié)的分布信息，保持了較高的一致性;(c)與(a)的一致性較差。從而顯示圖5(b)的建模方法是穩(wěn)定可靠的。

5 結(jié)論

1)進(jìn)一步驗(yàn)證了多面函數(shù)法比較適用于局域海洋地磁場(chǎng)建模，在地磁變化復(fù)雜的海域，先進(jìn)行趨勢(shì)面擬合，再對(duì)剩余的誤差進(jìn)行多面函數(shù)擬合，可以突出局部變化，擬合效果要優(yōu)于直接使用多面函數(shù)進(jìn)行擬合。

2)趨勢(shì)面擬合主要反映整體變化趨勢(shì)，當(dāng)局部磁力異常觀測(cè)值變化顯著時(shí)，采用抗差趨勢(shì)面法可以減小對(duì)趨勢(shì)面擬合的不良影響。在多面函數(shù)擬合前，采用抗差估計(jì)理論進(jìn)行質(zhì)量控制，可以更有效剔除異常觀測(cè)值對(duì)趨勢(shì)面擬合的影響。

3)使用正交最小二乘多面函數(shù)法得到的擬合精度優(yōu)于多面函數(shù)法，可以按照規(guī)則自適應(yīng)地選擇特征點(diǎn)，提高了人工選擇特征點(diǎn)的效率。對(duì)系數(shù)矩陣進(jìn)行正交化，不會(huì)影響最終的結(jié)果，因?yàn)檎换蟛]有參與多面函數(shù)系數(shù)的求解，僅是用于節(jié)點(diǎn)的選擇，因此可以作為一種自動(dòng)選擇節(jié)點(diǎn)的參考方法。

1 趙建虎，等.基于多面函數(shù)的局域地磁場(chǎng)建模方法研究［J］.海洋通報(bào)，2009，28(4):89-97.(Zhao Jianhu，et al.Study on construction of local marine geomagnetic field model based on multi-surface function［J］.Marine Science Bulletin，2009，28(4):89-97)

2 劉兆平，楊進(jìn)，武煒.地球物理數(shù)據(jù)網(wǎng)格化方法的選取［J］.物探與化探，2010，34(1):93-98.(Liu Zhaoping，Yang Jin and Wu Wei.The choice of gridding methods for geophysical data［J］.Geophysicalamp;Eochemical Exploration，2010，34(1):93-98)

3 張菊清，劉平芝.抗差趨勢(shì)面與正交多面函數(shù)結(jié)合擬合DEM數(shù)據(jù)［J］.測(cè)繪學(xué)報(bào)，2008，37(4):526-531.(Zhang Juqing and Liu Pingzhi.Combining fitting based on robust trend surface and orthogonal multiquadrics with application in DEM fitting［J］.Acta Geodaetica et Cartographica Sinica，2008，37(4):526-531)

4 楊元喜.抗差估計(jì)理論及其應(yīng)用［M］.北京:八一出版社，1993.(Yang Yuanxi.The theory and application of robust estimation［M］.Beijing:Bayi Publishing House，1993)

5 楊元喜.自適應(yīng)抗差最小二乘估計(jì)［J］.測(cè)繪學(xué)報(bào)，1996，25(3):206-211.(Yang Yuanxi.A daptively robust least squares estimation［J］.Acta Geodaetica et Cartographica Sinica，1996，25(3):206-211)

6 黃學(xué)功，等.地磁圖制備方法及其有效性評(píng)估［J］.北京航空航天大學(xué)學(xué)報(bào)，2009，35(7):891-895.(Huang Xuegong，et al.Geomagnetic mapping and validity estimation［J］.Journal of Beijing University of Aeronautics And Stronautics，2009，35(7):891-895)

STUDY ON CONSTRUCTION OF LOCAL MARINE GEOMAGNETIC FIELD MODEL BASED ON ORTHOGONAL LEAST SQUARE MULTI-SURFACE FUNCTION

Chang Yifeng，Chai Hongzhou and Wang Min
(Institute of Surveying and Mapping，Information Engineering University，Zhengzhou 450052)

Aim at the problem of quality control and node choice in modeling local marine geomagnetic field based on multi-surface function，a fitting method with combining robust trend surface and orthogonal least squares is proposed.In order to control the influence of some outstanding points to surface fitting，a robust fitting of the trend surface by using an equivalent weight is adopted and the adaptive node choosing method is proposed based on the effect of every node on the curve fitting calculated by using the orthogonal least squares as well.The results show that this method is more effective in accuracy and stability than normal method.

local marine geomagnetic field;robust estimate;orthogonal least square;multi-surface function;node

1671-5942(2012)03-0086-05

2011-11-09

常宜峰，男，1986年生，碩士研究生，主要從事測(cè)量數(shù)據(jù)處理理論與方法研究.E-mail:changyifeng1986@163.com

P207

A