一種基于加權(quán)多維標(biāo)度分析的多個非相關(guān)源TDOA/FDOA協(xié)同定位方法

王鼎，尹潔昕，*，張欣光，鄭娜娥

1.中國人民解放軍戰(zhàn)略支援部隊信息工程大學(xué) 信息系統(tǒng)工程學(xué)院，鄭州 450001

2.國家數(shù)字交換系統(tǒng)工程技術(shù)研究中心，鄭州 450002

3.北京航天長征飛行器研究所，北京 100076

4.中國人民解放軍戰(zhàn)略支援部隊信息工程大學(xué) 數(shù)據(jù)與目標(biāo)工程學(xué)院，鄭州 450001

眾所周知，輻射源定位技術(shù)已廣泛應(yīng)用于無線通信、目標(biāo)監(jiān)測、航天航空、安全管理等諸多領(lǐng)域［1-3］，其在工業(yè)生產(chǎn)和國防安全中均發(fā)揮重要作用。對于靜止輻射源，常使用的定位觀測量為到達(dá)時間差（Time Difference of Arrival，TDOA）［4-6］。對于運動輻射源，多普勒效應(yīng)會產(chǎn)生到達(dá)頻率差（Frequency Difference of Arrival，F(xiàn)DOA）［7-8］，此時聯(lián)合TDOA和FDOA兩種觀測量進(jìn)行定位可獲得更高精度。近些年來，國內(nèi)外學(xué)者提出若干TDOA/FDOA定位方法，主要分為迭代定位方法和閉式定位方法兩大類。

迭代定位方法包括泰勒（Taylor）級數(shù)迭代方法［9］、半正定松弛方法［10］、約束總體最小二乘（Constrained Total Least Squares，CTLS）估計方法［11］、約束加權(quán)最小二乘（Constrained Weighted Least Squares，CWLS）估計方法［12］、迭代約束加權(quán)最小二乘（Iterative Constrained Weighted Least Squares， ICWLS）估計方法［13］、Bi迭代方法［14］等。雖然這些方法在一定條件下能獲得漸近統(tǒng)計最優(yōu)的估計精度，但需要復(fù)雜計算量，并且大多對迭代初始值較敏感，有些還需要設(shè)置迭代步長，易導(dǎo)致發(fā)散和局部收斂。

閉式定位方法能直接給出輻射源位置和速度顯式表達(dá)式，雖然其需要對加權(quán)矩陣進(jìn)行更新，但不用選擇迭代初始值和設(shè)置迭代步長，通常對加權(quán)矩陣更新3次就能滿足性能要求，并不存在發(fā)散和局部收斂問題，一般具有更低復(fù)雜度。最經(jīng)典閉式定位方法是兩步加權(quán)最小二乘（Two Step Weighted Least Squares，TSWLS）估計方法及其改進(jìn)型方法［15-17］，該類方法通過引入中間變量建立2組偽線性觀測方程，2次使用線性加權(quán)最小二乘估計器獲得輻射源位置和速度估計值。經(jīng)典閉式定位方法雖能有效克服迭代定位方法中的缺點，但該類方法需要推導(dǎo)偽線性觀測方程，因此其普適性略差于迭代定位方法，更重要的是其產(chǎn)生“門限效應(yīng)”的誤差閾值相對較小，在大觀測誤差條件下的定位精度會出現(xiàn)快速下降現(xiàn)象。

另一類改進(jìn)型閉式定位方法是加權(quán)多維標(biāo)度定位方法［18-25］。多維標(biāo)度分析是一種將多維空間對象簡化到低維空間進(jìn)行處理，同時還保留對象間原始關(guān)系的數(shù)據(jù)分析方法。近年來，國內(nèi)外學(xué)者提出若干基于加權(quán)多維標(biāo)度分析的定位方法。具體而言，文獻(xiàn)［18-20］提出基于到達(dá)時間（Time of Arrival，TOA）的加權(quán)多維標(biāo)度定位方法，文獻(xiàn)［21］提出基于接收信號強(qiáng)度（Received Signal Strength，RSS）的加權(quán)多維標(biāo)度定位方法，文獻(xiàn)［22-23］針對無線傳感網(wǎng)節(jié)點定位問題，提出相應(yīng)加權(quán)多維標(biāo)度定位方法，文獻(xiàn)［24］提出基于TDOA的加權(quán)多維標(biāo)度定位方法，文獻(xiàn)［25］針對運動輻射源提出基于TDOA/FDOA的加權(quán)多維標(biāo)度定位方法。從這些文獻(xiàn)中不難發(fā)現(xiàn)，加權(quán)多維標(biāo)度定位方法具備2個優(yōu)勢。第1個優(yōu)勢是可以給出輻射源位置和速度顯式表達(dá)式；第2個優(yōu)勢是相比其他定位方法，其對大觀測誤差具有更強(qiáng)魯棒性和穩(wěn)健性，這是因為該類方法能充分利用標(biāo)量積矩陣的維度信息和特征結(jié)構(gòu)信息，并且標(biāo)量積矩陣包含與距離及其變化率相關(guān)的全部信息。

