CN104931968A - 基于FPGA的InSAR通道幅相误差估计方法 - Google Patents

Abstract

本发明公开了一种基于FPGA的通道间幅相误差特性频域加窗的估计方法，主要解决现有雷达系统多通道间幅相特性不一致对通道间相干性降低的问题。其实现步骤是：(1)输入通道数据；(2)对通道信号进行时频转换；(3)对通道信号的频谱进行频谱比对；(4)对频谱比对结果进行FFT变换；(5)对频谱比对的FFT变换结果进行加窗处理；(6)对加窗处理后的数据进行逆FFT操作，得到通道间的幅相误差特性估计。本发明充分利用了FPGA快速高效的FFT时频转换和时序逻辑控制功能，在较大杂散和噪声下可获得高质量的幅相误差特性估计结果，降低了杂散和噪声对估计结果的影响，可用于双通道和/或多通道的雷达实时信号处理系统。

Description

基于FPGA的InSAR通道幅相误差估计方法

技术领域

本发明属于信号处理技术领域，特别涉及一种干涉合成孔径雷达InSAR的幅相误差估计方法，可用于双通道和/或多通道的雷达实时信号处理系统。

背景技术

多通道的雷达信号处理系统由于制作工艺的差异，通道间不可避免的存在幅相误差，对于合成孔径雷达InSAR系统，幅相误差的存在会严重恶化信号的峰值旁瓣比PSLR、积分旁瓣比ISLR以及成像分辨率，使得InSAR干涉相位图产生偏差，降低图像的相干性，影响最终DEM的测高性能。研究实时高效的误差估计方法是解决通道间幅相失配的根本出发点，也是近年来现代雷达系统的一个研究热点。

干涉相位图是将雷达的主辅天线将各自通道接收到的信号经SAR成像后做干涉处理，得到一幅表示两个天线接收信号相位差的图像。基于FPGA的幅相误差估计是在实时处理过程中，实时估计出各通道相对于参考通道的幅相差异，并在各通道中进行校正，使各通道间幅相误差的差异最小，提高干涉相位图的相干性。现有的基于各通道回波信号进行频谱比对的方法，虽能够获得较为准确的幅相误差估计结果，但是当通道间幅相差异比较明显时，各种杂散噪声的影响会使得估计结果误差很大，使干涉相位图失配，甚至无法获得正确的干涉相位图。

何志华等在文章“干涉SAR模拟器通道幅相特性实时校正方法”(电子学报,2013,41(9):1710-1715)中采用时域相关加窗方法提取通道幅相特性。在数学运算中，时域相关的求解运用的是卷积的方法。根据傅里叶性质，时域上的卷积可以转化为频域上的相乘，从而减少了运算量，因此时域相关处理方法在频域内更容易实现。因而上述方法的不足之处是：在FPGA的实时处理过程中需要多次使用FFT进行时频转换，从而占用FPGA中大量的资源，运算量巨大。并且每次时频转换需由对应的时序控制，随着时频转换次数的增多对于控制时序的产生也提出更严格要求，容易出错。

发明内容

本发明的目的在于提出一种基于FPGA的InSAR通道幅相误差估计方法，以解决上述现有技术时序难以控制、运算量巨大和资源消耗过大的问题。

实现本发明目的的技术方案是：通过将各通道信号的频谱基于参考通道频谱做比对，作为通道间幅相误差初步估计的系统函数；通过对系统函数在频谱内加矩形窗进行重采样，并转回时域，得到最终估计的幅相误差特性函数。其主要步骤包括如下：

(1)输入通道数据：把双通道雷达实时信号处理系统中的一个通道作为主通道，另一个通道作为辅通道，并将主通道数据作为参考通道数据，辅通道数据作为待校准通道数据；

(2)分别对主、辅通道数据进行N₂点的FFT变换，得到主通道数据的频域形式为M(f)、辅通道数据的频域形式为S(f)；其中N₂＝2^nextpow2(N1)，nextpow2(N₁)为在大于等于N₁的二的整数次方中最小的数的幂；N₁为输入通道数据的复数据点数；

(3)根据步骤(2)的结果，获取主通道与辅通道之间的幅相误差估计值：

