王鼎，尹潔昕,*，高路，楊賓

1. 中國人民解放軍戰(zhàn)略支援部隊信息工程大學 信息系統(tǒng)工程學院，鄭州 450001 2. 國家數(shù)字交換系統(tǒng)工程技術研究中心，鄭州 450002 3. 北京航天長征飛行器研究所，北京 100076

眾所周知，輻射源定位在無線通信、目標監(jiān)測、電子對抗、導航遙測等工程領域中具有廣闊的應用前景。輻射源定位觀測量涵蓋空、時、頻、能量等多域參數(shù)，其中到達時間差(Time Difference of Arrival，TDOA)是應用較頻繁的一類觀測量。隨著現(xiàn)代通信技術和TDOA測量技術的不斷發(fā)展，基于多個傳感器的TDOA定位技術已經成為最為主流的輻射源定位手段之一，該類定位體制適用于寬帶信號，能夠獲得較高的定位精度。

近些年來，國內外學者提出了一系列TDOA定位方法，從計算方式上可分為參數(shù)搜索類方法、迭代類方法以及閉式類方法3大類。參數(shù)搜索類方法包括網格逐點搜索方法和重要性采樣方法等。迭代類方法包括凸優(yōu)化方法、泰勒(Taylor)級數(shù)迭代方法、約束總體最小二乘(Constrained Total Least Squares，CTLS)估計方法以及約束加權最小二乘(Constrained Weighted Least Squares，CWLS)估計方法等。雖然參數(shù)搜索類方法和迭代類方法的定位精度都能逼近克拉美羅下界(Cramér-Rao Lower Bound，CRLB)，但前者計算量比較龐大，后者大都對迭代初值較敏感，易導致局部收斂。相比較而言，閉式類方法具有更高的計算效率，其中包括球面內插(Spherical-Interpolation，SI)方法、總體最小二乘(Total Least Squares,TLS)估計方法、兩步加權最小二乘(Two Step Weighted Least Squares，TSWLS)估計方法以及加權多維標度(Weighted Multidimensional Scaling，WMDS)方法等。SI方法和TLS方法的定位性能難以達到CRLB。TSWLS方法的定位精度雖然能夠逼近CRLB，但是該方法產生“門限效應”的誤差閾值相對較小，在大觀測誤差條件下的定位誤差會大幅增加。WMDS方法利用了標量積矩陣的維度和特征結構信息，對大觀測誤差具有更強的魯棒性，能夠獲得漸近統(tǒng)計最優(yōu)的定位精度，是本文研究的重點。

上述TDOA定位方法均假定傳感器位置精確已知，當傳感器安放在艦載或機載平臺，又或是傳感器隨機布設時，其位置精確值通常無法獲知，其中會包含先驗觀測誤差，該誤差會嚴重影響TDOA定位精度。為了抑制此類誤差的影響，國內外學者提出了3種方法。第1種方法將傳感器位置先驗觀測誤差和TDOA觀測誤差同等看待，通過重新設置各類估計器中的加權矩陣來提高定位魯棒性。第2種方法將輻射源位置與傳感器位置進行聯(lián)合估計。第3種方法利用位置精確已知的標校源信息實現(xiàn)目標定位。研究結果表明，第3種方法的定位精度最高，但需要標校源TDOA信息，這在實際中難以獲得，而第2種方法不僅可以獲得更準確的傳感器位置估計值，還比第1種方法具有更高的誤差閾值。因此本文主要針對第2種方法展開研究。

除了TDOA觀測誤差和傳感器位置誤差外，影響TDOA定位精度的另一個因素是不同傳感器之間的同步時鐘偏差。在無線傳感網絡(Wireless Sensor Network，WSN)節(jié)點定位領域，國內外學者提出一系列聯(lián)合定位與時鐘同步的有效方法。在一個WSN中，錨節(jié)點(位置精確已知或者近似已知)與源節(jié)點(位置未知)之間進行單向或者雙向通信，基于此可以建立這2類傳感節(jié)點位置參數(shù)與時鐘同步參數(shù)之間的代數(shù)關系，從而實現(xiàn)定位與時鐘同步的聯(lián)合處理。

