專利名稱：采用總體最小二乘和均衡算法的到達時差定位方法
技術(shù)領(lǐng)域：
采用總體最小二乘和均衡算法的到達時差定位方法屬于空間目標(biāo)定位技術(shù)領(lǐng)域。
背景技術(shù)：
目標(biāo)的定位精度除了與目標(biāo)相對于各接收陣元的空間幾何位置有關(guān)外，還取決于定位參數(shù)的估計誤差以及陣元位置的測量誤差。在實際應(yīng)用中，受信噪比的影響和具體陣元位置測量技術(shù)(如全球定位系統(tǒng))的限制，這些誤差難以消除，且誤差方差也不盡相同。如果采用的定位算法不能將這些誤差因素考慮進去，則實際定位精度將受到嚴(yán)重影響。
TDOA定位技術(shù)具有較強的工程可實現(xiàn)性和較高的定位精度，因此被廣泛運用于雷達、聲納及輻射源定位系統(tǒng)中。
TDOA定位方法的原理1、常規(guī)的TDOA定位方法如圖1所示，設(shè)定三維空間中任意分布著M個信號接收陣元，不失一般性，假定陣元1處于坐標(biāo)原點，各陣元的位置分別為si＝[xiyizi]T(i＝1，...，M)，目標(biāo)的位置為u＝[xyz]T，這里[]T表示矩陣轉(zhuǎn)置。各陣元到原點的距離為Ri＝‖si‖，‖‖表示Euclidean范數(shù)。這樣，陣元i到目標(biāo)的距離就是ri＝‖si-u‖，陣元i、j到目標(biāo)的距離差(RD)就是ri，j＝ri-rj＝‖si-u‖-‖sj-u‖，(i，j＝1，...，M)(1)因為距離差除以信號傳播速度就是目標(biāo)信號到達兩陣元的時間差，因此不妨將RD作為觀測量建立觀測方程。令j＝1，(1)式可以變換為xix+yiy+ziz+ri，1r1＝-0.5(ri，12-Ri2)(i＝2，...，M) (2)將(2)式表示為矩陣的形式，得到A1u1＝b1(3)其中A1=-x2y2z2r2,1x3y3z3r3,1MMMMxMyMzMrM,1---b1=0.5r2,12-R22r3,12-R32MrM,12-RM2u1=xyzr1]]>當(dāng)系數(shù)矩陣A1的秩等于未知數(shù)個數(shù)n，且M≥n時，可按照經(jīng)典的最小二乘算法求出目標(biāo)的估計位置
u^1=(A1TA1)-1A1Tb1---(4)]]>當(dāng)只有觀測向量b1含有誤差，且各誤差分量相互獨立、方差相同，則最小二乘估計在方差最小的意義上是最佳的。但在實際應(yīng)用中，系數(shù)矩陣A1往往也有擾動，并且誤差的協(xié)方差矩陣不為單位陣(誤差方差不同或誤差分量具有相關(guān)性)，因此最小二乘估計法的實用性受到很大限制。
另一種常規(guī)方法是用泰勒展式(Taylor-series)將非線性觀測方程((1)式)轉(zhuǎn)化為線性方程，再通過迭代的方式估計未知向量。這種方法的主要缺點是需要設(shè)定初始值，而且因為迭代使得方法的實時性較差。
2、現(xiàn)有的TDOA定位方法主要有兩種方法一是用總體最小二乘(TLS)估計未知向量，這樣不僅考慮了觀測向量的擾動，也考慮了系數(shù)矩陣的擾動，使得定位精度得以提高。在這種方法中一般采用特征值分解(SVD)求解未知向量，因而又具有較高的數(shù)值穩(wěn)定性。在誤差分量的協(xié)方差陣為單位陣的情況下，TLS具有最佳的估計效果。但實際中，系數(shù)矩陣元素往往是與觀測向量不同量綱的各種測量值，因此這個條件很難滿足。方法二是改進的加權(quán)最小二乘法，此方法在構(gòu)造誤差矩陣時，不僅考慮了RD對觀測向量引起的擾動，同時也考慮了RD對系數(shù)矩陣引起的擾動，并將誤差矩陣的協(xié)方差陣進行求逆作為加權(quán)矩陣。不僅如此，該方法還考慮到未知向量u1中各元素的關(guān)系，即x2+y2+z2＝r12，并利用這一固有關(guān)系進一步提高了定位精度。方法二的最大缺點是沒有考慮陣元位置測量誤差，因為考慮該誤差時難以構(gòu)造出加權(quán)矩陣。在實際定位系統(tǒng)中大量存在的情況是不僅RD的估計具有誤差，陣元位置的測量也有誤差，并且誤差的方差也是不同的(由信噪比不同等因素引起)。

