紀元法, 朱亮亮, 孫希延, 嚴素清

(桂林電子科技大學廣西精密導航技術與應用重點實驗室, 廣西 桂林 541000)

0 引 言

病態(tài)問題廣泛存在于衛(wèi)星導航定位等工程測量中。在進行快速定位時,短時間內(nèi)只能觀測到幾個歷元,使得觀測信息不足,導致觀測方程病態(tài)。通過傳統(tǒng)的最小二乘法(least squares, LS)得到的模糊度浮點解與真值相差很大,嚴重影響定位精度[1]。解決觀測方程的病態(tài)性是快速提高定位精度的關鍵。

目前,針對病態(tài)問題的解算方法主要分為兩類:一類以降低法矩陣的條件數(shù)為目標,通過修改法矩陣的奇異值,降低其在求逆過程中對微小擾動的敏感性,從而得到更接近真值的結果。代表性的方法有嶺估計[2]、截斷奇異值分解法(truncated singular value decomposition,TSVD)[3]和Tikhonov正則化法[4]等;另一類是以遺傳算法(genetic algorithm,GA)為代表的智能優(yōu)化算法,這類方法不試圖改善法方程的病態(tài)性,而是通過最優(yōu)化方法得到適應度函數(shù)的近似全局最優(yōu)解,并將其作為病態(tài)方程的解[5]。兩類方法均存在不同程度的缺陷,前者中的TSVD方法直接舍棄小奇異值,在得到近似解的同時,損害了解估計的分辨率,影響解算效果[6],而嶺估計和Tikhonov正則化方法對正則化參數(shù)和正則化矩陣選取困難,且?guī)в兄饔^性[7];后者中的GA算法近乎遍歷的搜索使優(yōu)化速度過慢,且參數(shù)設置多,不利于實際工程應用[8]。

同時,若觀測值存在粗差,將會嚴重影響定位精度[9],已有的抗差部分嶺估計[10]、抗差嶺估計[11]等,基于抗差有偏估計的思想,有效地抑制了觀測值中的粗差,但解算精度不高。

本文提出差分進化(differential evolution,DE)算法結合Tikhonov正則化的病態(tài)方程解算方法。第一,變異過程中自適應地改變縮放因子,即自適應地改變當前個體在下一代新個體中所占的權重,從而在優(yōu)化的不同時期調(diào)節(jié)搜索范圍與尋優(yōu)速度之間的關系,提高病態(tài)方程解算精度。第二,結合Tikhonov正則化方法,將正則化項加入目標函數(shù)中,以減輕病態(tài)性,抑制噪聲和觀測誤差帶來的粗差影響[12],提高算法的穩(wěn)健性。

1 病態(tài)方程的數(shù)學模型

對于觀測方程，即

(1)

(2)

當觀測方程病態(tài)時,法矩陣ATA的條件數(shù)N很大,對其求逆不穩(wěn)定,此時若觀測向量中存在誤差,用最小二乘法直接求解會使誤差放大N倍[5],因此得到的結果與真值相差很大。

2 Tikhonov正則化原理

Tikhonov正則化方法是針對不適定問題提出來的[4],病態(tài)問題屬于不適定問題。對于式(1)的線性化觀測模型有

min=‖AX-L‖2+αXTRX

(3)

使式(3)最小化的參數(shù)X是所要求的式(1)的Tikhonov正則化解。其中α是正則化參數(shù),R是正則化矩陣,‖·‖2表示2范數(shù)。根據(jù)Tikhonov估計準則,可以取R=I,即正則化矩陣為單位陣。則式(3)的Tikhonov正則化解為

(4)

由式(4)可知,當正則化矩陣R確定后,Tikhonov正則化的關鍵是確定正則化參數(shù)α。相關學者已做了大量的研究,主要有嶺跡法、L-曲線法及直接確定參數(shù)法。其中L-曲線法是理論較為嚴密、也最有可能得到正確正則化參數(shù)的方法[1],因此本文選擇此方法。

3 自適應加權的DE算法

3.1 基本DE算法原理

DE算法由文獻[14]最先提出,具有易用性、尋優(yōu)速度快、穩(wěn)健性和強大的全局尋優(yōu)能力。其具體步驟如下:

步驟1設置種群大小NP、最大迭代次數(shù)MaxIter、加權因子F和交叉概率CR等參數(shù),在可行解空間內(nèi)對初始種群進行隨機初始化;

步驟2變異,種群內(nèi)個體的差分向量經(jīng)過加權,與種群內(nèi)相異的個體相加產(chǎn)生當前代變異個體vi,g=xr1,g+F×(xr2,g-xr3,g),其中i≠r1≠r2≠r3;

