CN101545976B - 反射层析激光雷达双折线反投影成像方法 - Google Patents

Abstract

一种反射层析激光雷达双折线反投影成像方法，实质是将反投影波阵面近似为双折线的方法，将激光雷达探测器所接收到的球面反射投影信号反投影到所对应的双折线上，通过这些双折线信号重建图像。本发明方法与球面直接反投影相比，双折线反投影算法计算量小，操作简单，易于实现。模拟仿真表明本发明方法可以替代球面反投影完成相应质量要求的图像重建。

Description

反射层析激光雷达双折线反投影成像方法

技术领域

本发明涉及反射层析激光雷达，特别是一种反射层析激光雷达双折线反投影成像方法，通过把反投影波面近似为双折线，将所测得的球面反射投影信号反投影到对应的双折线上，用双折线反投影信号完成相应距离和质量要求的图像重建。

背景技术

在医学计算机辅助层析(简称为CT)中，根据一维透射投影信息实现两维图像重建已经得到广泛的应用和发展。所得到的一维透射投影信息反应的是在多个不同角度下，物体内部不同部位透射系数的差异。反投影成像的内容为重建图像某一点的光强密度值可看作这一平面内所有经过该点的射线投影值之和。

如图1透射层析所示，f(x，y)为代建图像，L_r，φ为对应的直线r＝xcosφ+ysinφ，p(r，φ)为f(x，y)沿直线L_r，φ的透射特征系数的积分，即角度φ所对应的一维透射投影：

$p(r,\mathrm{\φ})={\∫}_{L_{r,\mathrm{\φ}}}f(x,y)\mathrm{ds}---\left(1\right)$

用反投影成像方法重建图像g(x，y)，则：

$g(x,y)=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}p(x\cos {\mathrm{\φ}}_i+y\sin {\mathrm{\φ}}_i,{\mathrm{\φ}}_i)\mathrm{\Δ\φ}---\left(2\right)$

其中，φ_i为第i次投影所对应的角度，Δφ为投影角度采样间隔，m为总的投影个数。CT反投影成像把取自有限物体空间的射线均匀地反投影到射线所及的无限空间的各点之上，使原先像素数值为零的点经反投影之后不再为零，产生星状伪迹。为了减少重建图像星状伪迹，引入了滤波反投影算法，对获得的投影进行适度的调整，调整之后的投影数据为：

$q(r,\mathrm{\φ})=F_1^{-1}\left(\right|k\left|F_1\left(p(r,\mathrm{\φ})\right)\right)---\left(3\right)$

其中，符号F₁、F₁ ^-1分别表示一维的傅里叶变换与傅里叶反变换，r为空间变量，k为频率变量，则经过滤波反投影重建图像g_FB(x，y)表达式如下：

$g_{\mathrm{FB}}(x,y)=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}q(x\cos {\mathrm{\φ}}_i+y\sin {\mathrm{\φ}}_i,{\mathrm{\φ}}_i)\mathrm{\Δ\φ}---\left(4\right)$

反射层析(基于距离)类似于上述CT滤波反投影重建图像的操作过程，其实质是将到源点某一距离的有限表面区域反射所得的投影信息，反投影到隔源点相同距离的所有区域，如图2反射层析激光雷达所示。所不同的是CT依据的透射特征系数，只需要0-π的透射投影信息，而反射层析依据的是反射特征系数，需要0-2π的反射投影信息。

