非線性模型精度評定的Bootstrap方法及其加權(quán)采樣改進方法

王樂洋，李志強

東華理工大學(xué)測繪工程學(xué)院，江西 南昌 330013

在大地測量數(shù)據(jù)處理領(lǐng)域，嚴(yán)格的線性模型并不多見[1-3]，而非線性特征一般更接近客觀事物的性質(zhì)，已滲透于現(xiàn)代空間測量技術(shù)的各個領(lǐng)域，如地球重力場[4]、攝影測量與遙感、海洋測繪等領(lǐng)域，所涉及的函數(shù)模型普遍是非線性模型。傳統(tǒng)的獲取非線性模型的參數(shù)估值與精度信息的方法是最小二乘(least squares，LS)法，其函數(shù)模型為高斯-馬爾可夫(Gauss Markov，GM)模型，為線性模型[5-8]。隨著測量手段的多樣化和智能化，對平差結(jié)果的精度要求也越高，若采用傳統(tǒng)線性化的理論方法處理具有較高精度的觀測數(shù)據(jù)，可能會產(chǎn)生大于觀測誤差的模型誤差，進而導(dǎo)致平差結(jié)果的精度損失與特征轉(zhuǎn)變。因此，線性近似方法可能無法滿足部分高精度的要求，對非線性平差參數(shù)估計與精度評定方法的研究成為測繪領(lǐng)域所研究的熱點問題之一[9-11]，也是提升成果質(zhì)量的重要需求。

已有能夠解決非線性模型中未知量最佳估值獲取問題的方法主要有線性化法、迭代解法及搜索算法3種[12-13]。線性化法是用線性模型的求解理論與方法近似求解，其弊端在于當(dāng)模型的非線性程度較強時，將嚴(yán)重影響結(jié)果質(zhì)量。迭代解法包含牛頓迭代法、高斯-牛頓法、信賴域法、擬牛頓法等，此類方法的局限性在于對迭代初值較為敏感，對較差的初值可能會出現(xiàn)迭代不收斂現(xiàn)象[9,13-14]。搜索算法包含粒子群算法、模擬退火法、遺傳算法等，該類方法無法獲取參數(shù)的解析表達式，通常以犧牲計算耗時來代替求導(dǎo)計算[15-17]。

已有能夠獲取參數(shù)估值二階精度信息的非線性精度評定方法主要有近似函數(shù)法和近似非線性函數(shù)的概率密度分布兩種方法。近似函數(shù)法利用泰勒公式對非線性模型展開并截取至二次項，根據(jù)誤差傳播定律獲取參數(shù)估值的方差-協(xié)方差陣。文獻[18]推導(dǎo)了非線性函數(shù)協(xié)方差傳播的一般公式和含有二次項的協(xié)方差傳播公式。文獻[19]將高階泰勒級數(shù)展開式方法用于GIS誤差傳播中，豐富了近似函數(shù)法在非線性精度評定的理論研究。文獻[20]將總體最小二乘(total least squares，TLS)視為一般的非線性模型，通過泰勒級數(shù)展開，導(dǎo)出了參數(shù)估值二階精度的協(xié)方差陣計算公式和改正數(shù)的偏差計算公式。近似函數(shù)法具有固定的方差陣解析表達式，但無法避免求導(dǎo)計算。當(dāng)遇到未知參數(shù)向量與觀測值向量之間存在復(fù)雜非線性關(guān)系或觀測數(shù)據(jù)量較大的多維非線性函數(shù)時，求導(dǎo)計算復(fù)雜且獲取困難。

近些年，通過采樣來代替求導(dǎo)計算的近似非線性函數(shù)概率密度分布的精度評定方法流行起來，包括蒙特卡羅(Monte Carlo，MC)法[21-23]、Jackknife方法[24]等。蒙特卡羅法通過重復(fù)模擬隨機樣本來近似非線性函數(shù)的概率密度分布，進而輸出待統(tǒng)計量的精度信息。文獻[22]采用蒙特卡羅法計算近似單位權(quán)方差的偏差，并將偏差改正后的單位權(quán)方差估值用于求解參數(shù)估值的均方誤差。針對蒙特卡羅法模擬次數(shù)的選擇具有主觀性且無法對精度進行直接控制的問題，文獻[23]將自適應(yīng)蒙特卡羅法用于非線性模型精度評定，采用對偶變量的思想，給出了獲取參數(shù)估值偏差的對偶自適應(yīng)蒙特卡羅法的計算流程。當(dāng)代計算機技術(shù)的發(fā)展使得蒙特卡羅方法的實現(xiàn)成為了可能，但該方法的模擬次數(shù)通常需要設(shè)置在105次以上，它往往以增加計算耗時來獲得精度信息，并且隨著采樣次數(shù)的變化其精度統(tǒng)計信息并不唯一[25-26]。Jackknife采樣方法按順序逐次減小特定數(shù)目的觀測值來生成一系列Jackknife樣本，從而獲取方差信息。文獻[24]將Jackknife方法用于TLS精度評定，給出了獲取參數(shù)估值精度信息的詳細采樣步驟，進一步擴展了非線性精度評定的近似函數(shù)概率分布法理論。但Jackknife方法需要保證足夠的觀測數(shù)據(jù)才能獲得合理的精度信息，并且隨著d值的增大，其計算量也會迅速增加[26]。