$H_1\left(f\right)=\frac{M\left(f\right)}{S\left(f\right)};$

(4)对主、辅通道间的幅相误差估计值H₁(f)进行FFT变换，得到主、辅通道间的频域幅相误差估计值H₂(f)；

(5)在频域幅相误差估计值H₂(f)的零频附近进行长度为64个单元的矩形窗截断，获取截断后的频谱并对频谱中的数据进行截位；

(6)对截断后的频谱进行IFFT变换，得到最终加窗后的幅相误差估计结果

本发明与现有的技术相比具有以下优点：

第一，本发明由于对通道间的频谱比对结果进行频域加窗操作，能够去除通道频谱比对结果中的高频分量，减少了各种杂散噪声的影响，使估计结果变得平滑，提高了干涉相位图的相干性。

第二，本发明由于仅使用一条主通道数据信号和一条辅通道数据信号，就可以得到比较理想的估计结果，因而具有运算量小，操作简便，易于实现的特点；

第三，在本发明中只使用了一次FFT和逆FFT变换，因而克服了时域相关加窗法由于多次时频转换而消耗大量资源的不足。

附图说明

图1为本发明的实现流程图；

图2为本发明仿真时在辅通道上添加的幅相误差特性曲线图；

图3为本发明仿真时的主、辅通道信号频谱图；

图4为用本发明仿真得到的通道间频谱比对的幅相特性估计结果图；

图5为用本发明仿真得到的频谱比对结果数据的频谱图；

图6为用本发明仿真得到的通道间幅相误差特性估计结果图。

具体实施方式

下面结合附图对本发明做进一步的描述。

参照图1，本发明在FPGA内实现的具体实施步骤如下：

步骤1，输入通道数据。

在进行通道间幅相误差特性估计时，将双通道雷达实时信号处理系统中的一个通道作为主通道，另一个通道作为辅通道；

把主通道和辅通道接收到的线性调频信号分别通过A/D转换和下变频操作，转换为16位定点的主通道实虚部数据和辅通道实虚部数据，这些数据包含了主通道与辅通道之间的幅相误差特性；

将主通道的16位定点实虚部数据作为参考通道数据，将辅通道的16位定点实虚部数据作为待校准通道数据。

步骤2，对通道数据进行时频转换。

对于雷达实时信号处理系统来说，由于其控制复杂，数据量大，所以系统主要是基于FPGA+DSP硬件架构来实现的。在这种结构中，FPGA主要用来采集数据、参数计算、距离向脉冲压缩、方位向预滤波、控制信号的产生与传输，DSP主要作用是对FPGA传输过来的雷达回波数据进行干涉处理。

在FPGA所承担的数据处理过程中，FFT变换是其最重要的应用之一，而FPGA内部的FFT变换是通过调用IP核实现，其输入数据点数必须为2ⁿ点，n为正整数。

通道数据的时频转换步骤如下：

(2a)判断数据输入点数N₁是否为2ⁿ点，若不是，则进行步骤(2b)；若是，则进行步骤(2c)。

(2b)在FPGA内通过计数器cnt_nr将输入数据的末端补零至N₂点，N₂为大于N₁的最小的二的整数次方数，并产生该N₂点数据的匹配使能信号；

(2c)在FPGA内通过计数器cnt_nr产生FFT变换过程所需的起始信号fft_sop和结束使能信号fft_sop：当cnt_nr＝1时，产生一个宽度为一个时钟周期的fft_sop高电平信号，作为FFT变换输入数据的开始信号；当cnt_nr＝N₂时，产生一个宽度为一个时钟周期的fft_eop高电平信号，作为FFT变换输入数据的结束信号；

(2d)设置FFT IP核的工作参数及模式：

FFT IP核的输入数据为：16位位宽的定点实虚部数据，FFT变换的点数N₂；

输出数据为：16位位宽的定点实虚部数据，6位位宽的指数项数据以及匹配的有效信号；

将工作模式设置为突发缓冲模式；将inverse信号置为1'b0，表示FFT IP核要执行FFT变换；将sink_error信号置为2’b00，表示无误差输入输出；将source_ready信号置为1’b0，表示输出准备完毕，可以输出。

步骤3，频谱比对。

