基于无控制点的双机载SAR图像目标定位的方法与流程

本发明属于sar定位技术领域，具体涉及一种基于无控制点的双机载sar图像目标定位的方法。

背景技术：

合成孔径雷达(syntheticapertureradar,sar)具有全天候、全天时的对地观测能力，sar在海洋监测方面得到了广泛的应用，并发挥了巨大的社会、经济和军事效益。因成像波段不同，sar系统可以获取不同于光学系统的地物目标信息。

在一般情况下，传统的sar图像定位是由一幅sar图像，辅以地球模型方程等，在定向参数解算的基础上计算相应地面点的三维信息。然而对于机载平台的雷达对地目标定位条件下，由于飞行高度的限制，地球曲率对机载雷达成像定位的影响可以忽略，则不满足地球模型方程。因此，对目标的三维坐标定位须再加上有限个控制点的信息，这在非合作区域条件下实现可能性较小。本文采用两幅sar图像信息结合定向参数等信息，结合距离-多普勒方程，实现对目标的三维坐标实现无控制点定位，此方法的应用场景更为广泛。

sar图像定位与三维信息提取技术可以对较大范围的区域实现全天时、全天候和高精度的目标定位。现有的双基机载sar技术常用的是两部雷达分别作为发射机和接收机的工作模式，这样的方法无法由目标在sar图像的像素位置信息，通过直接解析方程组获得目标的实际位置信息。

技术实现要素：

为解决现有技术中非合作区域目标定位的问题，本发明的目的在于提供一种能够在无控制点条件下、准确地检测和定位杂波背景sar图像中目标的方法。

实现本发明目的的技术方案为：一种基于无控制点的双机载sar图像单目标定位的方法，包括以下步骤：

对所获取的两幅sar图像分别采用目标连通域分割标记算法，划分出多个目标的像点区域，保存区域像点质心坐标信息；

结合基于sift特征的目标匹配算法，在复杂背景下将sar1、sar2图像进行目标匹配，输出同名目标分别在两幅图像中的像点坐标；

将像点坐标、载机飞行位置、速度信息和sar成像角度信息代入基于无控制点的双sar协同立体定位模型，对目标的实际三维坐标进行牛顿迭代解算。

进一步的，所述的目标连通域分割标记算法具体步骤为：

根据设定的阈值对获取的图像进行二值化，分离前景像素和背景像素，其中前景像素构成待检测连通区域；

对图像进行扫描，当在某一行中找到各待检测连通区域的行连通域时，分别统计这些行连通域的信息并保存，所述行连通域的信息包括行连通域开始的列序号、结束的列序号、连通像素个数、行连通域中所有像素的列序号之和及所有像素的行序号之和；

若在下一行扫描中找到位于各待检测连通区域的行连通域时，将改行的每个行连通域分别与上一行所有行连通域一一进行比较，比较从上一行的最后一个行连通域开始；如果没有一个满足八邻域连通的融合条件，则分配一个标记号给当前行连通域，并将该行连通域的信息保存；如果满足融合条件，则将满足融合条件的上下两个行连通域合并，并结合后形成的连通域分配一个标记号；

对于每一个分配到标记号的连通域，将该连通域的最小行序号、最大行序号、最小列序号、最大列序号、该连通域中所有像素列序号之和及所有像素行序号之和，以标记号为地址保存；

图像扫描完毕，所有行连通域融合完成，根据合并后的连通域信息，采用质心公式计算连通域质心坐标，保存图像各连通区域像素行列序号范围及对应质心坐标信息；

分别对两幅sar图像进行上述图像连通域分割及标记处理，分别输出两幅图像中各目标连通域质心坐标。

进一步的，基于sift特征的目标匹配算法过程具体为：

对两幅sar图像分别进行搜索特征点，并对提取的特征点的梯度方向及模值进行描述，对关键点进行匹配，剔除错配点；然后得到两幅图像中所匹配的同名点的对应关系，匹配点所属的目标连通域范围、质心坐标也可以得到；所述的同名点是指同一个目标分别在两幅图中的像素点位置；输出两幅图中同名目标的质心像素坐标。

进一步的，步骤2具体为：

2a)对于sar图像的二维图像i(x,y)，在不同尺度空间表示为：

