专利名称：采用总体最小二乘和均衡算法的到达时差定位方法
技术领域：
采用总体最小二乘和均衡算法的到达时差定位方法属于空间目标定位技术领域。
背景技术：
目标的定位精度除了与目标相对于各接收阵元的空间几何位置有关外，还取决于定位参数的估计误差以及阵元位置的测量误差。在实际应用中，受信噪比的影响和具体阵元位置测量技术(如全球定位系统)的限制，这些误差难以消除，且误差方差也不尽相同。如果采用的定位算法不能将这些误差因素考虑进去，则实际定位精度将受到严重影响。
TDOA定位技术具有较强的工程可实现性和较高的定位精度，因此被广泛运用于雷达、声纳及辐射源定位系统中。
TDOA定位方法的原理1、常规的TDOA定位方法如图1所示，设定三维空间中任意分布着M个信号接收阵元，不失一般性，假定阵元1处于坐标原点，各阵元的位置分别为si＝[xiyizi]T(i＝1，...，M)，目标的位置为u＝[xyz]T，这里[]T表示矩阵转置。各阵元到原点的距离为Ri＝‖si‖，‖‖表示Euclidean范数。这样，阵元i到目标的距离就是ri＝‖si-u‖，阵元i、j到目标的距离差(RD)就是ri，j＝ri-rj＝‖si-u‖-‖sj-u‖，(i，j＝1，...，M)(1)因为距离差除以信号传播速度就是目标信号到达两阵元的时间差，因此不妨将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]]>当系数矩阵A1的秩等于未知数个数n，且M≥n时，可按照经典的最小二乘算法求出目标的估计位置
u^1=(A1TA1)-1A1Tb1---(4)]]>当只有观测向量b1含有误差，且各误差分量相互独立、方差相同，则最小二乘估计在方差最小的意义上是最佳的。但在实际应用中，系数矩阵A1往往也有扰动，并且误差的协方差矩阵不为单位阵(误差方差不同或误差分量具有相关性)，因此最小二乘估计法的实用性受到很大限制。
另一种常规方法是用泰勒展式(Taylor-series)将非线性观测方程((1)式)转化为线性方程，再通过迭代的方式估计未知向量。这种方法的主要缺点是需要设定初始值，而且因为迭代使得方法的实时性较差。
2、现有的TDOA定位方法主要有两种方法一是用总体最小二乘(TLS)估计未知向量，这样不仅考虑了观测向量的扰动，也考虑了系数矩阵的扰动，使得定位精度得以提高。在这种方法中一般采用特征值分解(SVD)求解未知向量，因而又具有较高的数值稳定性。在误差分量的协方差阵为单位阵的情况下，TLS具有最佳的估计效果。但实际中，系数矩阵元素往往是与观测向量不同量纲的各种测量值，因此这个条件很难满足。方法二是改进的加权最小二乘法，此方法在构造误差矩阵时，不仅考虑了RD对观测向量引起的扰动，同时也考虑了RD对系数矩阵引起的扰动，并将误差矩阵的协方差阵进行求逆作为加权矩阵。不仅如此，该方法还考虑到未知向量u1中各元素的关系，即x2+y2+z2＝r12，并利用这一固有关系进一步提高了定位精度。方法二的最大缺点是没有考虑阵元位置测量误差，因为考虑该误差时难以构造出加权矩阵。在实际定位系统中大量存在的情况是不仅RD的估计具有误差，阵元位置的测量也有误差，并且误差的方差也是不同的(由信噪比不同等因素引起)。