随着激光的出现和发展，Charles L.Matson(参见文献R.M.Marino，R.N.Capes，W.E.Keicher et al..Tomographic image reconstruction from laser radarreflective projections[J]，SPIE Laser Radar III，1988，Vol.999：248-263)，Parker.J.K.(参见文献Jeffrey K.Parker，E.B.Cralg，D.I.Klick et al..Reflective tomography：image from range resolved laser radarmearsurements[J]，Applied Optics，1988，27(13)：2642-2643)，F.K.Knight(参见文献Frederick L.Knight，David Klick，Danette P.Ryan-Howard et al..Laserradar reflective tomography utilizing a steak camera for precise rangeresolution[J].Applied Optics，1989，28(12)：2196-2198)等很快的将反射层析算法引入到激光雷达中。根据激光短脉冲探测和非相干接收器得到基于距离的一维反射投影信息，激光雷达系统可以实现二维图像重建，具有很高的分辨率。当波阵面发生弯曲时，如图3反射层析激光雷达所示，直接的球面反投影滤波方式计算量大，需要的处理时间长。引入双折线的处理方法，在一定距离下，可以减少计算量，实现图像的重建。

发明内容

本发明的目的在于为反射层析激光成像雷达提供一种反射层析激光雷达双折线反投影成像方法，该方法的原理可靠，结构简单，易于实现。

本发明的具体技术解决方案如下：

一种反射层析激光雷达双折线反投影成像方法，其特点是包括以下步骤：

①由探测物体的目标尺寸x_max、以探测器和光脉冲二者中较低的分辨率作为系统的距离分辨率ΔR、激光脉冲宽度τ三个因素，根据奈奎斯特采样定律计算满足物体完全重建所需要的最小投影角度采样间隔Δφ和单个投影最少采样点数N_p：

$\mathrm{\Δ\φ}\≤\frac{360}{N_{\mathrm{\ψ}}}=\frac{180\mathrm{\ΔR}}{\mathrm{\π}{x}_{\max }},$

$N_p\≥{2x}_{\max }{u}_{\max }=2\frac{x_{\max }}{\mathrm{\ΔR}};$

②激光器发出光脉冲，用会聚透镜使光脉冲会聚后再次发散，通过调节透镜的焦距来改变光斑的大小，使到达探测目标的光斑外轮廓涵盖整个探测目标，获得需要光束最小发散角α＝arcsin(x_max/2L)，也就是所对应的双折线反投影的角度，其中L是激光器与待测目标之间的距离根据激光器发出光脉冲的时刻t₁，最早接收到反射信息的时刻t₂，最晚接收到反射信息的时刻t₃，探测器到待测目标的距离L为

其中c为光速；

③将初始待测目标所对应的角度定义为φ₀＝0°，激光器发出单脉冲照射到待测目标上，探测器记录在角度φ₀＝0°下所对应的经过目标反射加宽后的反射投影信息，即探测器所得到的光功率随时间的变化信息反映的是物体表面反射特征系数随物体深度的变化信息：

p₀(t，φ₀)＝e∫f₀(x，y)ds

其中，f₀(x，y)为物体表面的反射特征系数，e代表到达物体表面的均匀光强，p₀(t，φ₀)为角度φ₀下获得的反射投影信息，ds为物体表面积分面元；

④待测目标旋转角度Δφ后，探测角度变为φ₁＝φ₀+Δφ，激光器再次发出单个脉冲，探测器记录相应的反射投影信息p₁(t，φ₁)；

⑤判断投影角度φ₁是否大于360°，若小于360°，重复步骤④(i-1)次操作，待测物体旋转iΔφ，i为2以上的正整数，相应记录反射投影信息p_i(t，φ_i)，直至φ_i大于360°；

⑥滤波处理：利用公式r＝ct，c为光速，将所得到的随时间变化的投影信息p_i(t，φ_i)，转换成随距离变化的投影信息，对每个角度获得投影信息进行滤波处理，减少重建图像中的伪迹，具体滤波公式为： $q_i(r,{\mathrm{\φ}}_i)=F_1^{-1}\left(\right|u\left|F_1\left(p_i(r,{\mathrm{\φ}}_i)\right)\right),$ 其中u为频域变量，符号F₁、F₁ ^-1分别表示信号一维傅里叶变换和傅里叶反变换；

⑦建立重建图像矩阵g_FB1(x，y)，其像元坐标间隔为所对应的系统距离分辨率ΔR，为获得完整的目标重建图像，图像矩阵g_FB1(x，y)行数、列数要大于