l(x,y,σ)＝g(x,y,σ)i(x,y)，高斯核(x,y)代表点坐标，σ代表高斯正态分布的方差；

将不同尺度的高斯差分和图像进行卷积：

d(x,y,σ)＝(g(x,y,kσ)-g(x,y,σ))*i(x,y)＝l(x,y,kσ)-l(x,y,σ)

如果一个点在其26个邻域中是最大值或最小值，则该点判断为该尺度下的一个特征点由此得到图像中的特征点集c；

2b)利用特征点邻域的梯度方向分布特征，为每个特征点指定方向参数，使算子具备旋转不变性；

梯度的模值：

梯度的方向：θ(x,y)＝arctan{[l(x,y+1)-l(x,y-1)]/[l(x+1,y)-l(x-1,y)]}

2c)生成sift特征向量，将坐标轴旋转到关键点方向，以保持旋转不变性；每个特征点采用16种子点来描述，产生128个数据；

2d)特征点找出后寻找图像特征点的对应关系，采用最近邻法，即找每一个特征点在另一幅图像中的最近邻点；

假设两个特征点的特征向量分别为(a1,a2,...an)和(b1,b2,...bn)，则此两点之间的斜面距离可以表示为i∈(1,2,...n)，n为维数；比较该两点之间的最邻近umin和次邻近距离ul，当满足条件umin/ul＜r及距离比例阈值0＜r≤1时，判定为正确匹配点，否则为错误匹配。

2e)将两幅sar图像进行上述目标匹配处理，得到两幅图中相同目标的像点坐标对应组合，结合目标连通域坐标范围进行搜索，找到所匹配目标分别在sar1、sar2图中的标号区域，进一步得到相应质心t1(i^l,j^l)、t2(i^r,j^r)。

进一步的，基于无控制点的双sar协同立体定位模型具体为：

将两组同名点像点坐标代入距离多普勒成像模型，结合sar1、sar2的定向参数及sar系统参数，采用牛顿迭代算法最后迭代出最符合实际条件的目标三维坐标；其中定向参数包括机载sar的飞行速度和实时坐标，sar系统参数包括雷达发射波俯仰角。

进一步的，步骤3对目标的实际三维坐标进行牛顿迭代解算，具体为：

3a)由两幅sar图像中同名点的像点t1(i^l,j^l)、t2(i^r,j^r)代入距离公式及多普勒公式，得到同名像点坐标与对应地面点坐标(x,y,z)的关系由如下四个方程构成的方程组表示：

即：

其中，分别代表像点t1、像点t2成像瞬间sar1、sar2天线相位中心位置，分别代表像点t1、像点t2成像瞬间sar1、sar2天线相位中心速度，分别为sar1、sar2成像时的近距延迟，分别为sar1、sar2的斜距向采样间隔；

其中，sar1、sar2系统的多普勒频率分别为发射信号波长分别为λ^l、λ^r，发射信号俯仰角为α^l和α^r；

3b)计算同名点成像瞬间两天线相位中心位置、速度

若以二次多项式表示天线相位中心位置与成像时刻的关系，则由解算得到的定向参数，利用下式可以分别求得像点t1、像点t2成像瞬间天线相位中心位置和成像瞬间天线相位中心速度

式中，t'为每行间的时间间隔；t^l和t^r分别为像点t1、像点t2在左右图像上的成像时刻；分别为sar1、sar2的天线相位中心加速度初值；分别为sar1、sar2的天线相位中心速度初值；分别为sar1、sar2的天线相位中心位置初值。

3c)构建误差方程组

由r-d模型对地面点坐标的线性化形式，可知同名像点t1(i^lj^l)、t2(i^r,j^r)与所对应的地面点p(x,y,z)的线性化关系为：

c·δg-l＝0

其中，c是关于地面点坐标改正量的系数矩阵即：

式中，分别是关于像点t1成像时刻对应的sar1位置的函数；分别是关于像点t1成像时刻对应的sar1速度的函数；分别是关于像点t2成像时刻对应的sar2位置的函数；分别是关于像点t2成像时刻对应的sar2速度的函数；

系数阵c中的各元素分别为：

