CN105425231A - 一种基于分层投影和泰勒展开的多传感器多目标定位方法 - Google Patents

Abstract

本发明提供了一种基于分层投影和泰勒展开的多传感器多目标定位方法，它是基于双站距离空间投影定位方法，首先以较大的分辨率对观测区域进行栅格划分，为每个单元分配代表点，用代表点到传感器的距离提取回波并将回波累积到对应的像素单元，利用双站距离空间投影定位得到粗糙的目标位置，再选取目标位置附近一定区域作为投影区域，然后以更小的分辨率划分观测区域，只在投影区域进行投影定位，分多层完成定位，在得到粗略的目标位置基础上对目标附近区域进行高精度成像，大大减少了投影定位算法的运算量，最后采用泰勒展开算法，对估计的目标位置进行修正，进一步提高了多目标定位的精度。

Description

一种基于分层投影和泰勒展开的多传感器多目标定位方法

技术领域

本发明属于多基地雷达系统中的多目标定位的技术领域，它特别涉及到传感器网络下多目标定位的技术领域。

背景技术

与单双站雷达相比，传感器网络可以从多个方位获取目标的多个延时信息，在侦察对抗领域，尤其是对于弱目标的探测具有很大的优势。利用这些延时构建非线性距离方程组，然后解算方程组，就可以在一个脉冲周期得到目标的位置。但是在多目标定位中，必须从各个传感器的多个延时中找到每个目标对应的延时，即完成回波与目标的关联，从而建立正确的距离方程组，这是一个非常复杂的任务。传统的时差定位方法是建立和求解所有可能的方程组，存在运算量大，不容易收敛以及如何准确提取目标的问题。基于成像策略的多目标定位算法(如双站距离空间投影定位算法)，利用成像算法，对空中目标进行探测与定位，不需要建立相应的时差方程组，因此可以解决回波关联，方程组解算复杂以及极易不收敛等问题，但是由于成像算法需要遍历空间的每一个像素点，以及空间分辨率不高，使得存在运算量大及定位精度有待提高的问题(详见“ShiTongyan,ZhangXiaoling,ShiJun.Multi-targetPositioningforSensorNetworkBasedonImagingStrategy,ComputationalProblem-solving(ICCP),2013IEEEInternationalConference,Oct.2013,pp.171-174”)。

发明内容：

为了克服双站距离空间投影定位算法的运算量大和定位精度有待提高的问题，本发明提出了一种基于分层投影和泰勒展开的多传感器多目标定位方法，其特点是先以较大的分辨率对观测区域进行栅格划分，利用双站距离空间投影定位得到粗糙的目标位置，再选取目标位置附近一定区域作为投影区域，然后以更小的分辨率划分观测区域，只在投影区域进行投影定位，分多层完成定位，在得到粗略的目标位置基础上对目标附近区域进行高精度成像，大大减少了投影定位算法的运算量，最后采用泰勒展开算法，对估计的目标位置进行修正，进一步提高了多目标定位的精度。

为了方便描述本发明的内容，作以下术语定义：

定义1、双站延时

从发射机发射信号时开始计时，当信号经由目标散射，并由接收机接收到信号时终止计时，这段时间即称为该目标对相应接收机的双站延时。双站延时可由GPS等时间同步系统得到。

定义2、双站距离

发射机到目标的距离与该目标到接收机的距离之和，即双站延时乘以光速所得到的距离，称为双站距离。

定义3、双站距离方程

由原理图附图1可知，目标位置为x，接收机数目为N，接收机位置为r_i(i＝1,2,…,N)，发射机位置为T，则信号从发射机经由目标到达第i个接收机的距离为

‖T-x‖₂+‖x-r_i‖₂＝v*τ_i

式中，v为光速，*为乘号，τ_i为信号从发射机经由目标到达第i个接收机的双站延时，‖·‖₂表示向量的2范数，i＝1,2,…,N，上式称为关于该目标的双站距离方程。

定义4、距离单元编号

在定位系统中，将雷达探测的距离划分成若干较小的单元并将其编号，雷达根据目标回波所在单元的编号计算目标与雷达间的距离。定位系统根据接收机对各目标的双站距离，结合距离向分辨率，可以得到各目标回波的量化编号，称为距离单元编号，量化时会进行取整，下文用round(*)表示对*进行四舍五入操作。

定义5、接收距离波门

信号从发射到接收的最短时间长度，雷达在大于这一时间长度时才对数据进行接收，这一最短时间长度即为接收时间波门。与接收时间波门对应的为接收距离波门，即接收时间波门乘以光速所得到的距离。

定义6、观测空间

观测空间指目标、杂波及噪声等所在的实际的三维场景空间。

定义7、笛卡尔坐标系

笛卡儿坐标系(Cartesian坐标系)，是一种正交的直角坐标系。二维笛卡儿坐标系由两条相互垂直的数轴构成，交点为0点。三维笛卡儿坐标系在二维笛卡儿坐标系的基础上根据右手定则增加第三维坐标而形成。

定义8、回波关联

回波关联指在多目标多雷达的定位中，需要分清各雷达的回波信号哪些是来源于同一目标的，并把属于同一目标的回波信号组合起来。

定义9、信噪比

信号功率与环境噪声功率的比值称为信噪比。详见文献“丁鹭飞，耿富录.雷达原理(第三版).西安电子科技大学出版社.2009.8”。

定义10、单载频脉冲信号

单载频脉冲信号可以表示为其中f_c表示脉冲的载频，T_P表示脉冲宽度。详见“皮亦鸣,杨建宇,付毓生,杨晓波.合成孔径雷达成像原理.第一版.电子科技大学出版社.2007.3”。

定义11、回波概率信息

回波概率信息是对回波采样值的幅度进行加权，将幅度变为0到1之间，表征该回波数据为目标回波的概率。

定义12、概率积累

概率积累指在雷达系统将来自同一目标的回波概率信息进行累加以提高信噪比，增加信息量的过程。

定义13、零矩阵

零矩阵是指矩阵元素全为0的矩阵。

定义14、标准差

标准差是指总体各单位标准值与总体平均数离差平方的算术平均数的平方根。

本发明提出了一种基于分层投影和泰勒展开的多传感器多目标定位方法，它包括如下步骤：

步骤1、相关参数的初始化

初始化的参数均为已知，如下所示：电磁波传播速度，记为v；雷达空间距离分辨率，记为δ_r；雷达在距离向上的最大距离单元编号，记为整数L；以发射机所在位置为原点，建立三维笛卡尔坐标系；发射机位置，记为P_t；目标数目，记为Ntar；第i个目标的位置，记为Ptar(i)(i＝1,2,…,Ntar)，i为目标序号；目标散射截面积，记为σ；接收机数目，记为Nrec；第j个地面接收机位置，记为Prec(j)(j＝1,2,…,Nrec)，j为接收机序号；接收机的接收距离波门，均记为Rrec0；各接收机通道的信噪比，记为SNR；以单载频脉冲信号为雷达发射信号，脉冲载频记为f_c，脉冲宽度记为T_p；信号经发射机发射被第i个目标散射，由第j个接收机接收，再经去载频、离散采样后的回波信号记为Echo(i,j)，其表达式Echo(i,j)＝σ_i,j*exp(-pkR_i,j)，R_i,j表示第i个目标与第j个接收机的双站距离，表达式R_i,j＝‖P_t-Ptar(i)‖₂+‖Ptar(i)-Prec(j)‖₂，其中exp(·)是以自然指数为底的指数函数，σ_i,j为第i个目标相对于第j个接收机的目标散射截面积，p表示虚数单位即p²＝-1，‖·‖₂表示向量的L2范数，*表示乘号，π为圆周率，v表示光速。

步骤2、创建回波存储矩阵并求取回波概率信息矩阵

由发射机和接收机的时间同步系统获取空中各目标对不同接收机的回波信号Echo(i,j)以及相应的双站延时τ_i,j，其中i＝1,2,…,Ntar，j＝1,2,…,Nrec，τ_i,j表示第i个目标相对第j个接收机的双站延时，Echo(i,j)表示第j个接收机接收到的第i个目标散射的回波信号。

步骤2.1初始化回波存储矩阵

根据步骤1中雷达在距离向上的最大距离单元编号L以及地面接收机数目Nrec，定义一个L行Nrec列的回波存储矩阵，记为并将回波存储矩阵初始化为零矩阵。回波存储矩阵的每一元素，记为l＝1,2,…,L，c＝1,2,…,Nrec，其中l为回波存储矩阵的行数，c为回波存储矩阵的列数；

根据第j个接收机获得的第i个目标的双站延时τ_i,j及对应的目标回波信号Echo(i,j)，计算第i个目标对第j个接收机的双站距离为R_i,j＝v*τ_i,j，其中v表示光速，*表示乘号，i＝1,2,…Ntar；j＝1,2,…Nrec；

根据双站距离R_i,j，计算相应回波的距离单元编号ID_i,j，公式为：其中ID_i,j表示第i个目标的回波在第j个接收机中的距离单元编号，1≤ID_i,j≤L，且ID_i,j取值为正整数，L为最大的距离单元编号，Rrec0表示接收机的接收距离波门；

然后将经发射机发射被第i个目标散射，由第j个接收机接收的回波信号Echo(i,j)，存储在回波存储矩阵中的位置处，即由 ${\mathrm{Echo}}_L^{Nrec}(l,c)=Echo(i,j)$ (l＝ID_i,j,c＝j)，其中l＝ID_i,j，c＝j，i＝1,2,…Ntar，j＝1,2,…Nrec，得到回波存储矩阵

步骤2.2回波存储矩阵转换为回波概率信息矩阵