上述定位方法均假設(shè)傳感器位置和速度精確已知，當(dāng)傳感器安裝在機(jī)載或艦載平臺，又或是傳感器隨機(jī)布設(shè)時，傳感器位置和速度精確值可能無法獲知，僅能得其先驗觀測值，其中含有先驗觀測誤差（也稱模型誤差）［26］。為了抑制模型誤差的影響，國內(nèi)外學(xué)者提出兩類定位方法。第1類方法是將模型誤差與TDOA/FDOA觀測誤差等同看待，通過對各種估計器重新設(shè)置加權(quán)矩陣以提高對模型誤差的魯棒性。例如，文獻(xiàn)［27］提出基于TDOA/FDOA的魯棒總體最小二乘（Total Least Squares，TLS）定位方法，文獻(xiàn)［28］提出基于TDOA/FDOA的魯棒CTLS定位方法，文獻(xiàn)［29］提出基于TDOA/FDOA的魯棒CWLS定位方法，文獻(xiàn)［26］提出基于TDOA/FDOA的魯棒TSWLS定位方法，文獻(xiàn)［30］提出基于TDOA/FDOA的魯棒加權(quán)多維標(biāo)度定位方法，文獻(xiàn)［31］提出基于TDOA/FDOA/DDR（Differential Doppler Rate）的魯棒TSWLS定位方法。這些方法均可有效減弱傳感器位置和速度先驗觀測誤差對定位性能的影響。第2類方法是對輻射源與傳感器位置和速度進(jìn)行聯(lián)合估計，該類方法不僅可以給出更準(zhǔn)確的傳感器位置和速度估計值，還能比第1類方法具有更高誤差閾值。例如，文獻(xiàn)［32-33］提出基于Taylor級數(shù)迭代的輻射源與傳感器位置和速度聯(lián)合估計方法。

前面提到的定位方法僅考慮對單個輻射源進(jìn)行定位，在實際場景中有時需要對多個輻射源進(jìn)行定位，此時應(yīng)考慮多源協(xié)同定位，以提高整體定位精度。當(dāng)傳感器位置和速度存在先驗觀測誤差時，即便對非相關(guān)輻射源而言，通過協(xié)同定位也可有效提升每個輻射源的定位精度。事實上，協(xié)同定位本身就可以看成是抑制模型誤差的有效途徑。文獻(xiàn)［34］提出基于TSWLS的多輻射源TDOA/FDOA協(xié)同定位方法，文獻(xiàn)［35］在文獻(xiàn)［34］的基礎(chǔ)上提出基于TSWLS的改進(jìn)型多輻射源TDOA/FDOA協(xié)同定位方法，該方法實現(xiàn)了對多個非相關(guān)輻射源以及各個傳感器位置和速度的聯(lián)合估計，其在大觀測誤差條件下具有更高的定位性能。

基于上述討論，本文在傳感器位置和速度存在先驗觀測誤差情況下，研究多個非相關(guān)輻射源TDOA/FDOA協(xié)同定位方法。為了提高在大觀測誤差條件下定位性能，文中首次提出一種基于加權(quán)多維標(biāo)度分析的多個非相關(guān)源TDOA/FDOA協(xié)同定位方法。該方法包含2個計算階段。階段a基于多維標(biāo)度分析原理構(gòu)造兩組標(biāo)量積矩陣，并由此形成定位關(guān)系式，用于獲得多個非相關(guān)輻射源以及各個傳感器位置和速度的估計值；階段b建立階段a估計誤差的約束優(yōu)化模型，并由此得到階段a估計誤差的估計值，用于更新階段a的估計結(jié)果。此外，文中還對新方法進(jìn)行理論性能分析，從數(shù)學(xué)上證明其對輻射源以及傳感器位置和速度估計精度均能漸近逼近克拉美羅界（Cramér-Rao Bound，CRB）。與現(xiàn)有加權(quán)多維標(biāo)度定位方法相比，新方法具有以下4個優(yōu)勢：① 實現(xiàn)多輻射源協(xié)同定位；② 實現(xiàn)輻射源與傳感器位置和速度聯(lián)合估計；③ 利用2個階段確保最終定位結(jié)果的漸近統(tǒng)計最優(yōu)性；④ 利用矩陣正交變換技術(shù)解決定位關(guān)系式中誤差協(xié)方差矩陣秩虧損問題。