Bootstrap方法[27-29]稱自助法，是與Jackknife方法緊密相關(guān)的一種統(tǒng)計推斷方法，均用于估計或修正統(tǒng)計估計值的偏差或方差信息。與Jackknife方法一樣，Bootstrap方法通過檢查樣本數(shù)據(jù)內(nèi)的變化，而不是通過參數(shù)假設(shè)來估計一個統(tǒng)計量的變異性，但Jackknife方法并沒有Bootstrap方法那么普及[30]。文獻[27]首次提出了通過重復(fù)采樣得到的自助世界的經(jīng)驗分布來逼近總體真實分布的思想，并給出了其基本采樣策略和相關(guān)理論證明。由于Bootstrap方法通過重采樣代替求導(dǎo)計算，只需選擇采樣次數(shù)并結(jié)合有放回抽樣策略獲取自助樣本，最后解算每個樣本便可求得未知統(tǒng)計量的均值和協(xié)方差陣，因此歷年來有較多學(xué)者將該方法用于偏差分析和方差計算[31-34]。目前，Bootstrap方法的研究主要集中于該方法在統(tǒng)計學(xué)上的性質(zhì)，將該方法用于解決大地測量領(lǐng)域非線性模型精度評定的問題有待研究；考慮到Bootstrap方法的原始采樣策略為等概率采樣，對于不等概率的采樣方式的適用性還有待驗證。

針對以上問題，本文從獲取更為合理的精度信息出發(fā)，針對現(xiàn)有非線性精度評定方法的不足，將Bootstrap方法用于獲取非線性模型參數(shù)估值的精度信息；顧及Bootstrap方法原始的等概率采樣方式，提出加權(quán)采樣策略的改進方法，以至獲取更加合理的方差估值；最后通過算例驗證本文非線性模型精度評定的Bootstrap方法的有效性及本文改進方法的正確性和優(yōu)勢。

1 非線性精度評定及其二階近似方差-協(xié)方差

將非線性平差函數(shù)模型定義為

(1)

參照文獻[35]的推導(dǎo)思路，可得參數(shù)估值的二階近似方差-協(xié)方差陣為[35]

(2)

其中，J、H可以具體表示為

(3)

(4)

式中，

(5)

(6)

(7)

(8)

(9)

式中，ei=[0 0 … 1 … 0 0]T∈Rn×1，第i處為1，其余為零。

由式(2)可以看出，在非線性模型精度評定中，若采用近似函數(shù)法獲取非線性模型參數(shù)估值的二階近似協(xié)方差信息，其泰勒級數(shù)展開過程需要計算非線性模型的Jacobian矩陣和Hessian矩陣，且當(dāng)遇到數(shù)據(jù)量較大的多維非線性函數(shù)或非線性程度強的模型時，偏導(dǎo)數(shù)的計算十分困難。因此，為避免復(fù)雜的求導(dǎo)計算并獲取更加合理的精度信息，本文將Bootstrap采樣方法及其加權(quán)采樣的改進方法引入非線性模型精度評定中，試圖有效解決近似函數(shù)法難以對非線性函數(shù)求導(dǎo)的問題。

2 非線性模型精度評定的Bootstrap方法

非線性平差模型的參數(shù)估計與精度評定問題可以理解為，利用從總體中抽取的固定原始觀測樣本對非線性模型中的未知統(tǒng)計量及其精度信息進行統(tǒng)計推斷。Bootstrap方法是將原始樣本作為總體，利用有放回隨機抽樣法從總體分布函數(shù)中得到一系列獨立樣本(稱為自助樣本)，并通過計算每個自助樣本來獲取未知統(tǒng)計量抽樣分布的經(jīng)驗估計[27,30]，最后利用這個估計的抽樣分布來做總體推斷。該方法的優(yōu)勢性主要表現(xiàn)在它不需要對未知模型及分布做任何假設(shè)，也無須推導(dǎo)估計量的精確表達式，只需通過檢驗樣本內(nèi)統(tǒng)計量的變化便能夠估計未知參數(shù)的整個抽樣分布。

(10)

顧及非線性優(yōu)化算法解算得到的參數(shù)估值不再是無偏估計量，而Bootstrap方法經(jīng)有放回抽樣獲取的大量自助樣本能夠減小參數(shù)估值的偏差，在一定程度上能夠有效改善參數(shù)估值的質(zhì)量[32]。因此，本文將式(10)的結(jié)果作為未知參數(shù)的自助估計。

(11)

Bootstrap方法的重采樣要點是采樣這個過程的隨機項。因此，根據(jù)重采樣對象的不同，Bootstrap方法又可分為：重采樣觀測值的Bootstrap方法和重采樣殘差的Bootstrap方法[30]。Bootstrap方法通過重采樣隨機項獲取自助樣本的方式代替復(fù)雜的或難以實現(xiàn)的求導(dǎo)計算，可以用于求取非線性平差參數(shù)估值的方差或協(xié)方差陣。因此，本文將這兩種采樣策略引入非線性模型的精度評定問題中，給出了獲取參數(shù)方差的具體步驟。

2.1 重采樣觀測值的Bootstrap方法

重采樣觀測值的Bootstrap方法，通過對觀測值的有放回隨機抽樣將原始觀測數(shù)據(jù)采樣成多個自助樣本，原始樣本集中的每個樣本點被采樣到的概率相同，每個自助樣本的樣本容量與原始樣本相同。對于不等精度的觀測數(shù)據(jù)，需要保證自助樣本中的觀測值與其權(quán)值相對應(yīng)，因此在重采樣獲得觀測值自助樣本的過程中需要對觀測值的權(quán)Pi進行重采樣，兩者采樣方式一致，具體表現(xiàn)為以隨機數(shù)序列元素為自助樣本元素的下標(biāo)，對觀測值及其權(quán)值進行賦值，進而獲取自助樣本及對應(yīng)權(quán)陣。根據(jù)以上分析，總結(jié)得到重采樣觀測值的Bootstrap方法用于非線性精度評定的完整計算流程：

(1)假設(shè)觀測值N=(x1,x2,…,xn)為隨機項樣本X1，其權(quán)陣為P=diag(P1,P2,…,Pn)，將1/n的概率分別設(shè)置在(x1,x2,…,xn)各樣本值上。

(2)產(chǎn)生滿足均勻分布的n個隨機數(shù)(i1,i2,…,in)，其中1≤ij≤n,j=1,2,3,…,n。

(5)循環(huán)步驟(2)—步驟(4)M次，得到M組參數(shù)估計值。