步驟5將種群中保留下來的最優(yōu)個體作為下一代的初始種群,進入循環(huán)迭代過程,重復步驟2～步驟4,直到滿足最大迭代次數(shù)。

3.2 基本DE算法存在的問題

盡管在非線性優(yōu)化問題上DE算法具有強穩(wěn)健性,但它不可避免地存在陷入局部最優(yōu)解和搜索停滯的問題[15]。若參數(shù)設置不當,例如加權因子F和交叉率CR過小,會使算法前期收斂過快,造成整個種群過早匯集于某一局部最優(yōu)點,此時無論變異個體還是交叉生成的新個體,均與原種群中的個體沒有顯著差異,這就陷入了局部最優(yōu)解;當種群沒有收斂到極值點且具有多樣性時,即使交叉變異后產(chǎn)生了新個體,也找不到比當前種群更優(yōu)的候選解,造成搜索停滯的問題。自適應加權的DE算法可通過在搜索的不同階段調(diào)節(jié)尋優(yōu)速度與搜索范圍的關系,改善陷入局部最優(yōu)解和搜索停滯的問題。

3.3 自適應加權的DE算法

對于變異式vi,g=xr1,g+F×(xr2,g-xr3,g),更一般地可以寫成

vi,g=Ai,g×xr1,g+Fi,g×(xr2,g-xr3,g)

(5)

式中,i和r代表種群中的某個體;g是種群進化代數(shù);Ai,g是對個體xr1,g的加權因子;Fi,g是對差分向量xr2,g-xr3,g的加權因子;Ai,g和Fi,g分別為當前個體和差分向量的權重。較小的Ai,g和較大的Fi,g使算法能在更大范圍內(nèi)探索更多有潛力的解;反之,較大的Ai,g和較小的Fi,g使算法的搜索范圍減小,更快地收斂到局部最優(yōu)解。

基于上述分析,提出自適應改變加權因子Ai,g和Fi,g的DE算法改進形式為

(6)

Fi,g=0.5×(1.5-pi,g)

(7)

式(6)中,i=[1,2,…,NP],NP是種群大小;g是進化代數(shù);ni是第i個個體未更新的次數(shù);max是個體可容許的最大未更新次數(shù)。當個體未更新次數(shù)達到max時,令Ai,g+1=0.1,從而在下一代的進化過程中強制更新此個體,解決了搜索停滯問題。由于當前個體xr1,g的適應度值越小,其加權因子Ai,g越大,且差分向量xr2,g-xr3,g的加權因子Fi,g越小,根據(jù)多次實驗得出式(7)中的經(jīng)驗值為0.5和1.5。pi,g的值為

(8)

式中,minF(xg)是第g代中最優(yōu)個體的適應度值。可以看出,pi,g和加權因子Fi,g都是種群個體適應度值F(xi,g)的函數(shù),加權因子Fi,g的值隨著個體適應度值的增大而增大,這是符合邏輯的,因為較大的適應度值,代表了種群中性能較差的個體。當個體性能不佳時,增大的加權因子Fi,g可以使算法在更大范圍內(nèi)搜索潛在的解;而當個體性能較好時,減小的加權因子Fi,g使變異產(chǎn)生的個體與原個體差別減小,算法快速收斂于全局最優(yōu)解附近的值。通過上述機制,提升DE算法的性能,避免了陷入局部最優(yōu)解和搜索停滯的問題。

自適應加權的DE算法具體步驟如下:

步驟1設置種群大小NP、最大迭代次數(shù)MaxIter、初始加權因子Fi,0、Ai,0和交叉概率CR等參數(shù),在可行解空間內(nèi)對初始種群進行隨機初始化,初始化盡量均勻,使整個種群在可行解范圍內(nèi)均勻分布。

步驟2變異,種群內(nèi)個體的差分向量經(jīng)過加權,與相異的個體相加產(chǎn)生當前代的變異個體vi,g=Ai,g×xr1,g+Fi,g×(xr2,g-xr3,g),其中當前代變異個體的當前個體加權因子為

Fi,g=0.5×(1.5-pi,g)

步驟5將種群中保留下來的最優(yōu)個體作為下一代的初始種群,進入循環(huán)迭代過程,重復步驟2～步驟5,直到滿足最大迭代次數(shù)。

3.4 自適應加權DE算法結合Tikhonov正則化

構造目標函數(shù)是DE算法的必要步驟,病態(tài)問題中,需要根據(jù)誤差方程構造目標函數(shù)。式(1)寫成誤差方程的形式為