根据得到的存储矩阵求回波概率信息矩阵P，回波概率信息矩阵P的维数与相同，回波概率信息矩阵P计算公式为：P(l,c)＝1-exp(-o(l,c)²/2σ(c))，其中l＝1,2,…,L，c＝1,2,…,Nrec，o(l,c)是的幅度，σ(c)是矩阵第c列的幅度标准差。

步骤3、栅格划分

步骤3.1建立观测空间

在步骤1建立的笛卡尔坐标系中，设定一个三维长方体观测空间Ω_Sear，观测空间的长、宽、高分别平行于笛卡尔坐标系的x轴、y轴、z轴。观测空间长为M×δ_r×Nx,宽为M×δ_r×Ny,高为M×δ_r×Nz，其中Nx，Ny，Nz分别为沿观测空间x轴，y轴，z轴的采样点数，δ_r表示雷达空间距离分辨率，M取值为2的非负整数次方，表示栅格划分边长相对于雷达空间距离分辨率δ_r的倍数。

步骤3.2等间隔划分栅格

将步骤3.1观测空间Ω_Sear的中心记为P0，步骤1中的雷达空间距离分辨率δ_r作为最小单元，将观测空间划分为一组栅格，每一栅格单元大小为Mδ_r×Mδ_r×Mδ_r，然后为每一栅格单元标记一代表点，代表点坐标为P_xyz＝P0+[x-Nx/2,y-Ny/2,z-Nz/2]*M*δ_r，x,y,z分别是观测空间沿x轴，y轴，z轴的采样点序号，x＝1,2…Nx，y＝1,2…Ny，z＝1,2…Nz，P_xyz表示采样点序号为x,y,z的代表点的三维坐标，相邻代表点间的距离为Mδ_r。

步骤4、回波降分辨率处理

根据步骤2得到的回波概率信息矩阵P，采用降分辨率回波信息矩阵元素的计算公式 $P_M(k,c)=1-{\Π}_{i=1}^M(1-P(\left(k-1\right)\×M+i,c\left)\right),$ 得到降分辨率回波信息矩阵P_M，其中取值为小于等于其自身的最大正整数，c＝1,2,…,Nrec，其中M表示栅格划分边长相对于空间分辨率δ_r的倍数，L表示最大的距离单元编号，M取值为2的非负整数次方，表示栅格划分边长相对于雷达空间距离分辨率δ_r的倍数。

步骤5、投影处理

在步骤3建立的观测空间Ω_Sear中，寻找每一栅格代表点对应各接收机的回波信号，采用传统的投影方法对每一栅格代表点对应的各接收机的回波信号进行积累；具体步骤如下：

步骤5.1

首先创建一个三维图像矩阵Ima_Nx×Ny×Nz，矩阵的大小为Nx×Ny×Nz，并初始化为零矩阵，用于存储观测空间的成像信息；三维图像矩阵元素单元记为Ima_Nx×Ny×Nz(k1,k2,k3)，k1、k2和k3分别表示三维图像矩阵Ima_Nx×Ny×Nz的行号、列号、层号，其中k1＝1,2,…Nx，k2＝1,2,…Ny，k3＝1,2,…Nz，Nx，Ny，Nz分别是三维图像矩阵的行数，列数，层数；然后分别初始化观测空间沿x轴，y轴，z轴的采样点序号x＝1,y＝1,z＝1，接收机序号初始化为j＝1；转到步骤5.2；

步骤5.2

根据步骤3中栅格代表点的三维坐标P_xyz＝P0+[x-Nx/2,y-Ny/2,z-Nz/2]*M*δ_r和接收机j的位置Prec(j)，计算从发射机位置P_t经栅格代表点P_xyz到达第j个接收机的双站距离，记为R_xyz,j，其计算公式为：R_xyz,j＝‖P_t-P_xyz‖₂+‖P_xyz-Prec(j)‖₂；

由双站距离R_xyz,j计算栅格单元代表点P_xyz相对接收机j的距离单元编号，记为ID_xyz,j，其计算公式为：其中ID_xyz,j表示采样序号分别为x,y,z的栅格单元代表点P_xyz相对第j个接收机的距离单元编号，转到步骤5.3；

步骤5.3

将步骤4中降分辨率回波信息矩阵P_M的第ID_xyz行第j列数据P_M(ID_xyz,j,j)作为各方向采样点序号分别为x,y,z的代表点P_xyz对接收机j的回波信号，对回波信号进行积累，其表达式为：Ima_Nx×Ny×Nz(k1,k2,k3)＝Ima_Nx×Ny×Nz(k1,k2,k3)+P_M(ID_xyz,j,j)，其中k1＝x,k2＝y,k3＝z；令接收机序号j加1，转到步骤5.4；

步骤5.4

如果j≤Nrec，则转至步骤5.2；

如果j＞Nrec，则j置为1，x加1，转至步骤5.5；

步骤5.5

如果x≤Nx，则转至步骤5.2；

如果x＞Nx，则x置为1，y加1，转至步骤5.6；

步骤5.6

如果y≤Ny，则转至步骤5.2；

如果y＞Ny，则y置为1，z加1，转至步骤5.7；

步骤5.7

如果z≤Nz，则转至步骤5.2；

如果z＞Nz，遍历结束，得到一个关于成像空间的三维图像矩阵Ima_Nx×Ny×Nz，转到步骤6；

步骤6、提取目标位置

首先，根据步骤5获得的三维图像矩阵Ima_Nx×Ny×Nz，在三维图像矩阵Ima_Nx×Ny×Nz中，遍历三维图像矩阵Ima_Nx×Ny×Nz的数据，找出三维图像矩阵Ima_Nx×Ny×Nz的最大值，然后从三维图像矩阵Ima_Nx×Ny×Nz中剔除该最大值所在栅格单元对应的回波信号；在剔除最大值所在栅格单元对应的回波信号的三维图像矩阵Ima_Nx×Ny×Nz中，遍历三维图像矩阵Ima_Nx×Ny×Nz的数据，依次找出Ntar-1个最大值，……，同理找出Ntar个最大值，记录Ntar个最大值的序号为n(n＝1,2…,Ntar)；具体步骤如下：

步骤6.1

初始化序号n＝1，在三维图像矩阵Ima_Nx×Ny×Nz中，遍历三维图像矩阵Ima_Nx×Ny×Nz的数据，找出由步骤5得到的三维图像矩阵Ima_Nx×Ny×Nz元素的最大值，记做Max_n，转到步骤6.2；

步骤6.2

寻找第n个最大值Max_n对应三维图像矩阵中的位置，有Ima_Nx×Ny×Nz(k1_{Max_n},k2_{Max_n},k3_{Max_n})＝Max_n，其中，k1_{Max_n},k2_{Max_n},k3_{Max_n}分别表示三维图像矩阵Ima_Nx×Ny×Nz元素的第n个最大值Max_n对应的行号，列号和层号；三维图像矩阵元素第n个最大值对应的目标位置坐标记为P_n＝P0+[k1_{Max_n}-Nx/2,k2_{Max_n}-Ny/2,k3_{Max_n}-Nz/2]*M*δ_r，其中P_n表示获得的第n个目标的位置，令接收机序号初始化为j＝1，转到步骤6.3；

步骤6.3

首先，由第n个目标的位置P_n，计算从发射机位置P_t经P_n到达第j个接收机的双站距离R_n,j，R_n,j＝‖P_t-P_n‖₂+‖P_n-Prec(j)‖₂；

然后，由双站距离R_n,j计算第n个目标相对第j个接收机的距离单元编号ID_n,j， ${\mathrm{ID}}_{n,j}=round\left(\frac{R_{n,j}-Rrec0}{M*{\mathrm{\δ}}_r}\right);$ 转到步骤6.4；

步骤6.4

根据步骤4建立的降分辨率回波信息矩阵P_M和步骤6.3中计算的距离单元编号ID_n,j，在降分辨率回波信息矩阵P_M的第ID_n,j行第j列通过公式 找出目标n对第j个接收机的回波信号，记为P_M(l,c)(l＝ID_n,j,c＝j)；

然后，从三维图像矩阵Ima_Nx×Ny×Nz中剔除此回波信号，Ima_Nx×Ny×Nz(k1,k2,k3)＝Ima_Nx×Ny×Nz(k1,k2,k3)-P_M(ID_n,j,j)，其中k1＝k1_{Max_n},k2＝k2_{Max_n},k3＝k3_{Max_n}；令j加1，转到步骤6.5；

步骤6.5

如果j≤Nrec，则转至步骤6.3；

如果j＞Nrec，则得到剔除第n个目标信息后的三维图像矩阵Ima_Nx×Ny×Nz；令n加1，转到步骤6.6；

步骤6.6

如果n≤Ntar，则转至步骤6.1；

如果n＞Ntar，则转至步骤7；

步骤7目标附近区域投影定位

提取Ntar个目标在观测空间的位置P_n(n＝1,2,…,Ntar)后，利用这些位置信息得到目标可能存在的区域，然后在目标可能存在的区域进行更高精度的投影定位。具体步骤如下：

步骤7.1

如果M＞1，则转到步骤7.2；

如果M＝1，则转到步骤8；

步骤7.2

对观测空间重新划分栅格：令沿观测空间x轴的采样点数Nx＝2*Nx，沿观测空间y轴的采样点数Ny＝2*Ny，沿观测空间z轴的采样点数Nz＝2*Nz，M＝M/2，然后执行步骤3完成观测空间重新划分栅格。