這里給出文中使用的數(shù)學(xué)符號： ?表示矩陣Kronecker積；⊙表示矩陣點積；In表示n×n階單位矩陣，其中第m列向量記為；On×m表示n×m階全零矩陣；1n×m表示n×m階全1矩陣；=[On×1In]T；diag[·]表示由向量元素構(gòu)成的對角矩陣；blkdiag{·}表示由矩陣或向量作為對角元素構(gòu)成的塊狀對角矩陣；range{·}表示矩陣列空間；Λm-n表示滿足等式vec()=vec(Λm-nAnm)的0-1置換矩陣（其中Anm為任意n×m階矩陣）；Π[·]表示矩陣列空間的正交投影矩陣，Π⊥[·] 表示矩陣列補(bǔ)空間的正交投影矩陣；[·]?表示矩陣Moore-Penrose逆；A≥B表示A-B為半正定矩陣；＜a＞m表示向量a中的第m個元素。

1 TDOA/FDOA定位觀測模型

現(xiàn)有M個運動傳感器利用TDOA/FDOA觀測信息對N個非相關(guān)運動輻射源進(jìn)行定位。第m個傳感器位置和速度向量分別記為sm=并稱其為第m個傳感器位置-速度向量；第n個輻射源位置和速度向量分別為令并稱其為第n個輻射源位置-速度向量。在輻射源定位場景中，TDOA和FDOA可分別等價為距離差和距離差變化率，下文直接利用距離差和距離差變化率進(jìn)行建模和分析。

不失一般性，將第1個傳感器設(shè)為參考，將第n個輻射源與第m個傳感器之間的距離和距離變化率分別記為于是第n個輻射源與第m個和參考傳感器之間的距離差和距離差變化率可以分別表示為

實際中獲得的距離差和距離差變化率觀測量都含有誤差，它們可以分別表示為

式中：和表示距離差和距離差變化率觀測誤差。針對第n個運動輻射源，將式（2）中的2組等式寫成向量形式可得

現(xiàn)將式（3）中的2個等式合并可得

式中：

假設(shè)觀測誤差向量服從均值為零、協(xié)方差矩陣為的高斯分布。

雖然輻射源互不相關(guān)，但是不同輻射源對應(yīng)共同的傳感器位置-速度誤差，此時應(yīng)考慮對多個輻射源進(jìn)行協(xié)同定位，以獲得協(xié)同增益。為此需要將式（5）中的N個等式合并可得

式中：

觀測誤差向量服從均值為0、協(xié)方差矩陣為}的高斯分布。

在實際定位中，當(dāng)傳感器安裝在機(jī)載或艦載平臺，又或是傳感器隨機(jī)布設(shè)時，傳感器位置-速度向量可能無法精確獲得，僅能得到其先驗觀測值，即

注釋1在多源協(xié)同定位前需要首先估計輻射源個數(shù)，這屬于信號處理問題。如果是利用單天線（單通道）估計信號的TDOA/FDOA參數(shù)，則可以通過信號自相關(guān)矩陣的奇異值來確定信號個數(shù)［36］，或者構(gòu)造信號時頻圖像，并通過卷積神經(jīng)網(wǎng)絡(luò)檢測信號個數(shù)［37］。如果是基于陣列天線估計信號的TDOA/FDOA參數(shù)，則可以利用信息論準(zhǔn)則（例如AIC準(zhǔn)則、MDL準(zhǔn)則等）獲得信號個數(shù)［38-39］，或者利用深度學(xué)習(xí)技術(shù)確定信號個數(shù)［40］。

2 新方法的基本原理與計算步驟

2.1 新方法的總體描述

為了避免發(fā)散和局部收斂等問題，本節(jié)提出一種新的定位方法，該方法是由多維標(biāo)度分析衍生出，相比其他定位方法，其在大觀測誤差條件下具有更高定位性能。為獲得漸近統(tǒng)計最優(yōu)估計精度，新方法包含2個計算階段（記為階段a和階段b），每個階段的主要原理可描述如下：

1）階段a首先基于多維標(biāo)度分析構(gòu)造2組標(biāo)量積矩陣，并利用標(biāo)量積矩陣性質(zhì)構(gòu)建定位關(guān)系式，然后通過一階誤差分析推導(dǎo)定位關(guān)系式中的誤差協(xié)方差矩陣，并利用矩陣正交變換將其恢復(fù)為滿秩矩陣，從而確定加權(quán)矩陣，最后基于定位關(guān)系式和加權(quán)矩陣獲得多個非相關(guān)輻射源以及各個傳感器位置和速度估計值。