(6)利用式(10)獲取參數(shù)均值，利用式(11)獲取參數(shù)估值的方差。

重采樣觀測值的Bootstrap方法的采樣過程將原始觀測值采樣成了多個自助樣本，盡管每個自助樣本的樣本容量與原始樣本相同，但并不是改變觀測值的先后排列順序，有放回隨機抽樣過程使得每個自助樣本中可能存在重復(fù)的原始數(shù)據(jù)點，而另外一些樣本點沒有出現(xiàn)。因此，每個自助樣本將隨機地異于原始樣本，導(dǎo)致每個自助樣本獲得的參數(shù)估值存在細微差異。

2.2 重采樣殘差的Bootstrap方法

文獻[27]在提出Bootstrap方法的同時，對抽取獨立同分布樣本點的采樣方式進行了延伸和拓展，結(jié)合多元回歸分析給出了一種對殘差進行抽樣的采樣策略。具體做法是，首先使用所有樣本點建立模型并獲取殘差，然后對殘差進行采樣并根據(jù)模型的方程重構(gòu)觀測向量，最后計算參數(shù)并重復(fù)該步驟。因此，當(dāng)模型的自變量為固定常數(shù)，因變量為自變量和一個隨機誤差項的函數(shù)時，基于模型的誤差結(jié)構(gòu)，可將因變量的殘差設(shè)置為采樣過程的隨機項。

重采樣殘差的Bootstrap方法與重采樣觀測值的Bootstrap方法相比，主要的區(qū)別在于重采樣對象不同。采用重采樣殘差的Bootstrap方法解決精度評定問題的計算步驟為：

(7)重復(fù)步驟(4)—步驟(6)M次，共得到M組參數(shù)估值。

(8)通過式(10)獲取未知參數(shù)的自助估計值，利用式(11)計算參數(shù)估值的方差。

重采樣觀測值與重采樣殘差的樣本數(shù)據(jù)均來源于原始固定樣本，并未根據(jù)更多的觀測信息進行計算。因此，這兩種采樣方法均是利用相同的觀測數(shù)據(jù)最終得到未知參數(shù)估值及方差信息，僅是采樣過程中的隨機項不同，其采樣過程均不會產(chǎn)生額外的模型誤差，也不會改變模型態(tài)性；并且循環(huán)有放回隨機采樣過程獲得的大量自助樣本能夠減小迭代方法求解非線性參數(shù)估值的偏差，從而提高參數(shù)估值的精確度。

3 改進的Bootstrap方法

Bootstrap方法在重采樣隨機項元素構(gòu)建自助樣本過程中，每個元素被采樣到的概率相同且均為1/n?？梢钥闯觯珺ootstrap方法在采樣時將樣本數(shù)據(jù)視為等精度觀測，在構(gòu)建自助樣本過程中忽略了數(shù)據(jù)的權(quán)值信息，使得精度不等的觀測數(shù)據(jù)出現(xiàn)在自助樣本中的概率卻相同，因此可能無法較好地逼近總體分布函數(shù)，在一定程度上有損理論嚴(yán)密性。

為使得Bootstrap方法能夠更加充分利用總體包含在樣本中的先驗信息和數(shù)據(jù)性質(zhì)，使得自助樣本的經(jīng)驗分布更加逼近總體分布函數(shù)，若利用觀測值的分布信息來構(gòu)建經(jīng)驗分布函數(shù)，通常需要借助分布假設(shè)，并且所得分布函數(shù)將嚴(yán)重依賴于所作的假設(shè)。當(dāng)觀測樣本的分布假設(shè)與實際觀測樣本的分布存在較大偏差時，所得樣本經(jīng)驗分布函數(shù)并不能很好地逼近總體分布函數(shù)，導(dǎo)致統(tǒng)計推斷結(jié)果與實際存在一定的偏差[36]。顧及觀測值的權(quán)的獲取方式較觀測值的分布信息獲取方式更為簡便，并且利用觀測權(quán)來構(gòu)建經(jīng)驗分布函數(shù)可以有效避免對觀測樣本分布的假設(shè)。因此，權(quán)作為比較觀測值之間精度高低的一種先驗信息，可用于調(diào)整隨機項的抽樣概率，以至獲取更加逼近總體分布函數(shù)的樣本經(jīng)驗分布函數(shù)。因此，本文將原始觀測樣本的權(quán)值納入重采樣獲取自助樣本的過程。在該過程中，權(quán)值被用于構(gòu)建一個自助樣本獲取準(zhǔn)則，通過對隨機項元素進行采樣可能性的重新分配，使得每個元素被抽樣到的概率與其對應(yīng)的觀測精度相匹配。

首先，對原始樣本所有元素的權(quán)值進行歸一化預(yù)處理，使得歸一化權(quán)值之和為1

(12)

式中，NORMi表示第i個觀測元素的歸一化權(quán)值；Pi是原始觀測樣本中第i個元素的權(quán)；n為原始樣本的樣本容量。

將NORMi作為隨機項中各元素xi的采樣概率，其概率值的大小取決于各元素所對應(yīng)的權(quán)值與所有觀測權(quán)之和的比值大小。因此，可構(gòu)建Bootstrap采樣過程中隨機變量的分布律函數(shù)

proxi=NORM{Y=xi}

(13)

式中，xi為隨機項X中的樣本值；Y為隨機項X中的隨機變量，其可能的取值為xi(i=1,2,…,n)；proxi表示事件{Y=xi}的概率。

(14)

由于Bootstrap法在采樣過程中隨機項的選擇不同，Bootstrap法可區(qū)分為重采樣觀測值的Bootstrap方法和重采樣殘差的Bootstrap方法，將以上加權(quán)采樣策略應(yīng)用到該兩種方法中，分別給出將改進方法用于非線性模型精度評定中的詳細計算步驟。

3.1 改進的重采樣觀測值的Bootstrap方法