频谱比对是为了获取两个通道间幅相误差的差异，得到待校准通道相对于参考通道的幅相误差特性，本发明选择在FPGA内完成频谱比对过程。该过程的除法主要采用FPGA内浮点型除法IP核实现。

本步骤的具体实现如下：

(3a)设由步骤2获得的主通道频谱数据的实虚部为a,b；辅通道频谱数据的实虚部为c,d；则主通道频谱M(f)＝a+jb，辅通道频谱S(f)＝c+jd，j为纯虚数；计算频谱比对H_e(f)为：

$\begin{array}{c}{H}_1\left(f\right)=\frac{M\left(f\right)}{S\left(f\right)}=\frac{a+\mathrm{jb}}{c+\mathrm{jd}}=\frac{(\mathrm{ac}+\mathrm{bd})+j(\mathrm{bc}-\mathrm{ad})}{c^2+d^2}\\ =\frac{\mathrm{ac}+\mathrm{bd}}{c^2+d^2}+j\frac{\mathrm{bc}-\mathrm{ad}}{c^2+d^2}\end{array}$

(3b)在FPGA内直接通过乘加运算,得到32位位宽的定点数据E,F,M值，其中：

E＝ac+bd

F＝bc-ad；

M＝c²+d²

(3c)判断M是否为0，若为0值，则将该时刻对应的E,F,M值均置为常数1，若不为0值，则保持E,F,M为原值；

(3d)将32位位宽的E,F值左移9位并更新位宽为41位，M值保持不变；

(3e)使用FPGA内提供的ALTFP_CONVERT IP核完成步骤(3d)中数据E,F,M从定点类型到浮点类型的转换，在ALTFP_CONVERT IP核的设置中，Mode选项选择“An integer to floating point”，Integer data选项中Custom设置为41位，当对M进行从定点类型到浮点类型的转化时，Integer data选项中Custom设置为32位。输出result均选择为32位。该过程，数据的输出比输入延迟6个时钟周期；

(3f)使用FPGA内提供的ALTFP_DIV IP核完成步骤(3e)中浮点数据的除法过程，输入输出数据位宽均选择为32位。该过程数据的输出比输入延迟14个时钟周期；

(3g)使用FPGA内提供的ALTFP_CONVERT IP核完成步骤(3f)中的数据从浮点类型向定点类型的转换，在ALTFP_CONVERT IP核的设置中，Mode选项中选择“Afloating point to an integer”，Integer data选项中设置为32位，输入输出数据位宽均选择为32位，该过程数据的输出比输入延迟6个时钟周期；

(3h)选取步骤(3g)输出结果的数据最高位为符号位，低15位为有效数据位，重新组合成新的16位数据输出，作为该频谱比对模块的输出结果，并详细计算各IP核工作所需的匹配使能信号。

此时，由步骤3得到的频谱比对数据受到杂散噪声的影响，误差比较大，毛刺现象比较严重，为使估计结果更准确平滑，需要对频谱比对数据进行频域加窗处理。

步骤4，FFT变换。

由于步骤2进行了N₂点的FFT变换，所以步骤3进行处理的数据点数为N₂点，因此步骤3结果的数据点数为N₂点，可以直接调用FFT IP核进行FFT变换。

FFT变换的具体步骤如下：

(4a)FPGA内通过计数器cnt_fft产生FFT变换过程的起始信号fft_sop和结束使能信号fft_eop：当cnt_fft＝1时，产生一个宽度为一个时钟周期的fft_sop高电平信号，作为FFT变换输入数据的开始信号；当cnt_fft＝N₂时，产生一个宽度为一个时钟周期的fft_eop高电平信号，作为FFT变换输入数据的结束信号；

(4b)设置FFT IP核的工作参数：

将inverse信号置为1'b0，表示FFT IP核要执行FFT变换；将sink_error信号置为2’b00，表示无误差输入输出；将source_ready信号置为1’b0表示输出准备完毕，可以输出。

步骤5，加窗操作。

对步骤4中FFT变换之后的实虚部数据进行添加矩形窗截断，去除高频分量，保留低频分量，由于零频在频谱的两端，所以矩形窗添加在该频谱的两端。

加窗操作的具体步骤为：