2）階段b首先基于多維標(biāo)度分析中引入的中間變量，構(gòu)建關(guān)于階段a定位誤差的約束優(yōu)化模型，然后基于此模型對階段a的定位誤差進(jìn)行估計，最后利用該估計值對階段a的定位結(jié)果進(jìn)行更新，以獲得漸近統(tǒng)計最優(yōu)的估計性能。

2.2 階段a的計算原理

2.2.1 構(gòu)造標(biāo)量積矩陣

在多維標(biāo)度分析中需要首先構(gòu)造標(biāo)量積矩陣，由于本文考慮多源協(xié)同定位，因此需要依次對每個輻射源構(gòu)造標(biāo)量積矩陣。首先針對第n個輻射源定義如下4維復(fù)坐標(biāo)和復(fù)速度向量

式中：j表示虛數(shù)單位；dn1=rn1-rn1=0和=-=0?；谑剑?1）可定義如下復(fù)坐標(biāo)和復(fù)速度矩陣

2.2.2 構(gòu)建定位關(guān)系式

下面利用標(biāo)量積矩陣構(gòu)建新的定位關(guān)系式，與其它常規(guī)定位方法中的關(guān)系式不同，這里的關(guān)系式能充分利用標(biāo)量積矩陣的維度信息和特征結(jié)構(gòu)信息，從而提高在大觀測誤差條件下的定位性能。首先由第1組標(biāo)量積矩陣{Wn}1≤n≤N可以得到如下結(jié)論。

命題1針對第n個輻射源，分別定義如下實向量vn和實矩陣Vn

證明首先由Moore-Penrose逆矩陣的性質(zhì)可得

利用矩陣和求逆公式可知

然后將式（19）代入式（18）中可得

另一方面，將式（15）和式（16）代入式（12）中的第1式中可知

最后聯(lián)合式（13）、式（20）以及式（21）可得

證畢。

若定義矩陣

則由命題1可知

將式（24）兩邊對時間求導(dǎo)可得

式中

最后將式（24）和式（25）合并可知

若忽略誤差二階項，則由式（28）可得

式中：ΔWn、、ΔTn以及分別表示、、以及中的誤差一階項。

利用一階誤差分析可以將式（29）右邊第1項表示為關(guān)于觀測誤差向量和ε(s)的線性函數(shù)，即

式中：

利用一階誤差分析可以將式（29）右邊第2項表示為關(guān)于觀測誤差向量ε(d)n和ε(s)的線性函數(shù)，即

利用一階誤差分析可以將式（29）右邊第3項表示為關(guān)于觀測誤差向量ε(d)n和ε(s)的線性函數(shù)，即

式中：

其中

利用一階誤差分析可以將式（29）右邊第4項表示為關(guān)于觀測誤差向量ε(d)n和ε(s)的線性函數(shù)，即

式中

其中

將式（30）、式（34）、式（38）以及式（42）代入式（29）中可得

式中：由式（46）可知，誤差向量ξn漸近服從零均值高斯分布，并且其協(xié)方差矩陣為

式中：Qn1和Qn2均為列正交矩陣；Rn為上三角滿秩矩陣。為了得到協(xié)方差矩陣為滿秩的誤差向量，可以將矩陣左乘以誤差向量ξn，并且結(jié)合式（28）、式（46）以及式（48）可知

由式（49）可知，誤差向量δn的協(xié)方差矩陣為

由于誤差向量δn的維數(shù)為2(M-1)，其與TDOA/FDOA觀測量個數(shù)相等，因此容易驗證COV(δn)為滿秩矩陣。

2.2.4 估計準(zhǔn)則及其最優(yōu)解

首先將矩陣和按列分塊表示為=由式（49）可將誤差向量δn重新寫為

式中：

從式（49）中可以看出，N個誤差向量{δn}1≤n≤N中包含共同的模型誤差向量ε(s)，為了進(jìn)行多輻射源協(xié)同定位，需要將式（51）中的N個等式合并，如下式所示

式中：

結(jié)合式（49）和式（53）可知，誤差向量δc漸近服從零均值高斯分布，并且其協(xié)方差矩陣為

式中：

為了提高發(fā)生門限效應(yīng)時的誤差閾值，這里考慮對向量和ˉ進(jìn)行聯(lián)合估計，為此需要定義擴(kuò)維參數(shù)向量和擴(kuò)維誤差向量然后結(jié)合式（10）和式（53）可構(gòu)建階段a中的估計準(zhǔn)則，即

式中：

式（57）中的(COV(ζ))-1可以看成是加權(quán)矩陣，其作用在于抑制觀測誤差向量和模型誤差向量ε(s)的影響，其最優(yōu)解為

