王 珂，劉國林，付政慶，王路遙

1.山東理工大學建筑工程學院，山東 淄博 255000；2.山東科技大學測繪與空間信息學院，山東 青島 266590

非線性最小二乘法作為基本的非線性優(yōu)化方法，廣泛應(yīng)用于非線性函數(shù)模型的參數(shù)估計中[1-6]。然而，根據(jù)實際問題的來源，許多非線性函數(shù)模型具有特殊結(jié)構(gòu)，若函數(shù)模型中待求參數(shù)一部分變量是線性的，一部分變量是非線性的，并且函數(shù)模型為非線性函數(shù)的線性組合形式，數(shù)學上稱為可分離非線性最小二乘(separable nonlinear least squares，SNLS)問題[7-9]。測繪領(lǐng)域中存在許多這種模型結(jié)構(gòu)的參數(shù)估計問題，如激光雷達(light detection and ranging，LiDAR)全波形回波信號高斯函數(shù)分解模型中，高斯特征參數(shù)的求解[10-11]；空間直角坐標轉(zhuǎn)換模型中，旋轉(zhuǎn)參數(shù)、平移參數(shù)和尺度參數(shù)的求解[12-13]；神經(jīng)網(wǎng)絡(luò)預(yù)測模型中，激活函數(shù)的中心位置、寬度和權(quán)值參數(shù)的求解[14]等。

傳統(tǒng)求解可分離非線性最小二乘問題的方法是將所有的待求參數(shù)(非線性和線性)均看作非線性參數(shù)，然后進行泰勒級數(shù)展開，利用非線性迭代法或直接搜索法尋找問題的最優(yōu)解[15-16]。文獻[17]將LiDAR波形數(shù)據(jù)根據(jù)高斯分量的重要性進行篩選，然后利用非線性Levenberg-Marquardt(LM)優(yōu)化算法求解各波形分量非線性參數(shù)(波峰位置、波形寬度)和線性參數(shù)(振幅)的最優(yōu)值。文獻[18]對機載激光測深波形分解中LM與Expectation Maximization(EM)參數(shù)優(yōu)化方法進行了比較，就波形擬合效果而言，模擬數(shù)據(jù)分析和實測數(shù)據(jù)試驗均說明LM方法通?？傻玫奖菶M方法更接近原始波形的擬合結(jié)果，有利于進一步的波形分析，但受波形數(shù)據(jù)質(zhì)量的影響較大。文獻[19]針對三維坐標轉(zhuǎn)換模型參數(shù)求解中先驗精度未知和觀測數(shù)據(jù)質(zhì)量不佳的問題，以坐標轉(zhuǎn)換后的點位殘差加權(quán)平方和最小為原則，提出了基于Nelder-Mead單純形直接搜索的非線性參數(shù)求解方法，提高了三維坐標轉(zhuǎn)換參數(shù)的求解質(zhì)量。文獻[20]針對非線性最小二乘求解病態(tài)問題的不穩(wěn)定問題，提出了數(shù)值收斂解和收斂效率更優(yōu)的Frozen-Barycentre迭代法。文獻[21]針對參數(shù)迭代過程中雅可比矩陣計算復雜的問題，結(jié)合自動微分和隱式迭代預(yù)處理技術(shù)，采用最速下降法和信賴域法求解非線性最小二乘問題。文獻[22]依賴對精度水平的控制，當檢測到精度太低無法進行優(yōu)化時，對函數(shù)逼近方法進行改進，提出了一種基于動態(tài)精度評價函數(shù)和梯度函數(shù)的大規(guī)模非線性最小二乘的LM迭代算法，并證明了該算法的全局和局部收斂性。文獻[23]對基于變量投影算法參數(shù)分離的非線性最小二乘問題，提出采用LM優(yōu)化方法對非線性參數(shù)進行求解以及采用線性最小二乘或基于L-曲線的譜修正迭代方法對線性參數(shù)進行估計的混合算法，豐富了可分離非線性最小二乘的參數(shù)估計方法，但對基于變量投影算法參數(shù)分離的目標函數(shù)中殘差向量雅可比矩陣的計算方法并未進行深入研究。綜上分析可知，非線性參數(shù)估計的研究主要針對模型的構(gòu)建和參數(shù)優(yōu)化過程進行改進，取得了很好的結(jié)果[24-25]，但很少有考慮函數(shù)模型可分離的特殊結(jié)構(gòu)，并且將所有待求參數(shù)均作為非線性參數(shù)進行解算時，對參數(shù)初值的要求較高[26]。

