貝葉斯框架下基于凸優(yōu)化的系統(tǒng)偏差估計(jì)方法

周 林，劉先省，方擁軍，金 勇

(河南大學(xué)計(jì)算機(jī)與信息工程學(xué)院，河南 開(kāi)封 475004)

0 引言

多傳感器協(xié)同探測(cè)系統(tǒng)中，融合中心為了獲得被跟蹤目標(biāo)的“統(tǒng)一”信息(包括時(shí)間和空間上的“統(tǒng)一”)，必須首先要在各分布平臺(tái)上修正傳感器量測(cè)，即消除不“統(tǒng)一”的多傳感器量測(cè)誤差。量測(cè)誤差主要包含兩類：隨機(jī)誤差和系統(tǒng)偏差。前者可以通過(guò)濾波等方法消除，而后者則需利用配準(zhǔn)方法來(lái)改善[1]。系統(tǒng)偏差配準(zhǔn)的目的是對(duì)傳感器量測(cè)(諸如，雷達(dá)徑向距和方位角測(cè)量、時(shí)鐘誤差、傳感器位置等)偏差進(jìn)行估計(jì)和補(bǔ)償，以提高目標(biāo)定位和跟蹤精確性。

傳統(tǒng)上，研究系統(tǒng)偏差配準(zhǔn)問(wèn)題時(shí)通常假設(shè)系統(tǒng)偏差先驗(yàn)信息已知，即系統(tǒng)偏差在整個(gè)探測(cè)周期是某一未知固定值或呈現(xiàn)緩慢動(dòng)態(tài)演變，此時(shí)可利用最小二乘法、最大似然法、加權(quán)最小方差法、智能法、信息法、濾波等進(jìn)行解決[2-7]。但由于探測(cè)區(qū)域氣候、地形及照射光線各異、外來(lái)人為干擾增多、系統(tǒng)非線性和多模型等問(wèn)題，導(dǎo)致系統(tǒng)機(jī)動(dòng)以及系統(tǒng)偏差演化都難以建模，呈現(xiàn)隨機(jī)性、突發(fā)性等特性。傳統(tǒng)系統(tǒng)偏差估計(jì)方法不再適用解決隨機(jī)性的系統(tǒng)偏差配準(zhǔn)問(wèn)題。因此，考慮充分利用序貫量測(cè)信息，進(jìn)而準(zhǔn)確、實(shí)時(shí)估計(jì)隨機(jī)性系統(tǒng)偏差就尤為重要。

在利用多傳感器序貫量測(cè)信息，通過(guò)構(gòu)造目標(biāo)函數(shù)，進(jìn)而解決系統(tǒng)偏差估計(jì)的時(shí)候，無(wú)論目標(biāo)函數(shù)中的參數(shù)有無(wú)約束條件，擬牛頓法、共軛梯度法、可行方向法等典型策略都取得了較好效果[8-10]。但上述方法通常找到的是平穩(wěn)點(diǎn)或駐點(diǎn)(即局部最優(yōu)點(diǎn))，甚至找不到局部解，進(jìn)而也較難找到所有極小點(diǎn)集合中的最好點(diǎn)(即全局最優(yōu)點(diǎn))。盡管全局最優(yōu)算法的設(shè)計(jì)非必需，但它有助于解決優(yōu)化問(wèn)題，常見(jiàn)的優(yōu)化算法包括遺傳法、粒子群法、蟻群法、區(qū)間方法等[11-14]，它們具有全局最優(yōu)性和良好的約束條件適應(yīng)性，但實(shí)時(shí)性差、較難獲得易驗(yàn)證的全局最優(yōu)性條件。因此，考慮引入凸優(yōu)化方法來(lái)求解全局最優(yōu)解，該方法時(shí)間復(fù)雜度小，可用于可靠和迅速地優(yōu)化模型求解。