若將估計值（或者{}1≤m≤M）中的誤差向量記為（或者{}1≤n≤N）以及則有

2.2.5 理論性能分析

階段b是對階段a的估計結(jié)果進(jìn)行更新，其需要階段a估計結(jié)果的均方誤差矩陣，因此下面將對階段a的估計結(jié)果進(jìn)行理論性能分析。首先由線性加權(quán)最小二乘估計理論可知，誤差向量Δρa(bǔ)漸近服從零均值高斯分布，并且估計值的均方誤差矩陣為

式中：Pc=blkdiag{P1，P2， … ，PN}，其中

接著利用分塊矩陣求逆公式可知

將式（64）代入式（62）中可得

式中：

結(jié)合式（61）和式（65）可知，誤差向量漸近服從零均值高斯分布，并且估計值的均方誤差矩陣為

式中需要指出的是，雖然向量是漸近無偏估計值，但尚不是漸近統(tǒng)計最優(yōu)估計值，具體可見如下命題。

命題2若將未知參數(shù)的估計均方誤差的CRB記為則在一階誤差分析理論框架下滿足

命題2的證明見附錄A。階段a的估計值不具備漸近統(tǒng)計最優(yōu)性的根本原因在于，中間變量的引入導(dǎo)致向量ρ中的元素間存在相關(guān)性，這種相關(guān)性會使得誤差向量Δρa(bǔ)服從某個等式約束。階段b將利用該等式約束進(jìn)一步提高定位精度。

2.3 階段b的計算原理

首先推導(dǎo)誤差向量Δρa(bǔ)服從的等式約束。根據(jù)向量ρ第8n-4個元素和第8n個元素的定義可知

分別將式（68）和式（69）右邊第1項在點處進(jìn)行一階Taylor級數(shù)展開可得

式中：

將式（70）和式（71）中的2N個等式合并，可以得到如下關(guān)于誤差向量Δρa(bǔ)的線性等式

式中：

式（73）即為誤差向量Δρa(bǔ)滿足的等式約束。根據(jù)誤差向量Δρa(bǔ)服從的高斯分布特性，可以建立估計誤差向量Δρa(bǔ)的約束優(yōu)化模型，即

利用拉格朗日乘子技術(shù)可知［46］，式（75）的最優(yōu)解為

結(jié)合式（61）和式（76）可以得到誤差向量和的估計值分別為

于是輻射源位置-速度向量和傳感器位置-速度向量ˉ在階段b的估計值分別為

2.4 新方法的計算步驟與討論

綜合前面的討論，總結(jié)新方法的計算步驟，如圖1所示。針對上述新方法給出3點注釋：

圖1 本文新方法的計算流程圖
Fig.1 Flow chart of method proposed in this paper

3 新方法的理論性能分析

本節(jié)推導(dǎo)估計值的統(tǒng)計特性，主要推導(dǎo)其估計均方誤差，并與相應(yīng)的CRB進(jìn)行比較，從而證明其漸近統(tǒng)計最優(yōu)性，為此首先給出多個非相關(guān)源協(xié)同定位的CRB的表達(dá)式。

3.1 多個非相關(guān)源協(xié)同定位的CRB

基于文獻(xiàn)［34-35］中的結(jié)論可知，多個非相關(guān)源協(xié)同定位的CRB矩陣可以表示為

式中：均表示函數(shù)fc(ˉ)的Jacobian矩陣，它們可以分別表示為

其中：均表示函數(shù)f(ˉ)的Jacobian矩陣。由式（79）和式（80）可知，多個非相關(guān)源協(xié)同定位的CRB與輻射源位置緊密相關(guān)。

3.2 估計均方誤差

首先將式（75）中的等式約束代入式（76）中可得

式中由式（81）可知，估計值的均方誤差矩陣為

結(jié)合式（77）和式（82）可知，估計值的均方誤差矩陣為

3.3 漸近統(tǒng)計最優(yōu)性分析

下面證明估計值具有漸近統(tǒng)計最優(yōu)性，為此需要給出均方誤差矩陣的另一種表達(dá)式，具體可見如下命題。

命題3均方誤差矩陣可以表示為

證明：首先由式（83）可得

考慮矩陣這2個矩陣的列數(shù)之和等于8N+6M，并且滿足

于是有

由此可知

將式（88）代入式（85）中可得

不難驗證，將該式代入式（89）中可知式（84）成立。證畢。

命題4在一階誤差分析理論框架下，滿足

證明：首先根據(jù)向量ρ的定義可得

將式（90）代入式（84）中，并且利用式（65）可知

結(jié)合式（66）和式（A8）可得

將式（92）～式（94）代入式（91）中，并且結(jié)合式（79）可知命題4成立。證畢。

4 仿真實驗