該方法首先對觀測數(shù)據(jù)的權(quán)值進行歸一化處理，計算每個樣本元素被采樣到自助樣本的概率值；構(gòu)建采樣過程中的隨機變量的分布律函數(shù)；產(chǎn)生滿足該分布律函數(shù)的隨機數(shù)，通過賦值獲取自助樣本并計算參數(shù)；然后將生成自助樣本到獲取參數(shù)的整個過程循環(huán)M次；最后對未知統(tǒng)計量的均值和方差進行統(tǒng)計推斷。每個自助樣本的樣本容量與原始樣本相同，但原始樣本中的每個元素被采樣到的概率受其權(quán)值影響而不同，這取決于該樣本點的權(quán)值大小。根據(jù)以上分析，總結(jié)得到改進的重采樣觀測值的Bootstrap方法的計算流程：

(1)假設(shè)觀測值樣本N=(x1,x2,…,xn)為采樣過程中的隨機項X1，其權(quán)陣為P=diag(P1,P2,…,Pn)。

(2)利用式(12)對權(quán)陣進行歸一化處理，使得歸一化后的權(quán)值之和為1。

(4)通過式(13)、式(14)構(gòu)建隨機變量的經(jīng)驗分布函數(shù)。

(5)產(chǎn)生滿足該分布函數(shù)的n個隨機數(shù)(i1,i2,…,in)，其中1≤ij≤n,j=1,2,3,…,n。

(9)重復(fù)步驟(5)—步驟(8)M次，得到M組參數(shù)估計值。

(11)分別通過式(10)、式(11)獲取參數(shù)的均值和方差。

重采樣觀測值的Bootstrap方法及其改進方法均是將相同的原始觀測數(shù)據(jù)重采樣成了多個自助樣本，且獲取的自助樣本個數(shù)及每個自助樣本的樣本容量均相同，但改進方法由于采用了加權(quán)采樣的策略，將等概率重采樣轉(zhuǎn)變?yōu)椴坏雀怕什蓸臃绞?，使得?quán)值大的原始樣本元素出現(xiàn)在自助樣本中的概率增大，從而使得自助樣本能夠更加貼合原始樣本。

3.2 改進的重采樣殘差的Bootstrap方法

與重采樣殘差的Bootstrap方法一致，改進的重采樣殘差的Bootstrap方法仍將因變量的殘差作為采樣過程中的隨機項，主要區(qū)別在于改進方法在采樣前對權(quán)值進行了歸一化預(yù)處理，通過獲取隨機項中隨機變量的分布函數(shù)，使得對殘差項的等概率采樣轉(zhuǎn)化為加權(quán)采樣。改進的重采樣殘差的Bootstrap精度評定方法的計算步驟總結(jié)為：

(1)從總體中獲取固定的觀測向量L及其權(quán)陣P，作為擬合該模型的原始樣本。

(5)對權(quán)陣P的對角元素Pi進行歸一化處理。

(7)通過式(13)、式(14)構(gòu)建隨機變量的經(jīng)驗分布律函數(shù)。

(8)產(chǎn)生滿足該分布律的隨機數(shù)(i1,i2,…,in)，其中1≤ij≤n,j=1,2,3,…,n。

(12)循環(huán)步驟(8)—步驟(11)M次，得到M組參數(shù)估值。

(14)分別通過式(10)、式(11)獲取未知參數(shù)的自助估計值及方差。

重采樣殘差方法和改進的重采樣殘差方法在自助樣本獲取過程中，均未對原始觀測樣本進行采樣，因此均利用了原始樣本中的所有觀測值，但對隨機項的有放回抽樣使得每個自助樣本包含的元素不同，進而使得每個自助樣本得到的估值存在差異。

根據(jù)以上4種精度評定流程可以發(fā)現(xiàn)，不管是采樣觀測值，還是采樣殘差，抑或是改進方法，重復(fù)M次的Bootstrap采樣過程僅由有放回抽樣的隨機性和滿足的抽樣分布所決定，即第i次的自助樣本獲取過程并不會影響第j次的重采樣，這使得自助樣本的獲取過程是獨立的。此外，由于自助樣本中的所有元素均來源于原始樣本，未利用其他的觀測值進行計算，并且由于受自助法自身有放回抽樣性質(zhì)以及抽樣隨機性影響，使得自助樣本之間可能包含重復(fù)的隨機項元素，導(dǎo)致自助樣本存在一定的相關(guān)性。第i組參數(shù)估值僅僅由第i個自助樣本決定，因此每個自助樣本的解算過程是獨立的；但由于自助樣本i與自助樣本j可能包含相同元素，使得M組參數(shù)估值存在一定的相關(guān)性，使得各自助樣本的參數(shù)估值在自助均值左右浮動并且較為接近。但在采用矩法估計計算自助方差時，并不能簡單忽略M組參數(shù)估值的分布規(guī)律或者將其視為相同，這是因為自助重采樣的目的便是要得到一個較好的抽樣分布估計來對模型中的未知統(tǒng)計量進行統(tǒng)計推斷。

4 案例研究與分析

經(jīng)試驗，Bootstrap采樣方法適用于線性平差模型的參數(shù)估計與精度評定問題，通過有放回隨機抽樣能夠在一定程度上提升參數(shù)質(zhì)量，并且獲取合理的參數(shù)估值精度信息，表明在線性模型精度評定問題中，Bootstrap方法具有可靠性和優(yōu)勢性。與最小二乘法相比，盡管實施Bootstrap采樣會增加精度評定過程計算耗時，但自助法具有操作簡單、易于編碼實現(xiàn)的優(yōu)點，并且其結(jié)果存在一定程度上的質(zhì)量提升，不失為一種較好的精度評定方法。在非線性模型精度評定問題中，本文第一節(jié)已證明已有的基于泰勒級數(shù)展開的近似函數(shù)法，在精度評定過程中需要計算非線性函數(shù)的一階或二階偏導(dǎo)數(shù)，當(dāng)遇到高維且非線性程度較強的模型時，求導(dǎo)過程復(fù)雜難以獲取。而自助法用于解決非線性模型精度評定問題最大的優(yōu)勢在于，獲取參數(shù)精度信息過程中能夠避免導(dǎo)數(shù)計算，該方法是通過對原始樣本的重采樣過程來代替復(fù)雜的Jacobian矩陣或Hessian矩陣計算過程，最后利用抽樣分布估計對未知參數(shù)進行統(tǒng)計推斷。因此，為充分體現(xiàn)Bootstrap方法在測量平差精度評定問題中的優(yōu)勢性，本文考慮的是非線性模型情形。