(5a)通过计数器cnt_num保留步骤4的结果在与区间内的实虚部数据，其余范围内的数据均置0，其中为所选矩形窗长度M的一半，M选择32或64或128均可；

(5b)根据步骤4中FFT变换后指数项的值，对步骤(5a)中计算结果的数据进行截位，得到截位后的值result：

result＝data×2^-exp/2ⁿ，

其中，data为FFT IP核输出结果中数据的实部或虚部，exp为FFT IP核输出结果中的指数项数据，n为保证result值的最大模值不超过2¹⁵-1的一个量值，该公式的运算过程可理解为：当data为正数时，在右移(exp+n)位时高位补0；当data为负数时，在右移(exp+n)位时高位补1。

步骤6，逆FFT操作。

将步骤5中的输出结果做逆FFT操作，完成对该信号的重采样滤波。

逆FFT操作的具体步骤为：

(6a)FPGA内通过计数器cnt_ifft产生IFFT变换过程的起始和结束使能信号：当cnt_ifft＝1时，产生一个宽度为一个时钟周期的ifft_sop高电平信号，作为IFFT变换输入数据的开始信号；当cnt_ifft＝N₂时，产生一个宽度为一个时钟周期的ifft_eop高电平信号，作为IFFT变换输入数据的结束信号；

(6b)将inverse信号置为1'b1，表示FFT IP核要执行IFFT变换；将sink_error信号置为2’b00，表示无误差输入输出；将source_ready信号置为1’b0表示输出准备完毕，可以输出；

(6c)根据FFT IP核输出的指数项对输出的实虚部进行数据截位，保证最终数据为16位有效值，得到一组N₂点实、虚部输出数据以及匹配的控制使能信号，截位方法与步骤(5b)相同。

此时，由步骤6得到的N₂个实、虚部数据即为主辅通道间的幅相误差估计的实、虚数据。

下面结合仿真实验实例对本发明的效果做进一步的说明。

1、仿真条件：

参照下表的仿真参数，利用MATLAB产生两个通道的含有幅相误差特性的数据。假设主通道为无幅相误差特性的通道，即主通道输入为理想线性调频信号；辅通道为含有幅相误差特性的通道，即辅通道输入为含有幅相误差特性的线性调频函数信号；在频域添加的幅相误差参数如表1所述，幅相误差添加的噪声为高斯白噪声，其幅相特性如附图2所示。

表1

参数	符号	数值	单位
				带宽	Br	100	M
发射脉冲宽度	Tp	10	μs
				距离采样率	Fr	120	MHz
距离调频率	Kr	10000	GHz/s
				误差幅频主值		[-1.5～1.5]	dB
误差相频主值		[-2.5～2.5]	°

2、仿真过程

第一步，在进行仿真时，给辅通道的输入信号添加误差信号。

给辅通道输入信号在幅频上添加大小在[-1.5dB,+1.5dB]之间且包含有随机高斯白噪声的正弦误差信号，如图2(a)所示，其中，横坐标轴表示采样点数，纵坐标轴表示添加的幅频特性大小，单位为dB，曲线值的大小，反映了误差相对于该通道内信号频谱幅度的强弱；

给辅通道输入信号在相频上添加大小在[-2.5°,+2.5°]之间且包含有随机高斯白噪声的线性相位误差信号，如图2(b)所示，其中，横坐标轴表示采样点数，纵坐标为相位值，单位为°，曲线值的大小，反映了误差相对于该通道内信号频谱相位的变化大小。

第二步，输入主、辅通道信号。

主通道的输入信号为理想的线性调频信号，其频谱如图3(a)所示；

辅通道输入信号为含有幅相误差的线性调频信号，其频谱如图3(b)所示。

第三步，频谱比对。

对主辅通道输入信号的频谱进行比对，频谱比对结果的幅频误差特性估计如图4(a)所示，图4(a)的横坐标轴表示误差估计时使用的FFT点数，纵坐标轴表示误差估计的幅频特性曲线的幅度值；频谱比对结果的相频误差特性估计图4(b)所示，图4(b)的横坐标轴表示估计时使用的FFT点数，纵坐标轴表示估计的相频特性曲线的相位值；频谱比对结果用软件Quartus Ⅱ中的SignalTap进行实时采样，结果如图4(c)所示。

