CN115640705A - 一种水下光电成像的psf与bsf仿真模拟方法 - Google Patents

Abstract

本发明提供了一种水下光电成像的PSF与BSF仿真模拟方法，属于水下光电成像技术领域。首先搭建水下光电成像的仿真模型；PSF函数的模拟通过搭建一个余弦发射器模拟朗伯点光源作为光源和准直辐亮度计作为接收器；BSF函数的模拟通过搭建一个激光发射器光源作为准直光源和余弦探测器作为接收器；然后设定水体参数，基于MATLAB蒙特卡罗方法进行仿真模拟；最后，进行数据处理计算，分别获取PSF与BSF的函数值。本发明通过蒙特卡罗数值模拟的方式实现两种函数的获取，解决了函数模型计算推算误差大、机械实际测量成本高等问题，针对不同的水体环境，可用于评估水下系统性能。同时证明了两种函数的互易性，推动了实际应用。

Description

一种水下光电成像的PSF与BSF仿真模拟方法

技术领域

本发明属于水下光电成像技术领域，尤其涉及一种水下光电成像的PSF与BSF仿真模拟方法。

背景技术

为了描述水下光场，近年来诸多从事海洋光学的研究人员相继提出了一些可靠的水下光场的分布模型，试图通过这些模型来描述光场的能量分布，其中包括BSF（beamspread function光束扩散函数）和PSF（point spread function点扩散函数）函数。BSF和PSF用来描述较远距离时的光场能量分布，其中BSF光束扩散函数描述准直光束经过在介质后的光场的空间照度分布，而PSF点扩散函数则描述一个余弦发射体经过介质后，通过对周围光场的辐亮度的测量，然后将该变量对余弦发射体法线辐射强度归一化之后的参量。

戈登上世纪末通过辐射传输方程和光传播的互易定理在理论方面证明了二者在数值上的互易性。将BSF和PSF均看作关于空间位置的函数，PSF函数的计算用到一个朗伯余弦发射体，加一个接收视场很小的辐亮度计。而BSF的计算则用到一个准直光源外加一个照度计。

当前获取两种函数的主要方法，主要分为理论计算与实际测量两种。理论计算法通常使用光迹追踪及透镜理论点扩散函数模型。计算模型建立在光迹追踪的基础上，光迹追踪本是建立在大气环境中的一套理论推算方法。然而现实中的水体是一种浑浊介质，浑浊介质中各种粒子的存在会导致光迹追踪方法的失效。

实际测量法由于测量会随着物距和视场的变化而发生明显的变化，因此针对不同物距、不同的视场均需要获取扩散函数，需要设置不同的工装设备，会明显增加获取成本；同时，针对PSF，由于在不同物距下，光源的孔径也需要调整，以保证其为点光源，否则会导致最终拍摄过程当中出现未获取光源或光源孔径过大的问题。此外，由于该种方式受到拍摄过程中噪声、曝光等因素的影响，使得获取函数的精度较低，针对精度较高的任务，这种方法是不适用的。

发明内容

针对两种函数获取时，现有方式存在的不足（函数模型计算推算误差大，机械实际测量成本高等），本发明通过以蒙特卡罗算法为基础的数值模拟方式，建立了BSF和PSF的仿真模型，弥补计算和机械获取的不足。用MATLAB仿真拟出光源、水体吸收和散射以及接收器的模型，然后连接各个模块分别测出PSF函数和BSF函数的数值。

本发明提出了一种水下光电成像的PSF与BSF仿真模拟方法，并包括以下步骤：

步骤1，搭建水下光电成像的PSF与BSF的仿真模型；PSF函数的模拟通过搭建一个余弦发射器模拟朗伯点光源作为光源和准直辐亮度计作为接收器；BSF函数的模拟通过搭建一个激光发射器光源作为准直光源和余弦探测器作为接收器；

步骤2，设定水体参数，基于MATLAB蒙特卡罗方法进行仿真模拟；首先，两种模拟发射器分别发射光子包，然后，光子通过设定好参数的水体，不断通过随机步长和随机散射方向，更新每个光子的传播位置和能量，最后，两种模拟接收器分别接收光子；

步骤3，数据处理计算；对经由余弦发射器发射，准直辐亮度计获取的光亮度E₁，在获取光亮度后，就知道了接收平面上的光场分布，进而可以得到了PSF函数的值；对经由准直光源发射，余弦探测器获取的光照度E₂，可以计算得出BSF函数的值。

优选的，所述步骤1中余弦发射器模拟了一个朗伯点光源，基于朗伯定律，用公式可以表示为:

将坐标轴z轴规定为I₀的方向，每个光子与z轴方向的夹角记为天顶角θ，θ的取值范围为[0,90°)，同时每一个光子发射方向在x-O-y平面内的投影同x轴正向的夹角定义为方位角φ，φ的取值范围为[0,360°)，对天顶角θ和方位角φ分别进行等间距分割，模拟光子的发射方向。

优选的，所述步骤1中准直光源是模拟了一束激光光源，将坐标轴z轴规定为激光发射的方向，每个光子与z轴方向的夹角记为天顶角θ，拟定θ的取值范围为[0,0.01°），同时每一个光子发射方向在 x-O-y 平面内的投影同x轴正向的夹角定义为方位角φ，φ的取值范围为[0,360°)，在对天顶角θ和方位角φ分别进行等间距分割，模拟光子的发射方向。

优选的，所述步骤1中的准直辐亮度计累计来自小于亮度计最大接收角α方向上的落于接收平面上的光子的全部能量，具体模拟步骤为：

S1，接收平面是长和宽均为 20cm 的平面，该平面上分布着不同位置和发射方向的光子，而且记录了每一个光子的能量大小，规定接收器在该平面上接收光子；将该平面划分成 201×201 的像素点，即对每一个像素点赋予一个坐标，坐标的步进大小均为 1mm；

S2，对第一个像素点(-10,-10)，定义接收器的半径为 D，接收器的中心位于该像素点，即接收器中心坐标为(-10,-10)，如果该光子距离接收器中心的距离小于 D，就能够被接收器有效接收，并获得该光子的信息；

S3，对于(-10,-10)这一像素点，同样也是此刻接收器中心的位置，遍历接收平面的所有光子，统计能满足S2中能被接收器有效接收的光子，将不能被接收器接收的光子的能量设置为0，方便后续计算；

S4，对每一个光子进行分析，加一个判定条件，使光子的入射角度θ小于最大的接收角度ɑ，当光子的入射角度θ大于最大的接收角度ɑ将不会对像元的照度做出贡献，定义ɑ为某个角度，将大角度入射的光子的能量也初始为0，便于后续的计算，由于每一个接收到的光子均以接近垂直的效果被探测器接收到，因此统计光能量时，不需要乘以该光子入射的余弦值；

S5，对该平面所有处理好的光子进行能量的求和；未能被探测器接受到的光子或者被探测器所接受但并不能落在像面上的光子，其能量值已经设置为0，因此不会对接受到的总能量做出贡献，因此累加的总能量就是该次接收的总通量，记为此处的光通量；

S6，利用公式就将像面的照信息转换成本次测量的光源的亮度，公式如下：

其中，E₁为像面亮度，L为光源亮度，D为物镜孔径，f为物镜焦距。

优选的，所述步骤1中的余弦探测器累计来自各个方向上的落于接收平面上光子的能量的余弦值，具体模拟步骤如下：

S1，接收平面是长和宽均为 20cm 的平面，该平面上分布着不同位置和发射方向的光子，而且记录了每一个光子的能量大小，规定接收器在该平面上接收光子，将该平面划分成 201×201 的像素点，即对每一个像素点赋予一个坐标，坐标的步进大小均为 1mm；

S2，对第一个像素点(-10,-10)，测量该处的照度大小，定义接收器的半径为 D，接收器的中心位于该像素点，即接收器中心坐标为(-10,-10)，如果该光子距离接收器中心的距离小于D，就能够被接收器有效接收，并获得该光子的信息；

S3，对于(-10,-10)这一像素点，同样也是此刻接收器中心的位置，遍历接收平面的所有光子，统计能满足步骤二中能被接收器有效接收的光子，将不能被接收器接收的光子的能量设置为0，方便后续的计算；

S4，对每一个光子进行分析，显然，当光子的方向与该平面的法向所呈角为θ时，根据光通量的概念，即单位时间内通过某一曲面能量流，假设该光子的能量为e₀，实际被探测器接受到的能量值为：

S5，对该平面所有处理好的光子进行能量的求和，未能被探测器接受到的光子，其能量值已经设置为0，因此不会对接受到的总能量做出贡献，因此该步骤累加的总能量就是该次接收的总通量，记为此处的光通量；在规定完照度计口径之后，利用公式E= dΦ/dS就能通过该像素点的光通量转换成光照度的大小，其中Φ为照度计接收光子的总通量，S为光子投射的面积。

优选的，所述步骤2的具体过程为：

S21，两种模拟发射器分别发射光子包；对朗伯光源而言，以法线方向(0.0.1)为I₀方向，对朗伯发射体发射的光子能量进行初始化，规定方向发射的光子能量大小为W₀，则其他发射方向的光子能量为:

W₀ = W₀×cosθ

向规定的每个方向发射一定数量光子的光子包；

对激光发射器而言，对每一个发射的光子都赋予一个初始能量W₀，向规定的每个方向发射一定数量光子的光子包；

S22，光子通过设定好参数的水体，不断通过随机步长和随机散射方向，更新每个光子的传播位置和能量；

单个光子随机步长s与衰减系数c的关系是：

其中，s为单个光子的随机步长，τ为[0,1]间随机数，c为提前设定的水体衰减系数，可由散射系数与吸收系数相加得到；

在蒙特卡罗模拟过程中，光子会与介质中粒子发生碰撞，每次碰撞后必须根据散射相函数通过抽样确定出散射方向，采用的Henyey-Greenstein散射相函数的解析形式可表示为：

其中，θ为光子下次动作时与z轴方向的夹角；g为先前输入的不对称因子，当g为0时，水体各向同性；rnd为自动产生的(0,1)区间内随机数；

S23，两种模拟接收器分别接收光子；对准直辐亮度计的接收而言，可以测量光场分布，以1mm 为步进，先固定像素点的 y 坐标，依此将 x 坐标从-10cm 处遍历到 10cm处，再将 y 坐标增加值-9.9cm 处进行 x 坐标的遍历，直至遍历到最后一个像素点(10,10)，将得到每一个像素点测量得到的光源亮度的大小E₁(x,y)，至此测出了该平面所测的光源的亮度分布，也即该平面上的光场分布；

对余弦探测器而言，以1mm 为步进，先固定像素点的y坐标，依次将x坐标从-10cm处遍历到 10cm 处，再将y坐标增加值-9.9cm 处进行 x 坐标的遍历，直至遍历到最后一个像素点(10,10)，将得到每一个像素点测量得到的光照度的大小E₂(x,y)，至此测出了该平面的照度分布，也即该平面上的光场分布。

优选的，所述步骤3中对经由余弦发射器发射，准直辐亮度计获取的光亮度E₁，在获取光亮度后，就知道了该平面上的光场分布，进而可以得到了PSF 函数的值，计算公式为：

其中，I₀是光源法线方向的发光强度，由于光源是朗伯点光源，光源发出的总光通量的大小等于光源的最大出射发光强度I₀乘以π立体角：

因此，在已知光源的光通量的情况下可以求出朗伯光源法线方向的发光强度I₀，实际上，由于仿真过程中，光子是限定在一秒内全部发射，因此，此处的光源的光通量Φ就是光源的光功率P；同时，亮度计是通过像面照度来反应实际物体的亮度，因而可以利用测得的像面的光照度，进而推算出光源的亮度，当选定了物镜之后，焦距随之确定下来，像面的照度值与光源的亮度成正比关系，而且当物镜焦距f <<物距 m 时，通过公式可以推导简化出如下关系：

其中，L为光源亮度，D为物镜孔径，f为物镜焦距；

进而可以导出L(x,y)，L(x,y)表示平面上点(x,y)测出的光源光亮度，PSF(x,y)表示平面上点(x,y)的点扩散函数的值。

优选的，所述步骤3中对经由准直光源发射，余弦探测器获取的光照度E₂，可以计算得出BSF函数的值，可以由下式给出：

其中 P 是光源的总功率，功率的定义是单位时间内能量值的大小，因此只要对发射时的光子的能量进行求和就能算出总功率，E(x,y)表示平面上像素点(x,y)处的光照度大小，BSF(x,y)表示平面上点(x,y)的光束扩散函数的值。

与现有技术相比，本发明提出一种水下光电成像的PSF与BSF仿真模拟方法，并产生如下有益效果：

1.本发明通过蒙特卡罗数值模拟的方式实现两种函数的获取，解决了函数模型计算推算误差大，机械实际测量成本高等问题，针对不同的水体环境，调整相应的水体参数，通过模拟的方式获取PSF与BSF函数，可用于评估水下系统性能，怎样的水体环境，拥有怎样的光学成像效果，提前获取水下效果；

2.本发明建立的模型，既可应用于光学成像所需的PSF，又可应用于光通信所需的BSF。解决了光迹追踪方法所忽略的浑浊介质粒子对光子传输的干扰，提高了获取结果的可靠性；同时模拟了实际测量，减少了各种设备及人物力投入，降低了获取成本，减少了拍摄过程中噪声、曝光等因素的应影响，提高函数获取的精度；

3.根据已有的理论证明，本发明通过两种函数的数值模拟，证明了两种函数在一定范围内，数值上的一致性，即PSF与BSF的互易性。例如在BSF的测量中，若条件受限，可采用PSF模型替代，仿真模型可实现移植；在实际测量过程测量仪器受限时，亦可实现两函数的移植。

附图说明

图1为水下光电成像的PSF与BSF仿真模拟方法的流程示意图。

图2为PSF的模型示意图。

图3为BSF的模型示意图。

图4为模拟余弦发射器示意图。

图5模拟准直光源示意图。

图6为蒙特卡洛仿真发射系统流程图。

图7为蒙特卡洛仿真单个光子水下传播流程图。

图8为蒙特卡洛仿真探测系统流程图。

图9为相同衰减长度和水体参数下PSF和BSF函数测量值图。

图10为相同水体参数下不同衰减长度时BSF和PSF的峰值高度图。

具体实施方式

下面结合具体实施例对发明进行进一步说明。

本发明提出一种水下光电成像的PSF与BSF仿真模拟方法，建立了一种基于MATLAB蒙特卡罗方法的水下光电成像PSF与BSF的仿真模型，模拟了目标光源（朗伯点光源或准直光源）在水体中传播并被接收器（准直辐亮度计或余弦探测器）的接收。通过该方法的建立，有效预评估水下光学系统。并通过函数仿真获取，验证了PSF与BSF的互易性，使得在一种散射函数获取困难的时候，可通过另一种散射函数的获取，达到等价的效果。整体过程如图1所示：

步骤1，搭建水下光电成像的PSF与BSF的仿真模型；PSF函数的模拟通过搭建一个余弦发射器模拟朗伯点光源作为光源和准直辐亮度计作为接收器；BSF函数的模拟通过搭建一个激光发射器光源作为准直光源和余弦探测器作为接收器；

步骤2，设定水体参数，基于MATLAB蒙特卡罗方法进行仿真模拟；首先，两种模拟发射器分别发射光子包，然后，光子通过设定好参数的水体，不断通过随机步长和随机散射方向，更新每个光子的传播位置和能量，最后，两种模拟接收器分别接收光子；

步骤3，数据处理计算；对经由余弦发射器发射，准直辐亮度计获取的光亮度E₁，在获取光亮度后，就知道了该平面上的光场分布，进而可以得到了PSF函数的值；对经由准直光源发射，余弦探测器获取的光照度E₂，可以计算得出BSF函数的值。

一、模型建立：

PSF函数的模拟需要搭建一个余弦发射器（光源）和准直辐亮度计（接收器）；BSF函数的模拟需要搭建一个准直光源（光源）和余弦探测器（接收器）。模拟示意图如图2和图3所示。

1．发射器模拟：

余弦发射器模拟了一个朗伯点光源。想要模拟一个朗伯发射体，需要基于朗伯定律，所谓朗伯定律，是指一个亮度在各个方向上都相等的发光面，在某一方向上的发光强度等于这个面垂直法线方向上的发光强度I₀乘以方向角的余弦。用公式可以表示为如下:

将坐标轴z轴规定为I₀的方向，每个光子与z轴方向的夹角记为天顶角θ，θ的取值范围为[0,90°)，[]表示可以取两边的值，（）表示不可以取那个值，这个取值范围就是大于等于0，小于90°；同时每一个光子发射方向在x-O-y平面内的投影同x轴正向的夹角定义为方位角φ，φ的取值范围为[0,360°)，对天顶角θ和方位角φ分别进行等间距分割，模拟光子的发射方向。如图4所示。

模拟准直光源将坐标轴z轴规定为激光发射的方向，每个光子与z轴方向的夹角记为天顶角θ，因为此处模拟了一束激光光源，故拟定θ的取值范围为[0,0.01°），同时每一个光子发射方向在 x-O-y 平面内的投影同x轴正向的夹角定义为方位角φ，φ的取值范围为[0,360°)。在对天顶角θ和方位角φ分别进行等间距分割，模拟光子的发射方向。如图5所示。

2．接收器模拟：

准直辐亮度计累计来自小于亮度计最大接收角α方向上的落于接收平面上的光子的全部能量；余弦探测器（照度计）累计来自各个方向上的落于接收平面上光子的能量的余弦值。

具体亮度计模拟步骤如下：

（1）接受平面是长和宽均为 20cm 的平面，该平面上分布着不同位置和发射方向的光子，而且记录了每一个光子的能量大小，规定接收器在该平面上接收光子。将该平面划分成 201×201 的像素点，即对每一个像素点赋予一个坐标，坐标的步进大小均为 1mm。

（2）对第一个像素点(-10,-10)，定义接收器的半径为 D，接收器的中心位于该像素点，即接收器中心坐标为(-10,-10)，显然，对于该接受平面上的每一个光子，如果该光子距离接收器中心的距离大于 D，也就是说光子打在了接收器的外部，因此不能够被接收器所接受，所以该光子不会对接收器接受到的能量值有贡献。相反如果该光子距离接收器中心的距离小于 D，就能够被接收器有效接收，并获得该光子的能量等信息。

（3）对于(-10,-10)这一像素点，同样也是此刻接收器中心的位置，遍历接收平面的所有光子，统计能满足步骤二中能被接收器有效接收的光子，将不能被接收器接收的光子的能量设置为0。方便后续的计算。

（4）对每一个光子进行分析，显然，当光子的方向与该平面的法向所呈角为θ，此时区别于照度计的设计，由于亮度计设计时有物镜的存在，所以对测量视场也有要求，在微观条件下，即光子如果以与接收器法线方向呈大角度入射进物镜，实际上并不能落在像孔内，而是散落在像孔 H 的外部，因此我们需要加一个判定条件，使光子的入射角度θ小于最大的接收角度ɑ，当光子的入射角度θ大于最大的接收角度ɑ将不会对像元的照度做出贡献，实验中我们定义ɑ为 0.5°。我们将大角度入射的光子的能量也初始为0，便于后续的计算。由于每一个接收到的光子均以接近垂直的效果被探测器接收到，因此统计光能量时，不需要乘以该光子入射的余弦值。

（5）对该平面所有处理好的光子进行能量的求和，由上述可知，未能被探测器接受到的光子或者被探测器所接受但并不能落在像面上的光子，其能量值已经设置为0，因此不会对接受到的总能量做出贡献，因此该步骤累加的总能量就是该次接收的总通量，记为此处的光通量。

（6）利用公式就将像面的照信息转换成本次测量的光源的亮度，公式如下：

其中，E₁为像面亮度，L为光源亮度，D为物镜孔径，f为物镜焦距。至此我们就完成了亮度计的模型。

具体照度计模拟步骤如下：

（1）接收平面是长和宽均为 20cm 的平面，该平面上分布着不同位置和发射方向的光子，而且记录了每一个光子的能量大小，规定接收器在该平面上接收光子，将该平面划分成 201×201 的像素点，即对每一个像素点赋予一个坐标，坐标的步进大小均为 1mm。

（2）对第一个像素点(-10,-10)，测量该处的照度大小，定义接收器的半径为 D，接收器的中心位于该像素点，即接收器中心坐标为(-10,-10)，显然，对于该接受平面上的每一个光子，如果该光子距离接收器中心的距离大于D，也就是说光子打在了接收器的外部，因此不能够被接收器所接收，所以该光子不会对接收器接收到的能量值有贡献。相反如果该光子距离接收器中心的距离小于D，就能够被接收器有效接收，并获得该光子的能量等信息。

（3）对于(-10,-10)这一像素点，同样也是此刻接收器中心的位置，遍历接收平面的所有光子，统计能满足（2）中能被接收器有效接收的光子，将不能被接收器接收的光子的能量设置为0。方便后续的计算。

（4）对每一个光子进行分析，显然，当光子的方向与该平面的法向所呈角为θ时，根据光通量的概念，即单位时间内通过某一曲面能量流，假设该光子的能量为e₀。实际被探测器接受到的能量值为：

（5）对该平面所有处理好的光子进行能量的求和，由上述可知，未能被探测器接受到的光子，其能量值已经设置为 0，因此不会对接受到的总能量做出贡献，因此该步骤累加的总能量就是该次接收的总通量，记为此处的光通量。在规定完照度计口径之后，利用公式E= dΦ/dS（Φ为照度计接收光子的总通量，S为光子投射的面积）就能通过该像素点的光通量转换成光照度的大小。至此就完成了照度计的模拟。

二、数值仿真过程：

1．模拟发射器发射一个光子包。

对朗伯光源而言，以法线方向(0.0.1)为I₀方向，对朗伯发射体发射的光子能量进行初始化，规定方向发射的光子能量大小为W₀.则其他发射方向的光子能量为:

W ₀ = W₀×cosθ

向规定的每个方向发射一定数量光子的光子包。

对激光发射器而言，对每一个发射的光子都赋予一个初始能量W₀。向规定的每个方向发射一定数量光子的光子包。具体发射流程如图6所示。

2.模拟光束穿过水体。

光子通过设定好的水体（已知散射、吸收系数），不断通过随机步长和随机散射方向，更新每个光子的传播位置和能量。

单个光子随机步长s与衰减系数c的关系是：

其中，s为单个光子的随机步长，τ为[0,1]间随机数，c为提前设定的水体衰减系数，可由散射系数与吸收系数相加得到。

在蒙特卡罗模拟过程中，光子会与介质中粒子发生碰撞，每次碰撞后必须根据相函数通过抽样确定出散射方向。因此根据相函数进行的散射方向抽样具有重要意义。Henyey-Greenstein散射相函数的解析形式可表示为：

其中，θ为光子下次动作时与z轴方向的夹角；g为先前输入的不对称因子，当g为0时，水体各向同性；rnd为自动产生的(0,1)区间内随机数。单个光子水下传播流程如图7所示。

3.模拟接收器接收光子。

对准直辐亮度计的接收而言，可以测量光场分布。以1mm为步进，先固定像素点的y坐标，依此将x坐标从-10cm处遍历到10cm处，再将y坐标增加值-9.9cm处进行x坐标的遍历，直至遍历到最后一个像素点(10,10)。将得到每一个像素点测量得到的光源亮度的大小E(x,y)，至此测出了该平面的所测的光源的亮度分布，也即该平面上的光场分布。

对余弦探测器而言（照度计）而言，以1mm 为步进，先固定像素点的y坐标，依次将x坐标从-10cm 处遍历到10cm 处，再将y坐标增加值-9.9cm 处进行x坐标的遍历，直至遍历到最后一个像素点(10,10)。我们将得到每一个像素点测量得到的光照度的大小E(x,y)，至此我们测出了该平面的照度分布，也即该平面上的光场分布。蒙特卡洛仿真探测系统如图8所示。

三、数据处理计算：

1．对经由朗伯发射体发射，准直辐亮度计获取的光亮度E₁。在获取光亮度后，就知道了该平面上的光场分布，进而可以得到了PSF 函数， PSF 可以由下式给出：

其中，I₀是光源法线方向的发光强度，由于光源是朗伯点光源，光源发出的总光通量的大小等于光源的最大出射发光强度I₀乘以π立体角：

因此，在已知光源的光通量的情况下可以求出朗伯光源法线方向的发光强度I₀，实际上，由于仿真过程中，光子是限定在一秒内全部发射，因此，此处的光源的光通量Φ就是光源的光功率P；同时，亮度计是通过像面照度来反应实际物体的亮度，因而可以利用测得的像面的光照度，进而推算出光源的亮度的，当选定了物镜之后，焦距随之确定下来，像面的照度值与光源的亮度成正比关系，而且当物镜焦距f <<物距 m 时，通过公式可以推导简化出如下关系：

其中，E1为像面亮度，L为光源亮度，D为物镜孔径，f为物镜焦距。

进而可以导出L(x,y)，L(x,y)表示平面上点(x,y)测出的光源光亮度，PSF(x,y)表示平面上点(x,y)的点扩散函数的值。

2．对经由准直光源发射，余弦探测器（照度计）获取的光照度E₂。BSF(x,y)表示平面上点(x,y)的光束扩散函数的值，可以由下式给出：

其中 P 是光源的总功率，功率的定义是单位时间内能量值的大小，因此只要对发射时的光子的能量进行求和就能算出总功率，E(x,y)表示平面上像素点(x,y)处的光照度大小，BSF(x,y)表示平面上点(x,y)的光束扩散函数的值。

3．为严格证明，在同样的水体条件下测出了BSF和PSF的分布图并观察二者的侧视图，通过对比数值大小来验证二者的等值性，图9是一个典型结果，其中点光源发射32400个光子，激光光源发射36000个光子，并对二者进行了初始能量的归一化，衰减长度为 2AL，散射系数为 1.1m^-1,吸收系数为 0.4m^-1。

从二者的侧视图可以看出，BSF和PSF的数值拟合程度是相当好的，该仿真结果很好的展示了在相同水体参数和同样衰减长度下，PSF和BSF数值上具有等值性这一特点。

为了验证当衰减距离改变时，二者依旧满足数值上的互易性，在设置二者相同的水体条件下，做出二者不同衰减长度下，图像峰值的折线图（图10）进行对比。其中点光源发射的光子数 N1为32400，激光光源发射光子数N2为36000,亮度计接收口径为1.8mm，照度计接受口径为2cm，水体散射系数为 1.1m^-1，吸收系数为0.4m^-1。衰减长度从1AL，以0.5AL为步进增加到4AL依次测量。

从二者随衰减距离变化峰值折线图可以看出，二者的拟合度是相当好的，这也证明了二者在不同距离处的函数分布均有等值性这一特点。

以上所述仅为本申请的优选实施例而已，并不用于限制本申请，对于本领域的技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。

上述虽然对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。

Claims (8)

1.一种水下光电成像的PSF与BSF仿真模拟方法，其特征在于，包括以下步骤：

步骤1，搭建水下光电成像的PSF与BSF的仿真模型；PSF函数的模拟通过搭建一个余弦发射器模拟朗伯点光源作为光源和准直辐亮度计作为接收器；BSF函数的模拟通过搭建一个激光发射器光源作为准直光源和余弦探测器作为接收器；

步骤2，设定水体参数，基于MATLAB蒙特卡罗方法进行仿真模拟；首先，两种模拟发射器分别发射光子包，然后，光子通过设定好参数的水体，不断通过随机步长和随机散射方向，更新每个光子的传播位置和能量，最后，两种模拟接收器分别接收光子；

步骤3，数据处理计算；对经由余弦发射器发射，准直辐亮度计获取的光亮度E₁，在获取光亮度后，就知道了接收平面上的光场分布，进而可以得到了PSF 函数的值；对经由准直光源发射，余弦探测器获取的光照度E₂，可以计算得出BSF函数的值。

2.如权利要求1所述的一种水下光电成像的PSF与BSF仿真模拟方法，其特征在于：所述步骤1中余弦发射器模拟了一个朗伯点光源，基于朗伯定律，用公式可以表示为:

将坐标轴z轴规定为I₀的方向，每个光子与z轴方向的夹角记为天顶角θ，θ的取值范围为[0,90°)，同时每一个光子发射方向在x-O-y平面内的投影同x轴正向的夹角定义为方位角φ，φ的取值范围为[0,360°)，对天顶角θ和方位角φ分别进行等间距分割，模拟光子的发射方向。

3.如权利要求1所述的一种水下光电成像的PSF与BSF仿真模拟方法，其特征在于：所述步骤1中准直光源是模拟了一束激光光源，将坐标轴z轴规定为激光发射的方向，每个光子与z轴方向的夹角记为天顶角θ，拟定θ的取值范围为[0,0.01°），同时每一个光子发射方向在 x-O-y 平面内的投影同x轴正向的夹角定义为方位角φ，φ的取值范围为[0,360°)，在对天顶角θ和方位角φ分别进行等间距分割，模拟光子的发射方向。

4.如权利要求1所述的一种水下光电成像的PSF与BSF仿真模拟方法，其特征在于，所述步骤1中的准直辐亮度计累计来自小于亮度计最大接收角α方向上的落于接收平面上的光子的全部能量，具体模拟步骤为：

S1，接收平面是长和宽均为 20cm 的平面，该平面上分布着不同位置和发射方向的光子，而且记录了每一个光子的能量大小，规定接收器在该平面上接收光子；将该平面划分成201×201 的像素点，即对每一个像素点赋予一个坐标，坐标的步进大小均为 1mm；

S2，对第一个像素点(-10,-10)，定义接收器的半径为 D，接收器的中心位于该像素点，即接收器中心坐标为(-10,-10)，如果该光子距离接收器中心的距离小于 D，就能够被接收器有效接收，并获得该光子的信息；

S3，对于(-10,-10)这一像素点，同样也是此刻接收器中心的位置，遍历接收平面的所有光子，统计能满足S2中能被接收器有效接收的光子，将不能被接收器接收的光子的能量设置为0，方便后续计算；

S4，对每一个光子进行分析，加一个判定条件，使光子的入射角度θ小于最大的接收角度ɑ，当光子的入射角度θ大于最大的接收角度ɑ将不会对像元的照度做出贡献，定义ɑ为某个角度，将大角度入射的光子的能量也初始为0，便于后续的计算，由于每一个接收到的光子均以接近垂直的效果被探测器接收到，因此统计光能量时，不需要乘以该光子入射的余弦值；

S5，对该平面所有处理好的光子进行能量的求和；未能被探测器接受到的光子或者被探测器所接受但并不能落在像面上的光子，其能量值已经设置为0，因此不会对接受到的总能量做出贡献，因此累加的总能量就是该次接收的总通量，记为此处的光通量；

S6，利用公式就将像面的照信息转换成本次测量的光源的亮度，公式如下：

其中，E₁为像面亮度，L为光源亮度，D为物镜孔径，f为物镜焦距。

5.如权利要求1所述的一种水下光电成像的PSF与BSF仿真模拟方法，其特征在于，所述步骤1中的余弦探测器累计来自各个方向上的落于接收平面上光子的能量的余弦值，具体模拟步骤如下：

S1，接收平面是长和宽均为 20cm 的平面，该平面上分布着不同位置和发射方向的光子，而且记录了每一个光子的能量大小，规定接收器在该平面上接收光子，将该平面划分成201×201 的像素点，即对每一个像素点赋予一个坐标，坐标的步进大小均为 1mm；

S2，对第一个像素点(-10,-10)，测量该处的照度大小，定义接收器的半径为 D，接收器的中心位于该像素点，即接收器中心坐标为(-10,-10)，如果该光子距离接收器中心的距离小于D，就能够被接收器有效接收，并获得该光子的信息；

S3，对于(-10,-10)这一像素点，同样也是此刻接收器中心的位置，遍历接收平面的所有光子，统计能满足步骤二中能被接收器有效接收的光子，将不能被接收器接收的光子的能量设置为0，方便后续的计算；

S4，对每一个光子进行分析，显然，当光子的方向与该平面的法向所呈角为θ时，根据光通量的概念，即单位时间内通过某一曲面能量流，假设该光子的能量为e₀，实际被探测器接受到的能量值为：

S5，对该平面所有处理好的光子进行能量的求和，未能被探测器接受到的光子，其能量值已经设置为0，因此不会对接受到的总能量做出贡献，因此该步骤累加的总能量就是该次接收的总通量，记为此处的光通量；在规定完照度计口径之后，利用公式E= dΦ/dS就能通过该像素点的光通量转换成光照度的大小，其中Φ为照度计接收光子的总通量，S为光子投射的面积。

6.如权利要求1所述的一种水下光电成像的PSF与BSF仿真模拟方法，其特征在于，所述步骤2的具体过程为：

S21，两种模拟发射器分别发射光子包；对朗伯光源而言，以法线方向(0.0.1)为I₀方向，对朗伯发射体发射的光子能量进行初始化，规定方向发射的光子能量大小为W₀，则其他发射方向的光子能量为:

W₀ = W₀×cosθ

向规定的每个方向发射一定数量光子的光子包；

对激光发射器而言，对每一个发射的光子都赋予一个初始能量W₀，向规定的每个方向发射一定数量光子的光子包；

S22，光子通过设定好参数的水体，不断通过随机步长和随机散射方向，更新每个光子的传播位置和能量；

单个光子随机步长s与衰减系数c的关系是：

其中，s为单个光子的随机步长，τ为[0,1]间随机数，c为提前设定的水体衰减系数，可由散射系数与吸收系数相加得到；

在蒙特卡罗模拟过程中，光子会与介质中粒子发生碰撞，每次碰撞后必须根据散射相函数通过抽样确定出散射方向，采用的Henyey-Greenstein散射相函数的解析形式可表示为：

其中，θ为光子下次动作时与z轴方向的夹角；g为先前输入的不对称因子，当g为0时，水体各向同性；rnd为自动产生的(0,1)区间内随机数；

S23，两种模拟接收器分别接收光子；对准直辐亮度计的接收而言，可以测量光场分布，以1mm 为步进，先固定像素点的 y 坐标，依此将 x 坐标从-10cm 处遍历到 10cm 处，再将y 坐标增加值-9.9cm 处进行 x 坐标的遍历，直至遍历到最后一个像素点(10,10)，将得到每一个像素点测量得到的光源亮度的大小E₁(x,y)，至此测出了该平面所测的光源的亮度分布，也即该平面上的光场分布；

对余弦探测器而言，以1mm 为步进，先固定像素点的y坐标，依次将x坐标从-10cm 处遍历到 10cm 处，再将y坐标增加值-9.9cm 处进行 x 坐标的遍历，直至遍历到最后一个像素点(10,10)，将得到每一个像素点测量得到的光照度的大小E₂(x,y)，至此测出了该平面的照度分布，也即该平面上的光场分布。

7.如权利要求1所述的一种水下光电成像的PSF与BSF仿真模拟方法，其特征在于，所述步骤3中对经由余弦发射器发射，准直辐亮度计获取的光亮度E₁，在获取光亮度后，就知道了该平面上的光场分布，进而可以得到了PSF 函数的值，计算公式为：

其中，I₀是光源法线方向的发光强度，由于光源是朗伯点光源，光源发出的总光通量的大小等于光源的最大出射发光强度I₀乘以π立体角：

因此，在已知光源的光通量的情况下可以求出朗伯光源法线方向的发光强度I₀，实际上，由于仿真过程中，光子是限定在一秒内全部发射，因此，此处的光源的光通量Φ就是光源的光功率P；同时，亮度计是通过像面照度来反应实际物体的亮度，因而可以利用测得的像面的光照度，进而推算出光源的亮度，当选定了物镜之后，焦距随之确定下来，像面的照度值与光源的亮度成正比关系，而且当物镜焦距f<<物距 m 时，通过公式可以推导简化出如下关系：

其中，L为光源亮度，D为物镜孔径，f为物镜焦距；

进而可以导出L(x,y)，L(x,y)表示平面上点(x,y)测出的光源光亮度，PSF(x,y)表示平面上点(x,y)的点扩散函数的值。

8.如权利要求1所述的一种水下光电成像的PSF与BSF仿真模拟方法，其特征在于，所述步骤3中对经由准直光源发射，余弦探测器获取的光照度E₂，可以计算得出BSF函数的值，可以由下式给出：

其中 P 是光源的总功率，功率的定义是单位时间内能量值的大小，因此只要对发射时的光子的能量进行求和就能算出总功率，E(x,y)表示平面上像素点(x,y)处的光照度大小，BSF(x,y)表示平面上点(x,y)的光束扩散函数的值。

Images