并将矩阵所有像元初始值设为零；

⑧取图像矩阵g_FB1(x，y)第一个像元所对应的点坐标(x₁，y₁)，按照角度φ_i(i＝0，1，2…N_φ)进行旋转，旋转之后的新坐标变为：

x′₁＝x₁cos(φ_i)+y₁sin(φ_i)，y′₁＝-x₁sin(φ_i)+y₁cos(φ_i)判断新坐标点(x′₁，y′₁)到角度φ_i、r_j所对应的双折线的距离是否满足小于或等于像素间隔ΔR的

倍，若满足以上条件，则图像矩阵元g_FB1(x₁，y₁)累加上反射投影值q_i(r_j，φ_i)Δφ，直至像元点坐标(x₁，y₁)完成所有旋转角度φ_i(i＝0，1，2…N_φ)和距离r(i＝0，1，2…N_p)的计算。

具体重建公式如下：

新坐标点(x′₁，y′₁)到角度φ_i、r_j所对应的双折线的距离d：

$d=\frac{|{-x}^{\′}\sin \mathrm{\α}+|y^{\′}|(1-\cos \mathrm{\α})+(L-r_j)\sin \mathrm{\α}|}{\sqrt{2-2\cos \mathrm{\α}}}$

若 $d\≤\frac{\sqrt{2}}{2}\mathrm{\ΔR},$ 则 $g_{\mathrm{FB}1}(x_1,y_1)=\underset{i=1}{\overset{N_{\mathrm{\φ}}}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_p}{\mathrm{\Σ}}}{q}_i(r_j,{\mathrm{\φ}}_i)\mathrm{\Δ\φ}$

⑨将图像矩阵g_FB1(x，y)的所有其它像元，依次重复步骤⑧操作，最终将得到双折线滤波反投影重建图像g_FB1。

本发明的技术效果：

本文在充分考虑激光发散角，波阵面呈现球面的状态。为此我们对反射层析激光雷达成像采用双折线滤波反投影的方法，将球面反射投影信息反投影到双折线上，利用双折线信息实现图像的重建。改善了直接球面反投影计算量大、处理时间长的问题。

附图说明

图1透射层析示意图；

图2反射层析示意图；

图3反射层析波阵面弯曲示意图；

图4反射层析反射投影探测示意图；

图中：1激光器；2会聚透镜；3探测目标；4光纤耦合头；5光纤；6探测器。

图5原正三角形图样，R为物体所占区域半径；

图6双折线反投影原理示意图；

图7双折线反投影重建图像二值化；

具体实施方式

下面结合附图和实施例进一步详细说明本发明，但不应以此限制本发明的保护范围。

一种反射层析激光雷达双折线反投影成像方法，包括以下步骤：

(1)由探测物体的目标尺寸、探测器距离分辨率、激光脉冲宽度τ三个因素，根据奈奎斯特采样定律计算满足物体完全重建所需要的最小投影角度采样间隔Δφ、单个投影最少采样点数N_p；

在实际选择中，距离分辨率要以探测器和光脉冲二者中较低的分辨率作为系统距离分辨率，否则会造成信息的丢失。例如，所用探测器带宽为4GHz，采样率最大为20GHz的示波器，脉冲激光器发出的脉冲宽度为8ps，可以计算示波器的距离分辨率(7.5mm)低于光脉冲所对应的距离分辨率(1.2mm)，应以示波器分辨率作为系统距离分辨率。假设系统距离分辨率为ΔR，转换到频域： $\mathrm{FT}\left\{\mathrm{rect}\left(\frac{x}{\mathrm{\ΔR}}\right)\right\}=\sin c\left(\mathrm{\ΔRu}\right)$

