几年前，我使用自制的FMCW雷达的第二版进行了一些简单的合成孔径雷达（SAR）成像实验。从那时起，我制作了雷达的第三个改进版本，但由于需要大量的工作，因此未进行任何SAR测量。之后我确实有计划进行一些SAR实验，但是直到现在为止，它才有足够的时间和动力。
合成孔径雷达（SAR）成像是一种通过在已知路径上移动单个天线来合成非常大的天线阵列的方法。如果场景中没有移动目标，则一架雷达沿一条路径进行多次测量的结果与一头与移动路径一样长的可笑的大型雷达的结果相同。
单个目标的SAR成像。随着雷达的移动，测得的距离跟随抛物线。
如果我们在一条直线路径上移动，而雷达指向的方向与该路径的方向成90度，则该距离将测量到单个目标，因此，我们发现测得的距离遵循抛物线。这直接来自勾股定理。 SAR成像问题是从测得的距离数据中找出目标位置。当然，在真实场景中，我们有多个目标，解决方案并不像上图中那样寻找最接近的方法在哪里。
解决此问题的算法很少，但我要使用的算法称为Omega-k算法。这是一种利用FFT的快速成像算法，它也使得在GPU上高效计算成为可能。该推导主要遵循郭和董的论文。
我拥有的雷达是调频恒定波（FMCW）雷达。它传输短的频率扫描。传输波形可以建模为：
$$ s_t（\ tau）= \ exp（j 2 \ pi f_c \ tau + \ pi \ gamma \ tau ^ 2），\ quad -T_s / 2 <\ tau ，其中\（j = \ sqrt {-1} \），\（f_c = \）RF载波频率，\（\ tau = \）时间变量，\（\ gamma = B / T_s = \）扫描带宽/扫描长度\（= \）扫描速率。 所发射的波在一定距离处反射离开目标，并在时间\（t_d \）之后被接收。忽略振幅，接收到的波是发射信号的延时复制：\（s_r（\ tau）= s_t（\ tau-t_d）\）。来自多个目标的信号相加。 接收器将接收的信号与发送的信号混合。这种混合称为去杂散，可以去除高频RF分量。结果是一个低频信号，通常在几千赫兹到兆赫兹之间，并且易于通过低成本ADC进行数字化。对于复数信号，我们对传输的信号进行复数共轭以获得低通积，得到的混合积为： $$ s _ {\ text {IF}}（\ tau）= s_t（\ tau-t_d）s_t ^ *（\ tau）= \ exp（-j 2 \ pi f_c t_d- j 2 \ pi \ gamma t_d \ tau + j \ pi \ gamma \ tau ^ 2）$$ 在SAR测量期间，雷达以恒定速度沿直线路径重复进行此测量。雷达在路径上的位置为：\（x = v \ tau + x_n \），其中\（v \）是雷达平台的速度，\（x_n = v n T_p \）。 \（n \）是测量的索引，\（T_p \）是发射重复间隔。 如果雷达目标在位置\（（（x_0，y_0）\））到目标的距离可以写成： 我们将路径的y坐标设置为0，x位置限制为\（-L / 2 由于电磁波以光速传播并且雷达信号需要传播到目标并返回到雷达，因此我们获得了接收信号时间延迟\（t_d = 2R（x）/ c \）的表达式，其中\（c \）为光的速度。 $$ s _ {\ text {IF}}（\ tau，x} = \ exp \ left（-j \ frac {4 \ pi} {c}（f_c + \ gamma \ tau）R（x）\ right）\ exp \ left（j \ frac {4 \ pi \ gamma ^ 2} {c ^ 2} R ^ 2（x）\ right）$$ 上式中的最后一项称为剩余视频相位项，它是去皮操作的不良副产品。应该先乘以\（\ exp（-j \ frac {4 \ pi \ gamma ^ 2} {c ^ 2} R ^ 2（x））\），然后将其删除。但是，这种形式不方便，因为它取决于\（R（x）\）。利用\（R（x）= c t_d / 2 \）和\（t_d \）可以用IF信号的频率表示的事实：\（f = -2 \ gamma R（x）/ c =- \ gamma t_d \ Rightarrowt_d =-\ frac {f} {\ gamma} \）我们可以将校正项写为\（\ exp（-j \ pi f ^ 2 / \ gamma）\）。这种形式可以容易地应用于傅立叶变换的信号。 $$ s（\ tau，x_n）= \ exp \ left（-j \ frac {4 \ pi} {c}（f_c + \ gamma \ tau）\ sqrt {y_0 ^ 2 +（x_n-x_0 + v \ tau ）^ 2} \ right）$$ 理想情况下，我们希望信号的形式为\（\ exp（-j 2 \ pi f_y y_0）\ exp（-j2 \ pi f_x x_0）\），然后我们可以应用二维傅立叶逆变换以获得以中心为中心的增量函数在\（（x_0，y_0）\）处聚焦图像。目前，信号\（s（\ tau，x_n）\）不在这种形式下，傅立叶逆变换没有任何有趣的意义。我们需要找到一些处理步骤以应用于信号以使其达到所需的形式，以便可以应用傅立叶逆变换。之所以专门考虑这种形式，是因为FFT可以非常有效地执行。 首先，请注意\（\ gamma \）的单位为Hz / s，\（\ tau \）的单位为s。乘积\（\ gamma \ tau \）的单位为Hz，因此它是一个频率。该乘积实际上是瞬时扫描的调制频率。我们做替换\（\ gamma \ tau \ rightarrow f_ \ tau \）来摆脱时间变量。 \（\ tau \）范围是\（-T / 2 \ ldots T / 2 \），新的\（f_ \ tau \）范围是\（-B / 2 \ ldots B / 2 \）。 $$ S（f_ \ tau，x_n）= \ exp \ left（-j \ frac {4 \ pi} {c}（f_c + f_ \ tau）\ sqrt {y_0 ^ 2 +（x_n-x_0 + \ frac { v f_ \ tau} {\ gamma}）^ 2} \ right）$$ 同样，代替使用频率，数学更简洁，而当使用波数时，算法的实现更容易。我们定义范围波数\（K_r = K_ {rc} + \ Delta K_r \）。 \（K_ {rc} = \ frac {4 \ pi f_c} {c} \），\（\ DeltaK_r = \ frac {4 \ pi f_ \ tau} {c} =-\ frac {2 \ pi B} { c} \ ldots \ frac {2 \ pi B} {c} \）。 $$ S（K_r，x_n）= \ exp \ left（-j K_r \ sqrt {y_0 ^ 2 +（x_n-x_0 + \ frac {vc \ Delta K_r} {4 \ pi \ gamma}）^ 2} \ right ）$$ 下一步是在方位角方向（运动方向）上进行傅立叶变换，以将\（x_n \）变量也移至频域。 $$ S（K_r，K_x）= \ int _ {-\ infty} ^ \ infty S（K_r，x_n）\ exp（-j K_x x_n）\，dx_n = \ int _ {-\ infty} ^ \ infty \ exp（ j \ Phi（x_n））\，dx_n $$ \（K_x = 2 \ pi f_x \）是方位角方向的波数。这个积分没有精确的解决方案，但是有一种方法可以使用称为固定相位原理（PSOP）的方法来计算非常精确的近似值。集成功能的阶段可以写成： $$ \ Phi（x_n）= -K_r \ sqrt {y_0 ^ 2 + \ left（x_n-x_0 + \ frac {v c \ Delta K_r} {4 \ pi \ gamma} \ right）^ 2}-K_x x_n $$ 如果我们为某些实际值绘制相位\（\ Phi（x_n）\），我们将得到一个如下图所示的图： 在一个点上，相位的导数为零（平稳点），并且函数变化缓慢，但远离该点，函数具有很高的振荡性。当我们对函数进行积分时，远离固定点的振荡会被抵消，并且主要是固定点周围的区域会贡献积分结果。 我们可以围绕固定点\（\ frac {d} {dx_n} \ Phi（x_n）\ rvert_ {x_n = x_n ^ \ star} = 0 \）展开函数，如\（\ Phi（x_n）= \ Phi （x_n ^ \ star）+ 0 + \ frac {1} {2} \ Phi ^ {''}（x_n-x_n ^ \ star）^ 2 \）。 $$ \ begin {aligned} S（K_r，K_x）＆\ approx \ exp（j \ Phi（x_n ^ \ star））\ int _ {-\ infty} ^ \ infty \ exp \ left（j \ frac {1} {2} \ Phi ^ {''}（x_n ^ \ star）（x_n-x_n ^ \ star）^ 2 \ right）\，d x_n \\＆= \ exp（j \ Phi（x_n ^ \ star）） \ int _ {-\ infty} ^ \ infty \ exp \ left（j \ frac {1} {2} \ Phi ^ {''}（x_n ^ \ star）s ^ 2 \ right）\，ds \\＆= \ exp（j \ Phi（x_n ^ \ star））\ sqrt {\ frac {2 \ pi j} {\ Phi ^ {''}（x_n ^ \ star）}} \ end {aligned} $$ 因为\（\ Phi（x_n）\）是纯实函数，所以如果\（\ mu \）是\（\ Phi（x_n ^ \ star）\）的符号，则平方根项可以写为\（ \ sqrt {\ frac {2 \ pi} {| \ Phi ^ {''}（x_n ^ \ star）|}} exp（j \ pi \ mu / 4）\）。二阶导数贡献振幅项和恒定相位项，它们都不对聚焦图像重要，而聚焦图像主要取决于对准相位。自开始以来，我们就忽略了振幅，而振幅最终是缓慢变化的函数，因此我们将其近似化。 将固定点插入上面的\（S（K_r，K_x）\）方程，我们得到积分的解： $$ S（K_r，K_x）\约\ exp \ left（j（-y_0 \ sqrt {K_r ^ 2-K_x ^ 2}-K_x x_0 + \ frac {c \ Delta K_r K_xv} {4 \ pi \ gamma} ）\ right）$$ 最后一项是由于扫描过程中雷达的运动引起的相位偏移。可以通过在相反的相位乘以指数将其删除。 \（x_0 \）项已经采用了正确的格式，因为它仅与\（K_x \）相乘，但是\（y_0 \）项同时取决于\（K_r \）和\（K_x \）。可以通过替换\（\ sqrt {K_r ^ 2-K_x ^ 2} \ rightarrow K_y \）来固定\（K_r，K_x \）依赖性。此步骤称为Stolt插值，因为它是通过将数据插值到新网格来实现的。 进行2D逆傅立叶变换可得出聚焦图像，其delta函数的中心为\（（x_0，y_0）\）。 Omega-k算法主要是大型FFT和插值。两者都可以在GPU上很好地实现，而GPU需要程序具有很大的并行性。精心编写的GPU实现应比CPU实现快几倍。为了方便起见，我将使用Tensorflowlibrary实现该算法。尽管它最常用于训练神经网络，但也可以用于其他目的。 以上推导是使用连续信号完成的，但实际上，雷达使用ADC对信号进行采样，从而得到离散采样。上面的大多数推导仍然有效，例如，在确保采样网格足够小以避免混叠时，还需要一些其他考虑。 import numpy as np import tensorflow as tf＃从磁盘加载捕获的数据和参数。 data，settings = load_data（...）＃'data'包含2D数组中捕获的数据。 ＃第一维是路径上扫描的索引，第二维是来自ADC的扫描的原始值。 ＃测量期间平台的移动速度。 v =设置['v']＃数字化信号的采样率。 fs =设置['fs']＃扫描长度。 tsweep =设置['tsweep']＃扫描的带宽。 bw =设置['bw']＃扫描的RF中心频率。 fc =设置['f0'] + bw / 2＃两次扫描之间的时间。 tdelay = settings ['tdelay']＃扫描速率。 gamma = bw / tsweep＃捕获的扫描数。 sweep_samples = len（data [0]）＃捕获的扫描之间的位置差。 delta_x =（tsweep + tdelay）* v＃波数轴kx = np。 linspace（-np.pi / delta_x，np.pi / delta_x，len（数据））dkr = np。 linspace（（4 * np .pi / c）*（-bw / 2），（4 * np .pi / c）*（bw / 2），sweep_samples）kr =（4 * np .pi / c）* fc + dkr ky0 =（kr [0] ** 2-kx [0] ** 2）** 0.5 ky_delta = kr [1]-kr [0]＃与kr的间距相同，以避免插值时出现混叠。 ＃插补后的Ky轴。 ky_interp = np。范围（ky0，kr [-1]，ky_delta） 实际数据与推导之间的直接区别是，myradar没有IQ采样，并且捕获的信号是纯真的。但是，我们可以轻松生成所需的虚部。如果我们对捕获的信号进行FFT，则由于该信号是纯实数的正和负频率分量是复共轭。但是，复信号应仅具有负频率分量（负，因为接收混频器中的延迟RF信号的频率低于LO信号的频率）。如果我们将正分量归零，然后进行逆FFT，则结果是具有正确属性的复杂信号。这种转换称为希尔伯特转换。我们还可以在范围方向上应用窗口功能，并同时进行RVP项乘法。我们将此步骤作为使用numpy的预处理步骤： def hilbert_rvp（x，fs，gamma）：＃最后一个维度是范围维度y = np。英尺fft（x，轴=-1）y [：，：y。 shape [1] // 2 + 1] = 0＃正频率为零＃剩余视频相位项补偿f = np。 linspace（-fs / 2，fs / 2，y。shape [1]）y * = np。 exp（1j * np。pi * f ** 2 / gamma）返回np。英尺ifft（y，axis =-1）＃范围方向的开窗功能。 ＃减少FFT的旁瓣。 w = np。 hanning（data。shape [1]）＃fs = ADC的采样率。 ＃gamma =扫描带宽/扫描长度。数据= hilbert_rvp（w *数据，fs，伽玛） 下一步将是方位角FFT，但是在此之前，我们将首先在方位角方向对数据进行零填充，因为目标方位角位置可能在移动路径的端点之外。没有零填充，这些目标将别名到路径内的位置。 ＃在方位轴的两侧对称添加“ zpad”零。数据= np。 pad（data，（（（zpad // 2，zpad // 2），（0，0）），'常量'） 现在预处理已经完成，其余的都在GPU上完成了。第一步是方位角FFT，由于功能的限制，需要以稍微绕行的方式进行。 ＃从预处理数据创建张量。 img = tf。常量（data，dtype = tf。complex64）＃Tensorflow FFT没有选择要转换的轴的选项＃并且它总是在最后一个轴上计算FFT。 ＃换位以交换FFT之前和之后的轴，以计算＃在第一个轴上的FFT。 img = tf。转置（img）img = tf。信号。 fft（img）img = tf。 transpose（img）＃移位频率分量，以使零频率为中心。 img1，img2 = tf。拆分（img，[img。shape [0]。值-img。shape [0]。value // 2，img。shape [0]。value // 2]，轴= 0）img = tf。 concat（[img2，img1]，轴= 0） ＃匹配过滤以补偿扫描期间的移动。 mf = tf。 exp（j * tf。cast（tf。expand_dims（dkr，0）* tf。expand_dims（kx，-1），tf。complex64）* tf。cast（c * v /（4 * np。pi * gamma）， tf。complex64））img = img * mf 经过匹配滤波后，我们应该进行Stolt插值。数据当前是在\（K_x，K_r \）轴上定义的，我们需要将其插值到\（K_x，K_y \）轴。通常，新的\（K_y \）轴点并不完全对应于现有网格上的点，因此我们需要进行插值。问题在于在Tensorflow中没有简单的方法进行插值。线性插值之类的简单插值方法虽然易于实现，但由于它们会引起频域失真，因此效果不佳。理想的插值将是Sinc插值，但是该公式需要对信号中的每个样本相乘以计算一个输出点，从而导致非常慢的\（O（n ^ 2）\）算法。 Lanczosinterpolation是高效算法和最小频域失真之间的一个很好的折衷方案。而不是使用具有无限支持的\（\ text {sinc}（x）\）进行插值，而是使用具有无限支持的内核\（L（x）\），以便在插值中仅需考虑附近的样本： $$ L（x）= \ begin {cases} \ text {sinc}（x）\，\ text {sinc}（x / a）＆\ text {if} \ -a 据我所知，实际上没有任何简单的方法可以在Tensorflow中使用现有操作有效地实现它。对图像中的每个点应用公式都会给计算图增加太多的操作，以至于它永远无法完成。我最终用C ++为它编写了具有CPU和GPU实现的自定义操作。它的代码太长，无法在此处包含，但您可以在Gitrepository中找到它。编译完该操作并编写一些Python代码以与其接口之后，Stolt插值可以单行编写： from interp_op import stolt_interp＃定义数据的原始网格。 ky = tf。常量（（kr ** 2-（kx ** 2）[:, np。newaxis]）** 0.5，dtype = tf。float32）＃Stolt插值＃'interp_order'是Lanczos内核顺序'a'。 img = stolt_interp（ky，img，ky_interp，interp_order） 我将雷达安装在自行车的后架上，并以恒定速度将自行车踩在一条直线上以进行测量。知道路径非常重要，因为实际位置和图像形成过程中使用的位置之间的任何差异都会导致图像质量下降。位置应该在波长的一小部分之内，以避免任何散焦误差。我的雷达工作在6 GHz，在大约200 m长的路径上达到几厘米的精度要求。这可能不会发生，但是我们稍后将看到如何至少部分纠正引入的位置误差。 正在成像的场景的照片。成像路径位于右侧的人行道上，雷达指向左侧的停车场。 上图中是正在成像的场景。这是一个停车场，上面有许多类似杆子的目标，应该在生成的图像中清晰可见。 原始捕获的数据看起来并不多。扫描长度为1 ms，采样率为1 MHz，因此每个扫描有1000个点，总共有6444个扫描。FMCW雷达的输出信号是来自每个可见目标的正弦波的叠加。频率取决于到目标的距离，距离越近，频率越低。幅度取决于反射的功率大小，取决于大小，形状，材料和到目标的距离。 在范围方向上进行FFT可以使其更容易读取X轴范围内的格式。 上面是原始数据在范围方向上的FFT。范围FFT并不是图像处理的一部分，但这是处理非成像测量的方式.FFT允许将X轴从波数更改为到目标的距离。 生成的停车场的SAR图像。上面的相机照片是在（90，0）朝（0,0）方向拍摄的。 图像主要是聚焦的，有一些可见的散布，但至少可以识别大多数目标。我不确定有些弯曲的工件来自何处。 放大前景对象时，可以看到运动偏差带来的一些拖影。这些可以通过自动聚焦算法解决。 使用来自停车场的以上数据，在不同程序上的图像形成时间。不包括程序启动和数据预处理时间。 利用GPU绝对有巨大的速度。 Numpy实现使用FFT进行FFT，Numpy和Stolt插值是在没有向量化的情况下在Python中编码的.Numpy版本需要22分钟30 s的时间来形成上述数据的图像，其中大部分数据都用在插值例程中。 TensorFlow CPU实现可计算CPU上的所有内容，Stolt插值是用C ++编写的自定义操作代码。它比Numpy实现快大约16倍，而大多数速度来自更快的C ++插值函数。图像形成需要1分钟20 s。 TensorFlow GPU Python代码与CPU版本完全相同。它在GPU上使用nVidia的cufft库计算FFT，并且在GPU上也使用自定义CUDA内核进行插值。图像形成仅需80毫秒。加速绝对是巨大的，比TensorFlow CPU快1000倍以上，比Numpy版本快一千倍以上。 GPU的加速实际上并不需要那么多。我认为部分原因是使用nVidia库的GPU实现比TensorFlow使用的任何CPU实现都要优化得多。 CPU FFT似乎不是多线程的，因此只有多线程FFT才能以较小的速度加快它的速度。 在自行车上移动雷达时，不可避免地会出现一些运动错误。如果有轻微 ......