一种光在水下湍流中近距离传播的相位屏仿真方法与流程

本发明涉及一种光在水下湍流中近距离传播的相位屏仿真方法，属于水下目标探测领域。

背景技术：

水下目标物的探测是海洋探测的重要分支之一。对于运动的水下目标，其在运动过程中必然会产生湍流。将目标物的直接探测转化为其产生的湍流场的间接探测，是一种重要的探测手段。因此，对水下湍流场的特性研究是开展水下探测等研究的主要问题。在利用激光进行水下湍流的特征研究时，对于光传播而言，湍流对其主要的影响就是由于水介质的折射率波动导致的光波的相位随机起伏。因此，建立光在湍流中的传播模型是进行水下湍流研究，开展水下目标探测的前提和基础。

建立水下湍流多层相位屏的算法模型是开展水下湍流研究的主要问题。进行水下湍流特性的研究可以为水下湍流探测、水下光通信等研究打下基础。传统的相位屏生成方法有Zernike多项式法，分形相位屏法，几何近似法，协方差矩阵预测相位屏法。这些方法都能较好地模拟湍流的特性，但是不适用于近距离的光传播模拟。光的近距离传播模型需要考虑湍流的连续性带来的多层相位屏之间的相关性，而传统模型中认为多层相位屏之间是相互独立的。

因此现有技术中，虽然水下湍流的光传播算法模型有很多，并且都较为成熟，但是在考虑近距离传播的实际情况时，需要对这些方法提出适当的修正来符合真实情况。

技术实现要素：

本发明的目的是为了解决现有技术的多层相位屏没有考虑到相位屏之间的相关性，从而导致相位屏模型不能真实模拟湍流特性的缺点，而提出一种光在水下湍流中近距离传播的相位屏仿真方法。

一种光在水下湍流中近距离传播的相位屏仿真方法，包括如下步骤：

步骤1)输入初始图像，初始图像用于生成相互间具有相关性的相位屏；初始图像中像素点的像素值用于表征相位值；初始图像中各个像素点的像素值和相对差值均被标记为已知值；初始图像和相位屏中的特定数量和位置的像素点被选定为基准点，并且每一行和每一列均具有相同数量的基准点；

步骤2)使用所有被标记为已知值的基准点，计算第a个待生成相位屏中第b个基准点的像素值，并将其像素值标记为已知值；a为从1到N的正整数，b为从1到M的正整数，N为待生成的相位屏个数，M为每个相位屏中基准点的个数；

步骤3)使用所有被标记为已知值的基准点的位置信息输入至三维协方差矩阵预测模型，得到第a个相位屏的第b个基准点的概率分布值，将概率分布值乘以与位置相关的随机波动值，得到相对像素差值。并将其绝对像素差值标记为已知值；

步骤4)将像素值与相对像素差值相加，得到第a个待生成相位屏的第b个基准点的绝对值；

步骤5)将b值加1；

步骤6)重复步骤2)至步骤5)，直到第a个待生成相位屏的M个基准点都计算完毕；

步骤7)将a值加1；

步骤8)重复步骤2至7，直到所有的N个相位屏中的所有的基准点都计算完毕；

步骤9)对于第1个待生成相位屏，生成(K-P+1)×(K-P+1)个进程，每个进程含有1个框架，每个框架含有第1个待生成相位屏中的P×P个像素；每个进程用于计算其含有的框架中所有像素点的绝对值，各个进程中的框架互不相同，相互之间的关系为：第1个进程的框架左上角的点为当前相位屏第1行第1列的点，第2个进程的框架左上角的点为当前相位屏第1行第2列的点……第(K-P+1)×(K-P+1)个进程的框架左上角的点为当前相位屏的第K行第K列的点。其中K为相位屏一条边的基准点的个数。每个框架中包括Q个基准点；

步骤10)使用并行计算的方法同时计算框架内每个点的绝对值，具体为：

使用框架内的被标记为已知值的基准点计算任意一个未被标记为已知值的点的绝对值，直至框架内的所有点都计算完成；

步骤11)所有进程计算完成后，将每个进程对相同的点计算得到绝对值求平均，得到最终每个点的绝对值；