本節(jié)進(jìn)行仿真實驗，基礎(chǔ)實驗條件為：利用7個運動傳感器獲得TDOA/FDOA信息（亦即距離差/距離差變化率信息）對運動輻射源進(jìn)行定位，傳感器位置和速度見表1，距離差/距離差變化率觀測誤差向量服從均值為零、協(xié)方差矩陣為的高斯分布，傳感器位置-速度先驗觀測誤差向量ε(s)服從均值為零、協(xié)方差矩陣為E(s)=σ22(IM?blkdiag{I3，I3/10})的高斯分布。

表1 傳感器位置和速度數(shù)值
Table 1 Numerical value of sensor position and velocity

傳感器序號1234567 x(s)m/m 2 400 1 700 2 300-1 600 2 500-1 100-2 200 y(s)m/m 1 800 2 100-1 900 1 700-1 600-1 400 1 500 z(s)m/m 1 300-2 200 1 600 1 500-1 200 2 400-2 500 x?(s)m/（m·s-1）-10 10-18-10-10 18 15 y?(s)m/（m·s-1）16-15-10 10 18-15-10 z?(s)m/（m·s-1）15 16 18-18 15-16 16

仿真實驗1假設(shè)有2個非相關(guān)輻射源，輻射源1位置和速度向量分別為u1=[-3 300，-3 600，3 800]Tm和=[20，10，-5]Tm/s，輻射源2位置和速度向量分別為u2=[7 200，7 500，-7 700]Tm和u?2=[-10，15，15]Tm/s，令20lg(σ1)=20lg(σ2)=-5 dBm，利用本文方法對這2個輻射源進(jìn)行協(xié)同定位，并且進(jìn)行5 000次蒙特卡羅實驗，圖2和圖3給出了輻射源1和輻射源2的定位結(jié)果散點圖與定位誤差橢圓曲線。

圖2 輻射源1的定位結(jié)果散點圖與定位誤差橢圓曲線（x-y平面）
Fig.2 Scatter plot of positioning results and elliptic curve of positioning error for Emitter 1 （x-yplane）

圖3 輻射源2的定位結(jié)果散點圖與定位誤差橢圓曲線（y-z平面）
Fig.3 Scatter plot of positioning results and elliptic curve of positioning error for Emitter 2 （y-zplane）

從圖2和圖3中可以看出，本文方法的定位結(jié)果散點圖的形狀與定位誤差橢圓的形狀一致，并且大概率對應(yīng)大橢圓，小概率對應(yīng)小橢圓，從而驗證了該方法的有效性。

仿真實驗2假設(shè)有3個非相關(guān)輻射源，輻射源1位置和速度向量分別為u1=[4 600，5 600，5 800]Tm和=[10，-25，20]Tm/s，輻射源2位置和速度向量分別為u2=[6 200，7 200，7 400]Tm和=[-20，15，10]Tm/s，輻射源3位置和速度向量分別為u3=[8 800，9 600，9 200]Tm和=[15，-10，-10]Tm/s，下面利用本文方法對這3個輻射源進(jìn)行協(xié)同定位，并將其定位精度與一些主流TDOA/FDOA定位方法進(jìn)行比較，其中包括文獻(xiàn)［27］中的TLS方法、文獻(xiàn)［31］中的TSWLS方法、文獻(xiàn)［34］中的TSWLS方法、文獻(xiàn)［35］中的改進(jìn)型TSWLS方法以及文獻(xiàn)［13］中的ICWLS方法，其中前面4種是閉式定位方法，最后1種是迭代定位方法。需要指出的是，雖然文獻(xiàn)［27］中的TLS方法、文獻(xiàn)［31］中的TSWLS方法以及文獻(xiàn)［13］中的ICWLS方法都是在單輻射源存在的場景下提出的，但卻可以將它們推廣至多輻射源協(xié)同定位的場景。此外，文獻(xiàn)［31］中的TSWLS方法聯(lián)合利用了DDR信息進(jìn)行定位，但為了公平比較，這里的仿真并未利用此信息。另一方面，選擇ICWLS方法作為比較對象是因為它是一種性能優(yōu)越的迭代定位方法［13］，為了進(jìn)行更全面的比較，這里給出其初始值為隨機(jī)值和真實值2種情形下的定位性能。首先令20lg(σ2)=0 dBm，改變σ1的數(shù)值，圖4給出了多輻射源位置和速度估計均方根誤差隨20lg(σ1)的變化曲線；圖5給出了多傳感器位置和速度估計均方根誤差隨20lg(σ1)的變化曲線。然后令20lg(σ1)=-5 dBm，改變σ2的數(shù)值，圖6給出了多輻射源位置和速度估計均方根誤差隨20lg(σ2)的變化曲線；圖7給出了多傳感器位置和速度估計均方根誤差隨20lg(σ2)的變化曲線。