执行步骤4得到降分辨率回波信息矩阵P_M和每一栅格单元代表点的坐标P_xyz，P_xyz＝P0+[x-Nx/2,y-Ny/2,z-Nz/2]*M*δ_r，其中，x＝1,2…Nx，y＝1,2…Ny，z＝1,2…Nz；根据目标的位置P_n(n＝1,2,…,Ntar)对应的三维图像矩阵Ima_Nx×Ny×Nz中位置Pos_n(k1_{Max_n},k2_{Max_n},k3_{Max_n})，更新目标所在观测空间的采样点序号为 $X\_n=(k{1}_{{\mathrm{Max}}_n}-$ $1)*2,Y\_n=(k{2}_{{\mathrm{Max}}_n}-1)*2,Z\_n=(k{3}_{{\mathrm{Max}}_n}-1)*2;$ 以X_n为中心，构建向量P_x＝[X_n-h,X_n-h+1,…,X_n+h]，以Y_n为中心，构建向量P_y＝[Y_n-h,Y_n-h+1,…,Y_n+h]，以Z_n为中心，构建向量P_z＝[Z_n-h,Z_n-h+1,…,Z_n+h]，其中h是以第n个目标在观测空间中的位置为中心点；P_x、P_y、P_z为选取的三维区域各维对应的采样点序号矩阵。

构造矩阵P_{n_ty}，各行的存取元素为以第n个目标为中心所选取的三维区域中各栅格代表点所对应的三维采样点序号，P_{n_ty}的行向量表示为(P_x(w),P_y(p),P_z(q))，其中w＝1,2,…,2*h+1，p＝1,2,…,2*h+1，q＝1,2,…,2*h+1，P_x(w)表示矩阵P_x的第w个元素，P_y(p)表示矩阵P_y的第p个元素，P_z(q)表示矩阵P_z的第q个元素，矩阵P_{n_ty}表示第n个目标附近选取区域的代表点的集合。

构建投影矩阵P_tyall＝[P_{1_ty}；P_{2_ty}；…；P_{Ntar_ty}]，令P_ty＝unique(P_tyall)，P_ty的行数为H，unique(P_tyall)表示将矩阵P_tyall中代表相同栅格代表点的行元素消去，使各行所代表的栅格代表点各不相同。

步骤7.3构建目标可能存在的区域

根据步骤7.2得到的矩阵P_ty，计算选取的三维区域的代表点坐标P_aero(v,:)＝P0+[P_ty(v,1)-Nx/2,P_ty(v,2)-Ny/2,P_ty(v,3)-Nz/2]*M*δ_r，其中，v为代表点序号，v＝1,2,…,H，H为步骤7.2中定义的P_ty的行数。

步骤7.4投影处理

根据步骤7.3建立的目标可能存在区域的代表点，寻找每一代表点对应的各接收机的回波信号，采用投影方法对每一代表点对应的各个接收机回波信号进行积累；具体步骤如下：

步骤7.4.1

转至步骤5.1，创建三维图像矩阵Ima_Nx×Ny×Nz；然后令接收机序号初始化为j＝1，代表点序号v初始化为v＝1；由步骤5.1转至步骤7.4.2；

步骤7.4.2

根据步骤7.3中选取的区域代表点坐标P_aero(v,:)和第j个接收机位置Prec(j)，计算从发射机P_t经栅格单元代表点P_aero(v,:)到达第j个接收机的双站距离，记为R_v,j，计算公式：R_v,j＝‖P_t-P_aero(v,:)‖₂+‖P_aero(v,:)-Prec(j)‖₂，其中，P_aero(v,:)表示第v行所有列的数据，v＝1,2,…,H；

然后，由双站距离R_v,j计算代表点P_aero(v,:)对第j个接收机的距离单元编号，记为ID_v,j， ${\mathrm{ID}}_{v,j}=round\left(\frac{R_{v,j}-Rrec0}{M*{\mathrm{\δ}}_r}\right),$ 转到步骤7.4.3；

步骤7.4.3

将步骤7.2中所建立的降分辨率回波信息矩阵P_M的第j列第ID_v,j行中的数据作为栅格单元代表点P_aero(v,:)所对应的回波信号，记为P₀(v,j)，P₀(v,j)表示代表点P_aero(v,:)对第j个接收机的回波信号，并将回波信号P₀(v,j)在成像空间中积累，回波信号P₀(v,j)在成像空间积累的表达式为：Ima_Nx×Ny×Nz(k1,k2,k3)＝Ima_Nx×Ny×Nz(k1,k2,k3)+P₀(v,j)，其中k1＝x,k2＝y,k3＝z；令接收机序号j加1，并转到步骤7.4.4；

步骤7.4.4

如果j≤Nrec，则转至步骤7.4.2；

如果j＞Nrec，则j置为1，v加1，转至步骤7.4.5；

步骤7.4.5

如果v≤H，则转至步骤7.4.2；

如果v＞H，则投影结束，得到一个关于成像空间的三维图像矩阵Ima_Nx×Ny×Nz，转至步骤7.5；

步骤7.5提取目标位置

采用步骤6所提方法，获取三维图像矩阵Ima_Nx×Ny×Nz中的Ntar个最大值，并记录Ntar个最大值的序号为n(n＝1,2…,Ntar)；具体步骤如下：

步骤7.5.1

找出三维图像矩阵Ima_Nx×Ny×Nz中的最大值，记做Max_n，转到步骤7.5.2；

步骤7.5.2

转至步骤6.2，获得第n个目标的位置坐标P_n，令接收机序号初始化为j＝1；转至步骤7.5.3；

步骤7.5.3

转至步骤6.3，计算第n个目标相对第j个接收机的回波的距离单元编号ID_n,j；再由步骤6.3转到步骤7.5.4；

步骤7.5.4

根据步骤4建立的降分辨率回波信息矩阵P_M和第n个目标相对第j个接收机的回波的距离单元编号ID_n,j，在降分辨率回波信息矩阵P_M的第j列第ID_n,j行中找出第n个目标对第j个接收机的回波信号，记为P_M(l,c)(l＝ID_n,j,c＝j)；然后，从三维图像矩阵Ima_Nx×Ny×Nz中剔除此回波信号，Ima_Nx×Ny×Nz(k1,k2,k3)＝Ima_Nx×Ny×Nz(k1,k2,k3)-P_M(l,c)，其中k1＝k1_{Max_n},k2＝k2_{Max_n},k3＝k3_{Max_n}；令j加1，转到步骤7.5.5；

步骤7.5.5

如果j≤Nrec，则转到至骤7.5.3；

如果j＞Nrec，令n加1，转到步骤7.5.6；

步骤7.5.6

如果n≤Ntar，则转至步骤7.5.1；

如果n＞Ntar，则转至步骤7.1；

步骤8泰勒级数展开定位

根据步骤7求得的目标位置P_n(n＝1,2,…,Ntar)，找到目标位置P_n对应在多个接收机中的距离信息DI_n(k)(k＝1,2,…,NOR_n)，构建时差定位方程组，以P_n为初始值，利用泰勒级数展开定位算法求解目标位置，NOR_n表示满足一定条件的接收机个数，初始化接收机的序号j为1，具体步骤如下：

步骤8.1提取目标距离信息

选择第n(n＝1,2,…,Ntar)个目标位置P_n和第j个接收机，由双站距离公式计算信号从发射机P_t经P_n到达第j个接收机的双站距离，记为R_n,j，R_n,j＝‖P_t-P_n‖₂+‖P_n-Prec(j)‖₂。

步骤8.1.1

由R_n,j计算目标n的回波距离单元编号分别求出距离RC_n,j最近的、第2近的和第3近的整数RI_n,j(i)(i＝1,2,3)，初始化i的值为1，n的值为1，k的值为1，转到步骤8.1.2。

步骤8.1.2

取出步骤2中回波概率信息矩阵P的第RI_n,j(i)行第j列的数据，判断该数据是否大于门限ε。如果数据值大于门限ε，则令距离信息DI_n(k)＝RI_n,j(i)，EI_n(k)＝j，k值加1，j值加1，并转到步骤8.1.4；如果数值小于等于门限ε，则转到步骤8.1.3。

步骤8.1.3

如果i＜3，i的值加1，转到步骤8.1.2；

如果i＝3，j的值加1，转到步骤8.1.4。

步骤8.1.4

如果j≤Nrec，则转到步骤8.1；

如果j＞Nrec，n的值加1并转到步骤8.1.5；

步骤8.1.5

如果n≤Ntar，则令j＝1并转到步骤8.1；

如果n＞Ntar，则转到步骤8.2；

步骤8.2构建方程组

根据得到的目标在多个接收机中对应的距离信息，构建定位方程组，初始化目标序号n＝1，具体步骤如下：

步骤8.2.1

根据步骤8.1得到的矩阵DI_n，构造距离差矩阵RD_n，RD_n(m)＝[DI_n(m+1)-DI_n(1)]*M*δ_r+Rrec0，其中n＝1,2,…,Ntar，m＝1,2,…,NOR_n-1。取出步骤7求得的目标位置P_n，初始化p为1，初始化系数矩阵A_n为NOR_n-1行3列的零矩阵，常数矩阵B_n为NOR_n-1行1列的零矩阵，并转到步骤8.2.2。

步骤8.2.2

根据步骤8.1得到的矩阵EI_n，令接收机序号j＝EI_n(p)，计算P_n到接收机j的距离d_n(p)＝‖P_n-Prec(j)‖₂，转到步骤8.2.3。

步骤8.2.3

如果p＝1，p的值加1并转到步骤8.2.2；

如果p≠1，转到步骤8.2.4。

步骤8.2.4

令q＝p-1，系数矩阵A_n的第q行元素为常数矩阵B_n的第q行元素为b_nq＝RD_n(q)-(d_n(p)-d_n(1))。p的值加1并转到步骤8.2.5。

步骤8.2.5

如果p≤NOR_n，转到步骤8.2.2；

如果p＞NOR_n，转到步骤8.3。

步骤8.3求解方程组

根据步骤8.2得到的系数矩阵A_n和常数矩阵B_n，利用传统的最小二乘算法得到未知量δ，其中，n＝1,2,…,Ntar，δ是位置修正量，求得δ后可对目标位置P_n进行修正，判断δ与门限ε的大小，多次迭代得到最终的目标位置，具体步骤如下：