(9)

目標函數(shù)可寫成如下形式:

minF(X)=(AX-L)T(AX-L)

(10)

將正則化項與目標函數(shù)直接相加,得到自適應加權DE算法與Tikhonov正則化相結合的目標函數(shù),即

minF(X)=(AX-L)T(AX-L)+αXTX

(11)

式中,α是正則化參數(shù),用第2節(jié)介紹的L-曲線法確定。

4 算例與分析

算例采用文獻[1]中的例4.2模擬病態(tài)問題,病態(tài)方程設計矩陣為

法矩陣N=ATA的條件數(shù)為1.289 2×105,方程嚴重病態(tài),待求向量有5個未知量,真值為Xtrue=[1,1,1,1,1]T,觀測噪聲Δ～N(0,σ2I),σ=1,根據(jù)觀測噪聲、系數(shù)陣A和真值Xtrue,隨機數(shù)發(fā)生器可產(chǎn)生觀測向量為

L=[-9.844,10.486,2.249,12.934,14.779,

0.648,21.943,1.892,9.665,12.171]T

設置種群大小NP=20,最大迭代次數(shù)MaxIter=500,初始加權因子Ai,g=Fi,g=0.5,交叉概率CR=0.9。待求參數(shù)搜索區(qū)間設為前3個未知量真值±1,后兩個未知量真值±5。式(11)基于誤差方程構建了新方法的適應度函數(shù),其正則化參數(shù)α由L-曲線法求得,如圖1所示。

圖1中,曲線上點(0.358 7,0.191 7)處的曲率最大,此點對應的α值為0.3,根據(jù)第2節(jié)的L-曲線法原理,所求的正則化參數(shù)α=0.3。

圖1 L-曲線法求得的正則化參數(shù)Fig.1 Tikhonov parameter by L-curve method

4.1 解算精度分析

由于智能優(yōu)化算法都是隨機搜索算法,不能保證每次搜索結果都相同,因此取100次優(yōu)化結果的平均值作為新算法的最終結果，即

XmDE=[1.332 2,0.823 5,1.058 7,0.557 4,1.124 2]T

同樣利用GA算法對上述算例進行優(yōu)化,得到100次優(yōu)化的均值為

Xga=[1.284 3,0.547 0,0.979 3,0.586 9,1.235 8]T

采用LS、TSVD、Tikhonov正則化法和DE算法分別解算本文算例,將歐式范數(shù)‖ΔX‖=sqrt((X-Xtrue)T(X-Xtrue))作為解算精度的評價標準,與新算法的結果對比,結果如表1所示。

表1 不同方案解算結果對比

根據(jù)表1,利用新算法得到的病態(tài)方程解的精度最優(yōu),與基本DE算法、GA算法、Tikhonov正則化法和TSVD法相比,其解算精度分別高約1倍、1.5倍、2倍和5倍。且LS求解精度最差,嚴重偏離真實值,這種情況下將不能得到接近真值的位置結果。

4.2 收斂時間分析

在快速定位解算病態(tài)方程時,已有的智能優(yōu)化算法沒有考慮運算時間,而快速定位對時間要求嚴格,是不得不考慮的問題。這里,先以算法的迭代次數(shù)為標準,假設每次迭代所用的時間相同,對算法優(yōu)化過程中的迭代次數(shù)進行統(tǒng)計?；贛atlab平臺,以第100次優(yōu)化為例,得到GA算法、基本DE算法及本文所提算法的優(yōu)化過程如圖2所示。

在程序中重新設置各算法的迭代次數(shù),要求比其收斂時間稍大,以滿足搜索到最優(yōu)解的要求。其中,GA算法的迭代次數(shù)為350,DE算法為100,新算法為25,分別統(tǒng)計它們的搜索時間,其結果為:新算法約0.122 1 s,DE算法約0.549 2 s,GA算法約1.837 5 s?？梢?新算法的收斂時間為基本DE算法的22.23%、GA算法的6.64%。因此,新算法的迭代次數(shù)最少,收斂速度最快,更能夠滿足快速定位中病態(tài)問題解算的要求。

圖2 不同算法的優(yōu)化過程Fig.2 Optimization process of different algorithms

4.3 穩(wěn)健性分析

為驗證新算法解算病態(tài)方程的穩(wěn)健性,在得到的觀測向量基礎上人為附加部分粗差:第4、5、9、10個觀測值上附加20%的粗差,觀測向量變?yōu)?/p>

L1=[-11.813,10.486,2.249,15.514,14.779,

0.648,21.943,2.704,11.599,12.171]T