δg是地面点坐标的改正量向量，δg＝[δxδyδz]^t；

l是r-d模型方程组的初值向量，

3d)计算三维坐标改正量

c矩阵的法方程式可表示为：

c^tcδg-c^tl＝0

解法方程，可求得地面点三维坐标的改正量向量δg：

δg＝(c^tc)^-1c^tl

再在上一次迭代的基础上修正三维坐标初值：

3e)限差判断

判断改正量是否小于给定限差，若改正量大于限差则返回步骤(3d)，利用改正后的三维坐标重新组误差方程式计算其改正量；若改正量小于或等于限差则停止迭代，输出计算出的地面点三维坐标。

本发明与现有技术相比，具有以下优点：(1)本发明能够更加有效地对sar图像中目标进行检测及定位；采用基于sift特征的目标检测及匹配方法，可以对尺度不同、图像旋转等情况下的图像进行有效的匹配；采用两机载sar数据融合定位的三维解算方法，可以不受相干性的限制，得到高精度定位结果；(2)双sar协同立体定位方法是基于两sar分别进行收发信号及成像的技术，可以在得到目标在图像中的像点坐标后直接进行解算，从而大大减少计算量，提高实时性。

附图说明

图1是本发明中双机载sar对目标协同成像立体定位结构示意图。

图2是本发明中双机载sar对目标协同成像立体定位总流程图。

图3是本发明中目标检测算法流程图。

图4是本发明中地面点三维坐标解算算法流程图。

具体实施方式

一种基于无控制点的双机载sar图像目标定位的方法，具体步骤如下：

(1)对所获取的两幅sar图像分别采用目标连通域分割标记算法，划分出多个目标的像点区域，保存区域像点质心坐标信息；

(2)结合基于sift特征的目标匹配算法，在复杂背景下将sar1、sar2图像进行目标匹配，输出同名目标分别在两幅图像中的像点坐标；

(3)将像点坐标、载机飞行位置、速度信息和sar成像角度信息等代入基于无控制点的双sar协同立体定位模型，对目标的实际三维坐标进行牛顿迭代解算。

所述的目标连通域分割标记算法，具体步骤为：

首先对两幅获取的sar图像同时进行目标连通区域分割及标记处理，分别对两幅图像中一个或多个目标区域进行标号，并保存对应区域的像点坐标范围和目标质心像点坐标。

所述的基于sift特征的目标匹配算法，具体步骤为：

对两幅sar图像分别进行搜索特征点、并对提取的特征点的梯度方向及模值进行描述，对目标区域特征进行匹配，剔除错配点。然后得到两幅图像中所匹配的同名目标的对应关系，匹配点所属的目标连通域范围、质心坐标也可以得到。这里所述的同名目标指的是同一个目标分别在两幅图中的像素点位置。输出两幅图中同名目标的一组质心像素坐标。

所述的双机载sar对目标协同成像立体定位的方法，其特征在于，基于无控制点的双sar协同立体定位模型具体为：

将两组同名目标质心像点坐标代入距离多普勒成像模型，结合sar1、sar2的定向参数及sar系统参数，采用牛顿迭代算法最后迭代出最符合实际条件的目标三维坐标；其中定向参数包括机载sar的飞行速度和实时坐标，sar系统参数包括雷达发射波俯仰角。

所述的距离多普勒模型包括距离公式及多普勒频率方程：

距离公式：rs²＝(x-xs)²+(y-ys)²+(z-zs)²＝(r0+mslant·j)²

可记：f1＝(x-xs)²+(y-ys)²+(z-zs)²-(r0+mslant·j)²

其中，(x,y,z)表示地面点目标坐标，(xs,ys,zs)为成像瞬间天线相位中心位置，r0为近距延迟，mslant为斜距向采样间隔，j为像点的距离向坐标。

多普勒频率方程：

可记：

其中，(vx,vy,vz)表示成像瞬间天线相位中心速度，λ为雷达发射波波长，rs为地面点目标到雷达平台的瞬时位置距离，fdc为多普勒频移参数，式中，r表示目标相对于飞机的瞬时斜距，v表示目标相对于飞机的瞬时速度，表示目标相对飞机的位置矢量，表示目标相对飞机的速度矢量，α代表雷达发射信号俯仰角。