步骤8.3.1

利用传统的最小二乘算法求得位置修正量δ_n＝(A_n ^TA_n)^-1A_n ^TB_n，将目标位置P_n的值修改为P_n+δ_n，令sum(δ_n)表示对δ_n所有元素的求和运算，转到步骤8.3.2。

步骤8.3.2

如果Δ＞ε，转到步骤8.2.1；

如果Δ≤ε，n的值加1并转到步骤8.3.3。

步骤8.3.3

如果n≤Ntar，转到步骤8.2.1；

如果n＞Ntar，转到步骤8.4。

步骤8.4输出定位结果

目标定位结束，得到最终的Ntar个目标的三维位置坐标P_n(n＝1,2,…,Ntar)。

本发明的创新点在于基于传感器网络的多目标定位中，基于双站距离空间的投影定位算法需要遍历整个场景区域，运算量巨大，限制了该算法的使用，并且在一定程度上定位精度也有待提高，本发明在双站距离空间投影定位方法的基础之上，提供了一种基于分层投影和泰勒展开的多传感器多目标定位方法。该方法首先根据传感器接收到的回波信号计算目标存在的概率信息，然后将不同传感器的概率信息进行积累，用分层投影的方法，在得到粗略的目标位置基础上对目标附近区域进行高精度成像，大大减少了投影定位算法的运算量，最后采用泰勒展开算法，对估计的目标位置进行修正，进一步提高了多目标定位的精度。

本发明的实质是一种将双站距离空间投影定位算法采用分层投影，逐级缩小目标可能存在的区域，以提高定位精度，然后与时差定位相结合，采用泰勒展开对定位结果进行修正的多目标定位方法。双站距离空间投影定位算法首先将目标所在的观测区域进行栅格单元划分，并且为每个单元分配代表点。然后，用代表点到传感器的距离提取回波并将回波累积到对应的像素单元，此方法需要对整个观测区域进行遍历，从而带来了巨大的运算量；然而，由于目标在空中的分布是极其稀疏的，对整个区域进行投影是不必要的，利用分层投影逐级缩小搜寻目标的区域，可以有效减少投影定位的运算量；最后，以投影定位结果作为目标位置初始值构建时差方程组，利用泰勒展开对定位结果进行修正，提高了多目标定位的精度。基于分层投影和泰勒展开的多目标定位的系统示意图如图1所示。

本发明的优点：本发明从基于成像策略的双站距离空间投影定位方法的实时性和定位的精度两方面对其进行了改进。本发明首先，在投影成像过程中采用分层投影，逐级逼近目标所在区域，对目标进行定位；然后，在投影定位的基础之上，对定位结果进行修正。本发明不仅极大的减少了投影定位的运算量，而且有效地提高了多目标定位的精度。

附图说明

图1为基于分层投影和泰勒展开的多目标定位的系统示意图

图2为基于分层投影和泰勒展开的多目标定位中的投影算法示意图

其中，A和B分别表示不同目标的回波信息，F表示虚假目标或噪声的回波信息。N表示接收机的个数，t表示时间。对同一目标在不同接收机中的回波信息在三维图像空间中进行累加，得到三维图像矩阵，从中提取目标所在的位置。

图3为本发明方法的流程图

具体实施方式

本发明采用计算机仿真的方法进行验证，所有步骤、结论都在MATLAB-R2013a上验证正确。具体实施步骤如下：

步骤1、相关参数的初始化

初始化多目标定位方法的参数如下：电磁波传播速度v＝3×10⁸m/s；雷达空间距离分辨率δ_r＝10m；雷达在距离向上的最大距离单元编号L＝10000；以发射机所在位置为原点，建立三维笛卡尔坐标系；发射机位置P_t＝[0,0,0]；目标数目Ntar＝5，目标随机分布在以[5,5,10]×1000m为中心，边长为200m的立方体区域中，第i(i＝1,2,…,5)个目标的位置为Ptar(i)＝[5,5,10]×1000+[(rand(1)-0.5),(rand(1)-0.5),(rand(1)-0.5)]×200m，i为目标序号；目标散射截面积为σ＝(30×rand(1)+1)m²；接收机数目为Nrec＝40，接收机随机分布在以[0,0,0]为中心，边长30Km的正方形区域中，第j(j＝1,2,…,40)个接收机位置为Prec(j)＝[rand(1)-0.5,rand(1)-0.5,rand(1)-0.5]×30×1000m，j为接收机序号；接收机的接收距离波门为Rrec0＝-25Km；各接收机通道的信噪比为SNR＝3dB；以单载频脉冲信号为雷达发射信号，脉冲载频为f_c＝1GHz，脉冲宽度为T_p＝0.0667us；信号经发射机发射被第i个目标散射，由第j个接收机接收，再经去载频、离散采样后的回波信号记为Echo(i,j)， $Echo(i,j)=(30\×rand(1)+1)\×\exp (-p\frac{2\mathrm{\π}\×1\×{10}^9}{3\×{10}^8}{R}_{i,j}),$ 其中，第i个目标与第j个接收机的双站距离R_i,j＝‖P_t-Ptar(i)‖₂+‖Ptar(i)-Prec(j)‖₂，i＝1,2,…,5，j＝1,2,…,40，exp(·)是以自然指数为底的指数函数，p表示虚数单位即p²＝-1，‖·‖₂表示向量的L2范数，π为圆周率，上述中rand(1)表示在0到1之间等概率的随机取一个数。

步骤2、创建回波存储矩阵并求取回波概率信息矩阵

由发射机和接收机的时间同步系统(GPS等)，收集空中各目标对不同接收机的回波信号Echo(i,j)以及相应的双站延时τ_i,j，其中i＝1,2,…,5，j＝1,2,…,40，τ_i,j表示第i个目标相对第j个接收机的双站延时，Echo(i,j)表示第j个接收机接收到的第i个目标散射的回波信号。

步骤2.1初始化回波存储矩阵

根据步骤1中雷达在距离向上的最大距离单元编号L＝10000以及地面接收机数目Nrec＝40，可以定义一个10000行40列的回波存储矩阵，记为并将其初始化为零矩阵。回波存储矩阵的每一元素，记为l＝1,2,…,10000，c＝1,2,…,40，其中l为回波存储矩阵的行数，c为回波存储矩阵的列数；根据第j个接收机获得的第i个目标的双站延时τ_i,j及对应的目标回波信号Echo(i,j)，计算第i个目标对第j个接收机的双站距离为R_i,j＝3×10⁸×τ_i,j，其中i＝1,2,…5；j＝1,2,…40；根据双站距离R_i,j，计算相应回波的距离单元编号其中ID_i,j表示第i个目标的回波在第j个接收机中的距离单元编号，1≤ID_i,j≤10000，且ID_i,j取值为正整数；然后将经发射机发射被第i个目标散射，由第j个接收机接收的回波信号Echo(i,j)，存储在回波存储矩阵中的位置处，其中l＝ID_i,j，c＝j，即(l＝ID_i,j,c＝j)；遍历目标和接收机，得到回波存储矩阵

步骤2.2回波存储矩阵转换为回波概率信息矩阵

根据得到的存储矩阵求回波概率信息矩阵P，回波概率信息矩阵P的维数与相同，其计算公式为：P(l,c)＝1-exp(-o(l,c)²/2σ(c))，其中l＝1,2,…,10000，c＝1,2,…,40，o(l,c)是的幅度，σ(c)是矩阵第c列的幅度标准差。

步骤3、栅格划分

步骤3.1建立观测空间

在步骤1建立的笛卡尔坐标系中，设定一个三维长方体观测空间Ω_Sear，使其可以包含所有目标，观测空间的长、宽、高分别平行于笛卡尔坐标系的x轴、y轴、z轴，沿x轴，y轴，z轴的采样点数均为Nx＝Ny＝Nz＝32，相对于雷达空间距离分辨率δ_r的倍数取值为M＝8，观测空间长为8×10×32，宽为8×10×32，高为8×10×32；

步骤3.2等间隔划分栅格

步骤3.1中观测空间Ω_Sear的中心坐标为P0＝[5,5,10]×1000，根据步骤1中的雷达空间分辨率δ_r＝10m以及相对于空间分辨率δ_r的倍数M＝8，将观测空间划分为一组栅格，每一栅格单元大小为(8×10)×(8×10)×(8×10)，然后为每一栅格单元标记一代表点，各代表点坐标表示为P_xyz＝[5,5,10]×1000+[x-32/2,y-32/2,z-32/2]×8×10，x,y,z分别是观测空间沿x轴，y轴，z轴的采样点序号，x＝1,2…32，y＝1,2…32，z＝1,2…32，P_xyz表示采样点序号为x,y,z的代表点的三维坐标，相邻代表点间的距隔为8×10m。

步骤4、回波降分辨率处理

根据步骤2得到的回波概率信息矩阵P，通过整合相邻单元的矩阵元素值完成回波降分辨率处理，得到降分辨率回波信息矩阵P_M，降分辨率回波信息矩阵元素的计算公式为： $P_M(k,c)=1-{\Π}_{i=1}^8(1-P(\left(k-1\right)\×8+i,c\left)\right),$ k＝1,2,…,1250，c＝1,2,…,40。

步骤5、投影处理

在步骤3建立的观测空间Ω_Sear中，寻找每一栅格代表点对应各接收机的回波信号，采用投影方法对每一栅格代表点对应的各接收机的回波信号进行积累；具体步骤如下：

步骤5.1