在上述研究中，一些方法假設傳感器時鐘偏差完全一致(即完全同步)，另一些方法假設各個傳感器時鐘偏差互不相同(即完全不同步)。然而，這2種建模方式可能與實際情況都不相符。事實上，對于空域中廣泛分布的傳感器，若2個傳感器間距較遠，在采集信號時參考的局部時鐘易出現(xiàn)差異，但是若2個傳感器間距較近，時鐘同步則容易得到滿足。對此，文獻[27]提出一種新的定位觀測模型，其依據(jù)傳感器之間的距離大小將傳感器分成若干組，每組傳感器共用相同時鐘偏差，但不同組之間的時鐘偏差互不相同。基于該定位觀測模型，文獻[27]還提出一種基于TSWLS原理的TDOA定位方法，該方法基于最小二乘估計框架依次給出輻射源位置、傳感器位置以及時鐘偏差的估計值。為了提高該方法在大觀測誤差條件下的定位精度，文獻[28]提出基于Taylor級數(shù)迭代的輻射源位置、傳感器位置以及時鐘偏差聯(lián)合估計方法，該方法發(fā)生“門限效應”的誤差閾值更高，但需要更多計算量，并且易受迭代初始值的影響，可能出現(xiàn)迭代發(fā)散和局部收斂等問題。

為了進一步提高定位性能，本文在同步時鐘偏差和傳感器位置誤差同時存在下提出一種基于WMDS原理的TDOA定位方法。該方法的TDOA觀測模型借鑒文獻[27]中的觀測模型，但是定位方法與文獻[27]中的方法截然不同，其并不是基于最小二乘估計框架衍生出的，而是根據(jù)多維標度分析原理提出的，可以提高對大觀測誤差的魯棒性。新方法首先通過構造消元矩陣消除同步時鐘偏差的影響；然后基于WMDS原理構建定位關系式，由此獲得輻射源位置與傳感器位置的估計值；最后利用最大似然估計準則得到同步時鐘偏差的估計值。所提方法實現(xiàn)了時鐘偏差參數(shù)與位置參數(shù)的解耦合估計，可以有效抑制時鐘偏差和傳感器位置誤差的影響。最后，本文對新方法的定位精度給出理論性能分析，結果表明其對全部參數(shù)的估計精度均能逼近CRLB。仿真實驗驗證新方法的優(yōu)越性和理論性能分析的有效性。

表1 本文使用的矩陣等式和不等式Table 1 Matrix equality and inequality used in this paper

1 定位觀測模型與問題形成

關于輻射源RDOA的觀測模型可以表示為

(1)

(2)

式中：

(3)

其中：觀測誤差服從均值為0、協(xié)方差矩陣為的高斯分布；表示同步距離偏差向量。

由于主站位于第1組，所有RDOA觀測量均以主站為參考，因此式(1)中第1組等式對應的同步時鐘偏差等于0。

(4)

式中：表示先驗觀測誤差，其服從均值為0、協(xié)方差矩陣為的高斯分布，并且與誤差向量統(tǒng)計獨立。為了便于表述，對先驗觀測誤差進行分塊，即

(5)

在上述建模過程中，由于RDOA和傳感器位置均來自觀測，因此應將RDOA觀測誤差和傳感器位置誤差建模成隨機型誤差，而同步時鐘偏差并沒有先驗觀測，其是完全未知的參數(shù)，因此應按照固定的系統(tǒng)偏差對其進行建模?；谶@2類誤差模型的差異性，本文將向量和劃分成一組參數(shù)，將向量單獨作為另一組參數(shù)，并且實現(xiàn)這2組參數(shù)的解耦合估計。

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

2.1 新方法的原理概述

針對觀測模型式(2)的代數(shù)特點，新方法的求解原理可概述為以下4點：