所述的牛顿迭代解算三维坐标的流程为计算同名点成像瞬间天线相位中心位置及速度，代入距离多普勒基本方程并构建误差方程。

下面参照附图，以杂波海洋环境下定位舰船目标为例，对本发明的示例性实施方式进行详细描述。

实施例

图1为本发明提出的一种基于无控制点的双机载sar图像目标定位的结构示意图，主要由两架固定翼无人机作为平台的两部sar雷达组成，两架无人机以相同方向在目标区域两侧飞行，同时对扫描区域进行成像。

图2为本发明提出的一种基于无控制点的双机载sar图像目标定位的总流程图，具体实施步骤为：

第一步，获取距离多普勒成像图像，对两个机载平台sar(分别为sar1、sar2)接收到的海面扫描区域的回波数据进行高分辨率成像，分辨率达到1米。采用目标连通域分割标记算法，对获取的两幅sar图像进行相应处理，划分为多个目标连通区域，保存区域像点质心坐标信息。

所述的目标连通域分割标记算法，具体步骤为：

1a)根据设定的阈值对获取的图像进行二值化，分离前景像素和背景像素，其中前景像素构成待检测连通区域；对图像进行扫描，当在某一行中找到各待检测连通区域的行连通域时，分别统计这些行连通域的信息并保存，所述行连通域的信息包括行连通域开始的列序号、结束的列序号、连通像素个数、行连通域中所有像素的列序号之和及所有像素的行序号之和；

1b)若在下一行扫描中找到位于各待检测连通区域的行连通域时，将改行的每个行连通域分别与上一行所有行连通域一一进行比较，比较从上一行的最后一个行连通域开始；如果没有一个满足八邻域连通的融合条件，则分配一个标记号给当前行连通域，并将该行连通域的信息保存；如果满足融合条件，则将满足融合条件的上下两个行连通域合并，并结合后形成的连通域分配一个标记号；

1c)对于每一个分配到标记号的连通域，将该连通域的最小行序号、最大行序号、最小列序号、最大列序号、该连通域中所有像素列序号之和及所有像素行序号之和，以标记号为地址保存；

1d)图像扫描完毕，所有行连通域融合完成，根据合并后的连通域信息，采用质心公式计算连通域质心坐标，保存图像各连通区域像素行列序号范围及对应质心坐标信息。

分别对两幅sar图像进行上述图像连通域分割及标记处理，分别输出sar1图像中的目标连通域像点范围矩阵il1、il2......ilm，对应区域质心坐标tl1(il1,jl1)、tl2(il2,jl2)……tlm(ilm,jlm)，sar2图像中的目标连通域像点范围矩阵ir1、ir2......irm，对应区域质心坐标tr1(ir1,jr1)、tr2(ir2,jr2)……trm(irm,jrm)。其中，(l1,l2......lm)指sar1图中目标区域的标号，(r1,r2......rm)指sar2图中目标区域的标号。

第二步，采用基于sift特征的目标匹配算法，流程图如图3所示，对sar1、sar2所成像进行杂波背景下舰船目标检测。所述的目标检测算法，具体步骤为：

2a)对于sar图像的二维图像i(x,y)，在不同尺度空间表示为：

l(x,y,σ)＝g(x,y,σ)i(x,y)，高斯核(x,y)代表点坐标，σ代表高斯正态分布的方差。

为了在尺度空间上检测到稳定的特征点，需要采用高斯差分尺度空间(dog)，即将不同尺度的高斯差分和图像进行卷积：

d(x,y,σ)＝(g(x,y,kσ)-g(x,y,σ))*i(x,y)＝l(x,y,kσ)-l(x,y,σ)

如果一个点在其26个邻域中是最大值或最小值，则该点判断为该尺度下的一个特征点由此得到图像中的特征点集c。

2b)利用特征点邻域的梯度方向分布特征，为每个特征点指定方向参数，使算子具备旋转不变性。

梯度的模值：

梯度的方向：θ(x,y)＝arctan{[l(x,y+1)-l(x,y-1)]/[l(x+1,y)-l(x-1,y)]}

