Loading [MathJax]/jax/element/mml/optable/BasicLatin.js
高级搜索

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

总变分图像复原方程的离散化方法

邹谋炎 刘小军

李东, 张成祥, 赵迪, 谭晓衡, 周喜川, 占木杨. 基于伪逆极坐标傅里叶变换的快速ISAR方位定标[J]. 电子与信息学报, 2019, 41(2): 262-269. doi: 10.11999/JEIT180299
引用本文: 邹谋炎, 刘小军. 总变分图像复原方程的离散化方法[J]. 电子与信息学报, 2001, 23(9): 861-867.
Dong LI, Chengxiang ZHANG, Di ZHAO, Xiaoheng TAN, Xichuan ZHOU, Muyang ZHAN. Fast Cross-range Scaling for ISAR Imaging Based on Pseudo Polar Fourier Transform[J]. Journal of Electronics & Information Technology, 2019, 41(2): 262-269. doi: 10.11999/JEIT180299
Citation: Zou Mouyan, Liu Xiaojun . A DISCRETIZATION METHOD FOR TOTAL VARIATION BASED IMAGE RESTORATION EQUATION[J]. Journal of Electronics & Information Technology, 2001, 23(9): 861-867.

总变分图像复原方程的离散化方法

A DISCRETIZATION METHOD FOR TOTAL VARIATION BASED IMAGE RESTORATION EQUATION

  • 摘要: 基于总变分的图像重建和复原导致一类最小化问题,它归结为解一个非线性椭圆型偏微分方程。为了得到数值解,必须将问题线性化和离散化。C.R.Vogel和M.E.Oman等人(1996)的定点迭代是一个良好的线性化方法。然而文献中已报道的离散化方法需要微分方程数值解的工具,比较繁杂。该文提出一种新的离散化方法,它只需要图像处理中的常规技术。图像反降晰和噪声抑制的实验结果表明该文的结果不亚于文献中报道的结果。
  • 逆合成孔径雷达(Inverse Synthetic Aperture Radar, ISAR)可以全天候、全天时、远距离获取非合作目标的高分辨图像,在军事和民用上都具有非常重要的应用价值[1, 2]。高效的距离多普勒(Range-Doppler, RD)算法或时频分析法可以得到非合作目标的高分辨ISAR图像;然而,所得的图像位于距离多普勒域,不能反映出目标的真实尺寸。为了后续更好地对目标进行识别,需要对目标ISAR图像进行方位定标(cross-range scaling);而估计转动角速度(Rotation Angle Velocity, RAV)是方位定标的一个难点和研究热点。

    针对该问题,近年来国内外学者先后提出多种ISAR方位定标算法[311]。文献[3]提出一种Radon-Ambiguity方法,通过在不同距离单元的平面上检测线性调频信号的调制率,进而估计出RAV;但在计算过程中需要2维搜索,计算量很大。文献[47]方法通过选择不同距离单元内的特显点,提取出其相位系数来估计RAV;然而在实际ISAR成像中很难找到理想的散射点,并且特显点的信噪比会影响RAV的估计精度,同时在不同距离单元估计出的RAV值起伏较大。文献[8]在转动速度、转动加速度的2维平面上以图像信息熵或对比度为准则进行搜索,当ISAR图像质量最好时,对应参数即为转动速度与转动加速度的估计值;然而,该方法需要进行多维搜索,无法兼顾算法的精度和运算效率。文献[9]提出一种基于压缩感知的方法,利用目标在成像空间内具有稀疏特性,实现目标高分辨成像和方位定标;但存在对噪声敏感且计算复杂度高等问题。文献[10]利用两幅RD图像之间的旋转相关特性估计RAV,此方法是在假设旋转中心已知的情况进行的,然而在实际图像处理中旋转中心很难获取,限制了此方法的实用性。文献[11]利用2维快速傅里叶变换(2-D FFT)以及极坐标映射估计RAV,该方法将ISAR图像的旋转变为极坐标域中沿极角的平移,进而估计RAV。虽然避免了旋转中心的估计,但是由2维频域向极坐标转换过程中的插值操作会带来RAV估计精度损失和计算复杂度高的问题。因此,需要进一步研究高效的ISAR图像方位定标方法。

    本文提出一种基于伪逆极坐标快速傅里叶变换(Pseudo Polar Fast Fourier Transform, PPFFT)的高效ISAR方位定标方法。首先利用高效的PPFFT把两幅不同时刻ISAR图像的旋转转化为伪极坐标域中沿极角方向的平移;然后定义了一种新的代价函数-积分相关函数,来粗估目标的RAV;最后,采用二分法快速得到最优的RAV,完成方位定标。相比于现有ISAR定标算法,本文算法避免了插值操作带来的精度损失和高计算复杂度问题。

    图1给出了经过平动补偿后ISAR成像的几何模型,其中目标以均匀角速度ω0旋转,雷达与目标旋转中心之间的距离为ra;则对旋转目标上任一散射点P(r0,θ0)tk时刻与雷达之间的瞬时斜距可写为

    图 1  ISAR成像的几何模型示意图
    r(tk)=r20+r2a2r0racos(θ0+ω0tk) (1)

    由于雷达到目标的距离远大于目标尺寸,在远场条件下照射到目标的电磁波可用平面波近似;因此,斜距表达式(1)可近似写为

    r(tk)rax0cos(ω0tk)+y0sin(ω0tk) (2)

    其中,x0=r0cosθ0;y0=r0sinθ0,根据距离时间延迟和多普勒频率变化特性,则该散射中心P在距离多普勒域上的坐标可表示为[10,11]

    X(tk)=fs2r(tk)c=X0+[x0cos(ω0tk)y0sin(ω0tk)]/ηr (3)
    Y(tk)=MPRF2λdr(tk)dtk=x0sin(ω0tk)+y0cos(ω0tk)ηa (4)

    其中,X0=ra/ηr为常数项,fs是距离向采样频率,PRF是脉冲重复频率,c是光速,λ是波长,M是累积脉冲数,ηr=c/(2fs)ηa=λPRF/(2Mω0)分别是距离向尺度因子和方位向尺度因子。

    把式(3)和式(4)表示为矩阵形式,可得散射点P在不同时刻tk1tk2得到的两幅RD图像I1I2中的位置具有式(5)关系[10,11]

    [Xc(tk2)Yc(tk2)]=SR(tk2)R1(tk1)S1[Xc(tk1)Yc(tk1)]=[cos(θd)Crsin(θd)sin(θd)/Crcos(θd)][Xc(tk1)Yc(tk1)] (5)

    式中,θd=ω0(tk2tk1)=ω0Δtk是两幅RD图像之间的旋转角度,Cr=ηa/ηr为相对尺度因子。从式(5)可知,第2幅RD图像I2是由第1幅RD图像I1经过旋转-拉伸变换得到的。下面就不同时刻得到的两幅RD图像在频域的关系进行分析;对给定的N×N大小的RD图像I(u,v),其2维傅里叶变换ˆI(ωx,ωy)可表示为

    ˆI(ωx,ωy)=N/21u,v=N/2I(u,v)2πiM(u ωx+v ωy) (6)

    根据傅里叶变换性质可知,对频率变量(ωx,ωy)在不同坐标下,可以得到不同的2维傅里叶变换形式;如2维频率变量在卡迪尔坐标采样网格下,(ωx,ωy)=(k,l), k,l=(M/2,M/21), M为采样点数,其2维傅里叶变换频域形式可写为

    ˆIc(k,l)I(k,l)=N/21u,v=N/2I(u,v)e2πiM(u ωx+v ωy)

    (7)

    从式(7)可知,由于在时域和频域都是均匀采样,因此可用2维快速傅里叶变换(2-D FFT)实现。如果2维频率变量(ωx,ωy)是在极坐标网格下采样,则在极坐标网格中:ωx=rkcosθl,ωy=rksinθl, rk=k,θk=2πl/L, k,l=(M/2,M/21);其中,ML分别表示沿极径rk和极角θl方向的采样点数,在极坐标网格下其2维傅里叶变换形式为

    ˆIp(ωx,wy)=Fp(rkcosθl,rksinθl)=++I(u,v)ej2πrk(ucosθl+vsinθl)dudv

    (8)

    根据式(5)和式(8),对不同时刻得到的两幅RD图像I1I2,在极坐标采样网格下做2维傅里叶变换,其在极坐标域的关系式为

    ˆI2(r,θ)=1Cr2ˆI1(r|Cr|,θ+θd) (9)

    从式(9)可知,两幅ISAR图像I1I2在RD域的旋转和拉伸分别对应于极坐标中极角的平移和极径的伸缩,然而当2维频率变量(ωx,ωy)处于极坐标网格时,其2维傅里叶变换没有相应的快速变换,计算复杂度急剧增加,给数据的实时处理带来困难。因此需要进一步研究高效的ISAR图像定标方法。

    根据第2节分析可知,不同时刻得到的两幅ISAR图像在RD域旋转对应为极坐标傅里叶变换域沿极角的平移;在此基础上,本节提出一种基于PPFFT的快速方位定标方法。

    伪逆极坐标傅里叶变换(PPFT)是一种在伪逆极坐标采样网络下的2维傅里叶变换,其2维频率变量在伪逆极坐标网格下采样形式为[12, 13]

    PP1P2 (10)

    式中,P1P2满足的条件为:P1{(2lNk,k)|N2lN2,NkN}, P2{(k,2lNk)|N2lN2,NkN},其中k表示“伪极径”,l表示“伪角”。图2(a)图2(b)分别表示式(10)中的伪逆极坐标网格下采样区域P1P2,完整的伪逆极坐标采样网格如图2(c)所示。根据对应的关系伪逆极坐标采样网格用(r,θ)表示可写为

    图 2  伪逆极坐标采样网格
    P1(k,l)=(r1k,θ1l), P2(k,l)=(r2k,θ2l) (11)

    其中,r1k,r2k,θ1l,θ2l满足的条件为:r1k=k4(lN)2+1,r2k=k4(lN)2+1, θ1l=π/2arctan(2lN),θ2l=arctan(2lN)

    类似于式(7)和式(8)定义的迪卡尔坐标和极坐标下2维傅里叶变换,ISAR图像在式(10)和图2表示的伪逆极坐标采样网格中,其2维傅里叶变换形式分别为[12,13]

    ˆI1PP(k,l)I(2lNk,k)=N/21u,v=N/2I(u,v)e2πiM(2lNku+kv)

    (12a)

    ˆI2PP(k,l)I(k,2lNk)=N/21u,v=N/2I(u,v)e2πiM(ku2lNkv) (12b)

    幸运的是伪逆极坐标采样网格下的2维傅里叶变换能够快速实现。图3描述了式(12a)快速实现的示意图;首先对原始图像沿y 轴方向补零,接着做快速傅里叶变换(FFT);然后,沿x 轴方向做Chirp-Z变换[13],得到了2维频率在伪逆极坐标采样网格下对应的区域P1;同理,可以得到对应的区域P2。将P1P2叠加就可得到整个图像的在伪逆极坐标采样网格下2维傅里叶变换结果。实现过程仅需进行2次1维FFT和2次1维Chirp-Z变换[12,13],避免了2维插值操作。因此PPFFT具有高的精度和较低的计算复杂度。

    图 3  伪逆极坐标采样区域P1的2维傅里叶变换快速实现示意图

    根据3.1节的PPFFT的原理,首先对不同时刻tk1tk2得到的两幅RD图像I1I2进行PPFFT,用Pk1Pk2分别表示两幅RD图像PPFFT后的幅度信息,则其具有式(13)关系:

    Pk2(rk,θl)=(ωx0,ωy0)+1Cr2Pk1(rk|Cr|,θl+Δθ) (13)

    其中,(ωx0,ωy0)为旋转中心,其位于2维频域零频处;式(13)表明两幅RD图像旋转对应伪逆极坐标系下沿伪极角θ的平移。因此,只要估计出tk2tk1时间段内两幅图像的旋转角度Δθ,就可以得到目标的RAV。为了估计两幅ISAR图像的旋转角度Δθ;本文定义了一种新的1维代价函数-径向积分相关函数,其表达式为

    Cor(k)=Ni=1φ1(θi)×φ2(θik) (14)

    其中,φm(θi)=Nj=1Pk,m(rj,θi),m=1,2为两幅RD图像经PPFFT后沿径向进行积分的结果;根据式(14),两幅RD图像的旋转角度Δθ可通过式(14)取最大值时的位置k获得。为了简化计算,利用傅里叶变换的卷积特性,最大值对应的位置k可由式(15)获得。

    k=argmax (15)

    式中,{F_\theta }\left( \cdot \right)表示傅里叶变换,{F_\theta }^{ - 1}\left( \cdot \right)表示傅里叶逆变换,下标\theta 表示沿极角方向。根据获得的两幅RD图像的旋转角度\Delta \hat \theta ,就可得到RAV的估计值\hat \omega = {{\Delta \theta } / {\left( {{t_2} - {t_1}} \right)}}以及相对尺度因子{C_{r0}}

    然而从式(5)可知,两幅图像之间除了旋转以外还存在尺度伸缩;尺度伸缩会对RAV的估计精度产生影响,所以为了更精确地估计RAV,需要粗估RAV对原始图像做旋转-拉伸补偿,从而精确估计出旋转角度。这里使用二分法求解最优的RAV;二分法是一种能够重新分配搜索区间并选择根所在的子区间进行2次搜索,进一步寻找真实根的经典优化算法[14]。本文利用二分法来得到最优的旋转角度,选择合适的搜索区间\left[ {a,b} \right],在小区间内实现迭代搜索,在完成多次迭代搜索就可得到目标精确RAV,从而完成ISAR图像方位定标。整个算法的实现过程如图4所示。

    图 4  提出的ISAR方位定标算法实现流程图

    本小节分析提出方法的计算复杂度,并与文献[11]提出的基于2维快速傅里叶变换和极坐标映射估计RAV的方法进行比较。假设距离维采样点数为{N_r},方位维采样点数为{N_a};且长度为N的信号做1次复数相乘需要的计算量为O\left( N \right),做1次FFT计算复杂度为O\left( {N{{\log }_2}N} \right)

    文献[11]提出的估计RAV方法主要有以下几步构成:首先,对不同时刻得到的两幅RD图像进行2维FFT操作变到2维频域;然后,利用插值操作从2维频域映射到极坐标域;最后,在极坐标域定义2维相关函数来粗估RAV,并利用优化算法来迭代求解最优RAV。其总的计算复杂度大约为

    \begin{align} {C_{[11]}} =& O\left[ \left( {Q + 2} \right)\left( {\left( {5{\rm{ + }}N_{\ker }^2} \right){N_r}{N_a}{\rm{ + }}{N_r}{N_a}{{\log }_2}{N_r}{N_a}} \right)\right. \\ & \left.+ \left( {Q + 1} \right)\left( {3{N_r}{N_a}{{\log }_2}\left( {{N_r}{N_a}} \right)} \right) \right] \end{align} (16)

    其中,{N_{\ker }}为插值核长度,Q为算法优化过程需要的迭代次数。

    根据图4提出的算法实现流程图,提出的ISAR方位定标算法主要有以下几步构成;首先对不同时刻得到的两幅RD图像进行PPFFT变换到伪逆极坐标域,N点Chirp-Z变换需要的计算复杂度为O\left( {N + 2N{{\log }_2}N} \right)[15];然后,在伪逆极坐标域定义1维积分相关函数来粗估RAV,并利用二分法来迭代求解最优RAV。提出的算法总的计算复杂度大约为

    \begin{align} {C_{{{本文}}}} =& O\left[ \left( {M + 2} \right)\left( {6{N_r}{N_a}{{\log }_2}\left( {{N_r}{N_a}} \right) + 4{N_r}{N_a}} \right) \right. \\ & \left. + \left( {M + 1} \right)\left( {3{N_a}{{\log }_2}{N_a}} \right) \right] \end{align} (17)

    其中,M为二分法求解过程所需要的迭代次数。

    根据式(16)和式(17)可知,文献[11]的方法由于需要使用插值操作和2维代价函数,导致其计算复杂度偏高。而本文提出的方位定标算法由于避免了插值操作,且定义了新的1维代价函数-积分相关函数;因此,其计算复杂度远低于文献[11]的方法。

    为了验证本文算法的有效性,现把所提方法应用到飞机模型的ISAR方位定标中。雷达工作参数和飞机目标运动参数见表1图5为经过平动补偿后在{t_{k1}}{t_{k2}}两个不同时刻进行RD成像所得到的ISAR图像,很明显第2幅RD图像是第1幅RD图像的旋转-拉伸结果。

    表 1  雷达参数和目标运动模型
    参数数值
    载波频率5.6 GHz
    波长0.0536 m
    传输信号带宽400 MHz
    距离采样频率512 MHz
    脉冲重复频率150 Hz
    有效回波脉冲600
    旋转角速度0.0436 rad/s
    下载: 导出CSV 
    | 显示表格
    图 5  两个不同时刻得到的RD图像

    {t_{k1}}{t_{k2}}两个时刻的RD图像进行PPFFT,得到图像在伪逆极坐标下的变换结果,将图像的旋转转化为沿伪极角方向的平移;接着对伪逆极坐标域的图像沿径向积分;图6(a)为利用式(14)定义的代价函数得到的相关曲线,根据峰值检测可得旋转目标的RAV为{\omega _0} = 0.0455 rad/s,最后根据式(5)对RD图像进行旋转-拉伸补偿并完成方位定标。

    图 6  粗估的1维相关曲线

    然而,由于受两幅RD图像之间拉伸的影响,一次粗估所得的RAV存在较大的误差,需要根据粗估得到的RAV补偿图像拉伸带来的影响,再进一步精细估计最优的RAV。采用二分法求最优RAV的迭代收敛曲线如图7所示,算法迭代了10次左右就趋于收敛,计算得到最优的RAV为{\omega _0} = 0.04358 rad/s。

    图 7  迭代收敛曲线

    为了更好地验证本文方法优势,本文和文献[11]方法进行比较。在本文仿真参数下,文献[11]方法得到最优的RAV为{\omega _0} = 0.04342 rad/s,定标后ISAR图像如图8(a)所示。而利用本文方法得到最优的RAV为{\omega _0} = 0.04358 rad/s,定标后ISAR图像如图8(b)所示。定标结果与目标真实尺寸基本一致;在RAV估计精度方面,文献[11]方法是0.41%,而本文方法是0.05%,说明本文方法能够对图像旋转角进行更为精确的估计,定标性能更好。

    图 8  方位定标后的ISAR图像

    为了更清晰直观比较计算复杂度,使用Intel Double-core处理器,CPU主频2.4 GHz,内存4 GB, Window 7操作系统,32位台式主机,在MATLAB 2014a环境下,分别运行了文献[11]方法和本文方法,运行时间如表2所示。由表2可知,本文提出的定标方法计算量远小于文献[11]方法。为了验证分析噪声环境下算法的稳健性,图9给出了本文方法和文献[11]方法在不同信噪比(SNR)下RAV估计的均方根误差(RMSE);其中,RMSE获得是通过100次独立的门特卡罗实验。从图9可知,在SNR较高时,相比文献[11]方法,本文方法具有更好的估计精度。随着SNR的降低,本文方法和文献[11]方法RAV估计精度趋于一致,但其计算复杂度远低于文献[11]的方法,与上述的仿真结果和理论分析相吻合。

    表 2  两种方法运行时间对比(s)
    方法名称粗估所用时间定标总时间
    文献[11]方法107.224958.659
    本文方法12.20396.712
    下载: 导出CSV 
    | 显示表格
    图 9  RAV估计的均方根误差随信噪比变化曲线

    Yak-42飞机模型如图10(a)所示,雷达载频{f_c}为5.52 GHz,脉冲重复频率{\rm PRF}为25 Hz,发射信号带宽{B_r}为400 MHz,方位向脉冲回波个数为128。对Yak-42实测时所得的回波数据进行旋转角度估计,计算得到最优的RAV为{\omega _0} = 0.0092 rad/s;根据估计所得出的RAV进行方位定标,得到定标后ISAR图像如图10(c)所示。同时,图10(b)给出了利用文献[11]方法方位定标后的ISAR图像。从图10可以看出,文献[11]方法定标得出的Yak-42的长度约为30.65 m,翼展约为31.86 m;本文定标得出的Yak-42的长度约为32.82 m,翼展约为31.24 m,更接近图10(a)中Yak-42飞机的实际尺寸;并且飞机轮廓清晰,几何尺寸明显,有利于后续对运动目标高效识别。

    图 10  实测数据图像

    本文提出一种基于PPFFT快速ISAR方位定标方法。该方法主要思想是通过PPFFT将RD图像的旋转运动映射为伪逆极坐标域中沿极角方向的平移运动,并通过定义了新的1维代价函数-积分相关函数和经典二分法估计得到精确的RAV。提出的方法避免了插值操作和2维代价函数,具有较低的计算复杂度,在工程实现中具有一定的优势。最后,计算机仿真和实测数据结果验证了本文方法的有效性。

  • L.I. Rudin, S. Osher, Feature-oriented image enhancement using shock filters, SIAM J. Num.Anal., 1990, 27, 919-940.[2]L.I. Rudin, S. Osher, E. Fatemi, Nonlinear total variation based noise removal algorithms.[J]. Phvsica D.1992,60:259-[3]T.F. Chan, H. M. Zhou, R. H. Chan, Continuation method for total variation denoising problems,In SPIE 1995, vol.2563, Advanced Signal Processing Algorithms, F. T. Luk, Ed., San Diego, CA.1995, [4]Y. Li, F. Santosa, A computational algorithm for minimizing total variation in image restoration,IEEE Trans. on Image Processing, 1996, 5(6), 987-995.[4]C.R. Vogel, M. E. Oman, Iterative methods for total variation denoising, SIAM J. Sci. Comput.,1996, 17(1), 227-238.[5]C.R. Vogel, M. E. Oman, Fast, robust total variation-based reconstruction of noisy, blurred images, IEEE Trans. on Image Processing, 1998, 7(6), 813 824.[6]R.E. Ewing, J. Shen, A multigrid algorithm for the cell-centered finite difference scheme, in Proc.6th Copper Mountain Conf. Multigid Methods. NASA Conf. Pub. 3224, April, 1993. [8]T.F. Chan, Chiu-Kwong, Total variation blind deconvolution, IEEE Trans. on Image Processing,1998, 17(3), 370-375.[7]Zou Mou-yan, R. Unbenhauen, On the computational model of a kind of deconvolution problems,IEEE Trans. on Image Processing, 1995, 4(10), 1464-1467.[8]E.H. Golub, C. F. Van Loan, Matrix Computations, Baltimore, MD: Johns Hopkins Univ. Press,1989. [11]C. Charalambous, F. K. Ghaddar, K. Kouris, Two iterative image restoration algorithms with application to nuclear medicine, IEEE Trans. on Medical Imaging, 1992, 11(1), 2-8.[9]Zou Mou-yan, Zou Xi, Global optimization: An auxiliary cost function approach, IEEE Trans.on Systems, Man, and Cybernetics, Part A, 2000, 30(3), 347-354.
  • 期刊类型引用(3)

    1. 孙培育,史悦. 基于离散傅里叶变换的数字图像空间域水印方法. 微型电脑应用. 2023(03): 132-134 . 百度学术
    2. 田丹. 基于MEMS传感器的室内外多源融合导航系统设计. 计算机测量与控制. 2022(08): 289-295 . 百度学术
    3. 刘也,叶钒,马岩,赵华. 基于高精度运动信息的ISAR分辨率评估方法. 电子与信息学报. 2021(09): 2728-2734 . 本站查看

    其他类型引用(1)

  • 加载中
计量
  • 文章访问数:  2433
  • HTML全文浏览量:  105
  • PDF下载量:  502
  • 被引次数: 4
出版历程
  • 收稿日期:  2000-01-21
  • 修回日期:  2000-09-14
  • 刊出日期:  2001-09-19

目录

/

返回文章
返回