1) 線性消元。由于式(2)是部分線性模型，因此可以通過矩陣變換構造矩陣的正交補子空間，從而消除同步距離偏差向量的影響，并形成新的降維觀測模型。

2) 多維標度分析。消元后的觀測模型仍然具備距離差結構形式，因此可引入多維標度分析方法，構造標量積矩陣，形成定位關系式，并獲得位置估計值。

3) 信息補償。由于在多維標度分析中引入了輔助變量，其與位置向量和相關，因此構建約束優(yōu)化模型對這部分信息進行有效補償，以獲得漸近統(tǒng)計最優(yōu)的定位結果。

4) 最大似然估計。在獲得位置向量和的估計結果之后，回到最初的最大似然估計模型中對線性參數(shù)進行估計，以實現(xiàn)多類型參數(shù)的解耦合估計。

基于上述原理概述，本文將新方法的計算過程劃分成3個階段(分別稱為階段a、階段b以及階段c)，圖1描述了其總體技術路線。

下面依次闡述每個階段的計算原理，并最終給出全部計算步驟。

圖1 本文新方法的總體技術路線示意圖Fig.1 Schematic diagram of overall technical route of new method proposed in this paper

2.2 階段a的計算原理

階段a首先構建消元矩陣以消除同步時鐘偏差的影響；接著針對每一組傳感器構造標量積矩陣，并利用標量積矩陣的性質形成定位關系式；然后通過一階誤差分析方法推導定位關系式中的誤差協(xié)方差矩陣，用于確定加權矩陣；最后基于定位關系式和加權矩陣同時獲得輻射源位置與傳感器位置的估計值。

2.2.1 構建消元矩陣

式(2)是關于向量的線性觀測模型，可利用消元矩陣消除其影響，根據(jù)矩陣的結構可將消元矩陣寫為

=blkdiag(,,…,)∈(-)×(-1)

(6)

(7)

式中：=(,)；=。顯然，觀測誤差服從均值為0、協(xié)方差矩陣為=的高斯分布。

(8)

式中:=()。式(8)的證明見附錄B，該式對于本文的理論性能分析至關重要。

為了便于描述，進行向量分塊

(9)

結合式(7)和式(9)可得

1≤≤

(10)

式中:

其中：

由于向量仍具有距離差的結構形式，因此利用多維標度分析方法對向量和進行聯(lián)合估計。

2.2.2 構造標量積矩陣及定位關系式

在多維標度分析中需要構造標量積矩陣，文中的定位場景包含組傳感器，因此應構造個標量積矩陣。借鑒文獻[11,17]中的WMDS定位方法，這里針對第組傳感器，首先定義復坐標向量

(11)

式中：j表示虛數(shù)單位，其滿足j=-1;1=0。利用復坐標向量進一步定義復坐標矩陣

(12)

(13)

由式(13)可知，針對組傳感器構造了個標量積矩陣{}1≤≤，利用它們可以建立定位關系式，首先有結論1。

分別定義實向量和實矩陣

(14)

(15)

×11≤≤

(16)

結論1的證明見附錄C?；诮Y論1不妨定義矩陣

(17)

將式(17)代入式(16)中，并且由向量的定義可知

1≤≤

(18)

將式(18)中的個等式進行合并可得

(19)

式中：

(20)

式(19)為最終獲得的定位關系式，利用式(19)可以構建輻射源定位準則。

注意到式(7)中包含-個觀測方程，但是式(19)中卻包含個觀測方程，這意味著式(19)是有冗余的，這種冗余性易導致誤差協(xié)方差矩陣出現(xiàn)秩虧損現(xiàn)象，文獻[11]提出利用正則化技術使其恢復為滿秩矩陣，但這僅僅是一種近似處理方式，本文提出利用矩陣正交變換使誤差協(xié)方差矩陣恢復為滿秩矩陣，并且通過理論性能分析證明這種處理方式能夠獲得漸近統(tǒng)計最優(yōu)的參數(shù)估計精度。

2.2.3 一階誤差擾動分析

(21)

(22)

(23)

(24)

基于式(19)可以定義誤差向量為