第四步，计算频谱比的频谱。

对频谱比对结果进行FFT变换，获得频谱比对结果的频谱，并使用fftshift函数将零频移至图像中心，如图5(a)所示，并对图5(a)的零频附近进行局部放大，结果如图5(b)所示。

第五步，加窗处理。

以图5(a)所示的零频位置为中心，向左右两个方向各保留长度为32个单元的数据，作为加窗处理结果。

加窗后的幅频特性的估计结果如图6(a)所示，图6(a)的横坐标表示估计时使用的FFT点数，纵坐标表示估计的幅度值；

加窗后相频特性的估计结果如图6(b)所示，图6(b)的横坐标表示估计时使用的FFT点数，纵坐标表示估计的相位。

加窗后幅相特性估计结果的实、虚部数据的实时采样结果如图6(c)所示。

三.仿真结果分析：

对比图2和图4，附图4估计的幅相误差特性具有较大的毛刺现象，尤其是在零频位置附近，频谱比对得到的结果含有较大的误差，不能准确的用于通道幅相误差校正。因此有必要对频谱比对结果做进一步的处理，使得估计结果能够更加准确地表示通道间的幅相误差特性。

对比图4和图6，图6估计的结果毛刺现象明显减少，各种杂散和噪声相对较少，曲线相对比较平滑，估计结果明显优于图4的结果，从而表明，本发明对通道间幅相误差估计具有很好的实用性和高效性。

综合上述分析，本发明在处理通道间幅相误差特性估计时，得到的结果要明显优于频谱比对的结果，频域加窗处理对估计结果的平滑性起到了很大的帮助，并且该方法的实时处理性能很高。

Claims (3)

1.基于FPGA的InSAR通道幅相误差估计方法，包括如下步骤：

(1)输入通道数据：把双通道雷达实时信号处理系统中的一个通道作为主通道，另一个通道作为辅通道，并将主通道数据作为参考通道数据，辅通道数据作为待校准通道数据；

(2)分别对主、辅通道数据进行N₂点的FFT变换，得到主通道数据的频域形式为M(f)、辅通道数据的频域形式为S(f)；其中nextpow2(N₁)为在大于等于N₁的二的整数次方中最小的数的幂；N₁为输入通道数据的复数据点数；

(3)根据步骤(2)的结果，获取主通道与辅通道之间的幅相误差估计值： $H_1\left(f\right)=\frac{M\left(f\right)}{S\left(f\right)};$

(4)对主、辅通道间的幅相误差估计值H₁(f)进行FFT变换，得到主、辅通道间的频域幅相误差估计值H₂(f)；

(5)在频域幅相误差估计值H₂(f)的零频附近进行长度为64个单元的矩形窗截断，获取截断后的频谱并对频谱中的数据进行截位；

(6)对截断后的频谱进行IFFT变换，得到最终加窗后的幅相误差估计结果

2.根据权利要求1所述的基于FPGA的InSAR通道幅相误差估计方法，其特征在于：所述步骤(2)中对通道数据进行N₂点的FFT变换，按如下步骤进行：

(2a)判断输入信号点数N₁是否满足N₁＝N₂，其中如果不满足则执行步骤(2b)，如果满足，则执行步骤(2c)。

(2b)在输入数据的末端补零至N₂点，并产生该N₂点数据的匹配使能信号；

(2c)使用计数器产生FFT变换的起始和结束使能信号；

(2d)将FFT IP核的inverse信号置为1'b0，表示FFT IP核要执行FFT变换；将sink_error信号置为2’b00，表示无误差输入输出；将source_ready信号置为1’b0表示输出准备完毕，输出相应的频域数据。

3.根据权利要求1所述的基于FPGA的InSAR通道幅相误差估计方法，其中所述步骤(5)中对频谱中的数据进行截位，通过如下公式进行：

result＝data×2^-exp/2ⁿ，

其中，result为截位后的值，data为FFT IP核输出结果中数据的实部或虚部，exp为FFT IP核输出结果中的指数项数据，n为保证result值的最大模值不超过2¹⁵-1的一个量值，该公式的运算过程可理解为：当data为正数时，在右移(exp+n)位时高位补0；当data为负数时，在右移(exp+n)位时高位补1。