研究發(fā)現(xiàn)Metropolis-Hastings(簡(jiǎn)稱M-H)算法是一種近似貝葉斯方法[15-16]，屬于有約束的后驗(yàn)概率最大問(wèn)題，新樣本的選擇誤差會(huì)不斷增加，從而導(dǎo)致魯棒性得不到保證。因此，本文以文獻(xiàn)[16]中的系統(tǒng)偏差最大似然函數(shù)為基礎(chǔ)，構(gòu)造系統(tǒng)偏差估計(jì)的拉格朗日乘子二次函數(shù)并討論構(gòu)造函數(shù)的凸性，進(jìn)而提出基于凸優(yōu)化的系統(tǒng)偏差估計(jì)方法。

1 問(wèn)題的描述

正如前面所述，當(dāng)復(fù)雜環(huán)境或未知因素等出現(xiàn)時(shí)，導(dǎo)致系統(tǒng)以及系統(tǒng)偏差動(dòng)態(tài)模型難以構(gòu)建,傳統(tǒng)方法不再適用。假設(shè)某一傳感器探測(cè)系統(tǒng)中，考慮兩個(gè)問(wèn)題：1)機(jī)動(dòng)目標(biāo)的狀態(tài)難以建模；2)考慮傳感器量測(cè)受外界隨機(jī)干擾。因此，可只對(duì)極坐標(biāo)系下系統(tǒng)的非線性量測(cè)建模：

(1)

為討論方便，將每個(gè)傳感器極坐標(biāo)系下量測(cè)模型轉(zhuǎn)換為笛卡兒坐標(biāo)系下，其模型如下：

(2)

由于機(jī)動(dòng)目標(biāo)狀態(tài)模型未知，若想獲取估計(jì)系統(tǒng)偏差和目標(biāo)狀態(tài)信息，可利用zk={zk(i);i=1,…,n}及似然函數(shù)最大化方法估計(jì)傳感器系統(tǒng)偏差bk和目標(biāo)狀態(tài)xk，即：

(3)

對(duì)n部傳感器，假設(shè)量測(cè)具有白噪聲，進(jìn)而可得量測(cè)似然函數(shù)：

(4)

(5)

式(5)中,K2是常數(shù)K1的對(duì)數(shù)。忽略不相關(guān)常數(shù)項(xiàng)，通過(guò)優(yōu)化得出目標(biāo)最大似然估計(jì)：

清代乾隆年代[12]，江河商運(yùn)發(fā)展，博羅老城、惠東梁化、惠陽(yáng)的淡水、多祝又逐漸成為惠州區(qū)域主要的商業(yè)圩鎮(zhèn)，陸續(xù)發(fā)展興盛。宋、清惠州兩次人口的暴增[13]，導(dǎo)致當(dāng)?shù)孛耖g信仰祭祀場(chǎng)所對(duì)應(yīng)惠州人口的聚居點(diǎn)而增加，羅浮山道教也兩次借機(jī)復(fù)興。同時(shí)，惠州地貌的特征使得臺(tái)地、丘陵、平原階地沿東江相間分布，造就治水民間信仰隨人口聚居沿東江分布。

(6)

利用貝葉斯估計(jì)理論，得出多目標(biāo)最大似然估計(jì)函數(shù)如下[16]：

(7)

2 基于凸優(yōu)化的系統(tǒng)偏差估計(jì)方法

2.1 拉格朗日乘子二次函數(shù)

凸優(yōu)化模型可對(duì)一般非線性優(yōu)化模型進(jìn)行局部逼近，是研究非線性規(guī)劃問(wèn)題的主要途徑。本節(jié)將式(7)的系統(tǒng)偏差的最大似然估計(jì)函數(shù)轉(zhuǎn)化為具有目標(biāo)函數(shù)、約束函數(shù)的優(yōu)化問(wèn)題，其中目標(biāo)函數(shù)核心是構(gòu)造拉格朗日乘子二次函數(shù)形式，約束函數(shù)既包含求解參數(shù)空間邊界不等式約束條件也包含量測(cè)差形式的等式約束條件。同時(shí)，對(duì)上述優(yōu)化問(wèn)題的凸性進(jìn)行判斷分析。