圖4 多輻射源位置和速度估計均方根誤差隨20lg（σ1）的變化曲線
Fig.4 RMSE of multiple-source position and velocity estimation versus 20lg（σ1）

圖5 多傳感器位置和速度估計均方根誤差隨20lg（σ1）的變化曲線
Fig.5 RMSE of multiple-sensor position and velocity estimation versus 20lg（σ1）

圖6 多輻射源位置和速度估計均方根誤差隨20lg（σ2）的變化曲線
Fig.6 RMSE of multiple-source position and velocity estimation versus 20lg（σ2）

圖7 多傳感器位置和速度估計均方根誤差隨20lg（σ2）的變化曲線
Fig.7 RMSE of multiple-sensor position and velocity estimation versus 20lg（σ2）

從圖4～圖7中可以看出：① 在大觀測誤差條件下本文方法的估計精度優(yōu)于TLS方法和3種TSWLS方法，新方法發(fā)生門限效應(yīng)的誤差閾值更高，這驗證了該方法相比其他閉式定位方法的優(yōu)越性；② 在大觀測誤差條件下，本文方法估計精度優(yōu)于初始值為隨機(jī)值的ICWLS方法，并且與初始值為真實值的ICWLS方法接近，注意到將真實值作為初始值在實際計算中難以實現(xiàn)，而本文方法無需初始值，這驗證了本文方法相比迭代定位方法的優(yōu)越性；③ 本文方法對輻射源定位性能漸近逼近相應(yīng)的CRB，這驗證了第3節(jié)理論性能分析的有效性；④ 本文方法可以提高傳感器位置和速度估計精度（相比其先驗觀測精度而言），并且對傳感器位置和速度估計均方根誤差均可以漸近逼近相應(yīng)的CRB，從而再次驗證了第3節(jié)理論性能分析的有效性。

仿真實驗3仍然考慮3個非相關(guān)輻射源存在的場景，輻射源1位置和速度向量為u1=[5cos(γ)cos(α)，5sin(γ)cos(α)，5sin(α)]Tkm和=[10，15，-10]Tm/s，輻射源2位置和速度向量為u2=[6cos(γ)cos(α)，6sin(γ)cos(α)，6sin(α)]Tkm和=[-10，5，-15]Tm/s，輻射源3位置和速度向量為u3=[7cos(γ)cos(α)，7sin(γ)cos(α)， 7sin(α)]Tkm和=[-10，15，15]Tm/s。由于本文方法可以通過簡化直接應(yīng)用于單輻射源存在場景（此時退化為非協(xié)同定位方法），于是下面不妨比較多輻射源協(xié)同定位和多輻射源非協(xié)同獨立定位2種方式下的定位精度，旨在驗證多源協(xié)同定位能產(chǎn)生協(xié)同增益。此外，還將文中方法與文獻(xiàn)［30］中的加權(quán)多維標(biāo)度定位方法進(jìn)行比較，這是因為該方法也考慮了傳感器位置和速度先驗觀測誤差，只是該方法是針對單輻射源的定位方法，目前尚沒有進(jìn)行多源TDOA/FDOA協(xié)同定位的加權(quán)多維標(biāo)度定位方法。另一方面，文獻(xiàn)［25］也是基于加權(quán)多維標(biāo)度的TDOA/FDOA定位方法，但是該方法沒有考慮傳感器位置和速度先驗觀測誤差，因此不適合與文中的新方法進(jìn)行性能比較。令20lg(σ1)=-10 dBm，20lg(σ2)=-5 dBm以及γ=12°，通過改變α的數(shù)值來改變輻射源位置。圖8給出了輻射源1的位置和速度估計均方根誤差隨α的變化曲線；圖9給出了輻射源2的位置和速度估計均方根誤差隨α的變化曲線；圖10給出了輻射源3的位置和速度估計均方根誤差隨α的變化曲線。

圖8 輻射源1的位置和速度估計均方根誤差隨α的變化曲線
Fig.8 RMSE of position and velocity estimation of Emitter 1 versusα

圖9 輻射源2的位置和速度估計均方根誤差隨α的變化曲線
Fig.9 RMSE of position and velocity estimation of Emitter 2 versusα

圖10 輻射源3的位置和速度估計均方根誤差隨α的變化曲線
Fig.10 RMSE of position and velocity estimation of Emitter 3 versusα