其中，符号FT表示傅里叶变换，符号rect表示矩形函数，符号sinc表示sinc函数。定义sinc函数第一个零点为图像频域截止频率u_max：u_max＝1/ΔR。假设目标物体最大尺寸为x_max，应用标准的奈奎斯特采样定律，要想完全重建物体图像，傅里叶频域采样间隔要小于1/x_max，可以得到反射层析旋转角度采样数：

$N_{\mathrm{\φ}}\≥2\mathrm{\π}{x}_{\max }{u}_{\max }=2\mathrm{\π}\frac{x_{\max }}{\mathrm{\ΔR}},$ 投影角度采样间隔： $\mathrm{\Δ\φ}\≤\frac{360}{N_{\mathrm{\ψ}}}=\frac{180\mathrm{\ΔR}}{\mathrm{\π}{x}_{\max }}.$

每个投影上采样点数：

$N_p\≥{2x}_{\max }{u}_{\max }=2\frac{x_{\max }}{\mathrm{\ΔR}};$

(2)如图4所示，脉冲激光器发出光脉冲，用会聚透镜使光脉冲会聚后再次发散，可以通过调节透镜焦距来改变光斑的大小，使到达目标物体的光斑外轮廓涵盖整个探测目标；若物体距离脉冲激光器足够远，也可以不用透镜，光束的发散角导致远场光斑覆盖整个物体即可；假若探测目标最大尺寸为x_max，双折线反投影所对应的角度为α＝arcsin(x_max/2L)，L为探测器到物体的距离，根据激光器发出光脉冲的时刻t₁，最早接收到反射信息的时刻t₂，最晚接收到反射信息的时刻t₃，探测器到物体的距离L为

其中c为光速；

(3)将初始测量目标所对应的角度定义为φ₀＝0°，激光器发出单脉冲照射到目标物体上，探测器记录在角度φ₀＝0°下，所对应的经过目标反射加宽后的反射投影信息。在距离比较远的情况下，到达物体表面波面光强是均匀的，那么反射回来的光功率信息的变化则是物体表面的反射特征系数引起的，即探测器所得到的光功率随时间的变化信息反映的是物体表面反射特征系数随物体深度的变化信息：

p₀(t，φ₀)＝e∫f₀(x，y)ds(6)

其中，f₀(x，y)为物体表面的反射特征系数，e代表到达物体表面的均匀光强，p₀(t，φ₀)为角度φ₀下获得的反射投影信息，ds为物体表面积分面元；

(4)当目标物体发生旋转角度Δφ后，探测角度变为φ₀+Δφ，激光器再次发出单个脉冲，在角度φ₀+Δφ下探测器记录相应的反射投影信息p₁(t，φ₁)，并更新当前角度φ₁＝φ₀+Δφ；

(5)判断投影角度φ₁是否大于360°，若小于360°，多次重复步骤(4)操作，当物体旋转iΔφ(i＝1，2，3……)时重新记录所对应的反射投影信息p_i(t，φ_i)，并更新角度φ_i＝φ₀+iΔφ，直至φ_i大于360°；

(6)利用公式r＝ct(c为光速)，将所得到的随时间变化的投影信息，转换成随距离变化的投影信息。对每个角度获得投影信息进行滤波处理，减少重建图像中的伪迹，具体滤波公式为： $q_i(r,{\mathrm{\φ}}_i)=F_1^{-1}\left(\right|u\left|F_1\left(p_i(r,{\mathrm{\φ}}_i)\right)\right),$ 其中u为频域变量，符号F₁、F₁ ^-1分别表示信号一维傅里叶变换及傅里叶反变换；

(7)建立图像g_FB1(x，y)，其像元坐标间隔为所对应的系统距离分辨率ΔR。为获得完整的目标重建图像，图像矩阵g_FB1(x，y)行数、列数要大于

(符号表示取整数上限，例如

)，并将矩阵所有像元初始值设为零；

(8)取图像矩阵g_FB1(x，y)第一个像元所对应的点坐标(x₁，y₁)，按照角度φ_i(i＝0，1，2…N_φ)进行旋转，旋转之后的新坐标变为：