步骤12)对于下一个相位屏，重复步骤9)至步骤11)，直至所有的N个相位屏中的所有点都计算完毕，并输出所述的N个相位屏。

本发明的有益效果为：在构建相位屏模型时，使用所有已知点的信息来构建新的相位屏，使得形成的相位屏之间具有相关性；在湍流折射率结构函数的模拟结果上较理论曲线有较好的符合度；相较于其他成熟的算法，本发明可以使得湍流在光的传播方向上具有一定的相关性。

附图说明

图1为本发明的光在水下湍流中近距离传播的相位屏仿真方法的流程图；

图2为本发明的光在水下湍流中近距离传播的相位屏仿真方法的一个实施例的流程图；

图3为本发明的光在水下湍流中近距离传播的相位屏仿真方法的一个实施例的示意图；

图4为随机入射光单元的示意图；

图5为光通过湍流时的相位畸变情况示意图；

图6为使用本发明的方法得到的相位屏的湍流折射率结构函数的模拟结果与理论曲线的对比图；图中星形离散点为仿真结果得到的点；图中的连续曲线为理论值对应的曲线；

图7为用本发明的方法得到的不同间隔的相位屏之间的相关系数示意图；图中横坐标为两个相位屏的间距；图中纵坐标为相关系数。

具体实施方式

具体实施方式一：本实施方式的光在水下湍流中近距离传播的相位屏仿真方法，如图1所示，包括如下步骤：

步骤1)输入初始图像，初始图像用于生成相互间具有相关性的相位屏；初始图像中像素点的像素值用于表征相位值；初始图像中各个像素点的像素值和相对差值均被标记为已知值；初始图像和相位屏中的特定数量和位置的像素点被选定为基准点，并且每一行和每一列均具有相同数量的基准点。

初始图像即为图3(a)，每个方格为一个像素点，初始图像也属于相位屏，图中每个三个点被标出的就是基准点，图3(a)中有9×9＝81个基准点。这些像素点的像素值和相对差值均被标记为已知值。相对差值是与像素点的位置相关的量，只要知道这个像素点是第几个相位屏中第几行第几列的点，就可以唯一确定一个相对差值。

从步骤2)开始，本发明将生成一系列的相位屏，并计算其中每一个点的像素值和相对差值。每一个待生成的相位屏中，像素点的个数和基准点的排布与初始图像是相同的，区别仅在于待生成的相位屏的每一个点的像素值和相对差值在一开始都是未知的，需要通过本发明的方法进行计算。

步骤2)使用所有被标记为已知值的基准点，计算第a个待生成相位屏中第b个基准点的像素值，并将其像素值标记为已知值；a为从1到N的正整数，b为从1到M的正整数，N为待生成的相位屏个数，M为每个相位屏中基准点的个数。

具体而言，就是使用已知的所有基准点来计算当前新的待生成的新相位屏中的未知点。例如，第一次计算是用初始图像中的81个基准点来计算第一个待生成的相位屏中的第一个(第一行第一列)基准点，具体计算的内容是计算像素值和相对差值。该点计算完毕后，使用初始图像的81个点和刚刚计算完的点，共82个点，来计算第一个待生成的相位屏中的第二个(第一行第五列)的基准点。以此类推，当前需要计算的每一个点都是使用所有先前计算过的已知点来确定的，这样做可以保证相位屏之间的相关性。

步骤3)使用所有被标记为已知值的基准点的位置信息输入至三维协方差矩阵预测模型，得到第a个相位屏的第b个基准点的概率分布值，将概率分布值乘以与位置相关的随机波动值，得到相对像素差值。并将其绝对像素差值标记为已知值；

三维协方差矩阵预测模型可以用如下公式表示：

和代表任意两个空间位置矢量；x₁,y₁,z₁为的三个坐标分量；x₂,y₂,z₂为的三个坐标分量；x_n,y_n,z_n为当前相位屏图像中第n个像素点的三个坐标分量；为除去所代表的点外，其余的所有点的空间位置矢量；为除去所代表的点外，其余的所有点的空间位置矢量；x′₁,y′₁,z′₁为的三个坐标分量；x′₂,y′₂,z′₂为的坐标分量；x′_n,y′_n,z′_n为当前相位屏图像中除去第n个像素点外，其余所有点的空间位置矢量对应的三个坐标分量；D_φ为相位结构函数；T为权重函数，其在离散区域内值为1，区域外为0；C为相位协方差函数。

