CN104035093A - 基于变标算法的星载sar地面加速运动目标成像方法 - Google Patents

基于变标算法的星载sar地面加速运动目标成像方法 Download PDF

Info

Publication number
CN104035093A
CN104035093A CN201410143273.1A CN201410143273A CN104035093A CN 104035093 A CN104035093 A CN 104035093A CN 201410143273 A CN201410143273 A CN 201410143273A CN 104035093 A CN104035093 A CN 104035093A
Authority
CN
China
Prior art keywords
range
exp
target
echo signal
correction
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.)
Granted
Application number
CN201410143273.1A
Other languages
English (en)
Other versions
CN104035093B (zh
Inventor
王彤
李永康
任倩倩
张颖
马南
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xidian University
Original Assignee
Xidian University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Xidian University filed Critical Xidian University
Priority to CN201410143273.1A priority Critical patent/CN104035093B/zh
Priority claimed from CN201410143273.1A external-priority patent/CN104035093B/zh
Publication of CN104035093A publication Critical patent/CN104035093A/zh
Application granted granted Critical
Publication of CN104035093B publication Critical patent/CN104035093B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Systems 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/88Radar or analogous systems specially adapted for specific applications
    • G01S13/89Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021SAR image post-processing techniques
    • G01S13/9029SAR image post-processing techniques specially adapted for moving target detection within a single SAR image or within multiple SAR images taken at the same time
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Systems 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/88Radar or analogous systems specially adapted for specific applications
    • G01S13/89Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9004SAR image acquisition techniques
    • G01S13/9011SAR image acquisition techniques with frequency domain processing of the SAR signals in azimuth
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Systems 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/02Systems using reflection of radio waves, e.g. primary radar systems; Analogous systems
    • G01S13/50Systems of measurement based on relative movement of target
    • G01S13/52Discriminating between fixed and moving objects or between objects moving at different speeds
    • G01S13/536Discriminating between fixed and moving objects or between objects moving at different speeds using transmission of continuous unmodulated waves, amplitude-, frequency-, or phase-modulated waves

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Signal Processing (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明属于雷达技术领域,公开了一种基于变标算法的星载SAR地面加速运动目标成像方法,其方法步骤为:1)对目标瞬时距离方程进行等效;2)将目标原始回波信号变到距离多普勒域;3)根据等效后的距离方程,在距离多普勒域进行补余RCMC;4)将补余RCMC后的信号变到二维频域;5)在二维频域进行距离压缩和一致RCMC;6)将一致RCMC后的信号变到距离多普勒域;7)在距离多普勒域进行相位校正和方位压缩;8)对相位校正和方位压缩后的信号进行方位向逆傅里叶变换,完成对目标的成像。本发明仅需知道等效速度ve,无需对目标参数进行高维搜索,计算效率高;利用变标实现距离徙动校正,不需要插值,便于工程实现,成像精确。

Description