發(fā)明內(nèi)容
本發(fā)明的目的在于克服已有技術(shù)的不足之處，針對更接近實際的定位系統(tǒng)誤差模型，提出一種采用對誤差矩陣進行預(yù)均衡并結(jié)合TLS算法的定位方法，利用均衡后的誤差矩陣的協(xié)方差陣為單位陣的特性，提高定位精度。
當(dāng)TDOA估計和陣元位置測量都存在誤差時，由(3)式可知，系數(shù)矩陣A1及觀測向量b1都存在擾動。因此，增廣矩陣[A1b1]可表示為由真值組成的矩陣A0與誤差矩陣E之和，即[A1b1]＝[A10b10]+[e′e]＝A0+E。設(shè)ri，1＝ri，10+Δri，1，si＝si0+Δsi＝[xi0yi0zi0]T+[ΔxiΔyiΔzi]T(i，j＝2，...，M)，并忽略掉二階誤差項，則誤差矩陣E為E=-Δx2-Δy2-Δz2-Δr2,1r2,10Δr2,1-Δx3-Δy3-Δz3-Δr3,1r3,10Δr3,1MMMMM-ΔxM-ΔyM-ΔzM-ΔrM,1rM,10ΔrM,1---(5)]]>
式中ri，10和si0分別表示RD和陣元位置的準(zhǔn)確值，Δri，1和[ΔxiΔyiΔzi]T分別表示RD和陣元位置的誤差，并假定為相互獨立的零均值高斯白噪聲隨機變量。若誤差矩陣不滿足獨立同分布結(jié)構(gòu)(IID)，例如各元素的方差不等或元素間不相互獨立等，則已有的TLS算法就不能達到最佳估計效果。因此，將傳統(tǒng)的TLS算法替代為求解未知向量使之滿足minimize‖DET‖F(xiàn)(6)subject to b1+e∈Range(A1+e′)式中D、T分別為起均衡誤差矩陣作用的左乘加權(quán)矩陣和右乘加權(quán)矩陣，‖‖F(xiàn)表示Frobenius范數(shù)。左乘加權(quán)矩陣D是對不同陣元的定位參數(shù)估計誤差進行加權(quán)均衡處理；右乘加權(quán)矩陣T是對不同誤差(如陣元位置測量誤差、距離差估計誤差)進行加權(quán)均衡處理。由此可見，此方法比WLS和TLS具有更大的靈活性。
本發(fā)明提出的采用總體最小二乘和均衡算法的到達時差定位方法，如圖2所示，其特征在于它用計算機來依次執(zhí)行以下步驟定位計算開始后，計算機從步驟2a同時轉(zhuǎn)移到并行執(zhí)行的步驟2b、2c和2d。其中步驟2b是根據(jù)前端TDOA估計和陣元位置測量分別提供的觀測量(距離差)ri，1(i＝2，...，M)和陣元位置矢量si(i＝1，...，M)構(gòu)造出系數(shù)矩陣A1及觀測向量b1，這里假定陣元數(shù)目為M。步驟2c是根據(jù)距離差估計均方根誤差σi1計算用于均衡誤差矩陣的左乘加權(quán)矩陣DD=1/σ210Λ001/σ31Λ0MMM00Λ1/σM1---(7)]]>矩陣D是對涉及不同陣元的陣元位置測量誤差、距離差估計誤差進行加權(quán)均衡處理。步驟2d是利用距離差估計方差σi12、陣元位置測量誤差方差σ2和距離差觀測量ri，1，計算用于均衡誤差矩陣的右乘加權(quán)矩陣T，該矩陣是對不同誤差(如陣元位置測量誤差、距離差估計誤差)進行加權(quán)均衡處理。在這一步驟中，首先得到矩陣E′＝DE的行向量的協(xié)方差矩陣為di=σ2/σi1200000σ2/σi1200000σ2/σi121-ri,10000-ri,10(ri,10)2(i=2,...,M)---(8)]]>于是可以由下式計算右加權(quán)矩陣TTTT=E[E′TE′]-1=(Σidi)-1---(9)]]>
從而使得均衡后的誤差矩陣滿足E[(DET)TDET]～I，這樣采用基于特征值分解的TLS算法就能夠得到未知向量的最佳估計，可以證明，這種估計是無偏的。由(9)式可知，矩陣T可以方便的通過Cholesky分解得到。由于ri，10是未知量，實際計算中可用觀測值ri，1代替。為了進一步減少運算量和提高算法的健壯性，可以將T近似為一對角陣，即將矩陣E[E′TE′]-1的非對角元素置0，通過仿真可知，由以上近似帶來的誤差很小，可以忽略不記。
完成以上步驟后，估計器轉(zhuǎn)移到步驟2e，對均衡后的增廣矩陣作特征值分解D[A1b1]T＝UΩVT(10)在(10)式中，U是(M-1)×(M-1)對角矩陣，V是5×5對角矩陣，Ω是(M-1)×5對角矩陣，其對角元素就是矩陣D[A1b1]T的特征值，并按降序排列。將矩陣V、T、Ω分塊為以下形式V=V11(4×4)V15(4×1)V51(1×4)V55(1×1),Ω=Ω1(5×5)0(M-6×5),T=T1(4×4)0(4×1)0(1×4)t5---(11)]]>式中分塊矩陣的下標(biāo)(i×j)表示該分塊矩陣是由原矩陣的1～i行及1～j列所對應(yīng)元素組成，這里Ω1＝diag(σ1σ2Λσ4σ5)，σi(i＝1...5)是矩陣D[A1b1]T的特征值。
隨后估計器轉(zhuǎn)移到步驟2f，比較Ω1中按降序排列的最后兩個特征值的大小，若前一個特征值σ4大于最小特征值σ5，則估計器轉(zhuǎn)移到步驟2g；反之，轉(zhuǎn)回到步驟2b、2c、2d，從新開始新一輪的定位計算。步驟2g利用特征值分解結(jié)果計算未知向量u1的估計值 即得到目標(biāo)位置的初始估計值u^1=-T1V15V55t5---(12)]]> 還可以表示為另一種等價形式u^1=(A1TA1-σ52I4)-1A1Tb1---(13)]]>其中I4表示4階單位陣。
接下來估計器轉(zhuǎn)移到步驟2h，根據(jù)u1的估計值構(gòu)成觀測向量b2=u^12,]]>并計算加權(quán)矩陣W2。顯然向量u1中的4個元素并非相互獨立的，x，y，z，和r1滿足關(guān)系式r12＝x2+y2+z2(14)這一先驗信息可以用來進一步提高目標(biāo)位置的估計精度。根據(jù)(14)式建立以下觀測方程組A2u2＝b2(15)
其中A2=100010001111---b2=u^1(1)2u^1(2)2u^1(3)2u^1(4)2---u2=x2y2z2]]>由于系數(shù)矩陣A2是常系數(shù)矩陣，因此誤差矩陣ε由估計值 的各誤差分量組成。根據(jù)誤差矩陣與加權(quán)矩陣的關(guān)系W2=E(ϵϵT)-1|u^1=u10]]>(u10為向量u1的真值)，并忽略二階誤差項，可以得到加權(quán)矩陣W2的表達式W2=(Ccov(u^1)CT)-1---(16)]]>其中C＝diag(x0y0z0r10)這里，[x0y0z0]表示目標(biāo)的真實位置。而(17)式中的 的協(xié)方差矩陣可由(13)式得到cov(u^1)=(A1TA1-σ52I4)-1A1Tcov(b1)((A1TA1-σ52I4)-1A1T)T---(17)]]>同樣忽略掉高階誤差項，觀測向量b1的協(xié)方差矩陣為cov(b1)＝4diag((r2，10)2σ2，12(r3，10)2σ3，12Λ(rM-1，10)2σM-1，12)(18)實際計算中，目標(biāo)的真實位置[x0y0z0]可由 中的對應(yīng)元素代替，而矩陣cov(b1)中的ri，10(i＝2，...,M)可由觀測值ri，1代替。
完成以上計算后，估計器轉(zhuǎn)移到步驟2i，應(yīng)用加權(quán)最小二乘計算估計值 u^2=(A2TW2A2)-1A2TW2b2----(19)]]>接著估計器轉(zhuǎn)到步驟2j，對前一步的輸出進行開方運算并乘上符號矩陣G，從而得到目標(biāo)位置的估計值 u^=G[u^2(1)0.5,u^2(2)0.5,u^2(3)0.5]---(20)]]>其中G用于確定 中各元素的符號，可由初始估計值 中相應(yīng)元素的符號得知，即G=diag(sign(u^1(1)),sign(u^1(2)),sign(u^1(3)))---(21)]]>式中sign表示符號函數(shù)(自變量為正取1，為負(fù)取-1)。最后估計器轉(zhuǎn)移到步驟2k，定位計算結(jié)束。
本發(fā)明的特征還在于它可以通過一個由若干乘法器組成的電路完成，并依次包含以下步驟
(1)，輸入以下測定的初始參數(shù)值陣元位置si＝[xiyizi]T，i＝1，...，M，M為陣元數(shù)，其中第1個陣元設(shè)在X-Y-Z坐標(biāo)系的原點，觀測量ri，1，i＝2，...，M，ri，1為陣元i、1到目標(biāo)的距離差，陣元i，1到目標(biāo)的距離差的誤差方差σi12，陣元位置測量誤差方差σ2；(2)，把下列數(shù)組通過各自的運算電路運算后，其結(jié)果送往相應(yīng)的寄存器組數(shù)組ri，12，σi，12經(jīng)過乘法器運算后，得到cov(b1)，送往cov(b1)寄存器組，si經(jīng)過平方和電路得到Ri2，再和ri，12經(jīng)過減法器得到b1，送往b1寄存器組，把si，ri，1送往A1寄存器組，σi，1經(jīng)過倒數(shù)電路得到1/σi，1，送往D寄存器組，σi12、σ2經(jīng)過除法電路得到σ2/σi12，再把σ2/σi12和ri，1一起送往di寄存器組，得到di，然后對di求和，經(jīng)求逆、分解電路后形成T，送往T寄存器組；(3)，把上述各寄存器組的輸出A1、b1、D、T送往矩陣乘法電路得到D[A1b1]T，再把D[A1b1]T經(jīng)過公知的特征值分解電路得到U、V、Ω，并送往U、V、Ω寄存器組；(4)，把U、V、Ω寄存器組中的V和T寄存器組中的T，經(jīng)過矩陣乘法電路，得到 和C，并送往 寄存器組；(5)，把從相應(yīng)的寄存器組得到的cov(b1)、A1、σ5和C送往矩陣乘法電路，求出W2，并送往W2寄存器組(6)，使 經(jīng)平方電路得到b2，把b2、W2送往矩陣乘法電路得到 再把 經(jīng)過開方電路得到 作為一個例子，設(shè)定陣元數(shù)M＝9，各陣元位置分別位于
，[-50 80 20]，[40 60-30]，[-20 40 40]，[60 30 30]，[-70 50-20]，[20 50 10]，[-40 20-40]，[30 30-10]，輻射源位置位于[3002400 200]，按RD估計及陣元位置測量精度，分別設(shè)RD的方差σi12(i＝2，...，9)等于10-5[149 16 25 36 49 64]，陣元位置測量誤差方差σ2等于[10-410-3.75…10-2.2510-2]。算法用MATLAB編程，并在清華同方X08-34930型計算機上運行。通過105次蒙特卡羅仿真，得到三種方法的定位均方根誤差(RMSE)仿真結(jié)果(如表1所示)，這三種方法分別是改進的WLS定位方法，TLS定位方法和本發(fā)明提出的定位方法。
表1隨陣元位置測量誤差方差變化的定位均方根誤差仿真結(jié)果 可見，本發(fā)明提出的定位方法具有三種方法中最高的定位精度。本發(fā)明的效果是，利用加權(quán)矩陣對誤差矩陣作均衡處理，使得均衡后的誤差矩陣的協(xié)方差矩陣近似為單位陣，并通過具有較高數(shù)值穩(wěn)定性的特征值分解方法得到目標(biāo)位置的估計值，從而大幅度提高了目標(biāo)的定位精度。所處理的誤差模型既包括有陣元位置測量誤差，又包括互不相等的距離差誤差，因而更符合實際應(yīng)用情況。比之于已有方法一，在距離差估計誤差不相等時，本方法具有更高的定位精度；比之于已有方法二，在存在陣元位置測量誤差時，本方法具有更好的定位性能。