解優(yōu)化問(wèn)題用到的目標(biāo)函數(shù)和約束函數(shù)即使是光滑的，有時(shí)候仍很難求解。為了更好描述后驗(yàn)概率函數(shù)并用優(yōu)化函數(shù)條更精確估計(jì)參數(shù)，目標(biāo)函數(shù)和約束函數(shù)之間的數(shù)學(xué)關(guān)系通過(guò)表達(dá)式描述，即用拉格朗日乘子二次函數(shù)表示，以第i個(gè)傳感器為例，其拉格朗日乘子二次函數(shù)為：

(8)

式(8)中，L表示系統(tǒng)偏差拉格朗日二次函數(shù)；γ表示包含待估系統(tǒng)偏差的表達(dá)式；λ1、λ2、λ3分別表示約束條件系數(shù)，它們根據(jù)實(shí)際問(wèn)題取值；f1、f2、f3分別表示約束條件。

傳統(tǒng)求解拉格朗日乘子二次函數(shù)最優(yōu)解的方法是用二次函數(shù)對(duì)待求參數(shù)求偏導(dǎo)并令偏導(dǎo)表達(dá)式等于零，進(jìn)而得到解析解。盡管該方法易于數(shù)學(xué)描述，但弊端是當(dāng)目標(biāo)函數(shù)非凸(或非凹)時(shí)，函數(shù)局部極值的存在會(huì)導(dǎo)致函數(shù)全局最優(yōu)解得不到確保。由于凸優(yōu)化模型是研究非線性規(guī)劃問(wèn)題的主要途徑，因此，可對(duì)上述拉格朗日乘子二次函數(shù)進(jìn)行凸性判斷。若函數(shù)呈現(xiàn)凸性，則其最優(yōu)解可利用凸優(yōu)化技術(shù)直接求解；若函數(shù)呈現(xiàn)非凸性，則可利用約束條件松弛、對(duì)偶等措施將非凸問(wèn)題轉(zhuǎn)換為凸問(wèn)題之后再求解。

2.2 二次函數(shù)的凸性判斷

判斷函數(shù)凸性的充要條件是函數(shù)的Hessian矩陣為半正定陣。求拉格朗日乘子二次函數(shù)L(γ,λ1,λ2,λ3)對(duì)系統(tǒng)偏差bk的Hessian矩陣，并將前述的A、γ等代入，可得如下Hessian矩陣表達(dá)式：

(9)

上述已用凸優(yōu)化方法證明拉格朗日乘子二次函數(shù)的凸性，結(jié)合拉格朗日乘子方法求解目標(biāo)函數(shù)，可將拉格朗日乘子二次函數(shù)最優(yōu)化問(wèn)題等價(jià)為有約束的凸優(yōu)化問(wèn)題，具體表達(dá)式如下：

(10)

L(γ,λ1,λ2,λ3)=

(11)

因此，根據(jù)拉格朗日乘子二次函數(shù)表達(dá)式，本文通過(guò)Matlab軟件附帶的專業(yè)凸優(yōu)化工具包來(lái)獲得每一個(gè)采樣時(shí)刻k使L最小的系統(tǒng)偏差bk。

2.3 算法描述

已有研究中，推導(dǎo)出的系統(tǒng)偏差最大似然估計(jì)函數(shù)，可利用M-H方法進(jìn)行尋優(yōu)求解，從而穩(wěn)態(tài)模擬參數(shù)貝葉斯估計(jì)方法[16]，但該方法存在兩個(gè)問(wèn)題：1)本質(zhì)上是一種迭代采樣方法，馬爾科夫鏈的鏈長(zhǎng)影響算法的實(shí)時(shí)性。2)通常較易找到局部最優(yōu)點(diǎn)而非全局最優(yōu)解。因此，為了提高估計(jì)實(shí)時(shí)性和全局最優(yōu)性，本文提出基于凸優(yōu)化的系統(tǒng)偏差估計(jì)方法，其流程如下圖1所示。

