基于TR‑MUSIC算法的光学扫描全息轴向定位方法与流程
本发明涉及光学扫描全息技术与空间定位领域，具体涉及基于tr-music算法的光学扫描全息轴向定位方法。
背景技术：
光学扫描全息技术，简称osh，是数字全息技术的一个重要分支。它利用光学扫描技术将物体的3维信息切片储存为2维信息，从而得到物体的全息图。该技术是1979年，poon和korpel在研究声光外差图像处理器的时候所提出。自该技术提出以来，已经在扫描全息显微镜、3d图像识别以及3d光学遥感等领域得到了广泛的应用。
在光学扫描全息技术领域中，由于全息图的重建过程需要使用到扫描镜到物体的距离，因此获得精确的轴向扫描距离一直以来都是国内外的研究热点，至今为止，相继有使用双波长激光源、自聚焦等技术应用到光学扫描全息技术当中，希望提高轴向分辨率，但其精度依然有待提高。
文献《opticalscanningholographywithmatlab》讲述了光学扫描全息技术的相关理论，利用光学外差的方法将三维物体信息以二维的方式存储，并给出了通过物体的复振幅函数与光学传递函数作卷积运算得到物体的切片全息图的方法。
文献《timereversalopticaltomographylocatingtargetsinahighlyscatteringturbidmedium》将时间反演用到了光学领域，并总结了利用tr-music算法实现对高散射浑浊介质中物体的定位，该方法具有很高的定位精度与准确性。
文献《depthresolutionenhancementinopticalscanningholographywithadual-wavelengthlasersource》分析了光学扫描全息技术当中影响轴向分辨率的一些相关因素，同时提出了一种具有较高分辨率的重建算法，但其轴向定位精度仍有可提升空间。
技术实现要素：
本发明要解决的技术问题是：提供一种基于tr-music算法的光学扫描全息轴向定位方法，解决现有技术光学扫描全息技术扫描轴向距离分辨率低的问题。
为实现上述目的，本发明采用的技术方案如下：
基于tr-music算法的光学扫描全息轴向定位方法，包括以下步骤：
步骤1、将激光采用第一偏振分束器分成两束，之后一束光依次通过第一光瞳和第一凸透镜后投射至第二偏振分束器，第二束光依次通过第二光瞳和第二凸透镜后投射至第二偏振分束器，第二偏振分束器将投射来的两束光聚光干涉形成菲涅尔波带板；
步骤2、首先采用步骤1中获得的菲涅尔波带板对物体进行扫描，然后利用光电探测器接受扫描后的透射光，经过解调后得到物体的由多个单层物体切片全息图组成的多层物体切片全息图，最后将所有全息图采用傅里叶变换后得到矩阵k；
步骤3、通过步骤2中k矩阵得到物体的多层切片对应的时间反演矩阵，并求得其特征值和特征向量，从而将全息图分解为信号子空间和噪声子空间；
步骤4、首先，利用有限元法将物体沿轴向等间距离散化；然后，对单层物体切片进行离散化处理，并将单层物体切片的各个单元作为测试目标，以求得各个单元作为测试目标时的全息图；最后，利用物体全息图噪声子空间与信号子空间的正交性得到所有单层物体切片的成像伪谱，并将所有单层物体切片的成像伪谱求和即得表征物体轴向位置的参量。
具体地说，步骤1中第一光瞳的函数为矩形1函数，第二光瞳的函数为狄拉克δ函数。
具体地说，步骤1中激光的光学传递函数如下：
将p1(x,y)＝1和p2(x,y)＝δ(x,y)代入式(1)中，则式(1)表示为下式：
则式(2)相应的空间冲击响应为：
其中，j表示虚数单位
具体地说，步骤2中物体的单层切片全息图的关系式如下：
其中，|γ(x,y；z)|2表示物体的复振幅函数，z0表示物体的轴向距离，h(x,y；z0)表示扫描位置在z0的点扩散函数，f和f-1分别表示傅里叶变换和傅里叶反变换，*表示卷积运算。
具体地说，步骤2中物体的双层切片全息图采用傅里叶变换得到矩阵k的公式如下：
k＝f{hc(x,y；z1)+hc(x,y；z2)}(5)
其中，x表示物体的横向坐标，y表示物体的纵向坐标，z1和z2分别表示切片的轴向距离，hc(x,y；z1)和hc(x,y；z2)分别表示物体在z1和z2处的全息图。
具体地说，步骤3的实现方法如下：
首先将式(5)中求得的矩阵k做奇异值分解即可得时间反演矩阵tx和ty，具体如下：
tx＝f-1{kkh}＝f-1{f[γ1*h1+γ2*h2]·{f[γ1*h1+γ2*h2]}h}
ty＝f-1{khk}＝f-1{{f[γ1*h1+γ2*h2]}h·f[γ1*h1+γ2*h2]}(6)
其中，h表示矩阵的hermite运算，khk表示k的共轭转置阵与k的乘积；kkh表示k与k的共轭转置阵的乘积，γ1＝|γ(x,y；z1)|2，γ2＝|γ(x,y；z2)|2，h1＝h(x,y；z1)，h2＝h(x,y；z2)，(h(x,y；z1))h＝h(x,y；-z1)，(h(x,y；z2))h＝h(x,y；-z2)，h(x,y；z1)表示菲涅尔波带在z1处的冲击响应，h(x,y；z2)表示菲涅尔波带在z2处的冲击响应，h(x,y；-z1)表示菲涅尔波带在-z1处的冲击响应，h(x,y；-z2)表示菲涅尔波带在-z2处的冲击响应，并且osh中全息图的重建过程为：
|γ(x,y；z0)|2＝(|γ(x,y；z0)|2)*h(x,y；z0)*h(x,y；-z0)(7)
因此，可以得到时间反演矩阵tx和ty如下：
其中，*表示共轭运算，
首先从式(8)可以看出时间反演阵tx和ty中的前两项只包含了在z1和z2处的物体的信息，由于后两项的值不为0，导致后时间反演阵包含了在z1和z2处物体信息的交叉项，这两项即为我们不需要的噪声项；然后，求得时间反演阵tx和ty的特征值λ和特征向量v1、v2，根据特征值的大小可以将原多层切片全息图分解为信号子空间和噪声子空间；由于式(8)中前两项为第一层全息图的自相关与第二层全息图的自相关之和，后两项理解为第一层全息图与第二层全息图的互相关，互相关的值小于自相关的值，因此在做特征值分解之时，根据特征值的大小可以将tx和ty中的前两项分到信号子空间，将(8)式中的后两项大部分信息分解到噪声子空间；又根据tx和ty中的前两项可以得到v1和v2分别表征物体x和y方向的位置信息，其中信号子空间与噪声子空间具有正交性：
＜v1(i＝1,...m),v1(i＝m+1,...n)＞＝0
＜v2(i＝1,...m),v2(i＝m+1,...n)＞＝0
其中，m(m≤n)表示非零特征值的个数，n表示特征值的总个数。
具体地说，步骤4中单层物体切片成像伪谱的方法如下：
用v1中噪声子空间对应的特征向量的共轭转置左乘k(xp)并求和得qx，用v2中噪声子空间对应的特征向量右乘k(xp)并求和qy：
其中xp表示测试目标的位置，k(xp)为测试目标xp通过式(4)求得的全息图，p表示单层切片离散化的第p个单元，p＝1,2,3...n2，v1(i)和v2(i)分别表示特征值从大到小排列时，tx和ty的第i个特征值对应的特征向量，h表示矩阵的hermite运算，m表示向量的第m个元素，当测试目标和探测目标的位置相同时，qx(xp)和qy(xp)的值约等于零，而测试目标不在探测目标位置时，qx(xp)和qy(xp)的值有限；求得测试目标xp的qx和qy之后，即可求得该点的x方向和y方向的成像伪谱为：
px(xp)＝||k(xp)||2/qx(xp)
py(xp)＝||k(xp)||2/qy(xp)(10)
其中px(xp)和py(xp)分别表示测试目标在x和y方向上的成像伪谱，式(10)中的px(xp)和py(xp)相乘可以得到测试目标xp的成像伪谱：p(xp)＝px(xp)py(xp)，当逐点求得p(xp)之后，可以得到单层切片的成像伪谱。
具体地说，步骤4中单层物体切片成像伪谱的方法如下：
首先，对测试目标xp运用式(4)求得的全息图k(xp)做奇异值分解，即求得k(xp)时间反演阵k(xp)(k(xp))h及(k(xp))hk(xp)的非零特征值对应的特征向量vx和vy；然后通过下式即可求得测试目标xp的qx和qy，即：
其中*表示共轭运算，t表示转置运算，v1(i)和v2(i)分别表示特征值从大到小排列时，tx和ty的第i个特征值对应的特征向量，vx和vy分别表示k(xp)[k(xp)]h及[k(xp)]hk(xp)的非零特征值对应的特征向量，k(xp)为测试目标xp通过式(4)求得的全息图，p表示单层切片离散化的第p个单元，p＝1,2,3...n2，h表示矩阵的共轭转置运算；当测试目标和探测目标的位置相同时，qx(xp)和qy(xp)的值约等于零，于是：
px(xp)＝||vx||2/qx(xp)
py(xp)＝||vy||2/qy(xp)
p(xp)＝px(xp)py(xp)(12)
其中px(xp)和py(xp)分别表示探测目标在x和y方向上的成像伪谱，p(xp)表示包含探测目标位置信息的成像伪谱。
具体地说，步骤4中根据单层物体切片的成像伪谱得到表征物体轴向距离的参数如下：
其中r(z)表示在z处的表征轴向位置的参量，p表示单层切片离散化的第p个单元，n表示特征值的总个数，其局部最大值表征了探测目标的轴向位置。
与现有技术相比，本发明具有以下有益效果：
(1)本发明首创性的使用tr-music算法全息图进行处理，获得物体切片的轴向超分辨率定位。
(2)本发明使用tr-music算法实现物体轴向超分辨率定位，具有很高的精度与分辨率，且受噪声的影响小，适用范围广。
(3)本发明使用了针对tr-music在光学扫描全息中的应用，给出了多层切片的相关理论及公式，即式(1)-式(13)，并且结合实际给出了两种获取单层物体切片成像伪谱的方法，最后利用所获取的单层物体切片成像伪谱进一步实现了获取高分辨率的轴向距离。
(4)本发明不仅实现方式简单、便于操作，同时具有很强的实用性，适合推广使用。
附图说明
图1为本发明流程示意图。
图2为本发明实施例采用的基本结构图。
图3a为本发明实施例1中仿真1的结果图。
图3b为本发明实施例1中仿真2的结果图。
图3c为本发明实施例1中仿真3的结果图。
图4为本发明实施例2中仿真4的结果图。
图5a为本发明实施例3中仿真5的结果图。
图5b为本发明实施例3中仿真6的结果图。
图5c为本发明实施例3中仿真7的结果图。
图6为本发明实施例4中仿真8的结果图。
图7a为本发明实施例5中z1＝5mm处的结果图。
图7b为本发明实施例5中z2＝8mm处的结果图。
图7c为本发明实施例5中其轴向位置分布情况图。
具体实施方式
下面结合附图说明和实施例对本发明作进一步说明，本发明的方式包括但不仅限于以下实施例。
如图1所示，本发明提供的基于tr-music算法的光学扫描全息轴向定位方法，采用tr-music算法对多层切片全息图进行处理，实现各层轴向的切片位置，具有很高的精度与分辨率；这种基于tr-music算法的光学扫描全息轴向超分辨率定位不仅实现方式简单、便于操作，同时具有很强的实用性，适合推广使用。本发明包括以下步骤：
步骤1、将激光采用第一偏振分束器分成两束，之后一束光依次通过第一光瞳和第一凸透镜后投射至第二偏振分束器，第二束光依次通过第二光瞳和第二凸透镜后投射至第二偏振分束器，第二偏振分束器将投射来的两束光聚光干涉形成菲涅尔波带板；
步骤2、首先采用步骤1中获得的菲涅尔波带板对物体进行扫描，然后利用光电探测器接受扫描后的透射光，经过解调后得到物体的由多个单层物体切片全息图组成的多层物体切片全息图，最后将所有全息图采用傅里叶变换后得到矩阵k；
步骤3、通过步骤2中k矩阵得到物体的多层切片对应的时间反演矩阵，并求得其特征值和特征向量，从而将全息图分解为信号子空间和噪声子空间；
步骤4、首先，利用有限元法将物体沿轴向等间距离散化；然后，对单层物体切片进行离散化处理，并将单层物体切片的各个单元作为测试目标，以求得各个单元作为测试目标时的全息图；最后，利用物体全息图噪声子空间与信号子空间的正交性得到所有单层物体切片的成像伪谱，并将所有单层物体切片的成像伪谱求和即得表征物体轴向位置的参量。
本发明采用如图2所示的结构来实现，激光光源laser发射的激光通过第一偏振分束器bs1后分成两束激光，一束激光通过第一反光镜m1后再依次通过第一光瞳p1(x,y)和第一凸透镜l1投射至第二偏振分束器bs2，另一束激光首先通过声光调制器aofs进行调制，调制后的激光束再依次通过第二反射镜m2、第二光瞳p2(x,y)和第二凸透镜l2后透射至第二偏振分束器bs2，两束激光经过第二偏振分束器bs2后干涉形成菲涅尔波带板进入扫描器，扫描器对物体扫描后，透射光被光电转换器接受，最后，经过解调传输至pc端，得到物体的全息图。本发明中he-ne激光器laser发射出激光的波长λ＝632.8nm，两个凸透镜(l1、l2)的焦距都为400mm，物体采用的切片尺寸为100μm×100μm。
为了帮助本领域技术人员能够更加直观的理解本发明，本发明人特提供以下实施例加以详细阐述。
实施例1
本发明使用tr-music算法实现物体轴向超分辨率定位，具有超高的精度与分辨率，现通过以下仿真研究加以阐述：
仿真1、z1＝34mm，z2＝34.4mm，z1和z2在(10，5)处放了一个探测目标时，如图3a所示。
仿真2、z1＝34mm，z2＝34.04mm，z1和z2在(10，5)处放了一个探测目标时，如图3b所示。
仿真3、z1＝100μm，z2＝100.004μm，z1和z2在(10，5)处放了一个探测目标时，如图3c所示。
根据图3中a和b可以看出对于轴向距离为34mm的物体，其分辨率可以达到0.4mm，而从图3c可知对于轴向距离为100μm的物体分辨率可以达到0.004μm。因此，认为该方法具有很高的分辨率。
上述仿真研究通过以下步骤实现：
步骤1、获得菲涅尔波带板，具体为：
首先，激光光源发出的激光经过第一偏振分束器bs1分束，再分别经过第一光瞳p1(x,y)和第二光瞳p2(x,y)后通过第二偏振分束器聚合后在待测物体前干涉形成菲涅尔波带板，其光学传递函数如下：
将p1(x,y)＝1和p2(x,y)＝δ(x,y)代入式(1)中，则式(1)表示为下式：
则式(2)相应的空间冲击响应为：
其中，j表示虚数单位
步骤2、用菲涅尔波带板扫描物体，得到物体的全息图，具体过程如下：
首先，用菲涅尔波带板对物体进行扫描；然后，利用光电探测器接受到透射光，经过解调之后，在pc端可以得到物体的多层切片全息图。该过程单层切片全息图的关系式如下：
其中，|γ(x,y；z0)|2表示物体的复振幅函数，f和f-1分别表示傅里叶变换和傅里叶反变换，z0表示物体的轴向距离，h(x,y；z0)表示扫描位置在z0的点扩散函数，*表示卷积运算。
根据物体的多层切片全息图傅里叶变换后得到矩阵k，从而得到物体的时间反演阵tx和ty。其中
k＝f{hc(x,y；z1)+hc(x,y；z2)}(5)
其中，x表示物体的横向坐标，y表示物体的纵向坐标，z1和z2分别表示切片的轴向距离，hc(x,y；z1)和hc(x,y；z2)分别表示物体在z1和z2处的全息图。
步骤3、通过k矩阵得到物体的多层切片对应的时间反演矩阵，并求得其特征值和特征向量，从而将全息图分解为信号子空间和噪声子空间。其实现过程如下：
对矩阵k做奇异值分解，根据式(5)可以求得矩阵k，然后，求得时间反演矩阵tx和ty。具体为：
tx＝f-1{kkh}＝f-1{f[γ1*h1+γ2*h2]·{f[γ1*h1+γ2*h2]}h}
ty＝f-1{khk}＝f-1{{f[γ1*h1+γ2*h2]}h·f[γ1*h1+γ2*h2]}(6)
其中，h表示矩阵的hermite运算，khk表示k的共轭转置阵与k的乘积；kkh表示k与k的共轭转置阵的乘积，γ1＝|γ(x,y；z1)|2，γ2＝|γ(x,y；z2)|2，h1＝h(x,y；z1)，h2＝h(x,y；z2)，(h(x,y；z1))h＝h(x,y；-z1)，(h(x,y；z2))h＝h(x,y；-z2)，h(x,y；z1)表示菲涅尔波带在z1处的冲击响应，h(x,y；z2)表示菲涅尔波带在z2处的冲击响应，h(x,y；-z1)表示菲涅尔波带在-z1处的冲击响应，h(x,y；-z2)表示菲涅尔波带在-z2处的冲击响应，并且osh中全息图的重建过程为：
|γ(x,y；z0)|2＝(|γ(x,y；z0)|2)*h(x,y；z0)*h(x,y；-z0)(7)
因此，可以得到：
其中，*表示共轭运算，
首先从式(8)可以看出时间反演阵tx和ty中的前两项只包含了在z1和z2处的物体的信息，由于后两项的值不为0，导致后两项包含了在z1和z2处物体信息的交叉项，这两项即为噪声项。
然后，求得时间反演阵tx和ty的特征值λ和特征向量v1、v2，根据特征值的大小可以将原多层切片全息图分解为信号子空间和噪声子空间。由于式(8)中前两项为第一层全息图的自相关与第二层全息图的自相关之和，后两项理解为第一层全息图与第二层全息图的互相关，互相关的值小于自相关的值，因此在做特征值分解之时，根据特征值的大小可以将tx和ty中的前两项分到信号子空间，将(8)式中的后两项大部分信息分解到噪声子空间。又根据tx和ty中的前两项可以得到v1和v2分别表征物体x和y方向的位置信息。其中信号子空间与噪声子空间具有正交性：
＜v1(i＝1,...m),v1(i＝m+1,...n)＞＝0
＜v2(i＝1,...m),v2(i＝m+1,...n)＞＝0
其中,m(m≤n)表示非零特征值的个数(探测目标个数)，n表示特征值的总个数。
步骤4、首先，利用有限元法将物体沿轴向等间距离散化；然后，对单层物体切片进行离散化处理，并将单层物体切片的各个单元作为测试目标，以求得各个单元作为测试目标时的全息图；最后，利用物体全息图噪声子空间与信号子空间的正交性得到所有单层物体切片的成像伪谱，并将所有单层物体切片的成像伪谱求和即得表征物体轴向位置的参量。其实现过程如下：
首先，利用有限元法将物体单层切片做等间距离散化，并将各个单元作为测试目标xp(n×n个点)；
然后，求得各个单元作为测试目标时的全息图k(xp)(p＝1,2,3...n2)。
最后，根据以下两种方法中的任意一种就可求得包含探测目标位置信息的成像伪谱：
方法1：用v1中噪声子空间对应的特征向量的共轭转置左乘k(xp)并求和得qx，用v2中噪声子空间对应的特征向量右乘k(xp)并求和qy：
其中xp表示测试目标的位置，k(xp)为测试目标xp通过式(4)求得的全息图，p表示单层切片离散化的第p个单元，p＝1,2,3...n2，v1(i)和v2(i)分别表示特征值从大到小排列时，tx和ty的第i个特征值对应的特征向量，h表示矩阵的hermite运算，m表示向量的第m个元素，当测试目标和探测目标的位置相同时，qx(xp)和qy(xp)的值约等于零。而测试目标不在探测目标位置时，qx(xp)和qy(xp)的值有限。求得测试目标xp的qx和qy之后，求得该点的x方向和y方向的成像伪谱为：
px(xp)＝||k(xp)||2/qx(xp)
py(xp)＝||k(xp)||2/qy(xp)(10)
其中px(xp)和py(xp)分别表示测试目标在x和y方向上的成像伪谱，式(10)中的px(xp)和py(xp)相乘可以得到测试目标xp的成像伪谱：p(xp)＝px(xp)py(xp)，当逐点求得p(xp)之后，可以得到单层切片的成像伪谱。
方法2：首先，对测试目标xp的全息图k(xp)做奇异值分解，即求得k(xp)时间反演阵k(xp)(k(xp))h及(k(xp))hk(xp)的非零特征值对应的特征向量vx和vy。然后通过下式即可求得测试目标xp的qx和qy，即：
其中*表示共轭运算，t表示转置运算，v1(i)和v2(i)分别表示特征值从大到小排列时，tx和ty的第i个特征值对应的特征向量，vx和vy分别表示k(xp)[k(xp)]h及[k(xp)]hk(xp)的非零特征值对应的特征向量，k(xp)为测试目标xp通过式(4)求得的全息图，p表示单层切片离散化的第p个单元，p＝1,2,3...n2，h表示矩阵的共轭转置运算，当测试目标和探测目标的位置相同时，qx(xp)和qy(xp)的值约等于零。于是：
px(xp)＝||vx||2/qx(xp)
py(xp)＝||vy||2/qy(xp)
p(xp)＝px(xp)py(xp)(12)
其中px(xp)和py(xp)分别表示探测目标在x和y方向上的成像伪谱，p(xp)表示包含探测目标位置信息的成像伪谱。
最终，根据上述两种方法中的任意一种求得的单层切片成像伪谱均可得到表征物体轴向位置的参数，具体如下：
其中r(z)表示在z处的表征轴向位置的参量，p表示单层切片离散化的第p个单元，n表示特征值的总个数，其局部最大值表征了探测目标的轴向位置，最终得出本仿真研究的结果如图3a、3b、3c所示。
实施例2
本发明使用tr-music算法实现物体轴向超分辨率定位，受噪声的影响小，适用范围广。现通过以下仿真研究加以阐述：
仿真4、z1＝100μm，z2＝100.008μm，z1和z2在(10，5)处放了一个探测目标，且加入60db的高斯白噪声时，如图4所示，可以看出本发明受噪声影响较小。
本仿真研究实现步骤与实施例1的仿真研究的实现步骤相同，此处不再赘述。
实施例3
本发明使用tr-music算法实现物体轴向超分辨率定位，能有效对单层中的多个点进行超分辨率轴向定位，具体采用仿真进行阐述为：
仿真5、z1＝100μm处放置4个点，z2＝101μm处放置4个点时，如图5a所示。
仿真6、z1＝100μm处放置4个点，z2＝101μm处放置3个点时，如图5b所示。
仿真7、z1＝100μm处放置4个点，z2＝101μm处放置2个点时，如图5c所示。
根据上述仿真可以看出本发明对于多点物体依然适用，且分辨率可以达到1μm。
本仿真研究实现步骤与实施例1的仿真研究的实现步骤相同，此处不再赘述。
实施例4
本发明使用tr-music算法实现物体轴向超分辨率定位，能有效对多层物体切片进行超分辨率轴向定位，具体采用仿真进行阐述为：
仿真8、z1＝100μm处放置4个点，z2＝101μm处放置3个点，z3＝102μm处放置2个点时，结果如图6所示，可以看出本发明对于三层物体依然适用。
本仿真研究实现步骤与实施例1的仿真研究的实现步骤相同，此处不再赘述。
实施例5
本发明使用tr-music算法实现物体轴向超分辨率定位，适用范围广，对于复杂图形依然适用，采用仿真9进行阐述，本仿真选用物体的切片尺寸为1mm×1mm，具体为：
如图7所示，a为在z1＝5mm图像，b为在z2＝8mm处的图像，c为其轴向位置分布情况。根据仿真结果，可以看出，对于复杂图形，本发明依然可以准确地得到物体的轴向位置。
本仿真研究实现步骤与实施例1的仿真研究的实现步骤相同，此处不再赘述。
综上所述，本发明首创性的使用tr-music算法全息图进行处理，采用本方法能有效获得物体切片的轴向超分辨率定位，本发明具有很高的精度与分辨率，且受噪声的影响小，适用范围广，适于在本技术领域或与本技术领域接近的技术领域大力推广。
上述实施例仅为本发明的优选实施方式之一，不应当用于限制本发明的保护范围，但凡在本发明的主体设计思想和精神上作出的毫无实质意义的改动或润色，其所解决的技术问题仍然与本发明一致的，均应当包含在本发明的保护范围之内。
- 还没有人留言评论。精彩留言会获得点赞！