本文從測繪領(lǐng)域中非線性函數(shù)模型的特殊結(jié)構(gòu)出發(fā)，首先，將線性參數(shù)通過變量投影算法用非線性函數(shù)表示，基于Moore-Penrose廣義逆矩陣微分和立體矩陣理論推導了殘差向量雅可比矩陣的計算過程；然后，對參數(shù)分離后的非線性最小二乘問題利用LM優(yōu)化方法進行求解，得到非線性參數(shù)的迭代方向和最優(yōu)估值，進而采用最小二乘方法計算線性參數(shù)的最優(yōu)解；最后，通過模擬試驗和LiDAR波形參數(shù)求解試驗驗證了基于Moore-Penrose廣義逆矩陣微分的可分離非線性最小二乘方法的有效性。

1 可分離非線性最小二乘法原理

設(shè)測量平差中的非線性函數(shù)模型和隨機模型為

(1)

式中，f(x)為非線性函數(shù)向量；x為待求參數(shù)向量；L為觀測值向量；σ2為單位權(quán)方差；矩陣Q和P分別為觀測向量協(xié)因數(shù)陣和權(quán)陣。在非線性最小二乘準則下建立目標函數(shù)為

(2)

(3)

在式(1)的非線性函數(shù)模型中，若待求參數(shù)x由線性參數(shù)向量a和非線性參數(shù)向量b兩部分組成，且fi(x)為非線性函數(shù)的線性組合，將其可以寫成以下形式

(4)

式中，p是線性參數(shù)a的維數(shù)；q是非線性參數(shù)b的維數(shù)；令n=p+q，n表示待求參數(shù)的總個數(shù)；m表示觀測次數(shù)。

令Φ(b)的列對應(yīng)的是與參數(shù)向量b相關(guān)的非線性函數(shù)φj(b),(j=1,2,…,p)。那么，非線性函數(shù)模型寫成矩陣形式為

V=f(a,b)-L=Φ(b)a-L

(5)

則式(5)的非線性最小二乘優(yōu)化問題為

(6)

若給定非線性參數(shù)b，則線性參數(shù)a的最小二乘最小范數(shù)解為

(7)

式中，Φ+(b)是Φ(b)的Moore-Penrose廣義逆[29]。

(8)

則原問題轉(zhuǎn)化為僅含有非線性參數(shù)的最小二乘問題?；谑?8)的算法也稱為變量投影算法[9]。

經(jīng)參數(shù)分離后的平差函數(shù)模型為

V(b)=Φ(b)Φ+(b)L-L

(9)

2 基于Moore-Penrose廣義逆矩陣微分的變量投影算法

根據(jù)最小二乘準則，待求參數(shù)的最優(yōu)值通過求極值方法得到[30]，其可分離非線性最小二乘的目標函數(shù)F(b)的一階導數(shù)為

(10)

(11)

又由于PΦ(b)Φ(b)=Φ(b)，因此

(12)

將式(12)代入式(11)，得變量投影算子的一階偏導為

(13)

(14)

(15)

根據(jù)立體矩陣的計算規(guī)則[31-32]，式(13)為

(16)

最終式(10)的目標函數(shù)一階導數(shù)為

F′(b)=

(17)

則殘差向量的Jacobian矩陣J(b)為

(18)

對目標函數(shù)F(b)用泰勒級數(shù)公式展開，并取至二階項，由以下優(yōu)化模型來得到迭代方向

(19)

由極值條件可得

(20)

式中

(21)

因此，式(20)迭代方向dk為

(22)

設(shè)第k次迭代的非線性參數(shù)向量為bk，從而基于變量投影的可分離非線性最小二乘LM算法的迭代方向dk為

(23)