圖1 算法流程圖Fig.1 The flowchart of algorithm

3 仿真結(jié)果與分析

為了驗(yàn)證所提算法可行性和有效性，仿真假設(shè)用配置在異地的兩部傳感器探測(cè)和跟蹤機(jī)動(dòng)目標(biāo)，且機(jī)動(dòng)目標(biāo)在2D坐標(biāo)系統(tǒng)中的運(yùn)動(dòng)模型如下：

同時(shí)，兩部傳感器都設(shè)定為二維極坐標(biāo)量測(cè)方式，即量測(cè)包含徑向距和方位角，其量測(cè)模型如下：

zk(i)=hi(xk)+bk(i)+wk(i)i=1,2

(12)

式(12)中，傳感器量測(cè)zk(i)=[ρk(i)θk(i)]T；量測(cè)噪聲分別為σρ(1)=150 m，σθ(1)=0.5°；σρ(2)=80 m，σθ(2)=1°，非線性測(cè)量函數(shù)hi(·)如下：

(13)

式(13)中,i=1,2, 且(ζi,ηi)為笛卡兒坐標(biāo)系下傳感器位置， (ζ1,η1)=(300 km,-200 km)，(ζ2,η2)=(80 km,75 km)。

假設(shè)系統(tǒng)偏差bk具有隨機(jī)性(即動(dòng)態(tài)演化模型未知)，其估計(jì)初始b0=[0 0 0 0]T。為了驗(yàn)證所提算法的可行性和有效性，將本文算法與先前研究的M-H算法和ML算法進(jìn)行比較。為了描述系統(tǒng)偏差的隨機(jī)性，本文凸優(yōu)化算法的待估系統(tǒng)偏差的在徑向距、方位角和俯仰角范圍rrange(i)、θrange(i)的上限和下限分別取整個(gè)仿真周期系統(tǒng)偏差模型的最大值和最小值，使之體現(xiàn)參數(shù)空間的廣域性和可行性。為了與M-H算法進(jìn)行比較，設(shè)置的相關(guān)參數(shù)為：閾值δ=0.4，Markov鏈長(zhǎng)L=1 000，舍去預(yù)迭代值的個(gè)數(shù)為m=200；仿真采樣100次；Monte Carlo次數(shù)為50次。算法運(yùn)行硬件條件為Intel Core i7-2600 CPU 4 GHz。

3.1 參數(shù)估計(jì)性能仿真

圖2為傳感器A和傳感器B在徑向距和方位角上隨機(jī)性系統(tǒng)偏差真實(shí)值且有突變發(fā)生。為了說(shuō)明算法的有效性和可行性，圖3和圖4給出一次仿真得出的三種算法和真實(shí)系統(tǒng)偏差比較關(guān)系曲線圖。由于隨機(jī)性系統(tǒng)偏差難以建模，考慮充分利用傳感器序貫量測(cè)的最大似然方法(簡(jiǎn)稱ML)、Metropolis-Hasting(簡(jiǎn)稱M-H)及本文提出的凸優(yōu)化方法(簡(jiǎn)稱Convex)進(jìn)行對(duì)比。從圖3和圖4可以看出，相比于其他兩種方法，Convex方法在徑向距系統(tǒng)偏差估計(jì)效果較好；在方位角系統(tǒng)偏差估計(jì)效果，Convex方法要優(yōu)于M-H方法且能估計(jì)出突變量，但是要遜于ML方法。原因在于系統(tǒng)偏差具有隨機(jī)性和突變性，M-H和Convex兩種方法在采樣、學(xué)習(xí)過(guò)程中隨機(jī)時(shí)變地在參數(shù)空間內(nèi)游走，導(dǎo)致噪聲存在時(shí)變方差，而ML方法的噪聲滿足某一近乎恒定方差，這在需要分辨率較高的方位角仿真圖中能得到驗(yàn)證。