基于变标算法的星载SAR地面加速运动目标成像方法
技术领域
本发明属于雷达技术领域,涉及一种基于变标算法的星载SAR地面加速运动目标成像方法,具体地说是一种基于变标算法(Chirp Scaling Algorithm,CSA)的地面加速运动目标成像方法,用于星载合成孔径雷达(Synthetic Aperture Radar,SAR)对地面加速度运动目标成像。
背景技术
星载SAR作为天基对地观测信息获取的重要手段,在地球观测、军事侦察、地质资源普查、灾情勘查、海洋观测和环境监测等方面具有重要的应用价值。自1978年世界上第一颗合成孔径雷达卫星(美国Seasat卫星)发射成功以来,星载SAR逐渐得到各国科研工作者的重视,并在近三十年内发展迅速。目前,星载SAR的研究和应用已成为一个国家综合国力和科技水平的标志之一。
与机载雷达系统相比,星载SAR系统具有运行高度高、覆盖面积广、安全性高等特点;与卫星光学成像系统相比,星载合成孔径雷达观测可以全天候、全天时地探测隐蔽目标,提供丰富的地表电磁散射特性信息。这些优点使得星载雷达被广泛应用于国民经济和国防领域。地面运动目标成像是星载合成孔径雷达的一个重要应用领域,它在战场侦察和交通监控中有重要作用,受到了广泛的关注。
SAR成像需要对目标进行精确的距离徙动校正和相位弯曲补偿。地面运动目标由于存在速度,其距离单元徙动和相位弯曲与静止目标不同。因此,如果直接将针对静止目标的成像方法用于对运动目标进行成像,由于距离徙动校正和相位弯曲补偿不准确,目标会出现散焦现象。为解决这个问题,各国学者提出了很多方法。但这些方法大都存在如下两个问题:1)需要知道目标位置参数和运动参数,或者需要对目标参数进行高维搜索,或者需要插值,计算量非常大,适用范围有限;2)提出的方法主要针对匀速运动的目标,而对于地面加速运动目标,没能提出一个精确的成像方法。
发明内容
本发明的主要目的在于提出一种基于变标算法的星载SAR地面加速运动目标成像方法,该方法无需对目标参数进行高维搜索,其成像过程中也无需插值,计算效率高,适用范围广。
本发明的基本思路是:1)对目标瞬时距离方程进行二次等效;2)将目标原始回波信号变到距离多普勒域;3)根据等效后的距离方程,在距离多普勒域进行补余RCMC;4)将补余RCMC后的信号变到二维频域;5)在二维频域进行距离压缩和一致RCMC;6)将一致RCMC后的信号变到距离多普勒域;7)在距离多普勒域进行相位校正和方位压缩;8)对相位校正和方位压缩后的信号进行方位向逆傅里叶变换,完成对目标的成像。
为了达到上述目的,本发明采用以下技术方案予以实现。
一种基于变标算法的星载SAR地面加速运动目标成像方法,其特征在于,包括以下步骤:
步骤1,对目标的瞬时距离方程进行二次等效,得到二次等效后的目标距离方程:
R ( t a ) = ( v e t a - x e ) 2 + y e 2
v e = ( v a - v x ) 2 + a y y 0 + v y 2
x e = y 0 v y ( v a - v x ) 2 + a y y 0 + v y 2
y e = ( v a - v x ) 2 + a y y 0 ( v a - v x ) 2 + a y y 0 + v y 2
其中,雷达工作在正侧视模式下,雷达平台速度为va,ta为慢时间;ta=0时,目标方位向速度、方位向加速度、距离向速度、距离向加速度分别为vx、ax、vy和ay,且此时,雷达位于坐标原点,目标位于(0,y0);R(ta)为ta时刻目标到雷达的瞬时距离,xe、ye为等效后目标的位置参数,ve为等效后雷达平台的速度。
目标的原始回波信号为:
s ( t a , t r ) = w a ( t a ) w r ( t r - 2 R ( t a ) / c ) exp { - j 4 π f c R ( t a ) c + jπ K r ( t r - 2 R ( t a ) c ) 2 }
其中,tr为快时间,c为光速,wa(ta)、wr(tr)分别为目标回波信号的方位包络和距离包络,fc为雷达发射信号的载频,Kr为雷达发射的线性调频信号的调频率;
步骤2,对地面加速度运动目标的原始回波信号进行方位向傅里叶变换,将目标的原始回波信号变换到距离多普勒域,得到距离多普勒域目标信号;
步骤3,根据二次等效后的目标距离方程,构造变标方程,对距离多普勒域目标信号进行补余距离徙动校正;
步骤4,对补余距离徙动校正后的目标信号进行距离向傅里叶变换,得到二维频域目标信号;
步骤5,在二维频域构造用于距离压缩和一致距离徙动校正的参考函数,将二维频域目标信号与构造的用于距离压缩和一致距离徙动校正的参考函数相乘,实现距离压缩和一致距离徙动校正;
步骤6,对距离压缩和一致距离徙动校正后的目标信号进行距离向逆傅里叶变换,将距离压缩和一致距离徙动校正后的目标信号变换到距离多普勒域,得到距离压缩和一致距离徙动校正后的距离多普勒域目标信号;
步骤7,在距离多普勒域,构造用于方位压缩和相位校正的参考函数,将距离压缩和一致距离徙动校正后的距离多普勒域目标信号与用于距离多普勒域方位压缩和相位校正的参考函数相乘,实现对距离压缩和一致距离徙动校正后的距离多普勒域目标信号的方位压缩和相位校正,得到方位压缩和相位校正后的目标信号;
步骤8,对方位压缩和相位校正后的目标信号进行方位向逆傅里叶变换,完成对地面加速运动目标的成像。
上述技术方案的特点和进一步改进在于:
(1)步骤2的具体子步骤为:
目标的原始回波信号表示为:
s ( t a , t r ) = w a ( t a ) w r ( t r - 2 R ( t a ) / c ) exp { - j 4 π f c R ( t a ) c + jπ K r ( t r - 2 R ( t a ) c ) 2 }
其中,tr为快时间,c为光速,wa(ta)、wr(tr)分别为目标回波信号的方位包络和距离包络,fc为雷达发射信号的载频,Kr为雷达发射的线性调频信号的调频率;
经方位向傅里叶变换后,距离多普勒域目标信号表示为:
s ( f a , t r ) = w a ( f a - f ac ) w r ( K m ( v e , f a ) K r ( t r - 2 y e cD ( f a , v e ) ) ) × exp { - jπ 4 y e D ( f a , v e ) f c c } exp { jπ K m ( v e , f a ) ( t r - 2 y e cD ( f a , v e ) ) 2 }
其中,Wa(fa)是目标信号方位频谱的包络,fa是方位频率,fac为目标多普勒中心频率, D ( f a , v e ) = 1 - c 2 f a 2 4 v e 2 f c 2 , K m ( v e , f a ) = K r + 2 v e 2 K r 2 f c 3 D 3 ( f a , v e ) cy e f a 2 .
(2)步骤3的具体子步骤为:
3a)根据二次等效后的目标距离方程和距离多普勒域回波目标信号的表达式
s ( f a , t r ) = w a ( f a - f ac ) w r ( K m ( v e , f a ) K r ( t r - 2 y e cD ( f a , v e ) ) ) × exp { - jπ 4 y e D ( f a , v e ) f c c } exp { jπ K m ( v e , f a ) ( t r - 2 y e cD ( f a , v e ) ) 2 }
目标的补余距离徙动表示为:其中Rref为场景中心距离;
距离多普勒域变标方程可构造为:
s sc ( f a , t r ) = exp { jπ K m ( v e , f a ) ( 1 D ( f a , v e ) - 1 ) ( t r - 2 R ref cD ( f a , v e ) ) 2 }
3b)将变换得到的距离多普勒域目标信号与构造的变标方程相乘,实现补余距离徙动校正:
s 1 ( f a , t r ) = s ( f a , t r ) s sc ( f a , t r ) = W a ( f a - f ac ) w r ( K m ( v e , f a ) K r ( t r - 2 y e cD ( f a , v e ) ) ) exp { - jπ 4 y e D ( f a , v e ) f c c } × exp { jπ K m ( v e , f a ) [ ( 1 D ( f a , v e ) - 1 ) ( t r - 2 R ref cD ( f a , v e ) ) 2 + ( t r - 2 y e cD ( f a , v e ) ) 2 ] }
其中,s1(fa,tr)为补余距离徙动校正后的目标信号。
(3)步骤5的具体子步骤为:
5a)根据二维频域目标信号表达式
S 1 ( f a , f r ) = W a ( f a - f ac ) W r ( f r ) exp { - j 4 πy e f c D ( f a , v e ) c } × exp { - jπ D ( f a , v e ) f r 2 K m ( v e , f a ) } exp { - j 4 πy e c f r } × exp { - j 4 π c ( 1 D ( f a , v e ) - 1 ) R ref f r } × exp { j 4 πK m ( v e , f a ) c 2 ( 1 - D ( f a , v e ) ) ( y e D ( f a , v e ) - R ref D ( f a , v e ) ) 2 }
目标距离调制和一致距离徙动的相位项为:
exp { - jπ D ( f a , v e ) f r 2 K m ( v e , f a ) - j 4 π c ( 1 D ( f a , v e ) - 1 ) R ref f r }
二维频域实现距离压缩和一致距离徙动校正的参考函数构造为:
H ( f a , f r ) = exp { jπ D ( f a , v e ) f r 2 K m ( v e , f a ) } exp { j 4 π c ( 1 D ( f a , v e ) - 1 ) R ref f r }
5b)将二维频域目标信号与距离压缩和一致距离徙动校正的参考函数相乘,即
S 2 ( f a , f r ) = S 1 ( f a , f r ) H ( f a , f r ) = W a ( f a - f ac ) W r ( f r ) exp { - j 4 πy e f c D ( f a , v e ) c } × exp { - j 4 πy e c f r } exp { j 4 π K m ( v e , f a ) c 2 ( 1 - D ( f a , v e ) ) ( y e D ( f a , v e ) - R ref D ( f a , v e ) ) 2 }
从而实现距离压缩和一致距离徙动校正。
(4)步骤7的具体子步骤为:
7a)根据距离压缩和一致距离徙动校正后的距离多普勒域目标信号的表达式
S 2 ( f a , t r ) = W a ( f a - f ac ) p r ( t r - 2 y e c ) exp { - j 4 πy e f c D ( f a , v e ) c } × exp { j 4 πK m ( v e , f a ) c 2 ( 1 - D ( f a , v e ) ) ( y e D ( f a , v e ) - R ref D ( f a , v e ) ) 2 }
目标方位调制和相位误差的相位项为:
exp { - j 4 π y e f c D ( f a , v e ) c + j 4 π K m ( v e , f a ) c 2 ( 1 - D ( f a , v e ) ) ( y e D ( f a , v e ) - R ref D ( f a , v e ) ) 2 } .
用于距离多普勒域方位压缩和相位校正的参考函数构造为:
H a ( f a ) = exp { j 4 πy e f c D ( f a , v e ) c } exp { - j 4 πK m K m ( v e , f a ) c 2 ( 1 - D ( f a , v e ) ) ( y e D ( f a , v e ) - R ref D ( f a , v e ) ) 2 }
7b)将距离压缩和一致距离徙动校正后的距离多普勒域目标信号与构造的用于距离多普勒域方位压缩和相位校正的参考函数相乘,即
S 3 ( f a , t r ) = S 2 ( f a , t r ) H a ( f a ) = W a ( f a - f ac ) p r ( t r - 2 y e c )
从而实现方位压缩和相位校正。
本发明与现有技术相比具有以下优点:a)本发明仅需知道等效速度ve,无需对目标参数进行高维搜索,计算效率高,适用范围广;b)本发明利用变标实现距离徙动校正,不需要插值,便于工程实现;c)本发明能实现对地面加速运动目标的精确成像。
附图说明
下面结合附图和具体实施方式对本发明做进一步详细说明。
图1是本发明的实现流程示意图;
图2斜距平面星载SAR系统观测几何图;其中横坐标表示方位向,纵坐标表示距离向;
图3是距离徙动校正前的目标轨迹图;其中横坐标表示距离单元,纵坐标表示多普勒单元;
图4是用本发明方法的距离徙动校正后的目标轨迹图;其中横坐标表示距离单元,纵坐标表示多普勒单元;
图5(a)是用本发明方法的方位压缩后的目标成像结果图;其中横坐标表示距离单元,纵坐标表示方位单元;
图5(b)是放大后的等高线图;其中横坐标表示距离向,纵坐标表示方位向。
具体实施方式
参照图1,进一步详细说明本发明的基于变标算法的星载SAR地面加速运动目标成像方法,其具体实现步骤如下:
步骤1,对目标的瞬时距离方程进行二次等效,得到二次等效后的目标距离方程:
R ( t a ) = ( v e t a - x e ) 2 + y e 2
v e = ( v a - v x ) 2 + a y y 0 + v y 2
x e = y 0 v y ( v a - v x ) 2 + a y y 0 + v y 2
y e = ( v a - v x ) 2 + a y y 0 ( v a - v x ) 2 + a y y 0 + v y 2
R(ta)为ta时刻目标到雷达的瞬时距离R(ta)。
1a)如图2所示,为斜距平面星载SAR系统观测几何图,其中,雷达工作在正侧视模式下,雷达平台速度为va,ta为慢时间。ta=0时,目标方位向速度、方位向加速度、距离向速度、距离向加速度分别为vx、ax、vy和ay,且此时,雷达位于坐标原点,目标位于(0,y0)。
因此,ta时刻目标到雷达的瞬时距离R(ta)方程可表示为:
R ( t a ) = ( y 0 + v y t a + 0.5 a y t a 2 ) 2 + ( v x t a + 0.5 a x t a 2 - v a t a ) 2
1b)对目标的瞬时距离方程两边进行平方,并进行数学上的等效变形可得:
R 2 ( t a ) = ( y 0 + v y t a + 0.5 a y t a 2 ) 2 + ( v x t a + 0.5 a x t a 2 - v a t a ) 2 = y 0 2 + 2 v y y 0 t a + [ ( v a - v x ) 2 + a y y 0 + v y 2 ] t a 2 + [ a y v y - a x ( v a - v x ) ] t a 3 + [ a x 2 4 + a y 2 4 ] t a 4
对于星载SAR系统,上式中的三次项和四次项与波长相比非常小,可以忽略。举个例子,考虑参数见表1的星载SAR系统,假设目标参数为:y0=1000km,vx=5m/s,vy=5m/s,ax=1m/s2,ay=1m/s2,经计算,忽略三次项和四次项带来的误差为0.00015m,远小于波长(0.0555m),可忽略。
因此,目标距离方程可写为:
R ( t a ) = y 0 2 + 2 v y y 0 t a + [ ( v a - v x ) 2 + a y y 0 + v y 2 ] t a 2
v e = ( v a - v x ) 2 + a y y 0 + v y 2
x e = - y 0 v y ( v a - v x ) 2 + a y y 0 + v y 2
y e = ( v a - v x ) 2 + a y y 0 ( v a - v x ) 2 + a y y 0 + v y 2
即目标距离方程可等效为:
R ( t a ) = ( v e t a - x e ) 2 + y e 2
1c)目标的原始回波信号可表示为:
s ( t a , t r ) = w a ( t a ) w r ( t r - 2 R ( t a ) / c ) exp { - j 4 π f c R ( t a ) c + jπ K r ( t r - 2 R ( t a ) c ) 2 }
其中,tr为快时间,c为光速,wa(ta)、wr(tr)分别为目标回波信号的方位包络和距离包络,fc为雷达发射信号的载频,Kr为雷达发射的线性调频信号的调频率。
步骤2,对地面加速度运动目标的原始回波信号进行方位向傅里叶变换,将目标的原始回波信号变换到距离多普勒域,得到距离多普勒域目标信号。
目标的原始回波信号可表示为:
s ( t a , t r ) = w a ( t a ) w r ( t r - 2 R ( t a ) / c ) exp { - j 4 π f c R ( t a ) c + jπ K r ( t r - 2 R ( t a ) c ) 2 }
其中,tr为快时间,c为光速,wa(ta)、wr(tr)分别为目标回波信号的方位包络和距离包络,fc为雷达发射信号的载频,Kr为雷达发射的线性调频信号的调频率。
经方位向傅里叶变换后,距离多普勒域目标信号可表示为:
s ( f a , t r ) = w a ( f a - f ac ) w r ( K m ( v e , f a ) K r ( t r - 2 y e cD ( f a , v e ) ) ) × exp { - jπ 4 y e D ( f a , v e ) f c c } exp { jπ K m ( v e , f a ) ( t r - 2 y e cD ( f a , v e ) ) 2 }
其中,Wa(fa)是目标信号方位频谱的包络,fa是方位频率,fac为目标多普勒中心频率, D ( f a , v e ) = 1 - c 2 f a 2 4 v e 2 f c 2 , K m ( v e , f a ) = K r + 2 v e 2 K r 2 f c 3 D 3 ( f a , v e ) cy e f a 2 .
步骤3,根据二次等效后的目标距离方程,构造变标方程,对距离多普勒域目标信号进行补余距离徙动校正。
3a)根据二次等效后的目标距离的等效方程和距离多普勒域目标信号的表达式,目标的补余距离徙动可表示为:
因此,距离多普勒域变标方程可构造为:
s sc ( f a , t r ) = exp { jπ K m ( v e , f a ) ( 1 D ( f a , v e ) - 1 ) ( t r - 2 R ref cD ( f a , v e ) ) 2 }
3b)将变换得到距离多普勒域目标信号与构造的变标方程相乘,实现补余距离徙动校正:
s 1 ( f a , t r ) = s ( f a , t r ) s sc ( f a , t r ) = W a ( f a - f ac ) w r ( K m ( v e , f a ) K r ( t r - 2 y e cD ( f a , v e ) ) ) exp { - jπ 4 y e D ( f a , v e ) f c c } × exp { jπ K m ( v e , f a ) [ ( 1 D ( f a , v e ) - 1 ) ( t r - 2 R ref cD ( f a , v e ) ) 2 + ( t r - 2 y e cD ( f a , v e ) ) 2 ] }
其中,s1(fa,tr)为补余距离徙动校正后的目标信号。
需要注意的是,在进行补余距离徙动校正的过程中,仅用了等效速度,而没有使用目标的位置参数、速度参数和加速度。此外,也没有使用插值。
步骤4,对补余距离徙动校正后的目标信号进行距离向傅里叶变换,得到二维频域目标信号。
补余距离徙动校正后,二维频域目标信号表达式为:
S 1 ( f a , f r ) = W a ( f a - f ac ) W r ( f r ) exp { - j 4 πy e f c D ( f a , v e ) c } × exp { - jπ D ( f a , v e ) f r 2 K m ( v e , f a ) } exp { - j 4 πy e c f r } × exp { - j 4 π c ( 1 D ( f a , v e ) - 1 ) R ref f r } × exp { j 4 πK m ( v e , f a ) c 2 ( 1 - D ( f a , v e ) ) ( y e D ( f a , v e ) - R ref D ( f a , v e ) ) 2 }
其中,Wr(fr)是目标信号距离频谱的包络。
步骤5,在二维频域构造用于距离压缩和一致距离徙动校正的参考函数,将二维频域目标信号与构造的用于距离压缩和一致距离徙动校正的参考函数相乘,实现距离压缩和一致距离徙动校正。
5a)由二维频域目标信号表达式可以看出,决定目标距离调制和一致距离徙动校正的相位项为: exp { - jπ D ( f a , v e ) f r 2 K m ( v e , f a ) - j 4 π c ( 1 D ( f a , v e ) - 1 ) R ref f r } . 因此,在二维频域实现距离压缩和一致距离徙动的参考函数可构造为:
H ( f a , f r ) = exp { jπ D ( f a , v e ) f r 2 K m ( v e , f a ) } exp { j 4 π c ( 1 D ( f a , v e ) - 1 ) R ref f r }
5b)将二维频域目标信号与距离压缩和一致距离徙动的参考函数相乘,实现距离压缩和一致距离徙动:
S 2 ( f a , f r ) = S 1 ( f a , f r ) H ( f a , f r ) = W a ( f a - f ac ) W r ( f r ) exp { - j 4 πy e f c D ( f a , v e ) c } × exp { - j 4 πy e c f r } exp { j 4 π K m ( v e , f a ) c 2 ( 1 - D ( f a , v e ) ) ( y e D ( f a , v e ) - R ref D ( f a , v e ) ) 2 }
步骤6,对距离压缩和一致距离徙动校正后的目标信号进行距离向逆傅里叶变换,将距离压缩和一致距离徙动校正后的目标信号变换到距离多普勒域,得到距离压缩和一致距离徙动校正后的距离多普勒域目标信号。
距离压缩和一致距离徙动校正后的距离多普勒域目标信号可表示为:
S 2 ( f a , t r ) = W a ( f a - f ac ) p r ( t r - 2 y e c ) exp { - j 4 πy e f c D ( f a , v e ) c } × exp { j 4 πK m ( v e , f a ) c 2 ( 1 - D ( f a , v e ) ) ( y e D ( f a , v e ) - R ref D ( f a , v e ) ) 2 }
其中,pr(tr)为距离冲激响应函数。
步骤7,在距离多普勒域,构造用于方位压缩和相位校正的参考函数,将距离压缩和一致距离徙动校正后的距离多普勒域目标信号与用于距离多普勒域方位压缩和相位校正的参考函数相乘,实现对距离压缩和一致距离徙动校正后的距离多普勒域目标信号进行方位压缩和相位校正,得到方位压缩和相位校正后的目标信号。
7a)由距离压缩和一致距离徙动校正后的距离多普勒域目标信号的表达式可以看出,决定目标方位调制和相位误差的相位项为:
exp { - j 4 π y e f c D ( f a , v e ) c + j 4 π K m ( v e , f a ) c 2 ( 1 - D ( f a , v e ) ) ( y e D ( f a , v e ) - R ref D ( f a , v e ) ) 2 } .
因此,用于距离多普勒域方位压缩和相位校正的参考函数可构造为:
H a ( f a ) = exp { j 4 πy e f c D ( f a , v e ) c } exp { - j 4 πK m K m ( v e , f a ) c 2 ( 1 - D ( f a , v e ) ) ( y e D ( f a , v e ) - R ref D ( f a , v e ) ) 2 }
7b)将变换得到的距离压缩和一致RCMC后的距离多普勒域目标信号与构造的用于距离多普勒域方位压缩和相位校正的参考函数相乘,实现方位压缩和相位校正:
S 3 ( f a , t r ) = S 2 ( f a , t r ) H a ( f a ) = W a ( f a - f ac ) p r ( t r - 2 y e c )
步骤8,对方位压缩和相位校正后的目标信号进行方位向逆傅里叶变换,完成对地面加速运动目标的成像。
对S3(fa,tr)进行方位向逆傅里叶变换,可得聚焦后的SAR图像域信号的表达式:
s ( t a , t r ) = p a ( t a - a 1 2 a 2 t a ) p r ( t r - 2 a 0 c ) × exp { j 2 π f ac t a }
其中,pa(ta)为方位冲激响应函数。
从上面的理论分析中可以看出,在整个成像过程中仅用了等效速度,而没有使用目标的位置参数、速度参数和加速度。此外,本发明仅通过复乘和快速傅里叶变换实现成像,而没有使用插值,计算效率高。
下面通过仿真实验对本发明的效果做进一步详细说明。
(1)距离压缩后的目标轨迹仿真。
星载SAR系统参数见表1,目标参数为:y0=1000km,vx=5m/s,vy=5m/s,ax=1m/s2,ay=1m/s2。在二维频域通过相位相乘实现距离压缩,然后进行距离向逆傅里叶变换获得距离压缩后的目标轨迹。
仿真结果见图3,由图3可以看到目标轨迹存在明显的距离徙动。
表1星载SAR系统参数
载频 5.4GHz 卫星速度 7545m/s
距离带宽 50MHz 场景中心距离 1000km
距离采样频率 80MHz 方位带宽 1244Hz
脉冲重复频率 3000Hz 脉宽 20μs
(2)用本发明进行距离徙动校正后的目标轨迹仿真。
本仿真中的参数设置与仿真(1)中的设置相同,仿真时仅用了等效速度,而没有使用目标的位置参数、速度参数和加速度,仿真时也没有使用插值,仿真结果见图4。由图4可以看出,目标的轨迹变成了一条直线,其距离徙动被很好地校正了。这个仿真实验说明本发明可以在仅知道等效速度的前提下实现对地面加速度运动目标的距离徙动校正。
(3)本发明方法的地面运动目标成像结果仿真。
本仿真中的参数设置与仿真(1)中的设置相同,仿真时仅用了等效速度,而没有使用目标的位置参数、速度参数和加速度,仿真结果见图5。图5(a)给出了方位压缩后的结果,图5(b)给出了放大后的等高线图。由图5(a)可以看出,目标被很好的聚焦了,由图5(b)可以看出,成像质量非常高。这个仿真实验证明了本发明可以在仅知道等效速度的前提下实现对地面加速运动目标的精确成像。