步骤4)将像素值与相对像素差值相加，得到第a个待生成相位屏的第b个基准点的绝对值。

步骤5)将b值加1。

步骤6)重复步骤2)至步骤5)，直到第a个待生成相位屏的M个基准点都计算完毕。

步骤7)将a值加1。

步骤8)重复步骤2至7，直到所有的N个相位屏中的所有的基准点都计算完毕。

步骤9)对于第1个待生成相位屏，生成(K-P+1)×(K-P+1)个进程，每个进程含有1个框架，每个框架含有第1个待生成相位屏中的P×P个像素；每个进程用于计算其含有的框架中所有像素点的绝对值，各个进程中的框架互不相同，相互之间的关系为：第1个进程的框架左上角的点为当前相位屏第1行第1列的点，第2个进程的框架左上角的点为当前相位屏第1行第2列的点……第(K-P+1)×(K-P+1)个进程的框架左上角的点为当前相位屏的第K行第K列的点。其中K为相位屏一条边的基准点的个数。每个框架中包括Q个基准点。

具体而言，步骤9)使用并行进程来计算每个框架中的像素值。框架的示意图如图3(b)所示，黑色粗线框即为前述的框架，图3(b)的框架中有3×3＝9个基准点。图3(b)的框架显示的是第1个进程的框架，第二个进程的框架是在第一个进程的框架向右平移4个像素单位形成的。以此类推。

步骤10)使用并行计算的方法同时计算框架内每个点的绝对值，具体为：

使用框架内的被标记为已知值的基准点计算任意一个未被标记为已知值的点的绝对值，直至框架内的所有点都计算完成。

需要说明的是步骤10)的计算指的是步骤9)中形成的所有(K-P+1)×(K-P+1)都同时进行本步骤的计算内容。这样多进程计算可以大大减少运算时间。本步骤计算结束后，即表示第1个待生成的相位屏中的所有点均计算完毕。

步骤11)所有进程计算完成后，将每个进程对相同的点计算得到绝对值求平均，得到最终每个点的绝对值。

即所有进程计算出来的点是有重叠的，步骤11)是将重复计算的点求一个平均值。

步骤12)对于下一个相位屏，重复步骤9)至步骤11)，直至所有的N个相位屏中的所有点都计算完毕，并输出的N个相位屏。

上述即为本发明的主要过程。简要概括起来是：首先利用分形方法获得符合湍流统计特性的单层二维随机相位屏。利用初始的相位屏，通过三维协方差矩阵预测方法获得具有一定相关性的多层框架。随后利用三维协方差矩阵预测方法，通过这种三维的预测方法生成的框架，其在层内的相位分布和层间的相位分布都能较好地符合水下湍流的特性。随后利用并行计算的方法，在保证结果的准确性的基础上，大幅度提升计算的速度。最后，将并行计算得到的结果作一定的加权平均计算，获得具有层间相关性的连续多层湍流相位屏。

图2是本发明的光在水下湍流中近距离传播的相位屏仿真方法的一个实施例的流程图。其中随机屏指的是相位屏，分配矩阵是指每个相位屏中基准点排列形成的矩阵，丰富矩阵指的是框架及其内部像素点形成的矩阵，从图2中可以看出利用当前已知的基准点计算未知基准点以及其余点的大致流程。

下面进一步解释本发明的设计思路以及显著效果：

三维协方差矩阵的使用使得在生成新的相位时，本发明由传统的仅考虑二维平面内的相位值拓展到了考虑三维空间中的相位值。而由于维度的扩大，使得计算量大幅度提升，因此在算法优化时，需要对用来参考的已知相位作一定要求的选择来达到减小计算量的目的。考虑到参考相位数量的减少，其计算得到的预测相位值会由于样本数量的不足而发生较理论值的偏差，因此，需要对参考相位的选择方法进行优化获得最佳的方案。