首先创建一个三维图像矩阵Ima_32×32×32，矩阵的大小为32×32×32，并初始化为零矩阵，用于存储观测空间的成像信息；三维图像矩阵元素单元记为Ima_32×32×32(k1,k2,k3)，k1、k2和k3分别表示三维图像矩阵Ima_32×32×32的行号、列号、层号，其中k1＝1,2,…32，k2＝1,2,…32，k3＝1,2,…32；然后分别初始化观测空间沿x轴，y轴，z轴的采样点序号x＝1,y＝1,z＝1，接收机序号初始化为j＝1；转到步骤5.2；

步骤5.2

根据步骤3中栅格代表点的三维坐标P_xyz和接收机j的位置Prec(j)，计算从发射机位置P_t经栅格代表点P_xyz到达第j个接收机的双站距离R_xyz,j，其计算公式为：R_xyz,j＝‖P_t-P_xyz‖₂+‖P_xyz-Prec(j)‖₂，由双站距离R_xyz,j计算栅格单元代表点P_xyz相对接收机j的距离单元编号其中ID_xyz,j表示采样序号分别为x,y,z的栅格单元代表点P_xyz相对第j个接收机的距离单元编号，转到步骤5.3；

步骤5.3

将步骤4中降分辨率回波信息矩阵P_M的第ID_xyz行第j列数据P_M(ID_xyz,j,j)作为各方向采样点序号分别为x,y,z的代表点P_xyz对接收机j的回波信号，对回波信号进行积累，其表达式为：Ima_32×32×32(k1,k2,k3)＝Ima_32×32×32(k1,k2,k3)+P_M(ID_xyz,j,j)，其中k1＝x,k2＝y,k3＝z；令接收机序号j加1，转到步骤5.4；

步骤5.4

如果j≤40，则转至步骤5.2；

如果j＞40，则j置为1，x加1，转至步骤5.5；

步骤5.5

如果x≤32，则转至步骤5.2；

如果x＞32，则x置为1，y加1，转至步骤5.6；

步骤5.6

如果y≤32，则转至步骤5.2；

如果y＞32，则y置为1，z加1，转至步骤5.7；

步骤5.7

如果z≤32，则转至步骤5.2；

如果z＞32，遍历结束，得到一个关于成像空间的三维图像矩阵Ima_32×32×32，转到步骤6；

步骤6、提取目标位置

首先，根据步骤5获得的三维图像矩阵Ima_32×32×32，找出三维图像矩阵Ima_32×32×32的最大值，然后从三维图像矩阵Ima_32×32×32中剔除该最大值所在栅格单元对应的回波信号；按此方法，依次找出5个最大值，并记录各最大值的序号为n(n＝1,2…,5)；具体步骤如下：

步骤6.1

初始化序号n＝1，找出由步骤5得到的三维图像矩阵Ima_32×32×32元素的最大值，记做Max_1，转到步骤6.2；

步骤6.2

寻找第n个最大值Max_n对应三维图像矩阵中的位置，有Ima_32×32×32(k1_{Max_n},k2_{Max_n},k3_{Max_n})＝Max_n，其中，k1_{Max_n},k2_{Max_n},k3_{Max_n}分别表示三维图像矩阵Ima_32×32×32元素的第n个最大值Max_n对应的行号，列号和层号；三维图像矩阵元素第n个最大值对应的目标位置坐标记为P_n＝[5,5,10]×1000+[k1_{Max_n}-32/2,k2_{Max_n}-32/2,k3_{Max_n}-32/2]*8*10，其中P_n表示获得的第n个目标的位置，令接收机序号初始化为j＝1，转到步骤6.3；

步骤6.3

首先，由第n个目标的位置P_n，计算从发射机位置P_t＝[0,0,0]经P_n到达第j个接收机的双站距离R_n,j，R_n,j＝‖P_t-P_n‖₂+‖P_n-Prec(j)‖₂，然后，由双站距离R_n,j计算第n个目标相对第j个接收机的距离单元编号ID_n,j，转到步骤6.4；

步骤6.4

如果0＜ID_n,j≤10000，则根据步骤4建立的降分辨率回波信息矩阵P_M和步骤6.3中计算的距离单元编号ID_n,j，在降分辨率回波信息矩阵P_M的第ID_n,j行第j列找出目标n对第j个接收机的回波信号，记为P_M(l,c)(l＝ID_n,j,c＝j)；然后，从三维图像矩阵Ima_32×32×32中剔除此回波信号，Ima_32×32×32(k1,k2,k3)＝Ima_32×32×32(k1,k2,k3)-P_M(ID_n,j,j)，其中k1＝k1_{Max_n},k2＝k2_{Max_n},k3＝k3_{Max_n}；令j加1，转到步骤6.5；

如果ID_n,j≤0或者ID_n,j＞10000，则直接令j加1，转到步骤6.5；

步骤6.5

如果j≤40，则转至步骤6.3；

如果j＞40，则得到剔除第n个目标信息后的三维图像矩阵Ima_32×32×32；令n加1，转到步骤6.6；

步骤6.6

如果n≤5，则转至步骤6.1；

如果n＞5，则转至步骤7；

步骤7目标附近区域投影定位

提取5个目标在观测空间的位置P_n(n＝1,2,…,5)后，就可以利用这些位置信息得到目标可能存在的区域，然后在目标可能存在的区域进行更高精度的投影定位。具体步骤如下：

步骤7.1

如果M＞1，则转到步骤7.2；

如果M＝1，则转到步骤8；

步骤7.2

首先令M＝M/2，Nx＝2*Nx，Ny＝2*Ny，Nz＝2*Nz，然后执行步骤3和步骤4，对观测空间重新划分栅格，并得到降分辨率回波信息矩阵P_M和每一栅格单元代表点的坐标P_xyz＝[5,5,10]×1000+[x-Nx/2,y-Ny/2,z-Nz/2]*M*10，其中，x＝1,2…Nx，y＝1,2…Ny，z＝1,2…Nz；根据目标的位置P_n(n＝1,2,…,5)对应的三维图像矩阵Ima_32×32×32中位置Pos_n(k1_{Max_n},k2_{Max_n},k3_{Max_n})，更新目标所在观测空间的采样点序号为 $X\_n=(k{1}_{{\mathrm{Max}}_n}-1)*2,Y\_n=(k{2}_{{\mathrm{Max}}_n}-1)*2,Z\_n=(k{3}_{{\mathrm{Max}}_n}-1)*2;$ 以X_n为中心，构建向量P_x＝[X_n-5,X_n-5+1,…,X_n+5]，以Y_n为中心，构建向量P_y＝[Y_n-5,Y_n-5+1,…,Y_n+5]，以Z_n为中心，构建向量P_z＝[Z_n-5,Z_n-5+1,…,Z_n+5]，其中5是以第n个目标在观测空间中的位置为中心点，选取的三维区域的每一维的采样点个数的一半；P_x、P_y、P_z为选取的三维区域各维对应的采样点序号。构造矩阵P_{n_ty}，存取元素为所选取的三维区域中各栅格代表点对应的采样点序号，P_{n_ty}各行的元素可表示为(P_x(w),P_y(p),P_z(q))，其中w、p、q可任意组合，分别取值为w＝1,2,…,2*5+1，p＝1,2,…,2*5+1，q＝1,2,…,2*5+1，矩阵P_{n_ty}表示第n个目标附近选取区域的代表点的集合。构建投影矩阵P_tyall＝[P_{1_ty}；P_{2_ty}；…；P_{5_ty}]，令P_ty＝unique(P_tyall)，P_ty的行数为H，unique(P_tyall)表示将矩阵P_tyall中具有相同元素的行消去，只保留一个，使各行元素各不相同。

步骤7.3构建目标可能存在的区域

根据步骤7.2得到的矩阵P_ty，计算选取的三维区域的代表点坐标P_aero(v,:)＝[5,5,10]×1000+[P_ty(v,1)-Nx/2,P_ty(v,2)-Ny/2,P_ty(v,3)-Nz/2]*M*10，其中，v为代表点序号，v＝1,2,…,H，H为步骤7.2中定义的P_ty的行数。

步骤7.4投影处理

根据步骤7.3建立的目标可能存在区域的代表点，寻找每一代表点对应的各接收机的回波信号，采用投影方法对每一代表点对应的各个接收机回波信号进行积累；具体步骤如下：

步骤7.4.1

转至步骤5.1，创建三维图像矩阵Ima_Nx×Ny×Nz；然后令接收机序号初始化为j＝1，代表点序号v初始化为v＝1；由步骤5.1转至步骤7.4.2；

步骤7.4.2

根据步骤7.3中选取的区域代表点坐标P_aero(v,:)和第j个接收机位置Prec(j)，计算从发射机P_t经栅格单元代表点P_aero(v,:)到达第j个接收机的双站距离R_v,j＝‖P_t-P_aero(v,:)‖₂+‖P_aero(v,:)-Prec(j)‖₂，其中，P_aero(v,:)表示第v行所有列的数据，v＝1,2,…,H；然后，由双站距离R_v,j计算代表点P_aero(v,:)对第j个接收机的距离单元编号 ${\mathrm{ID}}_{v,j}=round\left(\frac{R_{v,j}-Rrec0}{M*{\mathrm{\δ}}_r}\right),$ 转到步骤7.4.3；

步骤7.4.3

将步骤7.2中所建立的降分辨率回波信息矩阵P_M的第j列第ID_v,j行中的数据作为栅格单元代表点P_aero(v,:)所对应的回波信号，记为P₀(v,j)，P₀(v,j)表示代表点P_aero(v,:)对第j个接收机的回波信号，并将回波信号P₀(v,j)在成像空间中积累，其表达式为：Ima_Nx×Ny×Nz(k1,k2,k3)＝Ima_Nx×Ny×Nz(k1,k2,k3)+P₀(v,j)，其中k1＝x,k2＝y,k3＝z；令接收机序号j加1，并转到步骤7.4.4；

步骤7.4.4

如果j≤40，则转至步骤7.4.2；

如果j＞40，则j置为1，v加1，转至步骤7.4.5；