x′₁＝x₁cos(φ_i)+y₁sin(φ_i)，y′₁＝-x₁sin(φ_i)+y₁cos(φ_i)

判断新坐标点(x′₁，y′₁)到角度φ_i、r_j所对应的双折线的距离是否满足小于或等于像素间隔ΔR的

倍，若满足以上条件，则图像矩阵元g_FB1(x₁，y₁)累加上反射投影值q_i(r_j，φ_i)Δφ，直至像元点坐标(x₁，y₁)完成所有旋转角度φ_i(i＝0，1，2…N_φ)和距离r(i＝0，1，2…N_p)的计算。

具体重建公式如下：

新坐标点(x′₁，y′₁)到角度φ_i、r_j所对应的双折线的距离d：

$d=\frac{|{-x}^{\′}\sin \mathrm{\α}+|y^{\′}|(1-\cos \mathrm{\α})+(L-r_j)\sin \mathrm{\α}|}{\sqrt{2-2\cos \mathrm{\α}}}$

若 $d\≤\frac{\sqrt{2}}{2}\mathrm{\ΔR},$ 则 $g_{\mathrm{FB}1}(x_1,y_1)=\underset{i=1}{\overset{N_{\mathrm{\φ}}}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_p}{\mathrm{\Σ}}}{q}_i(r_j,{\mathrm{\φ}}_i)\mathrm{\Δ\φ}---\left(7\right)$

(9)将图像矩阵g_FB1(x，y)的所有其它像元，依次重复步骤(8)操作，最终将得到双折线滤波反投影重建图像g_FB1。

计算机仿真采用边长为3m的正三角形(如图5所示)，距离分辨率采样ΔR＝0.1m，则物体完全重建所需要的投影角度间隔Δψ≤3.3度，N_p≥34.64。模拟中Δψ＝1度，单投影采样点数N_p取满足采样条件的最小整数35，距离L＝10米。依次按照具体实施步骤进行计算机仿真，重建图像二值化结果如图7所示，可以证明本发明采用双折线滤波反投影的方法对反射层析激光雷达进行图像重建是一种行之有效的办法。计算量的减少可以直观的通过matlab运行时间来看，在同一外界条件下，仿真双折线反投影需要2分16秒，而球面法需要2分49秒。故双折线滤波反投影可以在减少计算量的前提下，实现相应距离、相应图像质量要求的图像重建。

Claims (1)

1.一种反射层析激光雷达双折线反投影成像方法，其特征在于包括以下步骤：

①由探测目标最大尺寸x_max、以探测器和光脉冲二者中较低的分辨率作为系统的距离分辨率ΔR、激光脉冲宽度τ三个因素，根据奈奎斯特采样定律计算满足探测目标完全重建所需要的最小探测角度采样间隔Δφ和单个投影最少采样点数N_p：

$\mathrm{\Δ\φ}\≤\frac{360}{N_{\mathrm{\ψ}}}=\frac{180\mathrm{\ΔR}}{{\mathrm{\πx}}_{\max }},$

$N_p\≥{2x}_{\max }{u}_{\max }=2\frac{x_{\max }}{\mathrm{\ΔR}};$

②激光器(1)发出光脉冲，用会聚透镜(2)使光脉冲会聚后再次发散，通过调节透镜的焦距来改变光斑的大小，使到达探测目标(3)的光斑外轮廓涵盖整个探测目标，获得需要光束最小发散角α＝arcsin(x_max/2L)，也就是所对应的双折线反投影的角度，其中L是探测器(6)与待测目标(3)之间的距离，根据激光器发出光脉冲的时刻t₁，最早接收到反射信息的时刻t₂，最晚接收到反射信息的时刻t₃，探测器(6)到探测目标(3)的距离L为

其中c为光速；

