一種基于加權多維標度分析的多個非相關源TDOA/FDOA協同定位方法

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

1.中國人民解放軍戰略支援部隊信息工程大學 信息系統工程學院，鄭州 450001

2.國家數字交換系統工程技術研究中心，鄭州 450002

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

4.中國人民解放軍戰略支援部隊信息工程大學 數據與目標工程學院，鄭州 450001

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

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

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

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

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

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

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

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

1 TDOA/FDOA定位觀測模型

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

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

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

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

現將式（3）中的2個等式合并可得

式中：

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

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

式中：

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

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

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

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

2.1 新方法的總體描述

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

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

2）階段b首先基于多維標度分析中引入的中間變量，構建關于階段a定位誤差的約束優化模型，然后基于此模型對階段a的定位誤差進行估計，最后利用該估計值對階段a的定位結果進行更新，以獲得漸近統計最優的估計性能。

2.2 階段a的計算原理

2.2.1 構造標量積矩陣

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

式中：j表示虛數單位；dn1=rn1-rn1=0和=-=0。基于式（11）可定義如下復坐標和復速度矩陣

2.2.2 構建定位關系式

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

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

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

利用矩陣和求逆公式可知

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

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

最后聯合式（13）、式（20）以及式（21）可得

證畢。

若定義矩陣

則由命題1可知

將式（24）兩邊對時間求導可得

式中

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

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

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

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

式中：

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

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

式中：

其中

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

式中

其中

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

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

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

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

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

2.2.4 估計準則及其最優解

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

式中：

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

式中：

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

式中：

為了提高發生門限效應時的誤差閾值，這里考慮對向量和ˉ進行聯合估計，為此需要定義擴維參數向量和擴維誤差向量然后結合式（10）和式（53）可構建階段a中的估計準則，即

式中：

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

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

2.2.5 理論性能分析

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

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

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

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

式中：

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

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

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

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

2.3 階段b的計算原理

首先推導誤差向量Δρa服從的等式約束。根據向量ρ第8n-4個元素和第8n個元素的定義可知

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

式中：

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

式中：

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

利用拉格朗日乘子技術可知［46］，式（75）的最優解為

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

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

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

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

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

3 新方法的理論性能分析

本節推導估計值的統計特性，主要推導其估計均方誤差，并與相應的CRB進行比較，從而證明其漸近統計最優性，為此首先給出多個非相關源協同定位的CRB的表達式。

3.1 多個非相關源協同定位的CRB

基于文獻［34-35］中的結論可知，多個非相關源協同定位的CRB矩陣可以表示為

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

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

3.2 估計均方誤差

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

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

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

3.3 漸近統計最優性分析

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

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

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

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

于是有

由此可知

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

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

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

證明：首先根據向量ρ的定義可得

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

結合式（66）和式（A8）可得

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

4 仿真實驗

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

表1 傳感器位置和速度數值
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假設有2個非相關輻射源，輻射源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個輻射源進行協同定位，并且進行5 000次蒙特卡羅實驗，圖2和圖3給出了輻射源1和輻射源2的定位結果散點圖與定位誤差橢圓曲線。

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

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

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

仿真實驗2假設有3個非相關輻射源，輻射源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個輻射源進行協同定位，并將其定位精度與一些主流TDOA/FDOA定位方法進行比較，其中包括文獻［27］中的TLS方法、文獻［31］中的TSWLS方法、文獻［34］中的TSWLS方法、文獻［35］中的改進型TSWLS方法以及文獻［13］中的ICWLS方法，其中前面4種是閉式定位方法，最后1種是迭代定位方法。需要指出的是，雖然文獻［27］中的TLS方法、文獻［31］中的TSWLS方法以及文獻［13］中的ICWLS方法都是在單輻射源存在的場景下提出的，但卻可以將它們推廣至多輻射源協同定位的場景。此外，文獻［31］中的TSWLS方法聯合利用了DDR信息進行定位，但為了公平比較，這里的仿真并未利用此信息。另一方面，選擇ICWLS方法作為比較對象是因為它是一種性能優越的迭代定位方法［13］，為了進行更全面的比較，這里給出其初始值為隨機值和真實值2種情形下的定位性能。首先令20lg(σ2)=0 dBm，改變σ1的數值，圖4給出了多輻射源位置和速度估計均方根誤差隨20lg(σ1)的變化曲線；圖5給出了多傳感器位置和速度估計均方根誤差隨20lg(σ1)的變化曲線。然后令20lg(σ1)=-5 dBm，改變σ2的數值，圖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中可以看出：① 在大觀測誤差條件下本文方法的估計精度優于TLS方法和3種TSWLS方法，新方法發生門限效應的誤差閾值更高，這驗證了該方法相比其他閉式定位方法的優越性；② 在大觀測誤差條件下，本文方法估計精度優于初始值為隨機值的ICWLS方法，并且與初始值為真實值的ICWLS方法接近，注意到將真實值作為初始值在實際計算中難以實現，而本文方法無需初始值，這驗證了本文方法相比迭代定位方法的優越性；③ 本文方法對輻射源定位性能漸近逼近相應的CRB，這驗證了第3節理論性能分析的有效性；④ 本文方法可以提高傳感器位置和速度估計精度（相比其先驗觀測精度而言），并且對傳感器位置和速度估計均方根誤差均可以漸近逼近相應的CRB，從而再次驗證了第3節理論性能分析的有效性。

仿真實驗3仍然考慮3個非相關輻射源存在的場景，輻射源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。由于本文方法可以通過簡化直接應用于單輻射源存在場景（此時退化為非協同定位方法），于是下面不妨比較多輻射源協同定位和多輻射源非協同獨立定位2種方式下的定位精度，旨在驗證多源協同定位能產生協同增益。此外，還將文中方法與文獻［30］中的加權多維標度定位方法進行比較，這是因為該方法也考慮了傳感器位置和速度先驗觀測誤差，只是該方法是針對單輻射源的定位方法，目前尚沒有進行多源TDOA/FDOA協同定位的加權多維標度定位方法。另一方面，文獻［25］也是基于加權多維標度的TDOA/FDOA定位方法，但是該方法沒有考慮傳感器位置和速度先驗觀測誤差，因此不適合與文中的新方法進行性能比較。令20lg(σ1)=-10 dBm，20lg(σ2)=-5 dBm以及γ=12°，通過改變α的數值來改變輻射源位置。圖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中可以看出：① 多源協同定位產生的性能增益是顯著的；② 本文方法（非協同定位版）的估計精度優于文獻［30］中的加權多維標度定位方法，這是因為后者尚未利用中間變量包含的有用信息；③ 隨著輻射源位置的改變，文獻［30］中的加權多維標度定位方法的估計性能并不總能達到相應的CRB，而本文方法包含2個計算階段，可以使其定位精度始終漸近逼近相應的CRB。

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

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

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

5 結論

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

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

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

3）本文方法的估計精度優于已有的加權多維標度定位方法。

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

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

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

式中：

于是有進一步由式（27）可得

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

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

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

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

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

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

附錄B

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

表 B1 新方法的計算復雜度
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)