Claims (5)

1.一种基于变标算法的星载SAR地面加速运动目标成像方法,其特征在于,包括以下步骤:
步骤1,对目标的瞬时距离方程进行二次等效,得到二次等效后的目标距离方程:
R ( t a ) = ( v e t a - x e ) 2 + y e 2
v e = ( v a - v x ) 2 + a y y 0 + v y 2
x e = y 0 v y ( v a - v x ) 2 + a y y 0 + v y 2
y e = ( v a - v x ) 2 + a y y 0 ( v a - v x ) 2 + a y y 0 + v y 2
其中,雷达工作在正侧视模式下,雷达平台速度为va,ta为慢时间;ta=0时,目标方位向速度、方位向加速度、距离向速度、距离向加速度分别为vx、ax、vy和ay,且此时,雷达位于坐标原点,目标位于(0,y0);R(ta)为ta时刻目标到雷达的瞬时距离,xe、ye为等效后目标的位置参数,ve为等效后雷达平台的速度;
目标的原始回波信号为:
s ( t a , t r ) = w a ( t a ) w r ( t r - 2 R ( t a ) / c ) exp { - j 4 π f c R ( t a ) c + jπ K r ( t r - 2 R ( t a ) c ) 2 }
其中,tr为快时间,c为光速,wa(ta)、wr(tr)分别为目标回波信号的方位包络和距离包络,fc为雷达发射信号的载频,Kr为雷达发射的线性调频信号的调频率;
步骤2,对地面加速度运动目标的原始回波信号进行方位向傅里叶变换,将目标的原始回波信号变换到距离多普勒域,得到距离多普勒域目标信号;
步骤3,根据二次等效后的目标距离方程,构造变标方程,对距离多普勒域目标信号进行补余距离徙动校正;
步骤4,对补余距离徙动校正后的目标信号进行距离向傅里叶变换,得到二维频域目标信号;
步骤5,在二维频域构造用于距离压缩和一致距离徙动校正的参考函数,将二维频域目标信号与构造的用于距离压缩和一致距离徙动校正的参考函数相乘,实现距离压缩和一致距离徙动校正;
步骤6,对距离压缩和一致距离徙动校正后的目标信号进行距离向逆傅里叶变换,将距离压缩和一致距离徙动校正后的目标信号变换到距离多普勒域,得到距离压缩和一致距离徙动校正后的距离多普勒域目标信号;
步骤7,在距离多普勒域,构造用于方位压缩和相位校正的参考函数,将距离压缩和一致距离徙动校正后的距离多普勒域目标信号与用于距离多普勒域方位压缩和相位校正的参考函数相乘,实现对距离压缩和一致距离徙动校正后的距离多普勒域目标信号的方位压缩和相位校正,得到方位压缩和相位校正后的目标信号;
步骤8,对方位压缩和相位校正后的目标信号进行方位向逆傅里叶变换,完成对地面加速运动目标的成像。
2.根据权利要求1所述的基于变标算法的星载SAR地面加速运动目标成像方法,其特征在于,步骤2的具体子步骤为:
目标的原始回波信号表示为:
s ( t a , t r ) = w a ( t a ) w r ( t r - 2 R ( t a ) / c ) exp { - j 4 π f c R ( t a ) c + jπ K r ( t r - 2 R ( t a ) c ) 2 }
经方位向傅里叶变换后,距离多普勒域目标信号表示为:
s ( f a , t r ) = w a ( f a - f ac ) w r ( K m ( v e , f a ) K r ( t r - 2 y e cD ( f a , v e ) ) ) × exp { - jπ 4 y e D ( f a , v e ) f c c } exp { jπ K m ( v e , f a ) ( t r - 2 y e cD ( f a , v e ) ) 2 }
其中,Wa(fa)是目标信号方位频谱的包络,fa是方位频率,fac为目标多普勒中心频率,
D ( f a , v e ) = 1 - c 2 f a 2 4 v e 2 f c 2 , K m ( v e , f a ) = K r + 2 v e 2 K r 2 f c 3 D 3 ( f a , v e ) cy e f a 2 .
3.根据权利要求1所述的基于变标算法的星载SAR地面加速运动目标成像方法,其特征在于,步骤3的具体子步骤为
3a)根据二次等效后的目标距离方程和距离多普勒域目标信号的表达式
s ( f a , t r ) = w a ( f a - f ac ) w r ( K m ( v e , f a ) K r ( t r - 2 y e cD ( f a , v e ) ) ) × exp { - jπ 4 y e D ( f a , v e ) f c c } exp { jπ K m ( v e , f a ) ( t r - 2 y e cD ( f a , v e ) ) 2 }
目标的补余距离徙动表示为:
距离多普勒域变标方程构造为:
s sc ( f a , t r ) = exp { jπ K m ( v e , f a ) ( 1 D ( f a , v e ) - 1 ) ( t r - 2 R ref cD ( f a , v e ) ) 2 }
3b)将变换得到的距离多普勒域目标信号与构造的变标方程相乘,实现补余距离徙动校正:
s 1 ( f a , t r ) = s ( f a , t r ) s sc ( f a , t r ) = W a ( f a - f ac ) w r ( K m ( v e , f a ) K r ( t r - 2 y e cD ( f a , v e ) ) ) exp { - jπ 4 y e D ( f a , v e ) f c c } × exp { jπ K m ( v e , f a ) [ ( 1 D ( f a , v e ) - 1 ) ( t r - 2 R ref cD ( f a , v e ) ) 2 + ( t r - 2 y e cD ( f a , v e ) ) 2 ] }
其中,s1(fa,tr)为补余距离徙动校正后的目标信号。
4.根据权利要求1所述的基于变标算法的星载SAR地面加速运动目标成像方法,其特征在于,步骤5的具体子步骤为:
5a)根据二维频域目标信号表达式
S 1 ( f a , f r ) = W a ( f a - f ac ) W r ( f r ) exp { - j 4 πy e f c D ( f a , v e ) c } × exp { - jπ D ( f a , v e ) f r 2 K m ( v e , f a ) } exp { - j 4 πy e c f r } × exp { - j 4 π c ( 1 D ( f a , v e ) - 1 ) R ref f r } × exp { j 4 πK m ( v e , f a ) c 2 ( 1 - D ( f a , v e ) ) ( y e D ( f a , v e ) - R ref D ( f a , v e ) ) 2 }
其中,Rref为场景中心距离。
目标距离调制和一致距离徙动的相位项为:
exp { - jπ D ( f a , v e ) f r 2 K m ( v e , f a ) - j 4 π c ( 1 D ( f a , v e ) - 1 ) R ref f r }
二维频域实现距离压缩和一致距离徙动校正的参考函数构造为:
H ( f a , f r ) = exp { jπ D ( f a , v e ) f r 2 K m ( v e , f a ) } exp { j 4 π c ( 1 D ( f a , v e ) - 1 ) R ref f r }
5b)将二维频域目标信号与距离压缩和一致距离徙动校正的参考函数相乘,即
S 2 ( f a , f r ) = S 1 ( f a , f r ) H ( f a , f r ) = W a ( f a - f ac ) W r ( f r ) exp { - j 4 πy e f c D ( f a , v e ) c } × exp { - j 4 πy e c f r } exp { j 4 π K m ( v e , f a ) c 2 ( 1 - D ( f a , v e ) ) ( y e D ( f a , v e ) - R ref D ( f a , v e ) ) 2 }
从而实现距离压缩和一致距离徙动校正。
5.根据权利要求1所述的基于变标算法的星载SAR地面加速运动目标成像方法,其特征在于,步骤7的具体子步骤为:
7a)根据距离压缩和一致距离徙动校正后的距离多普勒域目标信号的表达式
S 2 ( f a , t r ) = W a ( f a - f ac ) p r ( t r - 2 y e c ) exp { - j 4 πy e f c D ( f a , v e ) c } × exp { j 4 πK m ( v e , f a ) c 2 ( 1 - D ( f a , v e ) ) ( y e D ( f a , v e ) - R ref D ( f a , v e ) ) 2 }
目标方位调制和相位误差的相位项为:
exp { - j 4 π y e f c D ( f a , v e ) c + j 4 π K m ( v e , f a ) c 2 ( 1 - D ( f a , v e ) ) ( y e D ( f a , v e ) - R ref D ( f a , v e ) ) 2 } .
用于距离多普勒域方位压缩和相位校正的参考函数构造为:
H a ( f a ) = exp { j 4 πy e f c D ( f a , v e ) c } exp { - j 4 πK m K m ( v e , f a ) c 2 ( 1 - D ( f a , v e ) ) ( y e D ( f a , v e ) - R ref D ( f a , v e ) ) 2 }
7b)将距离压缩和一致距离徙动校正后的距离多普勒域目标信号与构造的用于距离多普勒域方位压缩和相位校正的参考函数相乘,即
S 3 ( f a , t r ) = S 2 ( f a , t r ) H a ( f a ) = W a ( f a - f ac ) p r ( t r - 2 y e c )
从而实现方位压缩和相位校正。
CN201410143273.1A 2014-04-10 基于变标算法的星载sar地面加速运动目标成像方法 Expired - Fee Related CN104035093B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201410143273.1A CN104035093B (zh) 2014-04-10 基于变标算法的星载sar地面加速运动目标成像方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201410143273.1A CN104035093B (zh) 2014-04-10 基于变标算法的星载sar地面加速运动目标成像方法