③将初始探测目标所对应的探测角度定义为φ₀＝0°，激光器(1)发出单脉冲照射到探测目标(3)上，探测器(6)记录在探测角度φ₀＝0°下所对应的经过探测目标反射加宽后的反射投影信息，即探测器所得到的光功率随时间的变化信息，反映的是探测目标表面反射特征系数随探测目标深度的变化信息：

p₀(t，φ₀)＝e∫f₀(x，y)ds

其中，f₀(x，y)为探测目标表面的反射特征系数，e代表到达探测目标表面的均匀光强，p₀(t，φ₀)为角度φ₀下获得的反射投影信息，ds为探测目标表面积分面元；

④探测目标(3)发生旋转角度的采样间隔Δφ后，探测角度变为φ₁＝φ₀+Δφ，激光器(1)再次发出单个脉冲，探测器(6)记录相应的反射投影信息p₁(t，φ₁)；

⑤判断探测角度φ₁是否大于360°，若小于360°，重复步骤④i-1次操作，探测目标旋转iΔφ，i为2以上的正整数，相应记录反射投影信息p_i(t，φ_i)，直至φ_i大于360°；

⑥滤波处理：利用公式r＝ct，c为光速，将所得到的随时间变化的投影信息p_i(t，φ_i)转换成随距离变化的投影信息，即距离反射投影p_i(r，φ_i)，其中r为距离变量，φ_i为第i次探测对应的探测角度，对每个探测角度获得投影信息进行滤波处理，减少重建图像中的伪迹，具体滤波公式为：

其中u为频域变量，q_i(r，φ_i)为滤波之后各个探测角度φ_i随距离r变化的反射投影信息，符号F₁、

分别表示信号一维傅里叶变换和傅里叶反变换；

⑦建立重建图像矩阵g_FB1(x，y)，其像元坐标间隔为所对应的系统距离分辨率ΔR，为获得完整的探测目标重建图像矩阵，该重建图像矩阵g_FB1(x，y)行数、列数要大于

并将重建图像矩阵所有像元初始值设为零；

⑧取所述的重建图像矩阵g_FB1(x，y)第一个像元所对应的点坐标(x₁，y₁)，按照探测角度φ_i(i＝0，1，2…N_φ)进行旋转，旋转之后的新坐标变为：

x′₁＝x₁cos(φ_i)+y₁sin(φ_i)，y′₁＝-x₁sin(φ_i)+y₁cos(φ_i)

判断新坐标点(x′₁，y′₁)到探测角度φ_i、距离r_j所对应的双折线的距离是否满足小于或等于像素间隔ΔR的

倍，若满足以上条件，则重建图像矩阵的矩阵元g_FB1(x₁，y₁)累加上距离反射投影值与旋转角度的采样间隔的乘积q_i(r_j，φ_i)Δφ，直至像元点坐标(x₁，y₁)完成所有旋转角度φ_i(i＝0，1，2…N_φ)和距离r_j(j＝0，1，2…N_p)的计算，其中q_i(r_j，φ_i)为探测目标第i次探测在探测角度φ_i、距离r_j所对应的滤波之后的反射投影值，

具体重建公式如下：

新坐标点(x′₁，y′₁)到角度φ_i、r_j所对应的双折线的距离d：

$d=\frac{|-x^{\′}\sin \mathrm{\α}+|y^{\′}|(1-\cos \mathrm{\α})+(L-r_j)\sin \mathrm{\α}|}{\sqrt{2-2\cos \mathrm{\α}}}$

若 $d\≤\frac{\sqrt{2}}{2}\mathrm{\ΔR},$ 则 $g_{\mathrm{FB}1}(x_1,y_1)=\underset{i=1}{\overset{N_{\mathrm{\φ}}}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_p}{\mathrm{\Σ}}}{q}_i(r_j,{\mathrm{\φ}}_i)\mathrm{\Δ\φ}$

⑨将重建图像矩阵g_FB1(x，y)的所有其它像元，依次重复步骤⑧操作，最终将得到双折线滤波反投影重建图像g_FB1。

Images