2c)生成sift特征向量，将坐标轴旋转到关键点方向，以保持旋转不变性。每个特征点采用16种子点来描述，则可以产生128个数据，此128维的特征描述向量，对光照、噪声、旋转和尺度都具有良好的不变性。

2d)特征点找出后寻找图像特征点的对应关系，采用最近邻法，即找每一个特征点在另一幅图像中的最近邻点。在理想状态下，两幅图像之间相同部分的特征点应该具有相同的特征点应该具有相同的特征描述向量，则距离最近。

假设两个特征点的特征向量分别为(a1,a2,...an)和(b1,b2,...bn)，则此两点之间的斜面距离可以表示为i∈(1,2,...n)，n为维数。比较该两点之间的最邻近umin和次邻近距离ul，当满足条件umin/ul＜r及距离比例阈值0＜r≤1时，判定为正确匹配点，否则为错误匹配。

2e)将两幅sar图像进行上述目标匹配处理，得到两幅图中相同目标的像点坐标对应组合，结合目标连通域坐标范围进行搜索，找到所匹配目标分别在sar1、sar2图中的标号区域，进一步得到相应质心t1(i^l,j^l)、t2(i^r,j^r)。

第三步，采用地面点三维坐标解算算法，对舰船目标的实际坐标进行解算。流程图如图4所示，具体步骤为：将两组同名点像点坐标代入距离多普勒成像模型，结合sar1、sar2的定向参数及sar系统参数，采用牛顿迭代算法解算三维坐标的流程为计算同名点成像瞬间天线相位中心位置及速度，代入距离多普勒基本方程并构建误差方程，最后迭代出最符合实际条件的目标三维坐标。

所述的距离多普勒模型导入牛顿迭代解算三维坐标，具体步骤为：

3a)由两幅sar图像中同名点的像点t1(i^l,j^l)、t2(i^r,j^r)代入距离公式及多普勒公式，得到同名像点坐标与对应地面点坐标(x,y,z)的关系由如下四个方程构成的方程组表示：

即：

其中，分别代表像点t1、像点t2成像瞬间sar1、sar2天线相位中心位置，分别代表像点t1、像点t2成像瞬间sar1、sar2天线相位中心速度，分别为sar1、sar2成像时的近距延迟，分别为sar1、sar2的斜距向采样间隔。

其中，sar1、sar2系统的多普勒频率分别为发射信号波长分别为λ^l、λ^r，发射信号俯仰角为α^l和α^r。

3b)计算同名点成像瞬间两天线相位中心位置、速度

若以二次多项式表示天线相位中心位置与成像时刻的关系，则由解算得到的定向参数，利用下式可以分别求得像点t1、像点t2成像瞬间天线相位中心位置和成像瞬间天线相位中心速度

式中，t'为每行间的时间间隔；t^l和t^r分别为像点t1、像点t2在左右图像上的成像时刻；分别为sar1、sar2的天线相位中心加速度初值；分别为sar1、sar2的天线相位中心速度初值；分别为sar1、sar2的天线相位中心位置初值。

3c)构建误差方程组

由r-d模型对地面点坐标的线性化形式，可知同名像点t1(i^lj^l)、t2(i^r,j^r)与所对应的地面点p(x,y,z)的线性化关系为：

c·δg-l＝0

其中，c是关于地面点坐标改正量的系数矩阵即：

式中，分别是关于像点t1成像时刻对应的sar1位置的函数；分别是关于像点t1成像时刻对应的sar1速度的函数；分别是关于像点t2成像时刻对应的sar2位置的函数；分别是关于像点t2成像时刻对应的sar2速度的函数。

系数阵c中的各元素分别为：

δg是地面点坐标的改正量向量，δg＝[δxδyδz]^t。

l是r-d模型方程组的初值向量，

3d)计算三维坐标改正量

c矩阵的法方程式可表示为：

c^tcδg-c^tl＝0

解法方程，可求得地面点三维坐标的改正量向量δg：

δg＝(c^tc)^-1c^tl

再在上一次迭代的基础上修正三维坐标初值：

3e)限差判断

判断改正量是否小于给定限差，若改正量大于限差则返回(3d)，利用改正后的三维坐标重新组误差方程式计算其改正量；若改正量小于或等于限差则停止迭代，输出计算出的地面点三维坐标。