圖1為定位系統(tǒng)目標(biāo)和陣元位置關(guān)系示意2為本發(fā)明的軟件實現(xiàn)流程圖。
圖3為按照本發(fā)明設(shè)計的定位估計器的硬件實現(xiàn)方框圖。
圖4為應(yīng)用本發(fā)明方法的目標(biāo)定位系統(tǒng)框圖。
具體實施例方式
下面根據(jù)附圖和兩個實施例更加詳細(xì)地解釋本發(fā)明實施例一本實施例采用軟件實現(xiàn)本發(fā)明提出的定位方法，如圖2所示。設(shè)陣元位置數(shù)為9，各陣元位置分別位于
，[-50 80 20]，[40 60-30]，[-20 40 40]，[60 30 30]，[-70 50-20]，[20 50 10]，[-40 20-40]，和[30 30-10]，輻射源位置為[300 2400 200]，按RD估計及陣元位置測量精度，設(shè)RD的方差σi12(i＝2，...，9)等于10-5[1 4 9 16 25 36 49 64]，陣元位置測量誤差方差σ2等于10-3。定位計算開始后，估計器從步驟2a同時轉(zhuǎn)移到并行執(zhí)行的步驟2b、2c和2d。其中步驟2b是根據(jù)前端TDOA估計和陣元位置測量分別提供的觀測量ri，1(i＝2，...，9)和陣元位置矢量si(i＝1，...，9)構(gòu)造出系數(shù)矩陣A1及觀測向量b1A1＝ b1＝49.9809 -79.9710 -19.9427 73.78831.0e+003*-39.9848 -59.9946 29.9165 61.320619.9417 -40.0340 -40.0427 39.9638-1.9276-60.0169 -29.9716 -30.0322 38.7308-1.169969.9681 -50.0313 19.9082 37.8466-1.0014-19.9884 -50.0261 -9.976952.6999-1.950040.0178 -20.0079 39.9515 10.8141-3.1838-30.0209 -29.9921 10.0189 32.3748-0.1114-1.7415-0.4259步驟2c是根據(jù)觀測量估計方差σi1用(7)式計算用于均衡誤差矩陣的左乘加權(quán)矩陣D，該矩陣是對不同陣元的定位參數(shù)估計誤差進行加權(quán)均衡處理，其值為D＝316.2278 0 0 0 0 0 000158.1139 0 0 0 0 000 0 105.4093 0 0 0 000 0 079.0569 0 0 000 0 0 063.2456 0 000 0 0 0 0 52.7046 000 0 0 0 0 0 45.175400 0 0 0 0 0 0 39.5285步驟2d是利用觀測量估計方差σi12、陣元位置測量誤差方差σ2和觀測量ri，1，用(8)、(9)式計算用于均衡誤差矩陣的右乘加權(quán)矩陣T，該矩陣是對不同誤差(如陣元位置測量誤差、觀測量估計誤差)進行加權(quán)均衡處理，其值為T＝0.0809 0 0 0 00 0.0809 0 0 00 00.0809 0 00 0 00.9256 00 0 0 00.0197完成以上步驟后，估計器轉(zhuǎn)移到步驟2e，實現(xiàn)矩陣的乘法運算并進而完成乘積矩陣的特征值分解D[A1b1]T＝U Ω VT，得到U＝-0.8964 0.0051 -0.3776 -0.0135 0.2301 -0.0161 0.0223 -0.0081-0.3470 -0.3794 0.7437 -0.3207 -0.1550 -0.1958 -0.0491 -0.1219-0.1602 -0.0218 -0.1130 0.4711 -0.8099 -0.2688 0.1035 0.0152-0.1423 0.3549 0.5092 0.6879 0.3141 0.1270 0.0728 -0.0353-0.1402 0.7028 0.1214 -0.3189 -0.3189 0.1991 -0.4735 0.0676-0.0833 -0.3417 0.0132 0.0863 -0.2269 0.9040 -0.0054 -0.0069-0.0417 0.3333 0.0869 -0.3045 -0.1200 0.1251 0.8699 0.0140-0.0431 -0.0889 0.0990 0.0044 0.0284 -0.0247 0.0152 0.9893Ω＝1.0e+004*2.7700 0 0 0 000.3402 0 0 00 00.1137 0 00 0 00.0534 00 0 0 0 0.00010 0 0 0 00 0 0 0 00 0 0 0 0V＝-0.0356 0.1172 -0.9081 -0.3812 0.12280.0810 0.0286 0.0818 0.1237 0.98520.0140 -0.0225 0.3921 -0.9159 0.0820-0.8699 -0.4851 -0.0165 -0.0006 0.08700.4849 -0.8658 -0.1212 -0.0233 -0.0017其中Ω包含有乘積矩陣的特征值，并按降序排列在對角線上。然后估計器轉(zhuǎn)移到步驟2f，比較按降序排列的最后兩個特征值的大小。由Ω的值可知σ4＝0.0534，σ5＝0.0001，顯然滿足σ4＞σ5，因此估計器轉(zhuǎn)移到步驟2g。步驟2g利用(12)式計算未知向量u1的估計值 即得到目標(biāo)位置的初始估計值u^1=1.0e+003*0.29152.33850.19462.3637]]>隨后，估計器轉(zhuǎn)移到步驟2h，根據(jù)u1的估計值構(gòu)成觀測向量b2=u^12,]]>并通過(16)、(17)和(18)式計算出加權(quán)矩陣W2W2＝0.0528-0.0030 0.00100.0020-0.0030 0.0009 0.0010 -0.00080.0010 0.0010 0.0614 -0.00140.0020-0.0008 -0.00140.0008完成這一步后，估計器轉(zhuǎn)移到步驟2i，應(yīng)用(19)式計算估計值 u^2=1.0e+006*0.09105.84350.0406]]>接著估計器轉(zhuǎn)到步驟2j，對前一步的輸出進行開方運算并乘上符號矩陣G，從而得到目標(biāo)位置的估計值 u^=1.0e+003*2.30162.41730.2015]]>這里G用于確定 中各元素的符號，可由初始估計值 中相應(yīng)元素的符號得知。最后估計器轉(zhuǎn)移到步驟2k，定位計算結(jié)束。
實施例二本實施例為采用硬件的方法實現(xiàn)本發(fā)明提出的定位方法。仍然以陣元位置數(shù)為9，觀測量估計誤差方差σi12(i＝2，...，9)互不相等，且陣元位置測量也存在誤差的情況為例，說明該硬件電路的定位過程。如圖3所示，由前端TDOA估計和陣元位置測量提供的第1到第8路估計參數(shù)輸入31，經(jīng)過簡單的運算電路32～39(圖中|·|2表示平方和電路，(·)-1表示求倒數(shù)電路)，得到系數(shù)矩陣313(A1)、觀測矩陣314(b1)、左加權(quán)矩陣315(D)和右加權(quán)矩陣316(T)的各組成元素。其中各路輸入的距離差的平方值ri，12與距離差誤差方差σi，12的乘積被存儲在一個寄存器組中，為后級提供觀測向量的協(xié)方差矩陣312(cov(b1))；各陣元位置矢量si(i＝1，...，9)和距離差ri，1被存儲到一個寄存器組，為后級提供系數(shù)矩陣313；求出各陣元到原點的距離Ri2并與距離差的平方值相減，從而得到觀測向量314的各分量值，并被存儲在一個寄存器組中；距離差誤差方差σi，1被取倒數(shù)后，作為左加權(quán)矩陣的對角元素存儲在一個寄存器組中；各路輸入的距離差ri，12及其誤差方差σi12，以及陣元位置測量誤差方差σ2各自構(gòu)成矩陣di的元素，并通過310電路——包括矩陣相加、求逆和Cholesky分解器，得到右加權(quán)矩陣316，并被存儲在一個寄存器組中。以上5個寄存器組311提供的矩陣被送到矩陣乘法電路317中，求出乘積矩陣318(D[A1b1]T)的值。隨后對該矩陣進行特征值分解319，并將酉陣U、V和特征值組成的對角陣320存儲在寄存器組321中，其中矩陣322(V)的分塊矩陣V15和V55在后級矩陣乘法電路324中與右加權(quán)矩陣316的分塊矩陣T1和t5作乘法運算，得到陣元位置的初始估計值 和由其相應(yīng)元素組成的矩陣328(C)。其中特征值分解電路319由黃浴等人提出的Systolic陣列實現(xiàn)(“一個實現(xiàn)狀態(tài)空間法的Systolic陣列系統(tǒng)設(shè)計”，通信學(xué)報，黃浴等，Vol.16，No.6，pp.25～32，1995)。矩陣328和前面得到的協(xié)方差矩陣312、系數(shù)矩陣313以及最小特征值σ5一起被送到矩陣乘法電路325中，計算出用于加權(quán)最小二乘運算的加權(quán)矩陣332(W2)。初始估計值327經(jīng)過平方電路329得到觀測矩陣330(b2)，并與加權(quán)矩陣332一起送入到矩陣乘法電路333，求出估計值 再經(jīng)過平方根電路335輸出目標(biāo)位置估計值 從上面的實施例可見，本定位方法與其前面介紹過的兩種現(xiàn)有的定位方法相比，具有更高的定位精度以及更強的實用性，而且其硬件實現(xiàn)也并不復(fù)雜。下面再通過一個例子來說明本發(fā)明在輻射源定位系統(tǒng)中的作用。
參照圖4，采用總體最小二乘及均衡算法的目標(biāo)定位系統(tǒng)包括接收信號的天線陣41、接收前端42、模數(shù)轉(zhuǎn)換(A/D)、TDOA(或RD)估計47、陣元位置測量410、以及圖3所示的定位估計器49。在本例中，天線陣元接收的各路信號42被分別送往各自的接收前端43，并轉(zhuǎn)化為中頻信號44。模數(shù)轉(zhuǎn)換45將模擬中頻信號變?yōu)閿?shù)字信號流46，輸給下一級TDOA估計47。TDOA估計47由此計算出各路信號之間的相對延遲，從而得到觀測量48用于定位估計。陣元位置測量410輸出各陣元位置411，與觀測量48一起輸入到定位估計器49，得到目標(biāo)估計位置412。
權(quán)利要求
1.采用總體最小二乘和均衡算法的到達時差定位方法，其特征在于它用計算機來依次執(zhí)行以下步驟(1)，初始化，向計算機輸入以下參數(shù)值陣元位置si＝[xiyizi]T，i＝1，...，M，M為陣元數(shù)，其中第1個陣元設(shè)在X-Y-Z坐標(biāo)系的原點，觀測量ri，1，i＝2，...，M，ri，1為陣元i、1到目標(biāo)的距離差，陣元i，1到目標(biāo)的距離差的誤差方差σi12，陣元位置測量誤差方差σ2；(2)，計算機同時執(zhí)行以下步驟(2.1)，計算機根據(jù)輸入?yún)?shù)值得到的距離差觀測量ri，1，i＝2，...，M，以及陣元位置si，i＝1，...，M，構(gòu)造出系數(shù)矩陣A1和觀測向量b1，其中，A1=x2y2z2r2,1x3y3z3r3,1MMMMxMyMzMrM,1,]]>b1=0.5r2,12-R22r3,12-R32MrM,12-RM2;]]>(2.2)，計算機根據(jù)距離差的均方根誤差σi1，i＝2，...，M，按下式計算均衡誤差矩陣用的左乘加權(quán)矩陣D，其中，D=1/σ210Λ001/σ31Λ0MMM00Λ1/σM1;]]>(2.3)，計算機利用距離差誤差方差σi12、陣元位置測量誤差方差σ2和距離差觀測量ri1，計算用于均衡誤差矩陣的右乘加權(quán)矩陣T，其中，TTT=(Σidi)-1,]]>di=σ2/σi1200000σ2/σi1200000σ2/σi12000001-ri,12000-ri,10(ri,10)2,]]>i＝2，...，M；上式中，ri，10是距離差準(zhǔn)確值，可用觀測值ri，1代替；(3)，計算機對經(jīng)過均衡處理的增廣矩陣D[A1b1]T作特征值分解D[A1b1]T＝UΩVT，其中，U是(M-1)×(M-1)對角矩陣，V是5×5對角矩陣，Ω是(M-1)×5對角矩陣，其對角元素是矩陣D[A1b1]T的特征值，將矩陣V、T、Ω分塊后得到V=V11(4×4)V15(4×1)V51(1×4)V55(1×1),]]>Ω=Ω1(5×5)0(M-6×5),]]>T=T1(4×4)0(4×1)0(1×4)t5,]]>這里，Ω1＝diag(σ1σ2Λσ4σ5)，σi(i＝1...5)是矩陣D[A1b1]T的特征值；(4)，判別σ4＞σ5？若σ4＞σ5，則執(zhí)行下一步驟，若σ4＞σ5不成立，則返回步驟(2)；(5)，計算機按總體最小二乘對向量u1＝[xyr1]T進行初始估計，求 u^1=-T1V15V55t5,]]>其中r1是第1個陣元到目標(biāo)的距離；(6)，計算機根據(jù)估計值 計算觀測向量b2，以構(gòu)成觀測方程組A2u2＝b2，并計算加權(quán)矩陣W2其中，b2=u^1(1)2u^1(2)2u^1(3)2u^1(4)2,]]>A2=100010000111,]]>u2=x2y2z2,]]>W2由以下各式得出W2=(C cov(u^1)CT)-1,]]>上式中，C＝diag(x0y0z0r10)，cov(u^1)=(A1TA1-σ52I4)-1A1Tcov(b1)((A1TA1-σ52I4)-1A1T)T,]]>cov(b1)＝4diag((r2，10)2σ2，12(r3，10)2σ3，12Λ(rM-1，10)2σM-1，12)，這里，I4表示4階單位陣，[x0y0z0]表示目標(biāo)的真實位置，r10表示陣元1到目標(biāo)的準(zhǔn)確距離，它們可分別由 中元素 和 代替，而矩陣cov(b1)中的ri，10，i＝2，...，M可由觀測值ri，1代替；(8)，計算機按加權(quán)最小二乘對向量u2進行估計，求 u^2=(A2TW2A2)-1A2TW2b2;]]>(9)，計算機按下式輸出目標(biāo)位置u＝[xyz]T的估計值 u^=G[u^2(1)0.5,u^2(2)0.5,u^2(3)0.5],]]>其中，G=diag(sign(u^1(1)),sign(u^1(2)),sign(u^1(3))).]]>
2.采用總體最小二乘和均衡算法的到達時差定位方法，其特征在于，它依次含有以下步驟(1)，輸入以下測定的初始參數(shù)值陣元位置si＝[xiyizi]T，i＝1，...，M，M為陣元數(shù)，其中第1個陣元設(shè)在X-Y-Z坐標(biāo)系的原點，觀測量ri，1，i＝2，...，M，ri，1為陣元i、1到目標(biāo)的距離差，陣元i，1到目標(biāo)的距離差的誤差方差σi12，陣元位置測量誤差方差σ2；(2)，把下列數(shù)組通過各自的運算電路運算后，其結(jié)果送往相應(yīng)的寄存器組數(shù)組ri，12，σi，12經(jīng)過乘法器運算后，得到cov(b1)，送往cov(b1)寄存器組，si經(jīng)過平方和電路得到Ri2，再和ri，12經(jīng)過減法器得到b1，送往b1寄存器組，si，ri，1送往A1寄存器組，σi，1經(jīng)過倒數(shù)電路得到1/σi，1，送往D寄存器組，σi12、σ2經(jīng)過除法電路得到σ2/σi12，再把σ2/σi12和ri，1一起送往di寄存器組，得到di，然后對di求和，經(jīng)求逆、分解電路后形成T，送往T寄存器組；(3)，把上述各寄存器組的輸出A1、b1、D、T送往矩陣乘法電路得到D[A1b1]T，再把D[A1b1]T經(jīng)過公知的特征值分解電路得到U、V、Ω，并送往U、V、Ω寄存器組；(4)，把U、V、Ω寄存器組中的V和T寄存器組中的T，經(jīng)過矩陣乘法電路，得到 和C，并送往 寄存器組；(5)，把從相應(yīng)的寄存器組得到的cov(b1)、A1、σ5和C送往矩陣乘法電路，求出W2，并送往W2寄存器組；(6)，使 經(jīng)平方電路得到b2，把b2、W2送往矩陣乘法電路得到 再把 經(jīng)過開方電路得到
全文摘要
采用總體最小二乘和均衡算法的到達時差定位方法，屬于空間目標(biāo)定位技術(shù)領(lǐng)域，其特征在于它采用對誤差矩陣進行預(yù)均衡并結(jié)合總體最小二乘法來定位，以便利用均衡后的誤差矩陣的協(xié)方差矩陣近似為單位陣的特性，來提高定位精度。它所處理的誤差模型既包括有陣元位置測量誤差，又包括有互不相等的距離差估計誤差，因而更符合實際情況。在距離差估計誤差不相等及存在陣元位置測量誤差時，相比于已有定位方法具有更高的定位精度。
文檔編號G01S7/00GK1553214SQ20031012170
公開日2004年12月8日 申請日期2003年12月19日 優(yōu)先權(quán)日2003年12月19日
發(fā)明者黃振, 陸建華, 黃 振 申請人:清華大學(xué)