保持原有參數(shù)不變的情況下,重新得到各方案的解算結果如表2所示。

表2 加入粗差后不同方案解算結果對比

對比表2與表1可知,加入粗差后,LS法的解算精度嚴重降低,TSVD法、GA法和DE法的解算精度均明顯降低,其‖ΔX‖值分別增加約0.55、0.54和0.49,新方法的解算精度仍保持最高,‖ΔX‖值增加約0.35。因此,新算法可在一定程度上克服粗差的影響,具有較好的穩(wěn)健性。

5 結 論

本文在基本DE算法的基礎上,提出自適應加權DE算法結合Tikhonov正則化解算病態(tài)問題的方法。通過仿真分析,無論與以降低方程病態(tài)性為目的的截斷奇異值法和Tikhonov正則化法相比,還是與以GA算法為代表的智能優(yōu)化算法相比,新算法都具有更高的解算精度;且新算法迭代次數(shù)最少,收斂時間最短,易實現(xiàn)快速定位;當觀測向量中混入粗差時,新算法的解算精度幾乎不變,仍保持最高,具有強穩(wěn)健性。另外,從算法過程來看,易應用于工程實際。

參考文獻：

[1] 王振杰. 大地測量中不適定問題的正則化解法研究[D]. 武漢:中國科學院測量與地球物理研究所, 2003.

WANG Z J. Research on the regularization solutions of ill-posed problems in geodesy[D]. Wuhan: Institute of Measurement and Geophysics, Chinese Academy of Sciences, 2003.

[2] ZHANG Y C, DUCHI J C, WAINWRIGHT M J. Divide and conquer kernel ridge regression: a distributed algorithm with minimax optimal rates[J].Journal of Machine Learning Research, 2013, 30(1):592-617.

[3] BOUHAMIDI A, JBILOU K, REICHEL L, et al. An extrapolated TSVD method for linear discrete ill-posed problems with Kronecker structure[J]. Linear Algebra & Its Applications,2011,434(7): 1677-1688.

[4] LIU C S. A dynamical Tikhonov regularization method for solving nonlinear ill-posed problems[J].Computer Modeling in Engineering & Sciences,2011,76(2): 109-132.

[5] GUO Q Y, HU ZH Q. Application of genetic algorithm to solve ill-posed equations for GPS rapid positioning[J]. Geomatics & Information Science of Wuhan University, 2009, 34(1):32-64.

[6] REICHEL L, RODRIGUEZ G. Old and new parameter choice rules for discrete ill-posed problems[M]. New York: Springer-Verlag, 2013.

[7] FAN Q, ZHANG N. Ill-posed problems robust solution of improved fruit fly optimization algorithm combining with Tikhonov regularization method[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(6):670-676.

[8] 馬永杰, 云文霞. 遺傳算法研究進展[J]. 計算機應用研究, 2012, 29(4):1201-1206.

MA Y J, YUN W X. Research progress of genetic algorithm[J]. Application Research of Computers, 2012, 29(4):1201-1206.

[9] WANG Q X, XU T H, XU G C. Application of combining method of outlier detection and robust estimation to GPS kinematic relative positioning[J]. Geomatics & Information Science of Wuhan University, 2011, 36(4):476-480.

[10] 歸慶明,韓松輝,隋立芬,等.抗差部分嶺估計及其在GPS快速定位中的應用[J].大地測量與地球動力學,2006,26(2):62-65.

GUI Q M, HAN S H, SUI L F, et al. Robust partial ordinary ridge estimator and its applications in GPS positioning[J]. Journal of Geodesy and Geodynamics, 2006, 26(2): 62-65.

[11] MARONNA R A. Robust ridge regression for high-dimensional data[J]. Technometrics, 2011, 53(1):44-53.

[12] HE R, HUANG S X, ZHOU C T, et al. Genetic algorithm with regularization method to retrieve ocean atmosphere duct[J]. Acta Physica Sinica, 2012, 61(4):273-335.

[13] XU Y B, PEI Y, DOND F. An extended L-curve method for choosing a regularization parameter in electrical resistance tomography[J]. Measurement Science & Technology, 2016, 27(11):114002.

[14] STORN R, PRICE K. Differential evolution-a simple and efficient heuristic for global optimization over continuous spaces[J]. Journal of Global Optimization, 1997, 11(4): 341-359.

[15] ZHANG H F, ZHOU J Z, ZHANG Y C, et al. Short term hydrothermal scheduling using multi-objective differential evolution with three chaotic sequences[J]. International Journal of Electrical Power & Energy Systems, 2013, 47(1): 85-99.