圖5和圖6為經(jīng)過(guò)50次Monte Carlo仿真后三種算法關(guān)于傳感器A和傳感器B關(guān)于徑向距、方位角的均方根誤差比較圖。從圖5和圖6可以看出，Convex方法在傳感器A和B徑向距的RMSE要好于M-H方法和ML方法；但傳感器A和B方位角的RMSE，Convex方法基本上處于M-H方法和ML方法之間，分析可知Convex方法在方位角估計(jì)精度并沒(méi)有定性于一定比M-H方法或ML方法中的任一種方法好。

圖2 傳感器A和B的系統(tǒng)偏差Fig.2 Systematic biases of the sensor A and B

圖3 傳感器A的系統(tǒng)偏差估計(jì)比較圖 Fig.3 Comparison of systematic bias estimation on sensor A

圖4 傳感器B的系統(tǒng)偏差估計(jì)比較圖Fig.4 Comparison of systematic bias estimation on sensor B

圖5 傳感器A的RMSE比較圖Fig.5 Comparison of RMSE on sensor A

圖6 傳感器B的RMSE比較圖Fig.6 Comparison of RMSE on sensor B

為了更直觀地描述整個(gè)仿真過(guò)程的估計(jì)性能，表1給出經(jīng)50次Monte Carlo仿真后三種算法估計(jì)出的傳感器A和B在徑向距和方位角系統(tǒng)偏差的平均均方誤差。

表1 三種算法的平均均方根誤差

3.2 復(fù)雜度分析

為了驗(yàn)證Convex方法的實(shí)時(shí)性，經(jīng)50次Monte Carlo仿真可知，三種算法的平均時(shí)間復(fù)雜度如表2所示。

分析原因，可依據(jù)ML、M-H及Convex算法的時(shí)間復(fù)雜度計(jì)算，其分別為O(n2·k)，O(L·n)，O(max{n3,n2·m,F})(其中，n表示系統(tǒng)狀態(tài)維數(shù)，k表示第k次采樣，L表示Markov鏈長(zhǎng)，m表示約束不等式約束條件數(shù)目，F(xiàn)表示對(duì)所有不等式約束條件及其一階、二階導(dǎo)的計(jì)算代價(jià))[17]。

表2 估計(jì)方法的平均運(yùn)行時(shí)間

綜合仿真曲線圖和數(shù)據(jù)可知，在參數(shù)空間尋優(yōu)求解范疇，Convex方法整體上比M-H方法耗時(shí)少、精度高；在整個(gè)仿真采樣時(shí)間段的估計(jì)精度范疇，Convex方法比ML方法精度高。

4 結(jié)論

在貝葉斯估計(jì)框架下，本文充分利用多傳感器的序貫量測(cè)信息，首先從量測(cè)最大似然函數(shù)入手，推算得到與狀態(tài)參數(shù)無(wú)關(guān)的多傳感器量測(cè)最大似然函數(shù)；然后將機(jī)動(dòng)目標(biāo)的系統(tǒng)模型從狀態(tài)空間投影到系統(tǒng)偏差空間，并將投影表達(dá)式作為系統(tǒng)偏差二次函數(shù)的等式約束條件；隨后，將待估的參數(shù)空間范圍作為不等式約束條件，并聯(lián)合等式約束條件，將最大似然估計(jì)問(wèn)題轉(zhuǎn)化為具有目標(biāo)函數(shù)、約束函數(shù)的凸優(yōu)化問(wèn)題；最后利用凸優(yōu)化工具包實(shí)現(xiàn)多傳感器系統(tǒng)偏差的優(yōu)化求解。本文新方法解決了當(dāng)隨機(jī)性的系統(tǒng)偏差難以建模時(shí)，傳統(tǒng)估計(jì)方法帶來(lái)的非全局最優(yōu)解、非實(shí)時(shí)性等問(wèn)題。