式中，μk>0，J(bk)是殘差向量V(bk)(式(9))關(guān)于bk的一階偏導構(gòu)成的矩陣。

(24)

殘差向量的雅可比矩陣J(b)的每一行由相應(yīng)向量函數(shù)的梯度轉(zhuǎn)置得到，即

(25)

直接由式(13)可得

(26)

然后計算新的迭代點

JT(bk)V(bk)

(27)

(28)

令系數(shù)矩陣

(29)

(30)

結(jié)合以上參數(shù)求解過程，基于變量投影算法參數(shù)分離的LM優(yōu)化方法步驟如下。

(1)給定非線性參數(shù)初值b0以及閾值0<ε?1，并置迭代次數(shù)k=0。

(2)計算非線性函數(shù)矩陣Φ(bk)和雅可比矩陣J(bk)。

(3)根據(jù)式(23)和式(27)計算可分離非線性參數(shù)的迭代方向dk和第k+1次的參數(shù)值bk+1。

當非線性函數(shù)矩陣Φ(bk)秩虧或病態(tài)時，則需要通過解算不適定問題的方法進行求解[35-36]，參數(shù)收斂結(jié)果也會有所不同。

3 試驗分析

為了驗證基于Moore-Penrose廣義逆矩陣微分的可分離非線性最小二乘方法在參數(shù)估計中的有效性，采用指數(shù)函數(shù)擬合的模擬試驗及LiDAR全波形分解試驗與傳統(tǒng)的參數(shù)不分離的非線性最小二乘LM優(yōu)化方法進行了對比，試驗環(huán)境為Matlab 2016b，2.80 GHz PC，Windows10系統(tǒng)。

3.1 指數(shù)函數(shù)模型的擬合

模擬試驗源于放射性物質(zhì)的衰變過程，其函數(shù)模型表達式為

f(β,λ)=β1exp(-λ1ti)+β2exp(-λ2ti)+ξ

(31)

為了獲得待求參數(shù)的最優(yōu)估計，首先根據(jù)實際觀測值列出誤差方程，并通過變量投影算法將線性參數(shù)進行消除，得到基于變量投影算法的誤差方程，然后基于Moore-Penrose廣義逆矩陣的微分和立體矩陣理論計算目標函數(shù)的雅可比矩陣，采用LM算法對非線性參數(shù)λ=[λ1λ2]T進行估計，當目標函數(shù)(殘差平方和)的梯度小于10-6時，迭代終止。

結(jié)合衰變模型的結(jié)構(gòu)特點，根據(jù)參數(shù)初值選擇在真值附近或初值選擇常用的0和1得到4組不同的參數(shù)初值為[3 2 10 1]T、[0 1 10 1]T、[3 2 0 1]T和[0 1 0 1]T，然后采用本文提出的參數(shù)分離的LM優(yōu)化方法對非線性參數(shù)進行優(yōu)化，與參數(shù)不分離的LM優(yōu)化方法進行結(jié)果對比。給定4組參數(shù)初值，參數(shù)分離方法和參數(shù)不分離方法得到的衰減因子估值均為[3.30 2.68 9.54 1.32]T。由該參數(shù)估值得到的擬合曲線如圖1所示。

圖1 參數(shù)不分離方法和參數(shù)分離方法的擬合曲線對比Fig.1 Comparison of fitted curves between parameter non-separation method and parameter separation method

由圖1可知，參數(shù)分離方法和參數(shù)不分離方法得到的指數(shù)函數(shù)模型曲線與觀測值擬合較好。

從迭代次數(shù)、殘差平方和、均方根誤差、擬合優(yōu)度、相關(guān)系數(shù)和最大差值等評價指標對參數(shù)分離的LM優(yōu)化方法得到的結(jié)果進行分析，為了減少計算機性能等影響因素對時間的影響，將線性參數(shù)和非線性參數(shù)初值均在真值附近的試驗中參數(shù)不分離方法的計算時間設(shè)為單位1，其余試驗的時間對其進行比值，評價指標結(jié)果見表1。