步骤7.4.5

如果v≤H，则转至步骤7.4.2；

如果v＞H，则投影结束，得到一个关于成像空间的三维图像矩阵Ima_Nx×Ny×Nz，转至步骤7.5；

步骤7.5提取目标位置

根据步骤6所提方法，获取三维图像矩阵Ima_Nx×Ny×Nz中的5个最大值，并记录5个最大值的序号为n(n＝1,2…,5)；具体步骤如下：

步骤7.5.1

找出三维图像矩阵Ima_Nx×Ny×Nz中的最大值，记做Max_n，转到步骤7.5.2；

步骤7.5.2

转至步骤6.2，获得第n个目标的位置坐标P_n，令接收机序号初始化为j＝1；转至步骤7.5.3；

步骤7.5.3

转至步骤6.3，计算第n个目标相对第j个接收机的回波的距离单元编号ID_n,j；再由步骤6.3转到步骤7.5.4；

步骤7.5.4

根据降分辨率回波信息矩阵P_M和第n个目标相对第j个接收机的回波的距离单元编号ID_n,j，在降分辨率回波信息矩阵P_M的第j列第ID_n,j行中找出第n个目标对第j个接收机的回波信号，记为P_M(l,c)(l＝ID_n,j,c＝j)；然后，从三维图像矩阵Ima_Nx×Ny×Nz中剔除此回波信号，Ima_Nx×Ny×Nz(k1,k2,k3)＝Ima_Nx×Ny×Nz(k1,k2,k3)-P_M(ID_n,j,j)，其中k1＝k1_{Max_n},k2＝k2_{Max_n},k3＝k3_{Max_n}；令j加1，转到步骤7.5.5；

步骤7.5.5

如果j≤40，则转到至骤7.5.3；

如果j＞40，令n加1，转到步骤7.5.6；

步骤7.5.6

如果n≤5，则转至步骤7.5.1；

如果n＞5，则转至步骤7.1；

步骤8泰勒级数展开定位

根据步骤7求得的目标位置P_n(n＝1,2,…,5)，找到P_n对应在多个接收机中的距离信息DI_n(k)(k＝1,2,…,NOR_n)，构建时差定位方程组，以P_n为初始值，就可以利用泰勒级数展开定位算法求解目标位置，NOR_n表示满足一定条件的接收机个数，初始化接收机的序号j为1，具体步骤如下：

步骤8.1提取目标距离信息

选择第n(n＝1,2,…,5)个目标位置P_n和第j个接收机，由双站距离公式计算信号从发射机P_t经P_n到达第j个接收机的双站距离R_n,j，R_n,j＝‖P_t-P_n‖₂+‖P_n-Prec(j)‖₂。

步骤8.1.1

由R_n,j计算目标n的回波距离单元编号分别求出距离RC_n,j最近的、第2近的和第3近的整数RI_n,j(i)(i＝1,2,3)，初始化i的值为1，n的值为1，k的值为1，转到步骤8.1.2。

步骤8.1.2

取出步骤2中回波概率信息矩阵P的第RI_n,j(i)行第j列的数据，判断该数据是否大于门限ε＝0.5。如果数值大于门限0.5，则令DI_n(k)＝RI_n,j(i)，EI_n(k)＝j，k值加1，j值加1，并转到步骤8.1.4；如果数值小于等于门限0.5，则转到步骤8.1.3。

步骤8.1.3

如果i＜3，i的值加1，转到步骤8.1.2；

如果i＝3，j的值加1，转到步骤8.1.4。

步骤8.1.4

如果j≤40，则转到步骤8.1；

如果j＞40，n的值加1并转到步骤8.1.5；

步骤8.1.5

如果n≤5，则令j＝1并转到步骤8.1；

如果n＞5，则转到步骤8.2；

步骤8.2构建方程组

根据得到的目标在多个接收机中对应的距离信息，就可以构建时差定位方程组，初始化目标序号n＝1，具体步骤如下：

步骤8.2.1

根据步骤8.1得到的矩阵DI_n，构造距离差矩阵RD_n，RD_n(m)＝[DI_n(m+1)-DI_n(1)]*M*10+(-25Km)，其中n＝1,2,…,5，m＝1,2,…,NOR_n-1。取出步骤7求得的目标位置P_n，初始化p为1，初始化系数矩阵A_n为NOR_n-1行3列的零矩阵，常数矩阵B_n为NOR_n-1行1列的零矩阵，并转到步骤8.2.2。

步骤8.2.2

根据步骤8.1得到的矩阵EI_n，令接收机序号j＝EI_n(p)，计算P_n到接收机j的距离d_n(p)＝‖P_n-Prec(j)‖₂，转到步骤8.2.3。

步骤8.2.3

如果p＝1，p的值加1并转到步骤8.2.2；

如果p≠1，转到步骤8.2.4。

步骤8.2.4

令q＝p-1，系数矩阵A_n的第q行元素为常数矩阵B_n的第q行元素为b_nq＝RD_n(q)-(d_n(p)-d_n(1))。p的值加1并转到步骤8.2.5。

步骤8.2.5

如果p≤NOR_n，转到步骤8.2.2；

如果p＞NOR_n，转到步骤8.3。

步骤8.3求解方程组

根据步骤8.2得到的系数矩阵A_n和常数矩阵B_n，可以利用最小二乘算法得到未知量δ，其中，n＝1,2,…,5，δ是位置修正量，求得δ后可对目标位置P_n进行修正，判断δ与门限ε＝0.5的大小，多次迭代得到最终的目标位置，具体步骤如下：

步骤8.3.1

利用最小二乘算法求得位置修正量δ_n＝(A_n ^TA_n)^-1A_n ^TB_n，将P_n的值修改为P_n+δ_n，令sum(δ_n)表示对δ_n所有元素的求和运算，转到步骤8.3.2。

步骤8.3.2

如果Δ＞0.5，转到步骤8.2.1；

如果Δ≤0.5，n的值加1并转到步骤8.3.3。

步骤8.3.3

如果n≤5，转到步骤8.2.1；

如果n＞5，转到步骤8.4。

步骤8.4输出定位结果

目标定位结束，得到最终的5个目标的三维位置坐标P_n(n＝1,2,…,5)。

Claims (1)

1.一种基于分层投影和泰勒展开的多传感器多目标定位方法，其特征是它包括以下步骤：

步骤1、相关参数的初始化

初始化的参数均为已知，如下所示：电磁波传播速度，记为v；雷达空间距离分辨率，记为δ_r；雷达在距离向上的最大距离单元编号，记为整数L；以发射机所在位置为原点，建立三维笛卡尔坐标系；发射机位置，记为P_t；目标数目，记为Ntar；第i个目标的位置，记为Ptar(i)(i＝1，2，…，Ntar)，i为目标序号；目标散射截面积，记为σ；接收机数目，记为Nrec；第j个地面接收机位置，记为Prec(j)(j＝1，2，…，Nrec)，j为接收机序号；接收机的接收距离波门，均记为Rrec0；各接收机通道的信噪比，记为SNR；以单载频脉冲信号为雷达发射信号，脉冲载频记为f_c，脉冲宽度记为T_p；信号经发射机发射被第i个目标散射，由第j个接收机接收，再经去载频、离散采样后的回波信号记为Echo(i，j)，其表达式Echo(i，j)＝σ_i，j*exp(-pkR_i，j)，R_i，j表示第i个目标与第j个接收机的双站距离，表达式R_i，j＝||P_t-Ptar(i)||₂+||Ptar(i)-Prec(j)||₂，其中exp(·)是以自然指数为底的指数函数，σ_i，j为第i个目标相对于第j个接收机的目标散射截面积，p表示虚数单位即p²＝-1，||·||₂表示向量的L2范数，*表示乘号，π为圆周率，v表示光速；

步骤2、创建回波存储矩阵并求取回波概率信息矩阵

由发射机和接收机的时间同步系统获取空中各目标对不同接收机的回波信号Echo(i，j)以及相应的双站延时τ_i，j，其中i＝1，2，…，Ntar，j＝1，2，…，Nrec，τ_i，j表示第i个目标相对第j个接收机的双站延时，Echo(i，j)表示第j个接收机接收到的第i个目标散射的回波信号；

步骤2.1初始化回波存储矩阵

根据步骤1中雷达在距离向上的最大距离单元编号L以及地面接收机数目Nrec，定义一个L行Nrec列的回波存储矩阵，记为并将回波存储矩阵初始化为零矩阵；回波存储矩阵的每一元素，记为l＝1，2，…，L，c＝1，2，…，Nrec，其中l为回波存储矩阵的行数，c为回波存储矩阵的列数；

根据第j个接收机获得的第i个目标的双站延时τ_i，j及对应的目标回波信号Echo(i，j)，计算第i个目标对第j个接收机的双站距离为R_i，j＝v*τ_i，j，其中v表示光速，*表示乘号，i＝1，2，…Ntar；j＝1，2，…Nrec；

根据双站距离R_i，j，计算相应回波的距离单元编号ID_i，j，公式为：其中ID_i，j表示第i个目标的回波在第j个接收机中的距离单元编号，1≤ID_i，j≤L，且ID_i，j取值为正整数，L为最大的距离单元编号，Rrec0表示接收机的接收距离波门；

然后将经发射机发射被第i个目标散射，由第j个接收机接收的回波信号Echo(i，j)，存储在回波存储矩阵中的位置处，即由其中l＝ID_i，j，c＝j，i＝1，2，…Ntar，j＝1，2，…Nrec，得到回波存储矩阵

步骤2.2回波存储矩阵转换为回波概率信息矩阵

根据得到的存储矩阵求回波概率信息矩阵P，回波概率信息矩阵P的维数与相同，回波概率信息矩阵P计算公式为：P(l，c)＝1-exp(-o(l，c)²/2σ(c))，其中l＝1，2，…，L，c＝1，2，…，Nrec，o(l，c)是的幅度，σ(c)是矩阵第c列的幅度标准差；