发明内容
本发明的目的在于克服已有技术的不足之处，针对更接近实际的定位系统误差模型，提出一种采用对误差矩阵进行预均衡并结合TLS算法的定位方法，利用均衡后的误差矩阵的协方差阵为单位阵的特性，提高定位精度。
当TDOA估计和阵元位置测量都存在误差时，由(3)式可知，系数矩阵A1及观测向量b1都存在扰动。因此，增广矩阵[A1b1]可表示为由真值组成的矩阵A0与误差矩阵E之和，即[A1b1]＝[A10b10]+[e′e]＝A0+E。设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和阵元位置的准确值，Δri，1和[ΔxiΔyiΔzi]T分别表示RD和阵元位置的误差，并假定为相互独立的零均值高斯白噪声随机变量。若误差矩阵不满足独立同分布结构(IID)，例如各元素的方差不等或元素间不相互独立等，则已有的TLS算法就不能达到最佳估计效果。因此，将传统的TLS算法替代为求解未知向量使之满足minimize‖DET‖F(6)subject to b1+e∈Range(A1+e′)式中D、T分别为起均衡误差矩阵作用的左乘加权矩阵和右乘加权矩阵，‖‖F表示Frobenius范数。左乘加权矩阵D是对不同阵元的定位参数估计误差进行加权均衡处理；右乘加权矩阵T是对不同误差(如阵元位置测量误差、距离差估计误差)进行加权均衡处理。由此可见，此方法比WLS和TLS具有更大的灵活性。
本发明提出的采用总体最小二乘和均衡算法的到达时差定位方法，如图2所示，其特征在于它用计算机来依次执行以下步骤定位计算开始后，计算机从步骤2a同时转移到并行执行的步骤2b、2c和2d。其中步骤2b是根据前端TDOA估计和阵元位置测量分别提供的观测量(距离差)ri，1(i＝2，...，M)和阵元位置矢量si(i＝1，...，M)构造出系数矩阵A1及观测向量b1，这里假定阵元数目为M。步骤2c是根据距离差估计均方根误差σi1计算用于均衡误差矩阵的左乘加权矩阵DD=1/σ210Λ001/σ31Λ0MMM00Λ1/σM1---(7)]]>矩阵D是对涉及不同阵元的阵元位置测量误差、距离差估计误差进行加权均衡处理。步骤2d是利用距离差估计方差σi12、阵元位置测量误差方差σ2和距离差观测量ri，1，计算用于均衡误差矩阵的右乘加权矩阵T，该矩阵是对不同误差(如阵元位置测量误差、距离差估计误差)进行加权均衡处理。在这一步骤中，首先得到矩阵E′＝DE的行向量的协方差矩阵为di=σ2/σi1200000σ2/σi1200000σ2/σi121-ri,10000-ri,10(ri,10)2(i=2,...,M)---(8)]]>于是可以由下式计算右加权矩阵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，通过仿真可知，由以上近似带来的误差很小，可以忽略不记。
完成以上步骤后，估计器转移到步骤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)]]>式中分块矩阵的下标(i×j)表示该分块矩阵是由原矩阵的1～i行及1～j列所对应元素组成，这里Ω1＝diag(σ1σ2Λσ4σ5)，σi(i＝1...5)是矩阵D[A1b1]T的特征值。
随后估计器转移到步骤2f，比较Ω1中按降序排列的最后两个特征值的大小，若前一个特征值σ4大于最小特征值σ5，则估计器转移到步骤2g；反之，转回到步骤2b、2c、2d，从新开始新一轮的定位计算。步骤2g利用特征值分解结果计算未知向量u1的估计值 即得到目标位置的初始估计值u^1=-T1V15V55t5---(12)]]> 还可以表示为另一种等价形式u^1=(A1TA1-σ52I4)-1A1Tb1---(13)]]>其中I4表示4阶单位阵。
接下来估计器转移到步骤2h，根据u1的估计值构成观测向量b2=u^12,]]>并计算加权矩阵W2。显然向量u1中的4个元素并非相互独立的，x，y，z，和r1满足关系式r12＝x2+y2+z2(14)这一先验信息可以用来进一步提高目标位置的估计精度。根据(14)式建立以下观测方程组A2u2＝b2(15)
其中A2=100010001111---b2=u^1(1)2u^1(2)2u^1(3)2u^1(4)2---u2=x2y2z2]]>由于系数矩阵A2是常系数矩阵，因此误差矩阵ε由估计值 的各误差分量组成。根据误差矩阵与加权矩阵的关系W2=E(ϵϵT)-1|u^1=u10]]>(u10为向量u1的真值)，并忽略二阶误差项，可以得到加权矩阵W2的表达式W2=(Ccov(u^1)CT)-1---(16)]]>其中C＝diag(x0y0z0r10)这里，[x0y0z0]表示目标的真实位置。而(17)式中的 的协方差矩阵可由(13)式得到cov(u^1)=(A1TA1-σ52I4)-1A1Tcov(b1)((A1TA1-σ52I4)-1A1T)T---(17)]]>同样忽略掉高阶误差项，观测向量b1的协方差矩阵为cov(b1)＝4diag((r2，10)2σ2，12(r3，10)2σ3，12Λ(rM-1，10)2σM-1，12)(18)实际计算中，目标的真实位置[x0y0z0]可由 中的对应元素代替，而矩阵cov(b1)中的ri，10(i＝2，...,M)可由观测值ri，1代替。
完成以上计算后，估计器转移到步骤2i，应用加权最小二乘计算估计值 u^2=(A2TW2A2)-1A2TW2b2----(19)]]>接着估计器转到步骤2j，对前一步的输出进行开方运算并乘上符号矩阵G，从而得到目标位置的估计值 u^=G[u^2(1)0.5,u^2(2)0.5,u^2(3)0.5]---(20)]]>其中G用于确定 中各元素的符号，可由初始估计值 中相应元素的符号得知，即G=diag(sign(u^1(1)),sign(u^1(2)),sign(u^1(3)))---(21)]]>式中sign表示符号函数(自变量为正取1，为负取-1)。最后估计器转移到步骤2k，定位计算结束。
本发明的特征还在于它可以通过一个由若干乘法器组成的电路完成，并依次包含以下步骤
(1)，输入以下测定的初始参数值阵元位置si＝[xiyizi]T，i＝1，...，M，M为阵元数，其中第1个阵元设在X-Y-Z坐标系的原点，观测量ri，1，i＝2，...，M，ri，1为阵元i、1到目标的距离差，阵元i，1到目标的距离差的误差方差σi12，阵元位置测量误差方差σ2；(2)，把下列数组通过各自的运算电路运算后，其结果送往相应的寄存器组数组ri，12，σi，12经过乘法器运算后，得到cov(b1)，送往cov(b1)寄存器组，si经过平方和电路得到Ri2，再和ri，12经过减法器得到b1，送往b1寄存器组，把si，ri，1送往A1寄存器组，σi，1经过倒数电路得到1/σi，1，送往D寄存器组，σi12、σ2经过除法电路得到σ2/σi12，再把σ2/σi12和ri，1一起送往di寄存器组，得到di，然后对di求和，经求逆、分解电路后形成T，送往T寄存器组；(3)，把上述各寄存器组的输出A1、b1、D、T送往矩阵乘法电路得到D[A1b1]T，再把D[A1b1]T经过公知的特征值分解电路得到U、V、Ω，并送往U、V、Ω寄存器组；(4)，把U、V、Ω寄存器组中的V和T寄存器组中的T，经过矩阵乘法电路，得到 和C，并送往 寄存器组；(5)，把从相应的寄存器组得到的cov(b1)、A1、σ5和C送往矩阵乘法电路，求出W2，并送往W2寄存器组(6)，使 经平方电路得到b2，把b2、W2送往矩阵乘法电路得到 再把 经过开方电路得到 作为一个例子，设定阵元数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估计及阵元位置测量精度，分别设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)仿真结果(如表1所示)，这三种方法分别是改进的WLS定位方法，TLS定位方法和本发明提出的定位方法。
表1随阵元位置测量误差方差变化的定位均方根误差仿真结果 可见，本发明提出的定位方法具有三种方法中最高的定位精度。本发明的效果是，利用加权矩阵对误差矩阵作均衡处理，使得均衡后的误差矩阵的协方差矩阵近似为单位阵，并通过具有较高数值稳定性的特征值分解方法得到目标位置的估计值，从而大幅度提高了目标的定位精度。所处理的误差模型既包括有阵元位置测量误差，又包括互不相等的距离差误差，因而更符合实际应用情况。比之于已有方法一，在距离差估计误差不相等时，本方法具有更高的定位精度；比之于已有方法二，在存在阵元位置测量误差时，本方法具有更好的定位性能。