表1 指數(shù)模型擬合試驗參數(shù)不分離LM方法和參數(shù)分離LM方法結(jié)果的評價指標對比Tab.1 Comparison of evaluation indexes between LM method without parameter separation and LM method with parameter separation in exponential model fitting experiment

結(jié)合參數(shù)最終估值結(jié)果和表2可知，4組試驗的評價指標除迭代次數(shù)和迭代時間外，其余指標結(jié)果相同，并且從擬合優(yōu)度和相關(guān)系數(shù)的數(shù)值可以看出，參數(shù)分離方法和參數(shù)不分離方法都能與觀測值擬合較好；迭代次數(shù)方面，當目標函數(shù)(殘差平方和)的梯度小于10-6時，由于第1組和第2組的非線性參數(shù)初值相同，參數(shù)分離的LM優(yōu)化方法不受線性參數(shù)初值的影響，因此迭代次數(shù)均為7，參數(shù)不分離方法由于線性參數(shù)的初值不同，迭代次數(shù)分別為23次和201次；第3組和第4組的非線性參數(shù)初值相同，參數(shù)分離的LM優(yōu)化方法的迭代次數(shù)均為15，參數(shù)不分離方法的迭代次數(shù)分別為1808次和1794次；迭代時間方面，第1組試驗中雖然參數(shù)分離的LM優(yōu)化方法迭代次數(shù)少，但是由于迭代過程中目標函數(shù)的一階導數(shù)計算復雜，導致計算時間比參數(shù)不分離的LM優(yōu)化方法多了約23%，但第2組、第3組和第4組試驗中，參數(shù)分離方法的計算時間均比參數(shù)不分離方法少，特別是第3組和第4組試驗，由于參數(shù)分離的LM優(yōu)化方法比參數(shù)不分離的LM優(yōu)化方法的迭代次數(shù)分別減少了1793次和1779次，因此參數(shù)分離方法的計算時間明顯減少，計算效率得到顯著提高。另外，第1組和第2組試驗參數(shù)分離方法的時間理論上應(yīng)該完全相同，但是由于計算機實時狀態(tài)的不同，所以計算時間有少許差別，第3組和第4組試驗的計算時間也是如此。

為了驗證參數(shù)分離方法對原問題病態(tài)性的改進，對基于變量投影算法參數(shù)分離的殘差向量雅可比矩陣的條件數(shù)進行分析，其4組初值對應(yīng)的參數(shù)分離方法和參數(shù)不分離方法的條件數(shù)結(jié)果見表2。

表2 不同參數(shù)初值情況下殘差向量雅克比矩陣條件數(shù)對比Tab.2 Comparison of condition number of Jacobian matrix of the residual vector with different initial parameters

由表2可知，同一組參數(shù)初值的參數(shù)分離方法雅克比矩陣的條件數(shù)比參數(shù)不分離方法更小，特別是第2組和第4組參數(shù)初值，原非線性最小二乘問題的雅可比矩陣的條件數(shù)均超過1018，存在明顯病態(tài)問題，而基于變量投影參數(shù)分離方法的雅可比矩陣的條件數(shù)分別為16.66和5.23，原問題的病態(tài)性得到明顯改善。

3.2 基于可分離非線性最小二乘的LiDAR波形分解

試驗采用海南某部分海域2012年12月測量得到的機載LiDAR測深數(shù)據(jù)進行波形分解，其波形數(shù)據(jù)采用Optech Aquarius機載測深系統(tǒng)測量得到，激光脈沖頻率為70 kHz，以1～2 GHz的高采樣率記錄各采樣點的振幅，采樣間隔為1 ns，采樣數(shù)量為287。由于在反射率低的情況下，電壓值會很小，因此，將觀測值進行量化，得到數(shù)字信號值(digital number)，即每一條波形由(ti,DNi),i=1,2,…，287組成。利用雙高斯模型對海面回波和海底回波建立模型，即

i=1,2,…,m

(32)

式中，a=[a0a1a2]T，b=[μ1μ2σ1σ2]T是待求線性參數(shù)和非線性參數(shù)；ξ是測量誤差。

圖2 不同參數(shù)初值的參數(shù)不分離方法和參數(shù)分離方法的LiDAR波形擬合曲線Fig.2 Comparison of LiDAR waveform fitted curves between parameter non-separation method and parameter separation method with different initial values of parameters

