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

权利要求书

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)。