图1为定位系统目标和阵元位置关系示意2为本发明的软件实现流程图。
图3为按照本发明设计的定位估计器的硬件实现方框图。
图4为应用本发明方法的目标定位系统框图。
具体实施例方式
下面根据附图和两个实施例更加详细地解释本发明实施例一本实施例采用软件实现本发明提出的定位方法，如图2所示。设阵元位置数为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估计及阵元位置测量精度，设RD的方差σi12(i＝2，...，9)等于10-5[1 4 9 16 25 36 49 64]，阵元位置测量误差方差σ2等于10-3。定位计算开始后，估计器从步骤2a同时转移到并行执行的步骤2b、2c和2d。其中步骤2b是根据前端TDOA估计和阵元位置测量分别提供的观测量ri，1(i＝2，...，9)和阵元位置矢量si(i＝1，...，9)构造出系数矩阵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是根据观测量估计方差σi1用(7)式计算用于均衡误差矩阵的左乘加权矩阵D，该矩阵是对不同阵元的定位参数估计误差进行加权均衡处理，其值为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)式计算用于均衡误差矩阵的右乘加权矩阵T，该矩阵是对不同误差(如阵元位置测量误差、观测量估计误差)进行加权均衡处理，其值为T＝0.0809 0 0 0 00 0.0809 0 0 00 00.0809 0 00 0 00.9256 00 0 0 00.0197完成以上步骤后，估计器转移到步骤2e，实现矩阵的乘法运算并进而完成乘积矩阵的特征值分解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其中Ω包含有乘积矩阵的特征值，并按降序排列在对角线上。然后估计器转移到步骤2f，比较按降序排列的最后两个特征值的大小。由Ω的值可知σ4＝0.0534，σ5＝0.0001，显然满足σ4＞σ5，因此估计器转移到步骤2g。步骤2g利用(12)式计算未知向量u1的估计值 即得到目标位置的初始估计值u^1=1.0e+003*0.29152.33850.19462.3637]]>随后，估计器转移到步骤2h，根据u1的估计值构成观测向量b2=u^12,]]>并通过(16)、(17)和(18)式计算出加权矩阵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完成这一步后，估计器转移到步骤2i，应用(19)式计算估计值 u^2=1.0e+006*0.09105.84350.0406]]>接着估计器转到步骤2j，对前一步的输出进行开方运算并乘上符号矩阵G，从而得到目标位置的估计值 u^=1.0e+003*2.30162.41730.2015]]>这里G用于确定 中各元素的符号，可由初始估计值 中相应元素的符号得知。最后估计器转移到步骤2k，定位计算结束。
实施例二本实施例为采用硬件的方法实现本发明提出的定位方法。仍然以阵元位置数为9，观测量估计误差方差σi12(i＝2，...，9)互不相等，且阵元位置测量也存在误差的情况为例，说明该硬件电路的定位过程。如图3所示，由前端TDOA估计和阵元位置测量提供的第1到第8路估计参数输入31，经过简单的运算电路32～39(图中|·|2表示平方和电路，(·)-1表示求倒数电路)，得到系数矩阵313(A1)、观测矩阵314(b1)、左加权矩阵315(D)和右加权矩阵316(T)的各组成元素。其中各路输入的距离差的平方值ri，12与距离差误差方差σi，12的乘积被存储在一个寄存器组中，为后级提供观测向量的协方差矩阵312(cov(b1))；各阵元位置矢量si(i＝1，...，9)和距离差ri，1被存储到一个寄存器组，为后级提供系数矩阵313；求出各阵元到原点的距离Ri2并与距离差的平方值相减，从而得到观测向量314的各分量值，并被存储在一个寄存器组中；距离差误差方差σi，1被取倒数后，作为左加权矩阵的对角元素存储在一个寄存器组中；各路输入的距离差ri，12及其误差方差σi12，以及阵元位置测量误差方差σ2各自构成矩阵di的元素，并通过310电路——包括矩阵相加、求逆和Cholesky分解器，得到右加权矩阵316，并被存储在一个寄存器组中。以上5个寄存器组311提供的矩阵被送到矩阵乘法电路317中，求出乘积矩阵318(D[A1b1]T)的值。随后对该矩阵进行特征值分解319，并将酉阵U、V和特征值组成的对角阵320存储在寄存器组321中，其中矩阵322(V)的分块矩阵V15和V55在后级矩阵乘法电路324中与右加权矩阵316的分块矩阵T1和t5作乘法运算，得到阵元位置的初始估计值 和由其相应元素组成的矩阵328(C)。其中特征值分解电路319由黄浴等人提出的Systolic阵列实现(“一个实现状态空间法的Systolic阵列系统设计”，通信学报，黄浴等，Vol.16，No.6，pp.25～32，1995)。矩阵328和前面得到的协方差矩阵312、系数矩阵313以及最小特征值σ5一起被送到矩阵乘法电路325中，计算出用于加权最小二乘运算的加权矩阵332(W2)。初始估计值327经过平方电路329得到观测矩阵330(b2)，并与加权矩阵332一起送入到矩阵乘法电路333，求出估计值 再经过平方根电路335输出目标位置估计值 从上面的实施例可见，本定位方法与其前面介绍过的两种现有的定位方法相比，具有更高的定位精度以及更强的实用性，而且其硬件实现也并不复杂。下面再通过一个例子来说明本发明在辐射源定位系统中的作用。
参照图4，采用总体最小二乘及均衡算法的目标定位系统包括接收信号的天线阵41、接收前端42、模数转换(A/D)、TDOA(或RD)估计47、阵元位置测量410、以及图3所示的定位估计器49。在本例中，天线阵元接收的各路信号42被分别送往各自的接收前端43，并转化为中频信号44。模数转换45将模拟中频信号变为数字信号流46，输给下一级TDOA估计47。TDOA估计47由此计算出各路信号之间的相对延迟，从而得到观测量48用于定位估计。阵元位置测量410输出各阵元位置411，与观测量48一起输入到定位估计器49，得到目标估计位置412。
权利要求
1.采用总体最小二乘和均衡算法的到达时差定位方法，其特征在于它用计算机来依次执行以下步骤(1)，初始化，向计算机输入以下参数值阵元位置si＝[xiyizi]T，i＝1，...，M，M为阵元数，其中第1个阵元设在X-Y-Z坐标系的原点，观测量ri，1，i＝2，...，M，ri，1为阵元i、1到目标的距离差，阵元i，1到目标的距离差的误差方差σi12，阵元位置测量误差方差σ2；(2)，计算机同时执行以下步骤(2.1)，计算机根据输入参数值得到的距离差观测量ri，1，i＝2，...，M，以及阵元位置si，i＝1，...，M，构造出系数矩阵A1和观测向量b1，其中，A1=x2y2z2r2,1x3y3z3r3,1MMMMxMyMzMrM,1,]]>b1=0.5r2,12-R22r3,12-R32MrM,12-RM2;]]>(2.2)，计算机根据距离差的均方根误差σi1，i＝2，...，M，按下式计算均衡误差矩阵用的左乘加权矩阵D，其中，D=1/σ210Λ001/σ31Λ0MMM00Λ1/σM1;]]>(2.3)，计算机利用距离差误差方差σi12、阵元位置测量误差方差σ2和距离差观测量ri1，计算用于均衡误差矩阵的右乘加权矩阵T，其中，TTT=(Σidi)-1,]]>di=σ2/σi1200000σ2/σi1200000σ2/σi12000001-ri,12000-ri,10(ri,10)2,]]>i＝2，...，M；上式中，ri，10是距离差准确值，可用观测值ri，1代替；(3)，计算机对经过均衡处理的增广矩阵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，则执行下一步骤，若σ4＞σ5不成立，则返回步骤(2)；(5)，计算机按总体最小二乘对向量u1＝[xyr1]T进行初始估计，求 u^1=-T1V15V55t5,]]>其中r1是第1个阵元到目标的距离；(6)，计算机根据估计值 计算观测向量b2，以构成观测方程组A2u2＝b2，并计算加权矩阵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]表示目标的真实位置，r10表示阵元1到目标的准确距离，它们可分别由 中元素 和 代替，而矩阵cov(b1)中的ri，10，i＝2，...，M可由观测值ri，1代替；(8)，计算机按加权最小二乘对向量u2进行估计，求 u^2=(A2TW2A2)-1A2TW2b2;]]>(9)，计算机按下式输出目标位置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)，输入以下测定的初始参数值阵元位置si＝[xiyizi]T，i＝1，...，M，M为阵元数，其中第1个阵元设在X-Y-Z坐标系的原点，观测量ri，1，i＝2，...，M，ri，1为阵元i、1到目标的距离差，阵元i，1到目标的距离差的误差方差σi12，阵元位置测量误差方差σ2；(2)，把下列数组通过各自的运算电路运算后，其结果送往相应的寄存器组数组ri，12，σi，12经过乘法器运算后，得到cov(b1)，送往cov(b1)寄存器组，si经过平方和电路得到Ri2，再和ri，12经过减法器得到b1，送往b1寄存器组，si，ri，1送往A1寄存器组，σi，1经过倒数电路得到1/σi，1，送往D寄存器组，σi12、σ2经过除法电路得到σ2/σi12，再把σ2/σi12和ri，1一起送往di寄存器组，得到di，然后对di求和，经求逆、分解电路后形成T，送往T寄存器组；(3)，把上述各寄存器组的输出A1、b1、D、T送往矩阵乘法电路得到D[A1b1]T，再把D[A1b1]T经过公知的特征值分解电路得到U、V、Ω，并送往U、V、Ω寄存器组；(4)，把U、V、Ω寄存器组中的V和T寄存器组中的T，经过矩阵乘法电路，得到 和C，并送往 寄存器组；(5)，把从相应的寄存器组得到的cov(b1)、A1、σ5和C送往矩阵乘法电路，求出W2，并送往W2寄存器组；(6)，使 经平方电路得到b2，把b2、W2送往矩阵乘法电路得到 再把 经过开方电路得到
全文摘要
采用总体最小二乘和均衡算法的到达时差定位方法，属于空间目标定位技术领域，其特征在于它采用对误差矩阵进行预均衡并结合总体最小二乘法来定位，以便利用均衡后的误差矩阵的协方差矩阵近似为单位阵的特性，来提高定位精度。它所处理的误差模型既包括有阵元位置测量误差，又包括有互不相等的距离差估计误差，因而更符合实际情况。在距离差估计误差不相等及存在阵元位置测量误差时，相比于已有定位方法具有更高的定位精度。
文档编号G01S7/00GK1553214SQ20031012170
公开日2004年12月8日 申请日期2003年12月19日 优先权日2003年12月19日
发明者黄振, 陆建华, 黄 振 申请人:清华大学