為評估Bootstrap方法和本文改進方法在獲取非線性模型參數(shù)估值精度信息方面的性能，以及探討在大地測量領(lǐng)域中的應(yīng)用，本文共設(shè)計了4個算例：兩個非線性強度不同的非線性回歸模型、邊角網(wǎng)平差模型、火山形變Mogi模型。通過算例1和算例2來驗證Bootstrap方法在非線性模型精度評定中的可行性及本文改進方法的正確性；通過算例3和算例4來展示Bootstrap方法在大地測量非線性精度評定中的應(yīng)用價值，并驗證本文改進方法在參數(shù)方差獲取方面的優(yōu)勢。對比試驗中，Monte Carlo方法的模擬次數(shù)取值為105。Bootstrap方法重采樣次數(shù)的選擇取決于統(tǒng)計量的具體研究問題并且依賴于需要做的檢驗[38]。文獻[38—41]指出，當(dāng)重采樣次數(shù)為50～200時，Bootstrap方法能夠得到較為合理的未知參數(shù)方差或中誤差；當(dāng)采樣次數(shù)超過200時，方差近似值的改進很小。經(jīng)試驗，通過增加Bootstrap重采樣次數(shù)在一定程度上能夠提升精度評定結(jié)果質(zhì)量，但隨著采樣次數(shù)的增加，中誤差結(jié)果趨于較為穩(wěn)定的值。這與已有文獻的結(jié)論較為一致。因此，本文依據(jù)已有研究的經(jīng)驗取值并顧及非線性平差精度評定問題，在保證較高的計算精度條件下，同時考慮較高的精度評定效率，將Bootstrap方法的采樣次數(shù)設(shè)置為103。案例研究中出現(xiàn)的字符及對應(yīng)含義列于表1。

表1 字母縮寫及對應(yīng)含義Tab.1 Abbreviations and their corresponding meanings

4.1 Bootstrap方法及本文改進方法的可行性及正確性驗證分析

針對非線性回歸模型的非線性強度不同，設(shè)計兩個仿真試驗，通過Bootstrap方法及加權(quán)采樣方法在參數(shù)估計中的應(yīng)用，以及與AF法、JK法、MC法在參數(shù)標(biāo)準(zhǔn)差估計中的性能評估對比，以清楚地揭示本文精度評定方法的可行性和有效性。

4.1.1 算例1

算例1采用的非線性回歸模型函數(shù)形式如下

(15)

式中，xi為n維自變量向量的元素；yi為對應(yīng)因變量值；εyi為因變量yi的隨機誤差；ξ1、ξ2為回歸參數(shù)。

4.1.2 算例2

為進一步探討AF法、JK法、MC法、Bootstrap方法及本文加權(quán)采樣方法在更復(fù)雜非線性情況下的方差估計性能和計算特性，采用比算例1非線性程度更強的回歸函數(shù)，其函數(shù)形式定義如下

(16)

式中，xi、yi分別為橫坐標(biāo)x和縱坐標(biāo)y的自變量；zi為n維自變量向量中的元素；εzi為zi的隨機誤差，ξ1、ξ2、ξ3、ξ4為回歸參數(shù)。

4.1.3 結(jié)果分析

在獲取模擬觀測數(shù)據(jù)后，分別采用LS法、GN法、BO法、MBO法、BR法及MBR法解算并獲取參數(shù)均值，其結(jié)果及差值二范數(shù)列于表2；再分別采用BO法、MBO法、BR法及MBR法對隨機項進行采樣，并通過高斯-牛頓解法解算自助樣本，獲取回歸系數(shù)的標(biāo)準(zhǔn)差及協(xié)方差信息，最后與FO法、SO法、JK法及MC法的結(jié)果進行對比，各方法的結(jié)果列于表3。

表2 非線性回歸模型參數(shù)估值及其與真值之間的二范數(shù)Tab.2 The estimated parameters and the norm of the nonlinear regression model

表3 非線性回歸模型參數(shù)估值的標(biāo)準(zhǔn)差和協(xié)方差計算結(jié)果Tab.3 The standard deviations and covariances of parameters estimations of the nonlinear regression model

由表2的參數(shù)估值及范數(shù)結(jié)果可以看出，在非線性強度不同的兩個仿真設(shè)置中，LS法的結(jié)果均為最差，表明LS法在將非線性模型線性近似的過程中，引入了較大的模型誤差，從而嚴(yán)重影響估計值的質(zhì)量。而相比于LS法，GN法得到的參數(shù)估值較優(yōu)，說明GN法能夠削弱非線性回歸模型線性近似所帶來的模型誤差的影響，其循環(huán)迭代過程可以不斷修正參數(shù)估值，在一定程度上能夠有效改善傳統(tǒng)LS法在對非線性模型線性化過程中所引起的精度損失，使其參數(shù)估值比傳統(tǒng)的線性近似方法更接近參數(shù)真值。BO法和BR法得到的參數(shù)估值與真值的范數(shù)均比LS法和GN法的結(jié)果更小，表明BO法和BR法通過實施Bootstrap并采用GN法解算自助樣本，在削弱線性化帶來的影響的基礎(chǔ)上，能夠在一定程度上減小參數(shù)估值的偏差，獲取更為精確的參數(shù)估值。與BO法相比，MBO法得到的參數(shù)估值與真值的范數(shù)更小，表明采用了加權(quán)采樣的重采樣觀測值方法在參數(shù)的精確度方面得到了較明顯的改善；而相比于BR法，MBR法也能夠提升參數(shù)估值的精確度。這表明MBO法及MBR法在獲取較高精度的參數(shù)估值方面比BO法和BR法更加有效。試驗結(jié)果表明，本文的加權(quán)采樣策略在非線性回歸模型的參數(shù)估計方面是可行且有效的，改進方法通過重采樣獲得的自助樣本能夠充分利用原始樣本的先驗信息，從而使得MBO法和MBR法獲得的估值較BO法和BR法更接近參數(shù)真值。

