CN113671497B - 基于圆柱对称模型的单通道sar目标三维坐标提取方法 - Google Patents
基于圆柱对称模型的单通道sar目标三维坐标提取方法 Download PDFInfo
- Publication number
- CN113671497B CN113671497B CN202110763425.8A CN202110763425A CN113671497B CN 113671497 B CN113671497 B CN 113671497B CN 202110763425 A CN202110763425 A CN 202110763425A CN 113671497 B CN113671497 B CN 113671497B
- Authority
- CN
- China
- Prior art keywords
- dimensional
- sub
- aperture
- image
- distance
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000605 extraction Methods 0.000 title claims abstract description 22
- 238000003384 imaging method Methods 0.000 claims abstract description 46
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 35
- 239000011159 matrix material Substances 0.000 claims abstract description 28
- 230000009466 transformation Effects 0.000 claims abstract description 21
- 238000000034 method Methods 0.000 claims abstract description 20
- 230000005012 migration Effects 0.000 claims abstract description 11
- 238000013508 migration Methods 0.000 claims abstract description 11
- 238000012545 processing Methods 0.000 claims abstract description 11
- 239000013598 vector Substances 0.000 claims description 57
- 238000001914 filtration Methods 0.000 claims description 10
- 230000006835 compression Effects 0.000 claims description 9
- 238000007906 compression Methods 0.000 claims description 9
- 230000035485 pulse pressure Effects 0.000 claims description 9
- 230000004044 response Effects 0.000 claims description 9
- 238000005070 sampling Methods 0.000 claims description 6
- 230000014509 gene expression Effects 0.000 claims description 5
- 230000001629 suppression Effects 0.000 claims description 5
- 230000000007 visual effect Effects 0.000 claims description 5
- 230000005540 biological transmission Effects 0.000 claims description 3
- 238000004364 calculation method Methods 0.000 claims description 3
- 238000002592 echocardiography Methods 0.000 claims description 3
- 230000001131 transforming effect Effects 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 abstract description 5
- 238000004088 simulation Methods 0.000 description 3
- 230000008859 change Effects 0.000 description 2
- 230000001427 coherent effect Effects 0.000 description 2
- 230000008878 coupling Effects 0.000 description 2
- 238000010168 coupling process Methods 0.000 description 2
- 238000005859 coupling reaction Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 208000004350 Strabismus Diseases 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
- G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
- G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
- G01S13/9021—SAR image post-processing techniques
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- Electromagnetism (AREA)
- Computer Networks & Wireless Communication (AREA)
- General Physics & Mathematics (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明属于雷达技术领域,公开了一种基于圆柱对称模型的单通道SAR目标三维坐标提取方法,在聚束模式下,将合成孔径分为两个子孔径,利用BPA算法和距离徙动算法分别对两个子孔径进行成像处理,得到聚焦于不同成像平面的同一场景的二维SAR图像对(包含主图像和辅助图像);利用SURF算法获得主图像与辅助图像之间的转换矩阵T;并利用转换矩阵T得到真实三维目标在两个子孔径图像中的二维坐标对;根据两幅子孔径SAR成像系统的几何特点,构建圆柱对称模型,利用得到的真实三维目标的二维坐标对,从SAR图像对中提取目标的三维坐标。该方法中模型求解简单,重构精度高,具有较强的实用性,适用于各种曲线路径的SAR平台。
Description
技术领域
本发明涉及雷达技术领域,具体涉及一种基于圆柱对称模型的单通道SAR目标三维坐标提取方法。
背景技术
合成孔径雷达(SAR)作为一种有源对地探测系统,通过雷达与装载平台之间的相对运动,在一定的积累时间内,将雷达在不同空间位置上接收到的宽带回波信号进行相干处理获得目标的二维图像,从而使人们真正看到目标的真实图像,具有全天时、全天候的工作特点。但是最终的SAR图像是二维的,不能得到具体的三维坐标,随着高分辨率合成孔径雷达(SAR)系统的快速发展,利用二维SAR成像的三维坐标提取技术近年来得到了广泛关注。常用的三维坐标提取方法是干涉SAR(InSAR)。然而,合成孔径雷达是一种基于相干相位的技术,存在着飞行轨迹数、相位展开误差等问题。
发明内容
针对现有技术中存在的问题,本发明的目的在于提供一种基于圆柱对称模型的单通道SAR目标三维坐标提取方法,该方法中模型求解简单,重构精度高,具有较强的实用性,适用于各种曲线路径的SAR平台。
本发明的主要技术思路为:首先,在聚束模式下,将合成孔径分为两个子孔径,利用反投影算法(BPA)和距离徙动算法分别对两个子孔径的图像进行成像处理,得到聚焦于不同成像平面的同一场景的二维SAR图像对(包含主图像和辅助图像)。其次,利用SURF算法获得主图像与辅助图像之间的转换矩阵T。接着,利用转换矩阵T,得到真实三维目标在两个不同视角(主图像和辅助图像)的子孔径图像中的二维坐标对。最后根据两幅子孔径SAR成像系统的几何特点,构建圆柱对称模型,利用得到的真实三维目标的二维坐标对,从SAR图像对中提取目标的三维坐标。
为了达到上述目的,本发明采用以下技术方案予以实现。
基于圆柱对称模型的单通道SAR目标三维坐标提取方法,包括以下步骤:
步骤1,获取全孔径雷达回波数据,根据雷达飞行轨迹将全孔径分为两个不同视角的子孔径k1、k2;
利用距离徙动算法对子孔径k1的回波信号在斜距平面进行成像处理,得到二维SAR辅助图像I1;利用反投影算法对子孔径k2的回波信号在地距平面进行成像处理,得到二维SAR主图像I2;
步骤2,利用快速稳健特征算法获得所述二维SAR辅助图像I1与所述二维SAR主图像I2之间的转换矩阵T;
步骤3,利用转换矩阵T,得到真实三维目标在两个子孔径图像中的二维坐标对;
步骤4,建立辅助图像圆柱对称模型和主图像圆柱对称模型,根据所述真实三维目标在两个子孔径图像中的二维坐标对,计算目标的真实三维坐标。
本发明技术方案的特点和进一步的改进为:
(1)步骤1中,所述距离徙动算法包含以下子步骤:
子步骤a),对子孔径k1的回波信号进行二维傅里叶变换,得到子孔径k1距离频域方位频域回波信号;
子步骤b),利用一致补偿函数,对所述子孔径k1距离频域方位频域回波信号进行一致补偿,得到子孔径k1一致粗补偿后二维频域回波信号;
子步骤c),对所述子孔径k1一致粗补偿后二维频域回波信号进行stold插值运算,得到子孔径k1位于斜距平面的二维SAR图像;
子步骤d),对所述子孔径k1位于斜距平面的二维SAR图像进行二维傅里叶逆变换,得到二维SAR辅助图像I1。
(2)步骤1中,所述反投影算法包含以下子步骤:
子步骤a),子孔径k2的回波信号进行距离向傅里叶变换,得到子孔径k2的距离频域方位时域回波信号;
子步骤b),利用距离脉压函数对所述子孔径k2的距离频域方位时域回波信号进行脉冲压缩处理,得到子孔径k2脉冲压缩后的回波信号;
子步骤c),对所述子孔径k2脉冲压缩后的回波信号进行距离向傅里叶逆变换,得到子孔径k2的二维时域回波信号;
子步骤d),将地距平面均匀的划分成网格,确定每一个网格到雷达平台的距离信息,根据该距离信息对所述子孔径k2的二维时域回波信号进行反投影处理,得到二维SAR主图像I2。
(3)所述傅里叶变换是利用以下公式完成的:
s1(f,t)=∫s0(τ,t)exp(-j2πfτ)dτ
其中,s1(f,t)表示距离频域为f、方位时间为t的距离频域回波信号,f的取值范围为fs表示SAR图像成像系统的采样频率,t的取值范围为(0,T),T表示SAR图像成像系统的运行时间,∫(·)dτ表示对距离时间τ进行积分操作,s0(τ,t)表示距离时域为τ、方位时间为t的子孔径回波信号,τ的取值范围为/>ro表示三维SAR图像成像区域中目标中心点的斜距,c表示光速,k表示SAR图像成像系统的采样点数,B表示SAR图像成像系统的带宽,exp(·)表示以自然常数e为底的指数操作,j表示虚数单位符号,π表示圆周率。
(4)所述距离脉压函数如下:
其中,s2(K,t)表示距离频域为K、方位时间为t的脉压后距离频域回波信号,K的取值范围为γ表示SAR图像成像系统发射信号的距离调频率。
(5)所述傅里叶逆变换是利用下式完成的:
s3(τ1,t)=∫s2(K,t)exp(j2πKτ)dK
其中,s3(τ1,t)表示距离时域为τ1、方位时间为t的子孔径脉冲压缩后距离时域信号,τ1的取值范围为∫(·)dK表示对距离频域K进行积分操作。
(6)所述一致补偿函数的相位为:
其中,f0为中心频率,fτ为距离向频率,fη为方位向频率,c为光速,v为平台速度,θ为斜视角,R0为中心斜距,Rref为参考斜距,θRFM为一致补偿相位;
所述stold插值运算是利用下式完成的:
其中,fτ'为stold插值后的新的域。
(7)步骤2中,所述快速稳健特征算法为:对主图像I2和辅助图像I1进行特征点提取、描述、匹配,通过两幅图像的特征点匹配对求出转换矩阵T,具体包含以下子步骤:
子步骤2.1,特征点提取:
(1)将待处理的图像变换成积分图像;
(2)建立基于箱式滤波器的特征点检测器作为卷积核,箱式滤波器矩阵表达式定义如下:
其中,Hopprox(x,σ)为利用箱式滤波器得到的过滤点,x为积分图像的坐标,σ为常数,Dxx(x,σ)、Dxy(x,σ)、Dyy(x,σ)分别为xx方向、yy方向、xy方向的偏导;
(3)改变卷积核的尺度大小,将不同尺度的卷积核与积分图像进行卷积,得到金字塔型的尺度空间,采用邻域非极大值抑制法提取所述金字塔型的尺度空间中的特征点;
子步骤2.2,特征描述符建立:
计算特征点周围的圆形邻域内所有点的haar小波响应,获取该特征点的主向量方向;
沿主向量方向构建正方形区域,将所述正方形区域均匀拆分为4个方形子区域,每个方形子区域都有其特定的四维描述符v0:
v0=(∑dx,∑dy,∑|dx|,∑|dx|)
其中,dx表示水平方向的Haar波响应;dy表示朝垂直方向的哈波响应;在所有4×4方形子区域中,描述符向量由64个元素组成;
子步骤2.3,特征点匹配并计算转换矩阵T:
描述符确认后,特征点的相似性可以通过Eculidean距离表示其中i表示描述符向量,xi表示主图像描述符,xi'表示辅助图像描述符,Eculidean距离最小的两个点形成一个特征点匹配对;计算从主图像到辅助图像的转换矩阵T:
其中,(xp,yp)、(yp-sla,Rp-sla)分别表示主图像和辅助图像中特征点的二维坐标;转换矩阵T包含六个未知参数,需要选择至少三组匹配度高的特征点对来计算转换矩阵T。
(8)所述领域非极大值抑制法为:将箱式滤波器得到的过滤点与金字塔型的尺度空间的其他点进行极大值对比,若过滤点的极大值最大,则该过滤点为提取尺度空间中的特征点。
(9)步骤4具体包含以下子步骤:
子步骤4.1,辅助图像圆柱对称模型的建立
散射体D、P、C位于三维观测场景中的同一方位角单元,其中C为位于地平面和散射体的参考散射体,D和P位于相同的斜距单元,则C、D、P点的瞬时倾斜范围RC(η)、RD(η)、RP(η)分别如下:
其中,Rc为参考散射体最小向量,v为平台速度,η为慢时间,y是三个散射体的方位角坐标,R0表示散射体D、P的最小倾斜范围;
由上式可知,从散射体D和散射体P接收到的回波是相同的,因此,辅助图像的圆柱对称模型可以表示为:
其中,H表示飞机飞行高度,x0、y0、z0表示三维观测场景的三维坐标点,Rp-sia表示二维观测场景距离向长度,yp-sia表示二维观测场景的方位向长度;
子步骤4.2,主图像圆柱对称模型的建立
散射体Pt=(x0,y0,z0)是三维观测空间中的任意散射体,根据反向投影算法的成像几何,主图像被聚焦于地平面,因此,基于圆柱对称特性,Pt将被聚焦在地平面中Pt'=(xp,yp,0)的位置,则从雷达平台到散射体的瞬时倾斜距离矢量可表示为:
|PrPt(η)|=|DPt-DPr(η)|
|PrPt'(η)|=|DPt'-DPr(η)|
其中,PrPt(η)为散射体Pt到飞机坐标Pr的距离矢量,DPt为散射体Pt到主图像子孔径基线AB的最小斜距向量,DPr(η)为载机飞行轨迹的瞬时矢量,PrPt'(η)为散射点Pt'到飞机坐标Pr的距离矢量,DPt'为散射点Pt'到主图像子孔径基线AB的最小斜距向量;Pt'为Pt的投影,随着慢时间η的变化,散射体Pt'与Pt的最小斜距向量必须相等,即:
|DPt|=|DPt′|
作为最小斜距向量,DPt'和DPt均与主图像子孔径基线AB正交,因此,三角形PtAB的面积与三角形Pt'AB的面积相等:
根据矢量原理,三角形的面积也可以表示为三角形两边的叉积:
其中,APt、APt'表示雷达平台初始位置到散射体Pt和Pt'的矢量;“×”算子是两个向量的叉乘运算;因此,根据三角形PtAB的面积与三角形Pt'AB的面积相等,则:
|APt×AB|=|APt'×AB|
借助雷达平台上的高精度INS和GPS系统,上式中的向量可以表示为坐标形式:
其中,i、j、k为三个互相垂直的单位向量,x0、y0、z0为Pt的三维点坐标,xA、yA、zA为A的三维点坐标,xB、yB、zB为B的三维点坐标,xP、yP、zP为Pt'的三维点坐标;将|APt×AB|=|APt'×AB|展开成三维坐标表示形式,即得主图像圆柱对称模型为:
其中,C1、C2表示常数,其他系数为:
α1=(zB-zA)2+(yB-yA)2
α2=(xB-xA)2+(zB-zA)2
α3=(xB-xA)2+(yB-yA)2
γ1=-2(yB-yA)(xB-xA)
γ2=-2(zB-zA)(yB-yA)
γ3=-2(zB-zA)(xB-xA)
子步骤4.3,结合所述辅助图像的圆柱对称模型和所述主图像圆柱对称模型,可以提取出目标的真实三维坐标(x0,y0,z0)。
与现有技术相比,本发明的有益效果为:
在聚束模式下,将完整的孔径分为两个不同视角的子孔径k1、k2,利用距离徙动算法对k1进行成像处理,获得二维SAR辅助图像I1,利用BPA算法对k2进行成像处理,得到二维SAR主图像I2;利用SAR成像系统的几何和特点,提出了一个圆柱形对称模型(CSM),通过二维SAR系统的圆柱对称性从SAR图像对中提取目标真正的三维坐标。与已发表的利用距离和多普勒模型(RDM)或有理多项式系数模型(RPC)的雷达测量工作相比,本发明的方法更实用,适用于各种曲线路径的SAR平台,且模型求解简单,重构精度高,具有较强的实用性。
附图说明
下面结合附图和具体实施例对本发明做进一步详细说明。
图1为基于圆柱对称模型的单通道SAR目标三维坐标提取方法的流程图;
图2(a)为SAR二维辅助图像成像几何图;
图2(b)为SAR二维主图像成像几何图;
图3(a)为斜距平面点目标SAR图;
图3(b)为地距平面点目标SAR图;
图4为利用本发明的方法获取的三维场景图。
具体实施方式
下面将结合实施例对本发明的实施方案进行详细描述,但是本领域的技术人员将会理解,下列实施例仅用于说明本发明,而不应视为限制本发明的范围。
基于圆柱对称模型的单通道SAR目标三维坐标提取方法,包括以下步骤:
步骤1,获取全孔径雷达回波数据,根据雷达飞行轨迹将全孔径分为两个不同视角的子孔径k1、k2;利用距离徙动算法对子孔径k1的回波信号在斜距平面进行成像处理,得到二维SAR辅助图像I1;利用反投影(BPA)算法对子孔径k2的回波信号在地距平面进行成像处理,得到二维SAR主图像I2;所述二维SAR辅助图像I1与所述二维SAR主图像I2构成二维SAR图像对。
具体的,反投影算法包含以下子步骤:
子步骤a),对SAR图像成像系统实时接收的曲线轨迹子孔径k2的回波信号进行距离向傅里叶变换,得到子孔径k2的距离频域方位时域回波信号;具体傅里叶变换是利用以下公式完成的:
s1(f,t)=∫s0(τ,t)exp(-j2πfτ)dτ (1)
其中,s1(f,t)表示距离频域为f、方位时间为t的距离频域回波信号,f的取值范围为fs表示SAR图像成像系统的采样频率,t的取值范围为(0,T),T表示SAR图像成像系统的运行时间,∫(·)dτ表示对距离时间τ进行积分操作,s0(τ,t)表示距离时域为τ、方位时间为t的子孔径回波信号,τ的取值范围为/>ro表示三维SAR图像成像区域中目标中心点的斜距,c表示光速,k表示SAR图像成像系统的采样点数,B表示SAR图像成像系统的带宽,exp(·)表示以自然常数e为底的指数操作,j表示虚数单位符号,π表示圆周率。
子步骤b),利用距离脉压函数对所述子孔径k2的距离频域方位时域回波信号进行脉冲压缩处理,得到子孔径k2脉冲压缩后的回波信号;具体距离脉压函数如下:
其中,s2(K,t)表示距离频域为K、方位时间为t的脉压后距离频域回波信号,K的取值范围为γ表示SAR图像成像系统发射信号的距离调频率。
子步骤c),对所述子孔径k2脉冲压缩后的回波信号进行距离向傅里叶逆变换,得到子孔径k2的二维时域回波信号;具体傅里叶逆变换是利用下式完成的:
s3(τ1,t)=∫s2(K,t)exp(j2πKτ)dK (3)
其中,s3(τ1,t)表示距离时域为τ1、方位时间为t的子孔径脉冲压缩后距离时域信号,τ1的取值范围为∫(·)dK表示对距离频域K进行积分操作。
子步骤d),将地距平面均匀的划分成网格,确定每一个网格到雷达平台的距离信息,根据该距离信息对所述子孔径k2的二维时域回波信号进行反投影处理,得到位于地距平面的二维SAR图像,即为二维SAR主图像I2。
具体的,距离徙动算法包含以下子步骤:
子步骤a),对SAR图像成像系统实时接收的曲线轨迹子孔径k1的回波信号进行二维傅里叶变换,得到子孔径k1距离频域方位频域回波信号;具体傅里叶变换的公式同反投影算法的公式(1)。
子步骤b),利用一致补偿函数,对所述子孔径k1距离频域方位频域回波信号进行一致补偿,得到子孔径k1一致粗补偿后二维频域回波信号;其中一致补偿函数的相位为下式:
其中,f0为中心频率,fτ为距离向频率,fη为方位向频率,c为光速,v为平台速度,θ为斜视角,R0为中心斜距,Rref为参考斜距,θRFM为一致补偿相位。
子步骤c),对所述子孔径k1一致粗补偿后二维频域回波信号进行stold插值运算,得到子孔径k1位于斜距平面的二维SAR图像;具体stold插值变换是利用下式完成的:
其中,f0为中心频率,fτ为距离向频率,fη为方位向频率,fτ'为stold插值后的新的域,c为光速,v为平台速度。左式表示未进行stold插值前的相位,其中距离频域与方位频域耦合十分严重;右式表示进行stold插值后的相位,可以看到插值后的信号未单频信号,完全消除了距离向与方位向的耦合。
子步骤d),对所述子孔径k1位于斜距平面的二维SAR图像进行二维傅里叶逆变换,得到子孔径k1二维时域信号,即为二维SAR辅助图像I1。具体傅里叶逆变换的公式同反投影算法的公式(3)。
步骤2,利用快速稳健特征(SURF)算法获得主图像与辅助图像之间的转换矩阵T。
具体的,SURF算法为:对主图像I2和辅助图像I1进行特征点提取、描述、匹配,通过两幅图像的特征点匹配对求出转换矩阵T;具体如下:
1、特征点提取:
(1)将待处理的图像变换成积分图像;
(2)建立基于箱式滤波器的特征点检测器作为卷积核,箱式滤波器矩阵表达式定义如下:
其中,x为积分图像的坐标,σ为1.2,Dxx(x,σ)、Dxy(x,σ)、Dyy(x,σ)分别为xx方向、yy方向、xy方向的偏导。Hopprox(x,σ)为利用箱式滤波器得到的过滤点。
(3)改变卷积核的尺度大小,将不同尺度的卷积核与积分图像进行卷积,得到金字塔型的尺度空间,在3×3×3尺度空间上采用邻域非极大值抑制法(将过滤点与3×3×3尺度空间的其他点进行极大值对比,若选取点的极大值最大,则为特征点)提取尺度空间中的特征点。
2、特征描述符建立:
计算特征点周围半径为6s的圆形邻域内所有点的haar小波响应,获取该特征点的主向量方向。
沿主向量方向构建正方形区域,将此区域均匀拆分为4个方形子区域。每个子区域都有其特定的四维描述符v0:
v0=(∑dx,∑dy,∑|dx|,∑|dx|) (7)
其中,dx表示水平方向的Haar波响应,此处的"水平"定义为所选特征点的方向。dy表示朝垂直方向的哈波响应。|dx|是反应|dy|的绝对值。因此,在所有4×4方形子区域中,描述符向量由64个元素组成。
3、特征点匹配并计算转换矩阵T:
描述符确认后,特征点的相似性可以通过Eculidean距离表示其中i表示描述符向量,xi表示主图像描述符,xi'表示辅助图像描述符,Eculidean距离最小的两个点形成一个特征点匹配对。此外,RANSAC算法用于消除特征点的误匹配。最后,可以计算从主图像到辅助图像的转换矩阵T:
其中,(xp,yp)、(yp-sla,Rp-sla)分别表示主图像和辅助图像中特征点的二维坐标。转换矩阵T包含六个未知参数。因此,需要选择至少三组匹配度高的特征点对来计算转换矩阵T。
步骤3,利用所述转换矩阵T,得到真实三维目标在两个不同视角(主图像和辅助图像)的子孔径图像中的二维坐标对。
步骤4,根据两幅子孔径SAR成像系统的几何特点,建立圆柱对称模型(CSM),利用步骤3中的二维坐标对,计算目标的真实三维坐标。
具体的,步骤4如下:
1、辅助图像圆柱对称模型的建立
子孔径k1系统的辅助成像几何如图2(a)所示,其中,ab表示子孔径基线。散射体D、P、C位于三维观测场景中的同一方位角单元,其中C为位于地平面和散射体的参考散射体,D和P位于相同的斜距单元。根据距离徙动算法的几何特性,二维成像平面是由两个正交向量L和Rc构成的斜距平面,其中L为雷达航迹向量,Rc为参考散射体最小向量。因此,瞬时斜视角可以通过D、P和C表示为:
其中,R0表示散射体D、P的最小倾斜范围,v为平台速度,y是三个散射体的方位角坐标,RC(η)、RD(η)、RP(η)分别表示C、D、P点的瞬时倾斜范围,|RC(η)|、|RD(η)|和|RP(η)|随着慢时间η的变化而变化。
因此,从散射体D和散射体P接收到的回波是相同的,这意味着这两个不同的散射体将被聚焦在斜距平面的相同像素点上,这种现象被称为圆柱对称,这也是单幅二维SAR图像不能用于三维坐标提取的原因。而圆柱对称是SAR成像系统特有的特性,可以用来表示观测场景中真实三维坐标与二维坐标之间的几何映射关系。辅助图像中三维观测场景与二维观测场景之间的CSM可以表示为:
其中,H表示飞机飞行高度,x0、y0、z0表示三维观测场景的三维坐标点,Rp-sia表示二维观测场景距离向长度,yp-sia表示二维观测场景的方位向长度。
2、主图像圆柱对称模型的建立
主图像子孔径k2系统的成像几何如图2(b)所示。AB表示主图像子孔径基线,Pr表示雷达平台。散射体Pt=(x0,y0,z0)是三维观测空间中的任意散射体,根据反向投影算法的成像几何,主图像被聚焦于地平面。因此,基于圆柱对称特性,Pt将被聚焦在地平面中Pt'=(xp,yp,0)的位置,由图2(b)可知,从雷达平台到散射体的瞬时倾斜距离矢量可表示为:
|PrPt(η)|=|DPt-DPr(η)|
|PrPt'(η)|=|DPt'-DPr(η)| (11)
其中,PrPt(η)为散射体Pt到飞机坐标Pr的距离矢量,DPt为散射体Pt到主图像子孔径基线AB的最小斜距向量,DPr(η)为载机飞行轨迹的瞬时矢量;PrPt'(η)为散射点Pt'到飞机坐标Pr的距离矢量,DPt'为散射点Pt'到主图像子孔径基线AB的最小斜距向量;Pt'为Pt的投影,随着慢时间η的变化,这两个散射体的瞬时斜距矢量必须与保持一致。因此,最小斜距向量必须相等:
|DPt|=|DPt′| (12)
作为最小斜距向量,DPt'和DPt均与主图像子孔径基线AB正交。因此,三角形PtAB的面积与三角形Pt'AB的面积相等:
根据矢量原理,三角形的面积也可以表示为三角形两边的叉积:
其中,APt、APt'表示雷达平台初始位置到散射体Pt和Pt'的矢量;“×”算子是两个向量的叉乘运算。因此,根据(13)和(14),由Pt与Pt'的几何关系及其投影可以表示为:
|APt×AB|=|APt'×AB| (15)
借助雷达平台上的高精度INS和GPS系统,式(15)中的向量可以表示为坐标形式:
其中,i、j、k为三个互相垂直的单位向量,x0、y0、z0为Pt的三维点坐标,xA、yA、zA为A的三维点坐标,xB、yB、zB为B的三维点坐标,xP、yP、zP为Pt'的三维点坐标。式(16)包含需要提取的三维坐标(x0,y0,z0);式(17)是常量,可以从主图像计算获得。根据(16)、(17)中的表达式,可以将式(15)展开成三维坐标表示形式:
其中,C1、C2表示常数,其他系数为:
α1=(zB-zA)2+(yB-yA)2
α2=(xB-xA)2+(zB-zA)2
α3=(xB-xA)2+(yB-yA)2 (19)
γ1=-2(yB-yA)(xB-xA)
γ2=-2(zB-zA)(yB-yA)
γ3=-2(zB-zA)(xB-xA) (21)
因此,结合式(10)和(18)可以建立如图2(b)所示的曲线SAR系统的CSM:
通过式(22)中的三个方程,可以提取出目标的真实三维坐标(x0,y0,z0)。
仿真实验
通过对以下实测数据的仿真实验来进一步验证本发明方法的有效性。
仿真实验利用距离徙动算法对斜距平面进行成像处理,获得SAR辅助图像I1如图3(a)所示,利用BPA算法对地距平面进行成像处理,获得SAR主图像I2如图3(b)所示。
采用本发明方法对图3中的两幅图像特征点提取,提取出两个不同的二维坐标。通过两个不同角度的二维点坐标可以通过本发明的方法得到三维点坐标。本发明方法更实用,它适用于各种各样的SAR平台与弯曲的路径,该方法适用于各种曲线路径的SAR平台,且模型求解简单,重构精度高,具有较强的实用性。
利用斜距平面和地距平面在同一点的二维坐标,再通过步骤4的几何关系可以获得最后的三维坐标如图4所示。
虽然,本说明书中已经用一般性说明及具体实施方案对本发明作了详尽的描述,但在本发明基础上,可以对之作一些修改或改进,这对本领域技术人员而言是显而易见的。因此,在不偏离本发明精神的基础上所做的这些修改或改进,均属于本发明要求保护的范围。
Claims (10)
1.基于圆柱对称模型的单通道SAR目标三维坐标提取方法,其特征在于,包括以下步骤:
步骤1,获取全孔径雷达回波数据,根据雷达飞行轨迹将全孔径分为两个不同视角的子孔径k1、k2;
利用距离徙动算法对子孔径k1的回波信号在斜距平面进行成像处理,得到二维SAR辅助图像I1;利用反投影算法对子孔径k2的回波信号在地距平面进行成像处理,得到二维SAR主图像I2;
步骤2,利用快速稳健特征算法获得所述二维SAR辅助图像I1与所述二维SAR主图像I2之间的转换矩阵T;
步骤3,利用转换矩阵T,得到真实三维目标在两个子孔径图像中的二维坐标对;
步骤4,建立辅助图像圆柱对称模型和主图像圆柱对称模型,根据所述真实三维目标在两个子孔径图像中的二维坐标对,计算目标的真实三维坐标。
2.根据权利要求1所述的基于圆柱对称模型的单通道SAR目标三维坐标提取方法,其特征在于,步骤1中,所述距离徙动算法包含以下子步骤:
子步骤a),对子孔径k1的回波信号进行二维傅里叶变换,得到子孔径k1距离频域方位频域回波信号;
子步骤b),利用一致补偿函数,对所述子孔径k1距离频域方位频域回波信号进行一致补偿,得到子孔径k1一致粗补偿后二维频域回波信号;
子步骤c),对所述子孔径k1一致粗补偿后二维频域回波信号进行stold插值运算,得到子孔径k1位于斜距平面的二维SAR图像;
子步骤d),对所述子孔径k1位于斜距平面的二维SAR图像进行二维傅里叶逆变换,得到二维SAR辅助图像I1。
3.根据权利要求1所述的基于圆柱对称模型的单通道SAR目标三维坐标提取方法,其特征在于,步骤1中,所述反投影算法包含以下子步骤:
子步骤a),子孔径k2的回波信号进行距离向傅里叶变换,得到子孔径k2的距离频域方位时域回波信号;
子步骤b),利用距离脉压函数对所述子孔径k2的距离频域方位时域回波信号进行脉冲压缩处理,得到子孔径k2脉冲压缩后的回波信号;
子步骤c),对所述子孔径k2脉冲压缩后的回波信号进行距离向傅里叶逆变换,得到子孔径k2的二维时域回波信号;
子步骤d),将地距平面均匀的划分成网格,确定每一个网格到雷达平台的距离信息,根据该距离信息对所述子孔径k2的二维时域回波信号进行反投影处理,得到二维SAR主图像I2。
4.根据权利要求2或3所述的基于圆柱对称模型的单通道SAR目标三维坐标提取方法,其特征在于,所述傅里叶变换是利用以下公式完成的:
s1(f,t)=∫s0(τ,t)exp(-j2πfτ)dτ
其中,s1(f,t)表示距离频域为f、方位时间为t的距离频域回波信号,f的取值范围为fs表示SAR图像成像系统的采样频率,t的取值范围为(0,T),T表示SAR图像成像系统的运行时间,∫(·)dτ表示对距离时间τ进行积分操作,s0(τ,t)表示距离时域为τ、方位时间为t的子孔径回波信号,τ的取值范围为/>ro表示三维SAR图像成像区域中目标中心点的斜距,c表示光速,k表示SAR图像成像系统的采样点数,B表示SAR图像成像系统的带宽,exp(·)表示以自然常数e为底的指数操作,j表示虚数单位符号,π表示圆周率。
5.根据权利要求4所述的基于圆柱对称模型的单通道SAR目标三维坐标提取方法,其特征在于,所述距离脉压函数如下:
其中,s2(K,t)表示距离频域为K、方位时间为t的脉压后距离频域回波信号,K的取值范围为γ表示SAR图像成像系统发射信号的距离调频率。
6.根据权利要求5所述的基于圆柱对称模型的单通道SAR目标三维坐标提取方法,其特征在于,所述傅里叶逆变换是利用下式完成的:
s3(τ1,t)=∫s2(K,t)exp(j2πKτ)dK
其中,s3(τ1,t)表示距离时域为τ1、方位时间为t的子孔径脉冲压缩后距离时域信号,τ1的取值范围为∫(·)dK表示对距离频域K进行积分操作。
7.根据权利要求2所述的基于圆柱对称模型的单通道SAR目标三维坐标提取方法,其特征在于,所述一致补偿函数的相位为:
其中,f0为中心频率,fτ为距离向频率,fη为方位向频率,c为光速,v为平台速度,θ为斜视角,R0为中心斜距,Rref为参考斜距,θRFM为一致补偿相位;
所述stold插值运算是利用下式完成的:
其中,fτ'为stold插值后的新的域。
8.根据权利要求1所述的基于圆柱对称模型的单通道SAR目标三维坐标提取方法,其特征在于,步骤2中,所述快速稳健特征算法为:对主图像I2和辅助图像I1进行特征点提取、描述、匹配,通过两幅图像的特征点匹配对求出转换矩阵T,具体包含以下子步骤:
子步骤2.1,特征点提取:
(1)将待处理的图像变换成积分图像;
(2)建立基于箱式滤波器的特征点检测器作为卷积核,箱式滤波器矩阵表达式定义如下:
其中,Hopprox(x,σ)为利用箱式滤波器得到的过滤点,x为积分图像的坐标,σ为常数,Dxx(x,σ)、Dxy(x,σ)、Dyy(x,σ)分别为xx方向、yy方向、xy方向的偏导;
(3)改变卷积核的尺度大小,将不同尺度的卷积核与积分图像进行卷积,得到金字塔型的尺度空间,采用邻域非极大值抑制法提取所述金字塔型的尺度空间中的特征点;
子步骤2.2,特征描述符建立:
计算特征点周围的圆形邻域内所有点的haar小波响应,获取该特征点的主向量方向;
沿主向量方向构建正方形区域,将所述正方形区域均匀拆分为4个方形子区域,每个方形子区域都有其特定的四维描述符v0:
v0=(∑dx,∑dy,∑|dx|,∑|dx|)
其中,dx表示水平方向的Haar波响应;dy表示朝垂直方向的哈波响应;在所有4×4方形子区域中,描述符向量由64个元素组成;
子步骤2.3,特征点匹配并计算转换矩阵T:
描述符确认后,特征点的相似性可以通过Eculidean距离表示其中i表示描述符向量,xi表示主图像描述符,x'i表示辅助图像描述符,Eculidean距离最小的两个点形成一个特征点匹配对;计算从主图像到辅助图像的转换矩阵T:
其中,(xp,yp)、(yp-sla,Rp-sla)分别表示主图像和辅助图像中特征点的二维坐标;转换矩阵T包含六个未知参数,需要选择至少三组匹配度高的特征点对来计算转换矩阵T。
9.根据权利要求8所述的基于圆柱对称模型的单通道SAR目标三维坐标提取方法,其特征在于,所述领域非极大值抑制法为:将箱式滤波器得到的过滤点与金字塔型的尺度空间的其他点进行极大值对比,若过滤点的极大值最大,则该过滤点为提取尺度空间中的特征点。
10.根据权利要求1所述的基于圆柱对称模型的单通道SAR目标三维坐标提取方法,其特征在于,步骤4具体包含以下子步骤:
子步骤4.1,辅助图像圆柱对称模型的建立
散射体D、P、C位于三维观测场景中的同一方位角单元,其中C为位于地平面和散射体的参考散射体,D和P位于相同的斜距单元,则C、D、P点的瞬时倾斜范围RC(η)、RD(η)、RP(η)分别如下:
其中,Rc为参考散射体最小向量,v为平台速度,η为慢时间,y是三个散射体的方位角坐标,R0表示散射体D、P的最小倾斜范围;
由上式可知,从散射体D和散射体P接收到的回波是相同的,因此,辅助图像的圆柱对称模型可以表示为:
其中,H表示飞机飞行高度,x0、y0、z0表示三维观测场景的三维坐标点,Rp-sia表示二维观测场景距离向长度,yp-sia表示二维观测场景的方位向长度;
子步骤4.2,主图像圆柱对称模型的建立
散射体Pt=(x0,y0,z0)是三维观测空间中的任意散射体,根据反向投影算法的成像几何,主图像被聚焦于地平面,因此,基于圆柱对称特性,Pt手被聚焦在地平面中Pt′=(xp,yp,0)的位置,则从雷达平台到散射体的瞬时倾斜距离矢量可表示为:
|PrPt(η)|=|DPt-DPr(η)|
|PrPt′(η)|=|DPt′-DPr(η)|
其中,PrPt(η)为散射体Pt到飞机坐标Pr的距离矢量,DPt为散射体Pt到主图像子孔径基线AB的最小斜距向量,DPr(η)为载机飞行轨迹的瞬时矢量,PrPt′(η)为散射点Pt′到飞机坐标Pr的距离矢量,DPt′为散射点Pt′到主图像子孔径基线AB的最小斜距向量;Pt′为Pt的投影,随着慢时间η的变化,散射体Pt′与Pt的最小斜距向量必须相等,即:
|DPt|=|DPt′|
作为最小斜距向量,DPt′和DPt均与主图像子孔径基线AB正交,因此,三角形PtAB的面积与三角形Pt′AB的面积相等:
根据矢量原理,三角形的面积也可以表示为三角形两边的叉积:
其中,APt、APt′表示雷达平台初始位置到散射体Pt和Pt′的矢量;“×”算子是两个向量的叉乘运算;因此,根据三角形PtAB的面积与三角形Pt′AB的面积相等,则:
|APt×AB|=|APt′×AB|
借助雷达平台上的高精度INS和GPS系统,上式中的向量可以表示为坐标形式:
其中,i、j、k为三个互相垂直的单位向量,x0、y0、z0为Pt的三维点坐标,xA、yA、zA为A的三维点坐标,xB、yB、zB为B的三维点坐标,xP、yP、zP为Pt′的三维点坐标;将|APt×AB|=|APt′×AB|展开成三维坐标表示形式,即得主图像圆柱对称模型为:
其中,C1、C2表示常数,其他系数为:
α1=(zB-zA)2+(yB-yA)2
α2=(xB-xA)2+(zB-zA)2
α3=(xB-xA)2+(yB-yA)2
γ1=-2(yB-yA)(xB-xA)
γ2=-2(zB-zA)(yB-yA)
γ3=-2(zB-zA)(xB-xA)
子步骤4.3,结合所述辅助图像的圆柱对称模型和所述主图像圆柱对称模型,可以提取出目标的真实三维坐标(x0,y0,z0)。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110763425.8A CN113671497B (zh) | 2021-07-06 | 2021-07-06 | 基于圆柱对称模型的单通道sar目标三维坐标提取方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110763425.8A CN113671497B (zh) | 2021-07-06 | 2021-07-06 | 基于圆柱对称模型的单通道sar目标三维坐标提取方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN113671497A CN113671497A (zh) | 2021-11-19 |
CN113671497B true CN113671497B (zh) | 2024-04-23 |
Family
ID=78538637
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110763425.8A Active CN113671497B (zh) | 2021-07-06 | 2021-07-06 | 基于圆柱对称模型的单通道sar目标三维坐标提取方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN113671497B (zh) |
Families Citing this family (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116087954B (zh) * | 2023-04-06 | 2023-06-30 | 中国科学院空天信息创新研究院 | 一种实孔径系统近距离偏焦合成孔径成像方法 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
EP2743727A2 (en) * | 2014-01-16 | 2014-06-18 | Institute of Electronics, Chinese Academy of Sciences | Method for implementing high-resolution wide-swath spaceborne SAR system |
CN109959932A (zh) * | 2019-04-08 | 2019-07-02 | 西安电子科技大学 | 基于下降段曲线轨迹的雷达前视三维成像方法 |
CN110673143A (zh) * | 2019-09-30 | 2020-01-10 | 西安电子科技大学 | 一种子孔径大斜视sar俯冲成像的两步处理方法 |
CN111208513A (zh) * | 2020-01-15 | 2020-05-29 | 西安电子科技大学 | 空间目标isar图像序列能量反向投影与三维重构方法 |
CN111273290A (zh) * | 2020-01-22 | 2020-06-12 | 西安电子科技大学 | 基于预成像曲线轨迹的三维sar成像方法 |
-
2021
- 2021-07-06 CN CN202110763425.8A patent/CN113671497B/zh active Active
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
EP2743727A2 (en) * | 2014-01-16 | 2014-06-18 | Institute of Electronics, Chinese Academy of Sciences | Method for implementing high-resolution wide-swath spaceborne SAR system |
CN109959932A (zh) * | 2019-04-08 | 2019-07-02 | 西安电子科技大学 | 基于下降段曲线轨迹的雷达前视三维成像方法 |
CN110673143A (zh) * | 2019-09-30 | 2020-01-10 | 西安电子科技大学 | 一种子孔径大斜视sar俯冲成像的两步处理方法 |
CN111208513A (zh) * | 2020-01-15 | 2020-05-29 | 西安电子科技大学 | 空间目标isar图像序列能量反向投影与三维重构方法 |
CN111273290A (zh) * | 2020-01-22 | 2020-06-12 | 西安电子科技大学 | 基于预成像曲线轨迹的三维sar成像方法 |
Non-Patent Citations (1)
Title |
---|
弧形阵列3D-SAR成像处理算法研究;黄增树;孙进平;谭维贤;黄平平;乞耀龙;;信号处理;20190425(第04期);全文 * |
Also Published As
Publication number | Publication date |
---|---|
CN113671497A (zh) | 2021-11-19 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN110488288B (zh) | 一种机载sar高分辨层析成像方法 | |
Feng et al. | An extended fast factorized back projection algorithm for missile-borne bistatic forward-looking SAR imaging | |
CN106772380A (zh) | 一种圆周合成孔径雷达成像方法 | |
CN108051809A (zh) | 基于Radon变换的运动目标成像方法、装置及电子设备 | |
CN111273290B (zh) | 基于预成像曲线轨迹的三维sar成像方法 | |
CN103454638B (zh) | 一种圆迹合成孔径雷达三维层析成像方法 | |
CN110673143A (zh) | 一种子孔径大斜视sar俯冲成像的两步处理方法 | |
CN112415515B (zh) | 一种机载圆迹sar对不同高度目标分离的方法 | |
CN103336278A (zh) | 多视角观测下前视三维sar成像方法 | |
CN104749570A (zh) | 一种移不变机载双基合成孔径雷达目标定位方法 | |
KR102151362B1 (ko) | 극좌표변환을 이용한 항공기기반 영상복원장치 및 이를 이용한 영상복원방법 | |
CN103308913A (zh) | 一种高速飞行器载双天线前视sar解模糊算法 | |
CN114545411B (zh) | 一种基于工程实现的极坐标格式多模高分辨sar成像方法 | |
CN103616688A (zh) | 一种提升三维干涉逆合成孔径雷达图像质量的方法 | |
Tang et al. | Focusing hypersonic vehicle-borne SAR data using radius/angle algorithm | |
CN108132466B (zh) | 一种机载阵列天线下视三维成像方法和系统 | |
Yang et al. | A novel three-step image formation scheme for unified focusing on spaceborne SAR data | |
Vu et al. | Phase error calculation for fast time-domain bistatic SAR algorithms | |
CN106154269A (zh) | 一种适用于圆周合成孔径雷达快速时域成像方法 | |
CN113671497B (zh) | 基于圆柱对称模型的单通道sar目标三维坐标提取方法 | |
CN115685200A (zh) | 一种高精度大前斜视sar成像运动补偿与几何校正方法 | |
CN109143236B (zh) | 适用于复杂飞行轨迹的双基聚束sar大场景成像方法 | |
CN112859018B (zh) | 一种基于图像几何校正的视频sar成像方法 | |
CN111880179A (zh) | 一种弹载弧线俯冲大斜视tops sar的成像方法 | |
CN110244267B (zh) | 弹载毫米波切向前视sar三维成像模型性能分析方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |