帶權最小二乘多項式擬合的泛克立格法

帶權最小二乘多項式擬合的泛克立格法

克立格法是一種基于隨機過程的統(tǒng)計預測方法。理論上，克立格估計具有最小估計散射差異的統(tǒng)計特征。此外，克立格應用中最重要的是半變異函數的確定，即協方差函數的確定。然而，由于觀測值較低，在實踐中難以應用更精確的協方差函數?？肆⒏窆浪闶欠窬哂凶钚』烙嫴町惖慕y(tǒng)計特征？從統(tǒng)計學的角度來看，在某些情況下，它是最小差分的最優(yōu)計算公式，但從數學的角度來看，它只是一種相對特殊的近似方法，即通過旋轉曲線生成的基函數。這已經被許多科學家證明。在這項工作中，我們將重點關注克立格估算的分析和解釋。1系數矩陣估計的可恢復算法設以節(jié)點(xi,yi)為中心生成的旋轉面為以此為基礎可建立其線性組合再設數據點為zi=z(xi,yi),i=1,2…,n,在數據點上可建立如下插值條件寫成矩陣形式為式中:適當地選擇節(jié)點和基函數可保證系數矩陣非奇異,因此由(3)式可得將(4)式代入(2)式可建立未觀測點(xp,yp)的插值函數為如將基函數取為各向同性的協方差函數則有它形式上與已知數學期望的簡單克立格估計相同.2克立格估計設z(x)是具有常數數學期望的平穩(wěn)隨機過程,我們的目的是通過觀測值zi=z(xi),i=1,2,…,n的加權線性組合最優(yōu)地確定z(xp),式中wi為待定的克立格權,利用條件極值可推得確定系數wi的方程組為式中:μ為拉格朗日乘數,引入矩陣記號(8)式可改寫為權系數與拉格朗日乘數可分別求解.由(9c)式有因為ITC-1I只是一個數,故有(13)式再回頭代入(11)式可得再將(14)式代入(8)式可得記M可理解為觀測值的廣義平均值,當觀測值獨立C=E或觀測值弱相關C≈E,則有即為簡單平均值.在(15)式兩端同時減去M可得(18)式與(7)式相對比可知帶一個未知均值的克立格估計可解釋為觀測值與觀測值的廣義均值之差被旋轉曲面生成的基函數所逼近.3拉格朗日乘數估計設平穩(wěn)隨機過程z(x)=m(x)+s(x)由二部分組成,第一部分m(x)為確定性漂移,第二部分s(x)為平穩(wěn)隨機信號.一般漂移可表示為多項式形式式中:fT(x)=[f0(x),f1(x),…,fk(x)]=[1,x,y,xy…],a=(a0,a1,…,ak)T.令zi,i=1,2,…,n為觀測值,z(xp)為待估值,并可由觀測值的加權線性組合求出由克立格估計準則應有限于篇幅,略去推導過程可得求解克立格權與拉格朗日乘數的方程組為式中:μ=(μ0,μ1,…,μk)T,fT=[f0(x),f1(x),…,fk(x)]T,F為由漂移部分組成的系數陣.C非奇異,由(22)式可得(23)式代入(22)式,易知然后將(24)式再次回代至(23)式,可知克立格權系數矢量為(25)式代入(20)式有記它就是(19)式的系統(tǒng)參數,因此(26)式可改寫為或(29)式寫成元素形式為z(xp)-fpTa=(30)式仍與(7)式類似,因此可以說在除去系統(tǒng)部分后,泛克立格與使用旋轉面作為基函數的函數逼近類似.此外與(22)式直接解泛克立格方程組相比,本文(27),(28)式先求系統(tǒng)漂移部分再求隨機信號部分,既降低了求逆矩陣階數,又提高了計算效率,并且也易于區(qū)分漂移部分與隨機部分的大小量級.4協方差參數值本節(jié)我們基于一定的數學分析觀點來分析協方差函數參數的選取.為討論方便以一維為例,(5)式略去下標p可改寫為(31)式表明預測值只是函數k(x)移位后的線性組合.取k(x)為最常用的Gauss函數Gauss函數衰減很快,它的局部支持區(qū)間受到了參數d的強烈控制.不同參數的Gauss函數可示于圖1,從圖1可看出,當d較小時,Gauss函數類似于不連續(xù)的脈沖函數δ(x).作為理想化的一個特例,我們考查ai≡1,xi=i時不同參數的Gauss函數之線性組合構成的函數形態(tài).d=1.0,d=0.5和d=0.1時Gauss函數移位后的線性組合可示于圖2.實際上ai≡1隱含著采樣值(或觀測值)為一水平直線的采樣,理想的逼近結果也應為一條直線.但事實并非如此,圖2c中的曲線類似于一把梳子,這是由于d=0.1時Gauss函數的支撐區(qū)間太窄所致.圖2b中的曲線已較圖2a好得多,但仍有許多連續(xù)的極大值和極小值,類似于一條正余弦曲線.只有圖2a達到了較理想的效果.這些事實說明協方差函數的逼近效果的確與其參數有很大關系,并且可以肯定地說合理的參數是存在的.首先我們假定信噪比或噪聲部分已被確定或已另外考慮,此處我們只考慮連續(xù)的信號部分.準則1協方差參數應使得證明:(1)d=min意味著Gauss函數具有較小的支撐區(qū)間,由此形成的矩陣帶寬將比較小.(2)在假定z(x)平穩(wěn)的情況下,半節(jié)點的預測值與節(jié)點上的預測值應有相同的數學期望,即在E{z(x)}=C(C為常數)或limC※0)的情況下,由于ai為觀測值的線性組合,有故(35)式代入(34)式有(36)式稍加變形也即(33)式,此處(33)式亦可從代數上理解為一條水平直線被核函數移位生成的基函數所逼近時,其半節(jié)點處的值應與節(jié)點處的值相等.證畢.為求出合理的協方差參數值,引入比值函數式中考慮到Gauss函數衰減很快,對區(qū)間[-4,4]以外的Gauss函數已視其為零略去不計.為更好地了解比值函數隨d的變化規(guī)律,圖3畫出了它的變化曲線.從圖3可得出兩條主要的結論:(1)合理的參數應位于區(qū)間[0.8,1.7],它可以使半節(jié)點與節(jié)點處的協方差函數線性組合大致相等;(2)d<0.5的參數應盡量避免使用,因為它使半節(jié)點處的協方差函數的線性組合遠小于節(jié)點處協方差函數的線性組合.準則2Gauss函數的參數可比照B樣條函數確定.說明:B樣條函數在函數擬合或逼近中具有許多優(yōu)點因而獲得了非常廣泛的應用.n次B樣條在數學上可方便地由矩形脈沖經n-1次卷積得到式中:*表示卷積,π(x)是矩形脈沖函數特別地,對3次B樣條有除有限支撐外,B樣條與Gauss函數比較相似,都為鐘形函數.Martheon和Dubrule曾討論過克立格估計和樣條間的等效性和相似性,特別是Unser等研究了Gauss函數與B樣條函數之間的相似關系,并建立了如下逼近公式取n=3,3次B樣條與Gauss函數已對比列于圖4.(41)式Gauss函數的參數轉化為常用的參數形式有不同階數的樣條對應的Gauss函數參數可列于表1.一般3～7次的樣條均可取得較好的逼近效果,否則階次過低或過高的逼近結果將不是太粗糙就是太平滑或支撐區(qū)間過寬.因此由類比可知Gauss函數的參數亦應與此相應,宜取為[0.82,1.16].這與我們前述的準則1的估計是相一致的.本節(jié)以上推導是以格網為單位長度導出的,如格網非單位長度時,可通過下列變換化格網為單位長度式中:Δx,Δy分別為格網長度,x,y和x′,y′分別為變換前后的坐標在實際應用時,數據一般是非格網散布在一定區(qū)域.設S為區(qū)域面積,N為域內數據個數,在假定x和y方向無異向性時等效的格網長度可直觀地由下式求出:試算證明用這種方法確定散亂數據的等效格網長度效果比較好.5gauss函數的海淡設理論模型為:f(x)=sin2πx,并令采樣間距Δx=0.0625,在0<x<1區(qū)間共有17個數據點,以這些數據點為基礎,用d=0.1,d=0.5和d=1.0(以Δx=0.0625為單位長度)時的Gauss函數去逼近該理論模型,圖5畫出了相應的逼近曲線示意圖.觀察圖5可以發(fā)現d=0.1和d=0.5時Gauss函數的逼近曲面產生了許多不應有的低谷和高峰,與原理論模型相差較大,這是由于d=0.1和d=0.5時Gauss函數支撐區(qū)間太窄,在半節(jié)點上達不到應有的預期結果,在半節(jié)點上出現了許多不應有的極小值.而d=1.0時Gauss函數產生了非常好的逼近效果,幾乎與理論模型一致.6協方差函數的函數應該說平穩(wěn)過程有著確定的譜S(w),通過Wiener＿Chintsch定理它又有著確定的協方差函數C(x),式中F和F-1分別為Fourier正演和反演.但是從另一方面來看,在一定條件下克立格可理解為由旋轉曲面生成的基函數進行的函數逼近,該基函數并不一定要取為被逼近過程的協方差函數,具有一定的可選擇性,并且,由于矩陣求逆的數值穩(wěn)定性等原因這種逼近還與觀測數據的取樣密度存在一定的內在聯系.對于同樣一個隨機過程,采樣密度變了,基函數(即“協方差函數”)中的參數亦應相應調整,否則將導致數值不穩(wěn)定性出現(如矩陣求逆失敗).:(1)在一定程度上與采樣數據的密度(或間隔)有一定聯系.(2)以一維Gauss函數為例,設Δx為數據間隔,則其最優(yōu)參數應取d=[0.8Δx,1.7Δ