由表3可知，在兩個不同的仿真設(shè)置下，二階近似函數(shù)法獲取的參數(shù)估值標(biāo)準(zhǔn)差在數(shù)值上略大于一階近似函數(shù)法得到的標(biāo)準(zhǔn)差結(jié)果，且更接近于蒙特卡羅方法獲取的參數(shù)標(biāo)準(zhǔn)差。這是因為二階近似函數(shù)法顧及了非線性平差模型的二階項，說明由一階近似函數(shù)法得到的結(jié)果可能高估了非線性平差模型參數(shù)估值的精度。但由于二階近似函數(shù)法也僅僅包含至二階項，忽略了高階項，因此也說明了蒙特卡羅方法獲取的精度信息更加合理。并且試驗證明，當(dāng)原始樣本量增大時，二階近似函數(shù)法所涉及的Hessian矩陣維數(shù)也將增大，也說明了無須求導(dǎo)計算的Bootstrap方法更加適用較大量的觀測數(shù)據(jù)。相比于蒙特卡羅方法的結(jié)果，Jackknife法獲得的標(biāo)準(zhǔn)差偏離較大，是因為Jackknife法易受原始樣本容量大小的影響且其重采樣過程獲取的Jackknife樣本較少，與文獻[24]中的結(jié)論一致。這表明Jackknife法往往會低估非線性模型參數(shù)估值的精度，因此該方法無法反映非線性模型參數(shù)估值的真實精度。

利用基于Bootstrap的方法對該非線性函數(shù)進行精度評定所獲得的標(biāo)準(zhǔn)差與近似函數(shù)法的結(jié)果在符號和量級上相同、在數(shù)值上相近，表明該4種方法均能夠?qū)Ψ蔷€性平差模型進行有效的精度評定。從計算結(jié)果上看，BO法和BR法獲取的參數(shù)的標(biāo)準(zhǔn)差小于Jackknife法的結(jié)果且更接近于蒙特卡羅方法的計算精度，表明BO法和BR法較Jackknife法更加精確。相比于BO法，MBO法獲取的參數(shù)估值標(biāo)準(zhǔn)差在數(shù)值上更小，且與蒙特卡羅方法的結(jié)果更為接近；MBR法計算得到的參數(shù)標(biāo)準(zhǔn)差比BR法的結(jié)果更小并且同樣更靠近于蒙特卡羅方法獲取的標(biāo)準(zhǔn)差，說明MBO法和MBR法能夠獲取更為合理的參數(shù)標(biāo)準(zhǔn)差。結(jié)合表2中的參數(shù)估值結(jié)果，說明本文加權(quán)采樣的獲取自助樣本方式不僅適用于重采樣觀測值，也適用于重采樣殘差情形。試驗結(jié)果表明，將Bootstrap方法用于解決非線性參數(shù)估計與精度評定問題是可行且有效的，本文加權(quán)采樣方法獲取的更為精確的參數(shù)估值和精度信息，也證明了本文改進方法更具可靠性和優(yōu)勢。

4.2 算例3：Bootstrap方法及本文改進方法在邊角網(wǎng)平差模型中的應(yīng)用

將文獻[42]中的例3-2-1進行改化，獨立不等精度觀測如圖1所示，邊角網(wǎng)有12個角度和6條邊長，其中A、B、C為已知點，P、D為待定點。起算數(shù)據(jù)和觀測值分別列于表4和表5。

圖1 邊角網(wǎng)Fig.1 The sketch of triangulateration network

表4 起算數(shù)據(jù)Tab.4 The data of known control points m

在獲取角度和邊長觀測數(shù)據(jù)后，首先采用余切公式計算得到待定點坐標(biāo)估值作為迭代初值，再分別采用BO法、MBO法、BR法及MBR法對隨機項采樣并通過高斯-牛頓迭代獲取未知點坐標(biāo)，其參數(shù)結(jié)果及差值二范數(shù)列于表6。通過算例1和算例2的結(jié)果可知，Jackknife法的結(jié)果不足以反映參數(shù)估值的真實精度；而一階近似方法的結(jié)果可能會高估參數(shù)估值的精度；并且根據(jù)邊長觀測方程及測角觀測方程可以看出，邊角網(wǎng)平差模型中的未知點坐標(biāo)與觀測值向量之間存在復(fù)雜的非線性關(guān)系，若采用泰勒級數(shù)展開并截取至二次項的方法獲取未知參數(shù)的方差信息，需要進行多次計算偏導(dǎo)，其求導(dǎo)過程十分復(fù)雜且較難獲得。針對以上問題，因此算例3不再使用一階近似函數(shù)法、二階近似函數(shù)法及Jackknife方法求解參數(shù)估值的精度信息，僅通過基于Bootstrap的4種方法與MC法進行對比，來驗證本文方法的正確性和優(yōu)勢。待定點坐標(biāo)標(biāo)準(zhǔn)差及協(xié)方差結(jié)果列于表7。

表5 角度和邊長觀測數(shù)據(jù)及其權(quán)值Tab.5 The data of angular and traverse measurement and corresponding weights

表6 邊角網(wǎng)平差模型坐標(biāo)估值及其與真值之間的二范數(shù)Tab.6 The estimated coordinate parameters and the norm of the triangulateration network model