(25)

(26)

式中:

(27)

(28)

利用誤差矩陣Δ的表達式可以將式(28)右邊第1項寫成關于觀測誤差和的線性函數(shù)，即

(29)

式中：

其中：

(30)

(31)

(32)

式(29)的證明見附錄E?；谡`差矩陣Δ的表達式可以將式(28)右邊第2項寫成關于觀測誤差和的線性函數(shù)，即

(33)

式中:

其中:

(34)

(35)

(36)

式(33)的證明見附錄F。

()=()≈()+

()

(37)

在注釋5中指出，由于定位關系式存在冗余性，易導致協(xié)方差矩陣()出現(xiàn)秩虧損現(xiàn)象，因而無法對該矩陣直接進行求逆運算，下面提出利用矩陣正交變換技術解決該問題。

首先對矩陣進行分解可得

(38)

(39)

由式(39)可知，誤差向量的協(xié)方差矩陣為

(40)

2.2.4 估計準則及其最優(yōu)估計值

(41)

式中:

(42)

基于式(41)可以獲得階段a中的估計準則，即

(43)

式中:

()=

(44)

式(43)中的(())可以看成是加權矩陣，其作用在于抑制觀測誤差和的影響，式(43) 的最優(yōu)估計值為

(45)

(46)

(47)

2.2.5 理論性能分析

(48)

(49)

將式(42)中的第2式和式(49)代入式(48)中可得

(50)

式中

(51)

(52)

2.3 階段b的計算原理

向量中的元素間存在約束關系，該約束關系會使得誤差向量Δ服從等式約束。階段b首先基于多維標度分析中引入的輔助變量，構建關于誤差向量Δ的約束優(yōu)化模型，然后基于此模型對誤差向量Δ進行估計，最后利用其估計值對階段a中的定位結果進行更新，以獲得漸近統(tǒng)計最優(yōu)的估計結果。

首先推導誤差向量Δ所服從的等式約束。由向量中的第3+個元素的定義可知

1≤≤

(53)

1≤≤

(54)

將式(54)中的個等式進行合并，可以得到關于誤差向量Δ的線性等式，即

(55)

式中:

(56)

式(55)為誤差向量Δ所應滿足的等式約束，其是線性等式約束。根據(jù)誤差向量Δ所服從的高斯分布特性，可以建立估計誤差向量Δ的約束優(yōu)化模型，即

(57)

利用拉格朗日乘子技術可知，式(57)的最優(yōu)估計值為

(58)

聯(lián)合式(47)和式(58)可以得到誤差向量Δ和Δ的估計值分別為

(59)

于是輻射源位置向量和傳感器位置向量在階段b的估計結果為

(60)

2.4 階段c的計算原理與新方法的計算步驟

(61)

式(61)是關于向量的二次優(yōu)化問題，其估計值可以表示為

(62)

結合階段a～階段c中的描述，下面總結新方法的計算步驟，如圖2所示。

圖2 本文新方法的計算流程圖Fig.2 Computational flow chart of new method proposed in this paper

下面針對圖2中描述的定位方法給出2點注釋：

階段a中的循環(huán)計算主要用于更新加權矩陣(())，以提高對大觀測誤差的穩(wěn)健性，其每次循環(huán)中無需設置迭代步長和收斂準則，并且循環(huán)次數(shù)是確定的。事實上，現(xiàn)有定位方法很多都需要通過循環(huán)計算來更新加權矩陣，以提高定位方法在大觀測誤差條件下的性能。

新方法的第1步通過消元矩陣消除了同步距離偏差向量的影響，從而實現(xiàn)了該向量與位置向量和的解耦合估計。采用這種求解方式的原因在于：① 只有首先在觀測模型中消除向量的影響，才能利用多維標度分析方法對輻射源位置和傳感器位置進行估計；② 解耦合估計方法比多參數(shù)聯(lián)合估計方法的計算復雜度更低，在求解過程中可以避免多參數(shù)之間的相互影響，具有更高的穩(wěn)健性。

3 新方法的理論性能分析

3.1 估計值的均方誤差及其漸近統(tǒng)計最優(yōu)性分析

3.1.1 估計均方誤差

首先將式(57)中的等式約束代入式(58)中可得

(63)

(64)

(65)

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

(66)

首先結合式(65)和表1第4個矩陣等式可得

(67)

(68)

于是有

(69)

由此可知

(70)

將式(70)代入式(67)中，并且利用表1第3個矩陣等式可得

(71)

基于結論3還可以進一步得到結論4。

首先結合式(50)和式(66)可得

(72)

將式(51)代入式(72)中的各個子矩陣塊中，并且利用式(G3)和式(G4)可知

((,))()(,)

(73)

((,))()(,)

(74)

((,))()(,)+()

(75)

將式(73)～式(75)代入式(72)中，并且結合式(8)可知結論4成立。證畢。

3.2 估計值的均方誤差及其漸近統(tǒng)計最優(yōu)性分析

3.2.1 估計均方誤差

首先將式(2)代入式(62)中可得

(-(,)Δ-(,)Δ)?Δ=

(76)

(())+(())()·

()(())--

(77)

式中：

=(())()·

()(())

(78)

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

在一階誤差分析理論框架下，滿足=(-1)×(-1)。

首先根據(jù)線性參數(shù)估計理論和式(48)可得

Δ≈-((()))(())=

(79)

然后結合式(60)和式(63)可知

(80)

將式(39)和式(79)代入式(80)中，并且利用關系式=可得

(81)

由式(81)可知

(82)

式中第2個等號利用了關系式=(-)×(-1)。最后將式(82)代入式(78)中可得=(-1)×(-1)。證畢。

將=(-1)×(-1)代入式(77)中，并且結合式(A3)和結論4可知

(())()·

()

(83)

由于最大似然估計器同樣具有漸近統(tǒng)計最優(yōu)性，因此這里有必要討論文中的新方法與最大似然估計方法之間的關系。首先，階段c中向量的估計值正是基于最大似然估計準則所獲得的。其次，向量和的估計值并不是直接由最大似然估計準則所產生，而是通過對原觀測模型進行變換處理后所獲得。根據(jù)統(tǒng)計信號處理理論可知，對原始觀測量進行代數(shù)變換后，在不損失觀測信息并且合理設置加權矩陣的條件下，能夠得到與最大似然等價的估計結果，這正是本文所采取的處理方式。因此，文中的新方法可以看成是最大似然估計方法的一種近似實現(xiàn)方式。另一方面，由于原觀測模型式(2)的非線性特征，2種方法的理論性能分析中都不可避免需要進行一階近似，當“門限效應”尚未產生時，一階近似誤差是可以忽略的，但是當“門限效應”產生時，一階近似誤差的影響就會逐漸顯現(xiàn)，但何時產生“門限效應”尚無完備的理論支撐，每種方法的主要區(qū)別在于產生“門限效應”時的誤差閾值大小。

4 仿真實驗

表2 傳感器3維位置坐標Table 2 3D position coordinate of sensors

將輻射源位置向量設為=[6 000, 6 000, 6 000]m，將同步距離偏差向量設為=[0, 50, 70]m，并令20lg()=0 dBm和20lg=5 dBm，下面利用文中的新方法對輻射源進行定位，并且進行5 000次蒙特卡羅實驗，圖3給出了輻射源定位結果散點圖與定位誤差橢圓曲線，圖中的橢圓曲線是利用文獻[31]中的方法計算所得。

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

圖3 輻射源定位結果散點圖與定位誤差橢圓曲線Fig.3 Scatter diagram of emitter location result and ellipse curve of location error

這里將本文新方法與文獻[17]中的WMDS定位方法、文獻[27]中的閉式定位方法以及文獻[28]中的迭代定位方法(實質上就是最大似然估計方法)進行比較，用于驗證新方法的優(yōu)越性。需要指出的是，文獻[17]中的WMDS定位方法并未考慮同步時鐘偏差的影響，而且該方法僅給出了輻射源位置估計結果，所以無法與其比較傳感器位置和同步距離偏差的估計精度。此外，與新方法相類似，文獻[27]中的閉式定位方法和文獻[28]中的迭代定位方法也具有漸近統(tǒng)計最優(yōu)性，因此需要在大觀測誤差條件下比較參數(shù)估計精度。將輻射源位置向量設為=[8 600, 9 200, 8 500]m，將同步距離偏差向量設為=[0, 150, 180]m，首先令20lg=20 dBm，圖4分別給出了輻射源位置、傳感器位置以及同步距離偏差估計均方根誤差隨著參數(shù)的變化曲線；然后令20lg=15 dBm，圖5分別給出了輻射源位置、傳感器位置以及同步距離偏差估計均方根誤差隨著參數(shù)的變化曲線。

從圖4和圖5中可以看出：① 新方法的定位精度明顯優(yōu)于文獻[17]中的WMDS定位方法，這是因為后者并未考慮同步時鐘偏差的影響，而同步時鐘偏差屬于系統(tǒng)偏差，當該偏差出現(xiàn)時會使得文獻[17]中的定位觀測模型失配，從而使其定位方法變成有偏估計方法，并最終導致估計性能急劇下降，甚至失效；② 在大觀測誤差情形下,新方法的估計精度高于文獻[27]中的閉式定位方法，前者產生“門限效應”的誤差閾值更高，這正是多維標度分析方法所帶來的性能增益；③ 當文獻[28]中的迭代定位方法取隨機值作為迭代初始值時，其在大觀測誤差情形下的估計誤差要明顯大于新方法，這里的隨機值是指在真實值的基礎上疊加一定的高斯隨機誤差；④ 當文獻[28]中的迭代定位方法取真實值作為迭代初始值時，其估計精度與新方法接近，但是將真實值作為迭代初始值在實際應用中難以實現(xiàn)，而新方法并不存在此問題，并且計算復雜度更低(見表3)；⑤ 新方法的參數(shù)估計均方根誤差能夠逼近相應的CRLB，從而驗證第4節(jié)理論性能分析的有效性。

圖4 輻射源位置、傳感器位置和同步距離偏差的估計均方根誤差隨參數(shù)σ1的變化曲線Fig.4 Curves of RMSE of emitter position, sensor position and synchronization distance bias estimation as a function of σ1

本節(jié)的最后將文中的新方法與文獻[27]中的閉式定位方法以及文獻[28]中的定位迭代方法的運行時間進行比較，用于間接說明3種定位方法的計算復雜度。仿真軟件為MATLAB R2020a，仿真程序在安裝有i7-3520 CPU的PC機上運行，令20lg=15 dBm和20lg=20 dBm，其余條件不變，并進行5 000次蒙特卡羅實驗。表3給出了3種定位方法的平均運行時間。從表3中可以看出，新方法的計算復雜度略高于文獻[27]中的閉式定位方法，但是明顯低于文獻[28]中的迭代定位方法。

圖5 輻射源位置、傳感器位置和同步距離偏差估計均方根誤差隨參數(shù)σ2的變化曲線Fig.5 Curves of RMSE of emitter position, sensor position and synchronization distance bias estimation as a function of σ2

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

5 結 論

本文在同步時鐘偏差和傳感器位置誤差同時存在的條件下研究了TDOA定位問題，主要結論包括：

1) 提出了一種基于加權多維標度分析的TDOA定位新方法，并利用一階誤差性能分析方法推導了新方法的參數(shù)估計性能，證明了其估計性能能夠漸近逼近CRLB。

2) 仿真實驗結果驗證了新方法的有效性和漸近統(tǒng)計最優(yōu)性。

3) 新方法的性能與初始值為真實值的迭代方法接近(以文獻[28]中的迭代定位方法為比較對象)，相比于文獻[27]中的閉式定位方法，其產生“門限效應”的誤差閾值更高，對于大觀測誤差具有更強的魯棒性。