步骤3、栅格划分

步骤3.1建立观测空间

在步骤1建立的笛卡尔坐标系中，设定一个三维长方体观测空间Ω_Sear，观测空间的长、宽、高分别平行于笛卡尔坐标系的x轴、y轴、z轴；观测空间长为M×δ_r×Nx，宽为M×δ_r×Ny，高为M×δ_r×Nz，其中Nx，Ny，Nz分别为沿观测空间x轴，y轴，z轴的采样点数，δ_r表示雷达空间距离分辨率，M取值为2的非负整数次方，表示栅格划分边长相对于雷达空间距离分辨率δ_r的倍数；

步骤3.2等间隔划分栅格

将步骤3.1观测空间Ω_Sear的中心记为P0，步骤1中的雷达空间距离分辨率δ_r作为最小单元，将观测空间划分为一组栅格，每一栅格单元大小为Mδ_r×Mδ_r×Mδ_r，然后为每一栅格单元标记一代表点，代表点坐标为P_xyz＝P0+[x-Nx/2，y-Ny/2，z-Nz/2]*M*δ_r，x，y，z分别是观测空间沿x轴，y轴，z轴的采样点序号，x＝1，2…Nx，y＝1，2…Ny，z＝1，2…Nz，P_xyz表示采样点序号为x，y，z的代表点的三维坐标，相邻代表点间的距离为Mδ_r；

步骤4、回波降分辨率处理

根据步骤2得到的回波概率信息矩阵P，采用降分辨率回波信息矩阵元素的计算公式得到降分辨率回波信息矩阵P_M，其中取值为小于等于其自身的最大正整数，c＝1，2，…，Nrec，其中M表示栅格划分边长相对于空间分辨率δ_r的倍数，L表示最大的距离单元编号，M取值为2的非负整数次方，表示栅格划分边长相对于雷达空间距离分辨率δ_r的倍数；

步骤5、投影处理

在步骤3建立的观测空间Ω_Sear中，寻找每一栅格代表点对应各接收机的回波信号，采用传统的投影方法对每一栅格代表点对应的各接收机的回波信号进行积累；具体步骤如下：

步骤5.1

首先创建一个三维图像矩阵Ima_Nx×Ny×Nz，矩阵的大小为Nx×Ny×Nz，并初始化为零矩阵，用于存储观测空间的成像信息；三维图像矩阵元素单元记为Ima_Nx×Ny×Nz(k1，k2，k3)，k1、k2和k3分别表示三维图像矩阵Ima_Nx×Ny×Nz的行号、列号、层号，其中k1＝1，2，…Nx，k2＝1，2，…Ny，k3＝1，2，…Nz，Nx，Ny，Nz分别是三维图像矩阵的行数，列数，层数；然后分别初始化观测空间沿x轴，y轴，z轴的采样点序号x＝1，y＝1，z＝1，接收机序号初始化为j＝1；转到步骤5.2；

步骤5.2

根据步骤3中栅格代表点的三维坐标P_xyz＝P0+[x-Nx/2，y-Ny/2，z-Nz/2]*M*δ_r和接收机j的位置Prec(j)，计算从发射机位置P_t经栅格代表点P_xyz到达第j个接收机的双站距离，记为R_xyz，j，其计算公式为：R_xyz，j＝||P_t-P_xyz||₂+||P_xyz-Prec(j)||₂；

由双站距离R_xyz，j计算栅格单元代表点P_xyz相对接收机j的距离单元编号，记为ID_xyz，j，其计算公式为：其中ID_xyz，j表示采样序号分别为x，y，z的栅格单元代表点P_xyz相对第j个接收机的距离单元编号，转到步骤5.3；

步骤5.3

将步骤4中降分辨率回波信息矩阵P_M的第ID_xyz行第j列数据P_M(ID_xyz，j，j)作为各方向采样点序号分别为x，y，z的代表点P_xyz对接收机j的回波信号，对回波信号进行积累，其表达式为：Ima_Nx×Ny×Nz(k1，k2，k3)＝Ima_Nx×Ny×Nz(k1，k2，k3)+P_M(ID_xyz，j，j)，其中k1＝x，k2＝y，k3＝z；令接收机序号j加1，转到步骤5.4；

步骤5.4

如果j≤Nrec，则转至步骤5.2；

如果j＞Nrec，则j置为1，x加1，转至步骤5.5；

步骤5.5

如果x≤Nx，则转至步骤5.2；

如果x＞Nx，则x置为1，y加1，转至步骤5.6；

步骤5.6

如果y≤Ny，则转至步骤5.2；

如果y＞Ny，则y置为1，z加1，转至步骤5.7；

步骤5.7

如果z≤Nz，则转至步骤5.2；

如果z＞Nz，遍历结束，得到一个关于成像空间的三维图像矩阵Ima_Nx×Ny×Nz，转到步骤6；

步骤6、提取目标位置

首先，根据步骤5获得的三维图像矩阵Ima_Nx×Ny×Nz，在三维图像矩阵Ima_Nx×Ny×Nz中，遍历三维图像矩阵Ima_Nx×Ny×Nz的数据，找出三维图像矩阵Ima_Nx×Ny×Nz的最大值，然后从三维图像矩阵Ima_Nx×Ny×Nz中剔除该最大值所在栅格单元对应的回波信号；在剔除最大值所在栅格单元对应的回波信号的三维图像矩阵Ima_Nx×Ny×Nz中，遍历三维图像矩阵Ima_Nx×Ny×Nz的数据，依次找出Ntar-1个最大值，……，同理找出Ntar个最大值，记录Ntar个最大值的序号为n(n＝1，2…，Ntar)；具体步骤如下：

步骤6.1

初始化序号n＝1，在三维图像矩阵Ima_Nx×Ny×Nz中，遍历三维图像矩阵Ima_Nx×Ny×Nz的数据，找出由步骤5得到的三维图像矩阵Ima_Nx×Ny×Nz元素的最大值，记做Max_n，转到步骤6.2；

步骤6.2

寻找第n个最大值Max_n对应三维图像矩阵中的位置，有Ima_Nx×Ny×Nz(k1_{Max_n}，k2_{Max_n}，k3_{Max_n})＝Max_n，其中，k1_{Max_n}，k2_{Max_n}，k3_{Max_n}分别表示三维图像矩阵Ima_Nx×Ny×Nz元素的第n个最大值Max_n对应的行号，列号和层号；三维图像矩阵元素第n个最大值对应的目标位置坐标记为P_n＝P0+[k1_{Max_n}-Nx/2，k2_{Max_n}-Ny/2，k3_{Max_n}-Nz/2]*M*δ_r，其中P_n表示获得的第n个目标的位置，令接收机序号初始化为j＝1，转到步骤6.3；

步骤6.3

首先，由第n个目标的位置P_n，计算从发射机位置P_t经P_n到达第j个接收机的双站距离R_n，j，R_n，j＝||P_t-P_n||₂+||P_n-Prec(j)||₂；

然后，由双站距离R_n，j计算第n个目标相对第j个接收机的距离单元编号ID_n，j，转到步骤6.4；

步骤6.4

根据步骤4建立的降分辨率回波信息矩阵P_M和步骤6.3中计算的距离单元编号ID_n，j，在降分辨率回波信息矩阵P_M的第ID_n，j行第j列通过公式找出目标n对第j个接收机的回波信号，记为P_M(l，c)(l＝ID_n，j，c＝j)；

然后，从三维图像矩阵Ima_Nx×Ny×Nz中剔除此回波信号，Ima_Nx×Ny×Nz(k1，k2，k3)＝Ima_Nx×Ny×Nz(k1，k2，k3)-P_M(ID_n，j，j)，其中k1＝k1_{Max_n}，k2＝k2_{Max_n}，k3＝k3_{Max_n}；令j加1，转到步骤6.5；

步骤6.5

如果j≤Nrec，则转至步骤6.3；

如果j＞Nrec，则得到剔除第n个目标信息后的三维图像矩阵Ima_Nx×Ny×Nz；令n加1，转到步骤6.6；

步骤6.6

如果n≤Ntar，则转至步骤6.1；

如果n＞Ntar，则转至步骤7；

步骤7目标附近区域投影定位

提取Ntar个目标在观测空间的位置P_n(n＝1，2，…，Ntar)后，利用这些位置信息得到目标可能存在的区域，然后在目标可能存在的区域进行更高精度的投影定位；具体步骤如下：

步骤7.1

如果M＞1，则转到步骤7.2；

如果M＝1，则转到步骤8；

步骤7.2

对观测空间重新划分栅格：令沿观测空间x轴的采样点数Nx＝2*Nx，沿观测空间y轴的采样点数Ny＝2*Ny，沿观测空间z轴的采样点数Nz＝2*Nz，M＝M/2，然后执行步骤3完成观测空间重新划分栅格；

执行步骤4得到降分辨率回波信息矩阵P_M和每一栅格单元代表点的坐标P_xyz，P_xyz＝P0+[x-Nx/2，y-Ny/2，z-Nz/2]*M*δ_r，其中，x＝1，2…Nx，y＝1，2…Ny，z＝1，2…Nz；根据目标的位置P_n(n＝1，2，…，Ntar)对应的三维图像矩阵Ima_Nx×Ny×Nz中位置Pos_n(k1_{Max_n}，k2_{Max_n}，k3_{Max_n})，更新目标所在观测空间的采样点序号为以X_n为中心，构建向量P_x＝[X_n-h，X_n-h+1，…，X_n+h]，以Y_n为中心，构建向量P_y＝[Y_n-h，Y_n-h+1，…，Y_n+h]，以Z_n为中心，构建向量P_z＝[Z_n-h，Z_n-h+1，…，Z_n+h]，其中h是以第n个目标在观测空间中的位置为中心点；P_x、P_y、P_z为选取的三维区域各维对应的采样点序号矩阵；