對比表6中BO法、MBO法、BR法及MBR法的結(jié)果可以看出，4種方法所得坐標(biāo)估值均與坐標(biāo)真值較為接近，再次驗證了Bootstrap方法及本文改進方法的正確性。同時，通過對比估值與真值的范數(shù)可知，MBO法所得估值較BO法的結(jié)果更接近于真值；與BR法的結(jié)果相比，MBR法所得結(jié)果的范數(shù)有較為明顯的改善。這表明MBO法及MBR法在獲取更高精度的邊角網(wǎng)坐標(biāo)估值方面比BO法和BR法更加有效。本算例的參數(shù)估值結(jié)果證明了Bootstrap方法在邊角網(wǎng)平差模型中依然適用，也再次證明了本文改進方法的有效性及加權(quán)采樣的必要性。

由表7可知，在邊角網(wǎng)平差模型中，BO法、BR法、MBO法及MBR法計算得到的待定點坐標(biāo)估值標(biāo)準(zhǔn)差在符號和量級上均與MC法的結(jié)果保持一致，同樣驗證了該4種方法在非線性平差精度評定問題中的正確性和有效性。由于邊角觀測方程的復(fù)雜非線性關(guān)系，使得二階近似函數(shù)法獲取參數(shù)中誤差較為困難，從側(cè)面體現(xiàn)了通過執(zhí)行Bootstrap采樣來避免復(fù)雜求導(dǎo)運算仍能獲取合理精度信息的優(yōu)勢，也說明了在面對高維復(fù)雜非線性模型的精度評定問題時，基于Bootstrap的方法具有更強的適用性。相比于BO法及BR法計算得到的標(biāo)準(zhǔn)差結(jié)果，MBO法和MBR法的結(jié)果在數(shù)值上更小，這表明通過充分利用觀測值先驗信息的加權(quán)采樣方法能夠在一定程度上提升參數(shù)精度的可靠性，也說明為更加合理反映參數(shù)估值的精度信息，采用加權(quán)采樣的抽樣方式更加有效。

表7 邊角網(wǎng)平差模型待定點坐標(biāo)標(biāo)準(zhǔn)差和協(xié)方差計算結(jié)果Tab.7 The standard deviations and covariances of the estimated coordinate parameters in the triangu-lateration network model

4.3 算例4：Bootstrap方法以及本文改進方法在火山形變Mogi模型反演中的應(yīng)用

為進一步驗證Bootstrap方法及本文改進方法在大地測量非線性精度評定中的廣泛適用性和應(yīng)用價值，算例4采用文獻[43]首次提出的Mogi模型。該模型是一種研究由火山爆炸源所引起的地表形變位移及重力變化的地球物理模型，它能夠有效地用于模擬和解釋大量的火山區(qū)地表形變。應(yīng)用火山區(qū)的地表形變觀測數(shù)據(jù)來反演巖漿壓力源參數(shù)，對火山的危險性評價有著重大的實際意義。地表位移又可分為垂直位移和水平位移，本文采用垂直位移單一反演的Mogi模型。

假設(shè)地面直角坐標(biāo)系為(x,y,z)，巖漿壓力源的中心坐標(biāo)表示為(x0,y0,D)，因此巖漿房膨脹所引起的地表垂直位移與等效壓力源參數(shù)之間的表達式可以表示為[44]

(17)

式中，Δh表示地表垂直位移；K表示為體積彈性模量；μ表示剪切模量；D表示為源的中心深度；r表示地表徑向距離且當(dāng)r=0時垂直位移取到最大值；ΔV為巖漿房體增量。

當(dāng)?shù)貧椴此山橘|(zhì)時，取泊松比v=0.25且顧及巖漿房體積的膨脹為半徑為R的等效球體，因此體積彈性模量K與剪切模量μ存在以下關(guān)系式

(18)

將式(18)代入式(17)，地表垂直位移可以化簡為

(19)

假設(shè)地表存在n個觀測點，則第i個監(jiān)測點到巖漿壓力源的地表徑向距離為

r2=(xi-x0)2+(yi-y0)2

(20)

式中，(x0,y0)為巖漿壓力源的中心在平面上的投影坐標(biāo)。

將式(20)代入式(19)，地表垂直位移可表示為一個多維非線性函數(shù)

Δhi+ei=fi(ξ)=

(21)

式中，i=1,2,…,n；ei為第i個監(jiān)測點的垂直位移Δhi的隨機誤差；f(·)表示未知參數(shù)向量ξ與觀測值之間的非線性映射關(guān)系；ξ=[ΔVDx0y0]T為待求的未知參數(shù)。

利用泰勒級數(shù)對式(21)在參數(shù)初值ξ0處展開并舍去二次及以上項，得到垂直位移反演壓力源參數(shù)的平差模型

(22)

通過仿真獲取某一火山區(qū)域的垂直位移觀測資料：觀測點的取值范圍為[x,y]∈[-5,5] km×[-5,5] km，相鄰垂直形變監(jiān)測點的間隔為0.5 km，每個觀測點的權(quán)值在[0,25]之間隨機取值。壓力源參數(shù)真值設(shè)置為：體積增量ΔV=6000×103m3、源的中心深度D=4 km、膨脹源的中心坐標(biāo)x0=0 km、y0=0 km；將地面監(jiān)測點坐標(biāo)及壓力源參數(shù)真值代入Mogi模型，正演得到地表垂直位移；最后對垂直位移模擬位移值施加均值為0、單位權(quán)方差為2 mm的隨機誤差，得到的垂直位移如圖2所示。

圖2 由火山Mogi模型正演得到的地表垂直形變Fig.2 The vertical displacement of forward the volcano Mogi model

由于垂直位移反演壓力源參數(shù)的平差模型是由泰勒級數(shù)對式(17)進行線性化處理后獲得，針對反演過程中可能出現(xiàn)的病態(tài)奇異情形[45]，算例4采用嶺估計法進行參數(shù)求解，分別采用BO法、MBO法、BR法和MBR法進行垂直位移單一反演，具體表現(xiàn)為在獲取自助樣本后，通過嶺估計方法來解決法矩陣求逆不穩(wěn)定的問題并獲取巖漿壓力源參數(shù)[46-47]。4種方法反演得到的參數(shù)估值結(jié)果及差值二范數(shù)列于表8。