由圖2可知，參數(shù)分離的LM優(yōu)化方法根據(jù)第1、2和3組的參數(shù)初值得到LiDAR波形曲線與原始觀測值擬合較好；參數(shù)不分離的LM優(yōu)化方法第1組和第2組試驗得到的LiDAR波形曲線與參數(shù)分離方法相比，在第1個波峰處相差較大，第2個波峰擬合較好，第3組試驗得到的LiDAR波形在兩個波峰處與參數(shù)分離方法得到的擬合曲線都存在差值；第4組試驗無論是參數(shù)分離方法還是參數(shù)不分離方法都未能與觀測值進行較好擬合，也說明了兩種方法都對參數(shù)初值具有一定的依賴性。

從迭代次數(shù)、計算時間，殘差平方和、均方根誤差、擬合優(yōu)度、相關(guān)系數(shù)和最大差值等評價指標與參數(shù)不分離的非線性最小二乘LM優(yōu)化方法進行對比，與模擬試驗一樣，并將第1組參數(shù)初值在近似值附近的參數(shù)不分離方法的計算時間設(shè)為單位1，結(jié)果見表3。

表3中加粗的值表示離近似值較遠的初值。由表3可知，當線性參數(shù)初值和非線性參數(shù)初值都選擇在近似值附近時(第1組)，參數(shù)分離方法和參數(shù)不分離方法的迭代次數(shù)都是10，迭代終止時，參數(shù)不分離方法的殘差平方和是6.31×104，參數(shù)分離方法的殘差平方和是6.13×104，從均方根誤差、擬合優(yōu)度、相關(guān)系數(shù)和最大偏差4個評價指標也可以看出，參數(shù)分離方法得到的參數(shù)估計值精度更高，但是相同的迭代次數(shù)參數(shù)分離方法計算所需的時間更多。

表3 LiDAR波形擬合試驗參數(shù)不分離的LM方法和參數(shù)分離的LM方法結(jié)果的評價指標對比Tab.3 Comparison of evaluation indexes between LM method without parameter separation and LM method with parameter separation in LiDAR waveform fitting experiment

當非線性參數(shù)初值選擇在最優(yōu)值附近，線性參數(shù)初值任意選擇時(第2組)，參數(shù)分離方法由于只與非線性參數(shù)有關(guān)，所以與第1組試驗結(jié)果相同；參數(shù)不分離方法與第1組試驗相比，迭代次數(shù)增加1，均方根誤差由14.82增大為16.67，結(jié)合其余評價指標可知，參數(shù)分離方法的擬合精度更高，但所需計算時間更多。

當線性參數(shù)初值選擇在最優(yōu)值附近，非線性的部分參數(shù)(波形寬度)初值任意選擇時(第3組)，參數(shù)分離方法仍然可以得到與第1組試驗相同的結(jié)果；而參數(shù)不分離方法均方根誤差由第1組的14.82增加為21.52，擬合優(yōu)度也由0.979 3降至0.956 4，由此可知，參數(shù)不分離方法更容易受參數(shù)初值的影響。

當線性參數(shù)初值選擇在最優(yōu)值附近，非線性的部分參數(shù)(波峰位置)初值任意選擇時，參數(shù)不分離方法和參數(shù)分離方法都未得到與第1組精度一致的參數(shù)估值，結(jié)合圖2可知，參數(shù)分離方法和參數(shù)不分離方法都對波峰位置參數(shù)的初值比較敏感。

對上述4組試驗中參數(shù)分離方法和參數(shù)不分離方法的殘差平方和變化過程進行對比，結(jié)果如圖3所示。

由圖3可知，第1、2和3組在迭代終止時，參數(shù)分離方法的殘差平方和更小，也說明了參數(shù)分離方法在求解可分離非線性最小二乘問題時的參數(shù)估計精度更高，第4組試驗在迭代終止時雖然參數(shù)不分離方法的殘差平方和更小，但結(jié)合表2可知，兩種方法得到均方根誤差(99.433 0和102.286 4)都與參數(shù)初值在近似值附近(第1組)的結(jié)果(14.616 5)相差較大，說明兩種方法均未得收斂到待估參數(shù)的最優(yōu)值。