Publications (2)

Publication Number Publication Date
CN104035093A true CN104035093A (zh) 2014-09-10
CN104035093B CN104035093B (zh) 2016-11-30

Family

ID=

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104833972A (zh) * 2015-05-08 2015-08-12 电子科技大学 一种双基地调频连续波合成孔径雷达频率变标成像方法
CN107102330A (zh) * 2017-04-24 2017-08-29 清华大学 一种面向星载地理参考条带sar的wnlcs成像方法
GB2552186A (en) * 2016-07-13 2018-01-17 Aptcore Ltd Signal Processing Apparatus for Generating a Range-Doppler Map
CN108061890A (zh) * 2016-11-08 2018-05-22 北京机电工程研究所 一种sar成像方法
CN110471060A (zh) * 2019-07-08 2019-11-19 中国空间技术研究院 一种高效的雷达动目标调频率估计方法
CN112505647A (zh) * 2020-09-25 2021-03-16 北京空间飞行器总体设计部 一种基于序贯子图像序列的动目标方位速度估计方法
CN113484861A (zh) * 2021-07-07 2021-10-08 中国科学院空天信息创新研究院 一种基于频率变标的THz步进频率SAR实时高分辨率成像方法
CN117518111A (zh) * 2023-09-15 2024-02-06 鹰飒科技(江苏)有限公司 回波信号处理方法、装置、电子设备及存储介质

Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101819274A (zh) * 2010-03-25 2010-09-01 北京航空航天大学 用于合成孔径雷达前斜视子孔径成像处理的扩展非线性变标方法
CN101915919A (zh) * 2010-07-12 2010-12-15 北京航空航天大学 利用北斗卫星信号的双基sar成像系统
CN102759734A (zh) * 2012-05-23 2012-10-31 北京理工大学 一种用于高分辨率星载sar的改进距离徙动校正算法的成像方法
JP2012242217A (ja) * 2011-05-18 2012-12-10 Mitsubishi Electric Corp 画像レーダ装置
CN103149554A (zh) * 2013-02-02 2013-06-12 西安电子科技大学 双基sar的变标逆傅里叶变换成像方法
CN103293520A (zh) * 2013-05-13 2013-09-11 西安电子科技大学 基于等效距离方程的sar地面运动目标成像方法
CN103323841A (zh) * 2012-09-03 2013-09-25 中国科学院电子学研究所 星机载移变式双基合成孔径雷达中的成像方法及装置
EP2650695A1 (en) * 2012-08-02 2013-10-16 Institute of Electronics, Chinese Academy of Sciences Imaging method for synthetic aperture radar in high squint mode
CN103412310A (zh) * 2013-08-26 2013-11-27 电子科技大学 双基地前视合成孔径雷达地面动目标检测方法与成像方法
CN103576150A (zh) * 2013-09-24 2014-02-12 西安电子科技大学 基于高超声速飞行器俯冲段的前斜视sar成像方法
CN103630896A (zh) * 2012-08-28 2014-03-12 中国科学院电子学研究所 非线性调频尺度成像的方法