构造矩阵P_{n_ty}，各行的存取元素为以第n个目标为中心所选取的三维区域中各栅格代表点所对应的三维采样点序号，P_{n_ty}的行向量表示为(P_x(w)，P_y(p)，P_z(q))，其中

w＝1，2，…，2*h+1，p＝1，2，…，2*h+1，q＝1，2，…，2*h+1，P_x(w)表示矩阵P_x的第w个元素，P_y(p)表示矩阵P_y的第p个元素，P_z(q)表示矩阵P_z的第q个元素，矩阵P_{n_ty}表示第n个目标附近选取区域的代表点的集合；

构建投影矩阵P_tyall＝[P_{1_ty}；P_{2_ty}；…；P_{Ntar_ty}]，令P_ty＝unique(P_tyall)，P_ty的行数为H，unique(P_tyall)表示将矩阵P_tyall中代表相同栅格代表点的行元素消去，使各行所代表的栅格代表点各不相同；

步骤7.3构建目标可能存在的区域

根据步骤7.2得到的矩阵P_ty，计算选取的三维区域的代表点坐标P_aero(v，：)＝P0+[P_ty(v，1)-Nx/2，P_ty(v，2)-Ny/2，P_ty(v，3)-Nz/2]*M*δ_r，其中，v为代表点序号，v＝1，2，…，H，H为步骤7.2中定义的P_ty的行数；

步骤7.4投影处理

根据步骤7.3建立的目标可能存在区域的代表点，寻找每一代表点对应的各接收机的回波信号，采用投影方法对每一代表点对应的各个接收机回波信号进行积累；具体步骤如下：

步骤7.4.1

转至步骤5.1，创建三维图像矩阵Ima_Nx×Ny×Nz；然后令接收机序号初始化为j＝1，代表点序号v初始化为v＝1；由步骤5.1转至步骤7.4.2；

步骤7.4.2

根据步骤7.3中选取的区域代表点坐标P_aero(v，：)和第j个接收机位置Prec(j)，计算从发射机P_t经栅格单元代表点P_aero(v，：)到达第j个接收机的双站距离，记为R_v，j，计算公式：R_v，j＝||P_t-P_aero(v，：)||₂+||P_aero(v，：)-Prec(j)||₂，其中，P_aero(v，：)表示第v行所有列的数据，v＝1，2，…，H；

然后，由双站距离R_v，j计算代表点P_aero(v，：)对第j个接收机的距离单元编号，记为ID_v，j，转到步骤7.4.3；

步骤7.4.3

将步骤7.2中所建立的降分辨率回波信息矩阵P_M的第j列第ID_v，j行中的数据作为栅格单元代表点P_aero(v，：)所对应的回波信号，记为P₀(v，j)，P₀(v，j)表示代表点P_aero(v，：)对第j个接收机的回波信号，并将回波信号P₀(v，j)在成像空间中积累，回波信号P₀(v，j)在成像空间积累的表达式为：Ima_Nx×Ny×Nz(k1，k2，k3)＝Ima_Nx×Ny×Nz(k1，k2，k3)+P₀(v，j)，其中k1＝x，k2＝y，k3＝z；令接收机序号j加1，并转到步骤7.4.4；

步骤7.4.4

如果j≤Nrec，则转至步骤7.4.2；

如果j＞Nrec，则j置为1，v加1，转至步骤7.4.5；

步骤7.4.5

如果v≤H，则转至步骤7.4.2；

如果v＞H，则投影结束，得到一个关于成像空间的三维图像矩阵Ima_Nx×Ny×Nz，转至步骤7.5；

步骤7.5提取目标位置

采用步骤6所提方法，获取三维图像矩阵Ima_Nx×Ny×Nz中的Ntar个最大值，并记录Ntar个最大值的序号为n(n＝1，2…，Ntar)；具体步骤如下：

步骤7.5.1

找出三维图像矩阵Ima_Nx×Ny×Nz中的最大值，记做Max_n，转到步骤7.5.2；

步骤7.5.2

转至步骤6.2，获得第n个目标的位置坐标P_n，令接收机序号初始化为j＝1；转至步骤7.5.3；

步骤7.5.3

转至步骤6.3，计算第n个目标相对第j个接收机的回波的距离单元编号ID_n，j；再由步骤6.3转到步骤7.5.4；

步骤7.5.4

根据步骤4建立的降分辨率回波信息矩阵P_M和第n个目标相对第j个接收机的回波的距离单元编号ID_n，j，在降分辨率回波信息矩阵P_M的第j列第ID_n，j行中找出第n个目标对第j个接收机的回波信号，记为P_M(l，c)(l＝ID_n，j，c＝j)；然后，从三维图像矩阵Ima_Nx×Ny×Nz中剔除此回波信号，Ima_Nx×Ny×Nz(k1，k2，k3)＝Ima_Nx×Ny×Nz(k1，k2，k3)-P_M(l，c)，其中k1＝k1_{Max_n}，k2＝k2_{Max_n}，k3＝k3_{Max_n}；令j加1，转到步骤7.5.5；

步骤7.5.5

如果j≤Nrec，则转到至骤7.5.3；

如果j＞Nrec，令n加1，转到步骤7.5.6；

步骤7.5.6

如果n≤Ntar，则转至步骤7.5.1；

如果n＞Ntar，则转至步骤7.1；

步骤8泰勒级数展开定位

根据步骤7求得的目标位置P_n(n＝1，2，…，Ntar)，找到目标位置P_n对应在多个接收机中的距离信息DI_n(k)(k＝1，2，…，NOR_n)，构建时差定位方程组，以P_n为初始值，利用泰勒级数展开定位算法求解目标位置，NOR_n表示满足一定条件的接收机个数，初始化接收机的序号j为1，具体步骤如下：

步骤8.1提取目标距离信息

选择第n(n＝1，2，…，Ntar)个目标位置P_n和第j个接收机，由双站距离公式计算信号从发射机P_t经P_n到达第j个接收机的双站距离，记为R_n，j，R_n，j＝||P_t-P_n||₂+||P_n-Prec(j)||₂；

步骤8.1.1

由R_n，j计算目标n的回波距离单元编号分别求出距离RC_n，j最近的、第2近的和第3近的整数RI_n，j(i)(i＝1，2，3)，初始化i的值为1，n的值为1，k的值为1，转到步骤8.1.2；

步骤8.1.2

取出步骤2中回波概率信息矩阵P的第RI_n，j(i)行第j列的数据，判断该数据是否大于门限ε；如果数据值大于门限ε，则令距离信息DI_n(k)＝RI_n，j(i)，EI_n(k)＝j，k值加1，j值加1，并转到步骤8.1.4；如果数值小于等于门限ε，则转到步骤8.1.3；

步骤8.1.3

如果i＜3，i的值加1，转到步骤8.1.2；

如果i＝3，j的值加1，转到步骤8.1.4；

步骤8.1.4

如果j≤Nrec，则转到步骤8.1；

如果j＞Nrec，n的值加1并转到步骤8.1.5；

步骤8.1.5

如果n≤Ntar，则令j＝1并转到步骤8.1；

如果n＞Ntar，则转到步骤8.2；

步骤8.2构建方程组

根据得到的目标在多个接收机中对应的距离信息，构建定位方程组，初始化目标序号n＝1，具体步骤如下：

步骤8.2.1

根据步骤8.1得到的矩阵DI_n，构造距离差矩阵RD_n，RD_n(m)＝[DI_n(m+1)-DI_n(1)]*M*δ_n+RrecO，其中n＝1，2，...，Ntar，m＝1，2，...，NOR_n-1；取出步骤7求得的目标位置P_n，初始化p为1，初始化系数矩阵A_n为NOR_n-1行3列的零矩阵，常数矩阵B_n为NOR_n-1行1列的零矩阵，并转到步骤8.2.2；

步骤8.2.2

根据步骤8.1得到的矩阵EI_n，令接收机序号j＝EI_n(p)，计算P_n到接收机j的距离d_n(p)＝||P_n-Prec(j)||₂，转到步骤8.2.3；

步骤8.2.3

如果p＝1，p的值加1并转到步骤8.2.2；

如果p≠1，转到步骤8.2.4；

步骤8.2.4

令q＝p-1，系数矩阵A_n的第q行元素为常数矩阵B_n的第q行元素为b_nq＝RD_n(q)-(d_n(p)-d_n(1))；p的值加1并转到步骤8.2.5；

步骤8.2.5

如果p≤NOR_n，转到步骤8.2.2；

如果p＞NOR_n，转到步骤8.3；

步骤8.3求解方程组

根据步骤8.2得到的系数矩阵A_n和常数矩阵B_n，利用传统的最小二乘算法得到未知量δ，其中，n＝1，2，...，Ntar，δ是位置修正量，求得δ后可对目标位置P_n进行修正，判断δ与门限ε的大小，多次迭代得到最终的目标位置，具体步骤如下：

步骤8.3.1

利用传统的最小二乘算法求得位置修正量δ_n＝(A_n ^TA_n)^-1A_n ^TB_n，将目标位置P_n的值修改为P_n+δ_n，令sum(δ_n)表示对δ_n所有元素的求和运算，转到步骤8.3.2；

步骤8.3.2

如果Δ＞ε，转到步骤8.2.1；

如果Δ≤ε，n的值加1并转到步骤8.3.3；

步骤8.3.3

如果n≤Ntar，转到步骤8.2.1；

如果n＞Ntar，转到步骤8.4；

步骤8.4输出定位结果

目标定位结束，得到最终的Ntar个目标的三维位置坐标P_n(n＝1，2，…，Ntar)。