為了更全面地對比參數(shù)不分離方法和參數(shù)分離方法的收斂結(jié)果，根據(jù)LiDAR波形分解中常用的峰值探測法得到波峰位置參數(shù)和振幅參數(shù)的近似初值，波峰寬度參數(shù)的近似初值由峰值一半處對應(yīng)的兩點位置得到，然后以近似初值為中心，初值的一半為鄰域設(shè)定線性參數(shù)和非線性參數(shù)的取值區(qū)間，在均勻分布區(qū)間內(nèi)隨機取200組線性參數(shù)和非線性參數(shù)的初值，進行迭代建模，進而分析200組試驗中兩種方法的均方根誤差、擬合優(yōu)度、相關(guān)系數(shù)和最大偏差等指標的分布情況，結(jié)果如圖4和表4所示。

表4 效果評價指標的統(tǒng)計結(jié)果Tab.4 The statistical results of effect evaluation indexes

由圖4和表4可知，參數(shù)分離方法的均方根誤差均值比參數(shù)不分離方法小34.721 3，說明整體上參數(shù)分離方法的精度更高；200次試驗中，擬合優(yōu)度大于0.9的參數(shù)分離方法占74.5%，參數(shù)不分離方法占32.5%，相關(guān)系數(shù)大于0.9的參數(shù)分離方法占76.5%，參數(shù)不分離方法占35%；由兩種方法相關(guān)系數(shù)的標準差數(shù)值可知，參數(shù)分離方法200次試驗的相關(guān)系數(shù)標準差更小，說明其算法穩(wěn)定性更好。

圖4 LiDAR波形擬合200次試驗的效果評價指標分布情況Fig.4 The distribution of evaluation indexes in 200 experiments of LIDAR waveform fitting

對200次試驗中非線性參數(shù)初值相同的參數(shù)分離方法的殘差向量雅可比矩陣的條件數(shù)與參數(shù)不分離方法進行對比，結(jié)果如圖5所示。

圖5 200次試驗中殘差向量雅可比矩陣的條件數(shù)對比Fig.5 The comparison of the condition number of the Jacobian matrix of residual vector in 200 experiments

綜上分析可知：①參數(shù)分離的LM優(yōu)化方法只對非線性參數(shù)初值具有一定的依賴性，比傳統(tǒng)參數(shù)不分離的LM方法對待求參數(shù)初值的要求更低；②參數(shù)分離的LM優(yōu)化方法在對參數(shù)優(yōu)化過程中，所需的迭代次數(shù)更少，計算效率更高；③若兩種方法的迭代次數(shù)相同，由于參數(shù)分離的LM優(yōu)化方法在迭代過程中需要對Moore-Penrose廣義逆矩陣的微分和立體矩陣進行計算，增加了計算復雜度，因此需要的時間更多。

4 結(jié) 論

基于Moore-Penrose廣義逆矩陣微分的可分離非線性最小二乘方法通過利用變量投影算法將原來含有兩類參數(shù)的優(yōu)化問題轉(zhuǎn)化為僅含有非線性參數(shù)的最小二乘問題，降低了待求參數(shù)的維數(shù)，迭代過程只依賴非線性參數(shù)的初值，避免了線性參數(shù)初值導致的病態(tài)問題。機載LiDAR測深全波形數(shù)據(jù)擬合試驗結(jié)果表明，參數(shù)分離后的非線性最小二乘方法在波形參數(shù)求解過程中，參數(shù)估計精度更高，算法穩(wěn)定性也更強，為機載LiDAR波形分解參數(shù)的求解提供了新的方法，也拓展了可分離非線性最小二乘方法在測繪領(lǐng)域的應(yīng)用，但基于Moore-Penrose廣義逆矩陣的微分和立體矩陣的雅可比矩陣的顯式表達增加了計算復雜度，導致得到相同參數(shù)估計精度時，計算時間更多，因此如何提高計算效率也是可分離非線性最小二乘解算需要進一步研究的內(nèi)容。