從圖8～圖10中可以看出：① 多源協(xié)同定位產(chǎn)生的性能增益是顯著的；② 本文方法（非協(xié)同定位版）的估計精度優(yōu)于文獻(xiàn)［30］中的加權(quán)多維標(biāo)度定位方法，這是因為后者尚未利用中間變量包含的有用信息；③ 隨著輻射源位置的改變，文獻(xiàn)［30］中的加權(quán)多維標(biāo)度定位方法的估計性能并不總能達(dá)到相應(yīng)的CRB，而本文方法包含2個計算階段，可以使其定位精度始終漸近逼近相應(yīng)的CRB。

仿真實驗4將本文方法與文獻(xiàn)［27］中的TLS方法、文獻(xiàn)［31］中的TSWLS方法、文獻(xiàn)［34］中的TSWLS方法、文獻(xiàn)［35］中的改進(jìn)型TSWLS方法以及文獻(xiàn)［13］中的ICWLS方法的運行時間進(jìn)行比較，以間接比較6種定位方法的計算復(fù)雜度。仿真軟件為MATLAB R2020a，仿真程序在安裝有i7-3520CPU的PC機(jī)上運行，令20lg(σ1)=10 dBm和20lg(σ2)=10 dBm，其余條件不變，并進(jìn)行5 000次蒙特卡羅實驗。表2給出了6種定位方法的平均運行時間。從表2中可以看出，新方法的計算復(fù)雜度高于4種閉式定位方法，但是低于ICWLS定位方法。

表2 6種定位方法的平均運行時間
Table 2 Average running time of six localization methods

定位方法本文方法文獻(xiàn)［27］方法文獻(xiàn)［31］方法文獻(xiàn)［34］方法文獻(xiàn)［35］方法文獻(xiàn)［13］方法平均運行時間/s 1.374 2 0.329 2 0.663 9 0.664 2 0.681 6 3.274 2

5 結(jié)論

為了提高在大觀測誤差條件下的定位性能，該文提出了一種基于加權(quán)多維標(biāo)度分析的多個非相關(guān)輻射源TDOA/FDOA協(xié)同定位方法，主要結(jié)論包括：

1）在大觀測誤差條件下，本文方法的估計精度優(yōu)于其他閉式定位方法，其具有更高的誤差閾值。

2）在大觀測誤差條件下，本文方法的估計精度與迭代初始值為真實值的迭代定位方法接近（以性能較優(yōu)的ICWLS定位方法為比較對象）。

3）本文方法的估計精度優(yōu)于已有的加權(quán)多維標(biāo)度定位方法。

4）本文方法可以進(jìn)一步提高傳感器位置和速度的估計精度（相比其先驗觀測精度而言）。

5）本文方法通過多源協(xié)同定位獲得了較高的協(xié)同增益。

附錄A首先定義如下向量和矩陣

式中：

于是有進(jìn)一步由式（27）可得

接著將定義式代入式（A3）中的矩陣Wn、、Tn以及中可以得到關(guān)于向量和ˉ的恒等式，將該恒等式兩邊分別對向量和求導(dǎo)可知

將式（48）代入式（A4）和式（A5）中可得

分別將式（A6）和式（A7）中的N個等式合并可知

將式（A8）代入式（79）中，并且利用式（66）可得

根據(jù)矩陣J的定義可以驗證JI6N+6M，由此可以將矩陣進(jìn)一步表示為

最后對比式（67）和式（A10），并且利用矩陣不等式A-1≥B(BTAB)-1BT（其中A為任意正定矩陣，B為任意列滿秩矩陣）可知命題2成立。

附錄B

附錄B中以實數(shù)乘法為參照給出了新方法的計算復(fù)雜度，具體見表B1。

表 B1 新方法的計算復(fù)雜度
Table B1 Computational complexity of proposed method

步驟步驟a-1步驟a-2步驟a-3步驟a-4步驟a-5步驟a-6步驟a-7步驟a-8步驟b-1步驟b-2步驟b-3步驟b-4各步驟計算量M(M+1)N(107M+O(53))N 15M2N 128MN3+16MN2+64N2(12M4+438M2+225M+225+O(53))N 16M(M-1)2N(16(M-1)3+24M2(M-1)+48M(M-1)2+O(8(M-1)3))N 32(M-1)2N3+128(M-1)N3+96M(M-1)N2+4M(M-1)N+(8N+6M)(6M+2(M-1)N)+O((6M+8N)3)32(M-1)2N3+128(M-1)N3+96M(M-1)N2+O((6M+8N)3)24N 12N 2(8N+6M)2N+(8N+6M)×(4N2+2N)+4N2+O(8N3)總計算量12M4N+64M2N3+192M2N2+256MN3+104M3N+56MN2+342M2N-96N3+68N2+36M2+489MN+309N+NO(8(M-1)3)+2O((6M+8N)3)+2NO(53)+O(8N3)