在算法的构建过程中，框架的选取和并行计算单元的选取主要影响计算量的大小，而权重因子的选取则影响了计算相位值与理论值得偏差程度。考虑权重因子的选取时，由于相邻两点的相关程度与二者之间的距离成反比，因此可以取距离的倒数作为权重因子来体现不同区域相位对计算区域的影响大小。并且由于距离较大的相位影响程度较小，则在选择参考相位时，可以相对减少距离较大位置的相位的数量。具体的权重因子需要具体情况具体分析，选择合适的表示形式来获得算法的最优方案。通过调整框架的大小，并行计算的计算单元的大小和加权平均的权重因子，可以得到不同统计特征的湍流模型，对于实际情况需要选取合适的参数。

具体实施方式二：本实施方式与具体实施方式一不同的是：三维协方差矩阵预测模型用于计算空间任意两点间的相位波动情况，其计算公式为：

相位协方差函数描述的是空间任意两点的相位波动情况，并且其表征的是相对波动情况，与真实的相位值无关。

r1，r2代表任意两个空间位置矢量；D_φ为相位结构函数；T为权重函数，其在离散区域内值为1，区域外为0；C为相位协方差函数。

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：

步骤2)中，第i个待生成相位屏中第j个基准点的像素值的计算公式为：

x_i,y_j,z_k分别为空间中任意一点在x方向、y方向以及z方向上的分量；其中A为预测矩阵，φ_m为已知的相位列矩阵，σ为插值标准差，R为方差为1均值为0的高斯随机变量。其中预测矩阵A的形式为：

插值标准差的形式为：

具体推导过程为：

对于一个k层i×j大小的随机屏组，可以将随机屏表示成k个i×j大小的矩阵，如第k层随机屏可以表示为：

对k层随机屏进行列向量化处理，可以得到一个一列i×j×k行的列矩阵：

φ_m＝(φ_1,1,1φ_2,1,1…φ_(i-1),j,kφ_i,j,k)^Τ

则可以利用已知的随机屏给出已知随机屏的协方差矩阵，其形式为：

当需要预测一个新的未知相位值位置的相位时，可以利用已知随机屏的协方差矩阵进行计算预测。如需要获得坐标为(x_i,y_j,z_k)的相位值，则计算公式为：

其中A为预测矩阵，φ_m为已知的相位列矩阵，σ为插值标准差，而R为方差为1均值为0的高斯随机变量。其中预测矩阵的形式和插值标准差的形式为：

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：初始图像及每个相位屏具有9×9个基准点，每个框架具有9×9个像素点，每个框架具有3×3个基准点。

就像前述提到的，框架如果选取的太大会导致计算每个框架内像素时较慢，框架选取太小会导致需要生成的进程过多，消耗更多的系统资源。因此本发明选取含有3×3个基准点的框架大小，可以达到最佳的计算效率。

其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：

计算任意四个已知点构成的最小矩形的几何中心点的绝对值，如果无法找到几何中心点，则按照像素点的顺序计算所有点的绝对值。

这个过程可以如图3(c)所示，先计算所有箭头指向的浅色点，这样的点有4个，这些点的特点为，都是框架内的基准点所能构成的最小矩形的几何中心点。依次类推，直至计算完框架内的所有点。本实施方式仅是提供一种较快计算完所有点的方法，本领域的技术人员可以毫无疑义地确定还可以使用其他方法计算这些点，例如按行列顺序计算，区别仅在于运算速度。

其它步骤及参数与具体实施方式一至四之一相同。

如图4所示，本发明在构建光在水下湍流中传播的模型时，将入射光看作由很多能量相等的光单元组成。光单元在表征湍流统计特性的多层相位屏中传播时，利用相位屏的梯度大小来确定光单元的偏折角度和传播路径。当计算的光单元数量达到足够多时，就可以近似得到最终的出射光的光强分布和相位分布，如图5所示，即最后的光波通过湍流后的结果。

从图6中可以看出，本发明在相位屏间距小于外尺度的范围内很接近与理论曲线。

从图7中可以看出本发明生成的相位屏之间具有明显的随相位屏间距变化的相关系数。

本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，本领域技术人员当可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。