表8 火山Mogi模型壓力源參數(shù)反演結(jié)果及其與真值之間的二范數(shù)Tab.8 Inversion results of pressure source parameters and the norm of the volcano Mogi model

由表8中的參數(shù)估值可以看出，本文提出的4種方法獲取的估值與真值的差異主要表現(xiàn)在巖漿房體積增量和巖漿壓力源的中心深度。MBO法所得結(jié)果優(yōu)于BO法，其參數(shù)估值與真值之間的范數(shù)減小了40.7%；與BR法相比，MBR法的范數(shù)減小了58.5%。結(jié)果表明，本文方法均能得到合理的參數(shù)估值，而采用了加權(quán)采樣的MBO法和MBR法在提升估值質(zhì)量方面更具優(yōu)勢。試驗結(jié)果說明，將Bootstrap方法用于法矩陣病態(tài)性嚴(yán)重的火山形變Mogi模型參數(shù)反演是可行且有效的，也說明了本文的加權(quán)采樣策略同樣適用于大地測量領(lǐng)域病態(tài)模型的參數(shù)估計問題。

由式(17)—式(22)可以看出，火山Mogi模型中的未知參數(shù)向量ξ與地表垂直位移觀測值向量Δh之間存在復(fù)雜的高維非線性關(guān)系，若采用對非線性模型泰勒級數(shù)展開的近似函數(shù)方法獲取未知參數(shù)的方差信息，需要進行多次計算偏導(dǎo)，使得精度評定過程十分復(fù)雜。因此，算例4僅采用本文提出的4種方法及MC法進行精度評定，各方法獲取的壓力源參數(shù)標(biāo)準(zhǔn)差結(jié)果列于表9。

由表9可知，BO法、BR法、MBO法及MBR法計算得到的參數(shù)估值標(biāo)準(zhǔn)差與MC法的結(jié)果在符號和量級上均保持一致，說明了本文Bootstrap方法精度評定的可行性和有效性。從估值的近似標(biāo)準(zhǔn)差結(jié)果來看，相比于BO法所獲得的參數(shù)標(biāo)準(zhǔn)差，MBO法計算得到的標(biāo)準(zhǔn)差比BO法的計算結(jié)果在數(shù)值上更小，且比BO法的結(jié)果更加接近于MC法所獲得的精度，說明MBO法能夠獲取比BO法更為合理的精度信息；與BR法的計算結(jié)果相比，MBR法的標(biāo)準(zhǔn)差有較為明顯的減小且更為靠近MC法的計算精度，結(jié)果同樣說明了MBR法計算得到的標(biāo)準(zhǔn)差比BR法更為精確。試驗結(jié)果表明，將Bootstrap方法用于獲取病態(tài)非線性平差模型的參數(shù)估值及精度信息依然是可靠且有效的，本文改進方法在集結(jié)Bootstrap方法的所有優(yōu)點外，通過加權(quán)采樣的抽樣方式還能夠更加充分地利用權(quán)值信息，從而能夠獲得更加精確的估值與精度信息值。

表9 火山Mogi模型壓力源參數(shù)的標(biāo)準(zhǔn)差和協(xié)方差計算結(jié)果Tab.9 The standard deviations and covariances of the pressure source parameters of the volcano Mogi model

5 結(jié) 論

為獲得統(tǒng)計意義上更為合理的精度信息，若采用泰勒展開的近似函數(shù)法，則無法避免地需要計算非線性函數(shù)的Jacobian矩陣和Hessian矩陣；若采用蒙特卡羅方法，則通常需要105以上次的模擬計算。因此，提出無需求導(dǎo)計算且計算效率相對較高的非線性精度評定方法顯得尤為重要。本文在闡述非線性模型精度評定原理及Bootstrap方法重采樣原理的基礎(chǔ)上，給出了非線性模型精度評定重采樣觀測值Bootstrap方法及重采樣殘差的Bootstrap方法，并分別給出了具體的采樣步驟，試圖通過這兩種方法在避免求導(dǎo)計算的情形下依然能獲取合理的精度評定結(jié)果。針對Bootstrap方法的重采樣過程是對模型隨機項的等概率采樣，本文通過對原始樣本的權(quán)重信息進行歸一化處理并獲取隨機變量的經(jīng)驗分布函數(shù)，將傳統(tǒng)自助法的等概率采樣轉(zhuǎn)化為加權(quán)采樣，嘗試能夠更加充分利用總體包含在原始樣本中的先驗信息和數(shù)據(jù)性質(zhì)，以至于自助樣本的經(jīng)驗分布更加逼近總體分布，進而獲取更為精確的參數(shù)方差估值。

對本文案例研究中的各方法進行對比分析可知，非線性模型精度評定的Bootstrap方法能夠獲取比近似函數(shù)法、Jackknife方法更為合理的參數(shù)標(biāo)準(zhǔn)差，這說明將Bootstrap方法用于解決非線性模型的精度評定問題是可行且有效的；同時，采用了本文給出的改進算法能夠得到比直接實施Bootstrap更加精確的精度信息，說明了本文加權(quán)采樣方法的可靠性與優(yōu)勢性。因此，建議在獲取非線性模型參數(shù)估值的方差-協(xié)方差時，采用本文的加權(quán)采樣方法。

本文方法適用于觀測值相互獨立的情況，Bootstrap方法在相關(guān)觀測數(shù)據(jù)的非線性精度評定問題，以及該方法的進一步應(yīng)用是接下來需要研究的內(nèi)容。隨著Bootstrap方法引入非線性模型精度評定中，拓展該方法在大地測量領(lǐng)域中的應(yīng)用，也將是非線性平差理論研究的一個重要補充。