紀(jì)元法, 朱亮亮, 孫希延, 嚴(yán)素清

(桂林電子科技大學(xué)廣西精密導(dǎo)航技術(shù)與應(yīng)用重點(diǎn)實(shí)驗(yàn)室, 廣西 桂林 541000)

0 引 言

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

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

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

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

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

對(duì)于觀測(cè)方程，即

(1)

(2)

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

2 Tikhonov正則化原理

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

min=‖AX-L‖2+αXTRX

(3)

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

(4)

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

3 自適應(yīng)加權(quán)的DE算法

3.1 基本DE算法原理

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

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

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

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

3.2 基本DE算法存在的問(wèn)題

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

3.3 自適應(yīng)加權(quán)的DE算法

對(duì)于變異式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è)體;g是種群進(jìn)化代數(shù);Ai,g是對(duì)個(gè)體xr1,g的加權(quán)因子;Fi,g是對(duì)差分向量xr2,g-xr3,g的加權(quán)因子;Ai,g和Fi,g分別為當(dāng)前個(gè)體和差分向量的權(quán)重。較小的Ai,g和較大的Fi,g使算法能在更大范圍內(nèi)探索更多有潛力的解;反之,較大的Ai,g和較小的Fi,g使算法的搜索范圍減小,更快地收斂到局部最優(yōu)解。

基于上述分析,提出自適應(yīng)改變加權(quán)因子Ai,g和Fi,g的DE算法改進(jìn)形式為

(6)

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

(7)

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

(8)

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

自適應(yīng)加權(quán)的DE算法具體步驟如下:

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

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

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

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

3.4 自適應(yīng)加權(quán)DE算法結(jié)合Tikhonov正則化

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

(9)

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

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

(10)

將正則化項(xiàng)與目標(biāo)函數(shù)直接相加,得到自適應(yīng)加權(quán)DE算法與Tikhonov正則化相結(jié)合的目標(biāo)函數(shù),即

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

(11)

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

4 算例與分析

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

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

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

0.648,21.943,1.892,9.665,12.171]T

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

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

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

4.1 解算精度分析

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

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

同樣利用GA算法對(duì)上述算例進(jìn)行優(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))作為解算精度的評(píng)價(jià)標(biāo)準(zhǔn),與新算法的結(jié)果對(duì)比,結(jié)果如表1所示。

表1 不同方案解算結(jié)果對(duì)比

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

4.2 收斂時(shí)間分析

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

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

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

4.3 穩(wěn)健性分析

為驗(yàn)證新算法解算病態(tài)方程的穩(wěn)健性,在得到的觀測(cè)向量基礎(chǔ)上人為附加部分粗差:第4、5、9、10個(gè)觀測(cè)值上附加20%的粗差,觀測(cè)向量變?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ù)不變的情況下,重新得到各方案的解算結(jié)果如表2所示。

表2 加入粗差后不同方案解算結(jié)果對(duì)比

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

5 結(jié) 論

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

參考文獻(xiàn)：

[1] 王振杰. 大地測(cè)量中不適定問(wèn)題的正則化解法研究[D]. 武漢:中國(guó)科學(xué)院測(cè)量與地球物理研究所, 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ìn)展[J]. 計(jì)算機(jī)應(yīng)用研究, 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] 歸慶明,韓松輝,隋立芬,等.抗差部分嶺估計(jì)及其在GPS快速定位中的應(yīng)用[J].大地測(cè)量與地球動(dòng)力學(xué),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.