Patent Citations (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101819274A (zh) * 2010-03-25 2010-09-01 北京航空航天大学 用于合成孔径雷达前斜视子孔径成像处理的扩展非线性变标方法
CN101915919A (zh) * 2010-07-12 2010-12-15 北京航空航天大学 利用北斗卫星信号的双基sar成像系统
JP2012242217A (ja) * 2011-05-18 2012-12-10 Mitsubishi Electric Corp 画像レーダ装置
CN102759734A (zh) * 2012-05-23 2012-10-31 北京理工大学 一种用于高分辨率星载sar的改进距离徙动校正算法的成像方法
EP2650695A1 (en) * 2012-08-02 2013-10-16 Institute of Electronics, Chinese Academy of Sciences Imaging method for synthetic aperture radar in high squint mode
CN103630896A (zh) * 2012-08-28 2014-03-12 中国科学院电子学研究所 非线性调频尺度成像的方法
CN103323841A (zh) * 2012-09-03 2013-09-25 中国科学院电子学研究所 星机载移变式双基合成孔径雷达中的成像方法及装置
CN103149554A (zh) * 2013-02-02 2013-06-12 西安电子科技大学 双基sar的变标逆傅里叶变换成像方法
CN103293520A (zh) * 2013-05-13 2013-09-11 西安电子科技大学 基于等效距离方程的sar地面运动目标成像方法
CN103412310A (zh) * 2013-08-26 2013-11-27 电子科技大学 双基地前视合成孔径雷达地面动目标检测方法与成像方法
CN103576150A (zh) * 2013-09-24 2014-02-12 西安电子科技大学 基于高超声速飞行器俯冲段的前斜视sar成像方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
张立峰等: "多通道SAR-GMTI通道均衡和动目标检测定位方法", 《西安电子科技大学学报(自然科学版)》 *
白雪茹等: "基于多普勒频谱相关的加速运动目标ISAR成像", 《系统工程与电子技术》 *

Cited By (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104833972A (zh) * 2015-05-08 2015-08-12 电子科技大学 一种双基地调频连续波合成孔径雷达频率变标成像方法
GB2552186A (en) * 2016-07-13 2018-01-17 Aptcore Ltd Signal Processing Apparatus for Generating a Range-Doppler Map
CN108061890A (zh) * 2016-11-08 2018-05-22 北京机电工程研究所 一种sar成像方法
CN108061890B (zh) * 2016-11-08 2021-07-13 北京机电工程研究所 一种sar成像方法
CN107102330A (zh) * 2017-04-24 2017-08-29 清华大学 一种面向星载地理参考条带sar的wnlcs成像方法
CN110471060A (zh) * 2019-07-08 2019-11-19 中国空间技术研究院 一种高效的雷达动目标调频率估计方法
CN110471060B (zh) * 2019-07-08 2021-03-26 中国空间技术研究院 一种高效的雷达动目标调频率估计方法
CN112505647A (zh) * 2020-09-25 2021-03-16 北京空间飞行器总体设计部 一种基于序贯子图像序列的动目标方位速度估计方法
CN112505647B (zh) * 2020-09-25 2023-11-10 北京空间飞行器总体设计部 一种基于序贯子图像序列的动目标方位速度估计方法
CN113484861A (zh) * 2021-07-07 2021-10-08 中国科学院空天信息创新研究院 一种基于频率变标的THz步进频率SAR实时高分辨率成像方法
CN113484861B (zh) * 2021-07-07 2023-07-25 中国科学院空天信息创新研究院 一种基于频率变标的THz步进频率SAR实时高分辨率成像方法
CN117518111A (zh) * 2023-09-15 2024-02-06 鹰飒科技(江苏)有限公司 回波信号处理方法、装置、电子设备及存储介质

Similar Documents

Publication Publication Date Title
CN102393518B (zh) 一种适用于大斜视角的机载sar成像方法
CN104007440B (zh) 一种加速分解后向投影聚束合成孔径雷达成像方法
CN102749621B (zh) 一种双基地合成孔径雷达频域成像方法
CN102226841B (zh) 基于高阶多项式距离方程的同步轨道sar成像方法
CN103207387B (zh) 一种机载相控阵pd雷达杂波的快速模拟方法
CN103913741B (zh) 一种合成孔径雷达高效自聚焦后向投影bp方法
CN106970386A (zh) 一种雷达多普勒波束锐化的优化方法
CN103412310B (zh) 双基地前视合成孔径雷达地面动目标检测方法与成像方法
CN102230964B (zh) 曲线轨迹模型下的geo sar调频率变标成像方法
CN104833972B (zh) 一种双基地调频连续波合成孔径雷达频率变标成像方法
CN103869311A (zh) 实波束扫描雷达超分辨成像方法
CN103293520B (zh) 基于等效距离方程的sar地面运动目标成像方法
CN105487074B (zh) 一种双基地合成孔径雷达数值距离多普勒成像方法
CN102819020B (zh) 一种俯冲模型方位向非线性变标的合成孔径雷达成像方法
CN104035095A (zh) 基于空时最优处理器的低空风切变风速估计方法
CN109270528B (zh) 基于全解析距离模型的一站固定式双站sar成像方法
CN103091674A (zh) 基于hrrp序列的空间目标高分辨成像方法
CN103235305B (zh) 星载超高分辨率滑动聚束sar成像方法
CN103336278A (zh) 多视角观测下前视三维sar成像方法
CN105093224A (zh) 一种大斜视合成孔径雷达成像处理方法
CN103308913A (zh) 一种高速飞行器载双天线前视sar解模糊算法
CN104749570A (zh) 一种移不变机载双基合成孔径雷达目标定位方法
CN103293521A (zh) 一种利用x波段雷达探测近海海域水深的方法
CN103983969A (zh) 基于二次拟合距离方程的地面加速运动目标成像方法
CN109188384B (zh) 空间目标回波动态观测的电磁仿真方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20161130

Termination date: 20210410