0 引言

水下光学无线通信^[1-3](Underwater Optical Wireless Communication, UOWC)尤其是无线激光通信技术，由于其保密性好、通信速率高等优点，近年来成为水下无线通信的重要研究方向.但由于水下环境的复杂性，以及不同水域之间的差异，激光在水下传输过程中，易受环境影响，且由于目标光源较小，很难与接收器实现精确对准，造成通信质量不佳或者无法实现通信^[4-6].因此需要对激光光束进行快速稳定的跟踪和对准.

在自由空间光通信中，对激光光束的捕获、对准和跟踪(Acquisition, Tracking and Pointing, ATP)系统的研究工作已经较为成熟^[7-9].在长距离的空地通信，以及星地通信间，ATP控制系统主要包括粗瞄准控制和精瞄准控制^[10].粗瞄准实现对目标的捕获以及粗跟踪，一般由较大的光学天线以及两轴伺服转台组成；精瞄准实现对目标的精确瞄准及跟踪，一般由快速倾斜镜、压电陶瓷、以及高帧率CCD相机组成^[11].但应用于卫星间的ATP系统设备结构复杂、成本高，且其接收端光学天线一般较大，很难直接应用于水下激光通信中.近年来，光电跟踪技术结合水下无线激光通信技术也应用于水下目标跟踪，可以实现接收端信号的稳定接收^[12].但在水下激光通信过程中，采用的光电探测器接收面积较小，激光在水下长距离传输过程中，由于光斑扩散，光斑直径远大于探测器直径，探测器的接收位置会对接收到的信号质量造成较大影响，因此在接收端需要更精确地对准控制来实现信号的接收.

研究者针对发光二极管(Light Emitting Diode，LED)光通信提出了一些对准解决方案.如PONTBRIAND C^[13]等通过增加接收器(光电倍增管)的面积来增加接收机的视野，在光通信过程中由于接收机可接收范围增大，避免了需要主动跟踪发射端.还有研究^[14]提出使用多个LED或多个光电二极管，系统通过增加发射器或接收器的数量实现通信链路的稳定，但这将导致更大的占用空间、更高的成本和更高的复杂性.Al-Rubaiai^[15-16]开发了一种新的基于LED的光通信系统，基于探测器在接收不同角度入射光束的接收效率不同的关系，采用扩展卡尔曼滤波算法(Extended Kalman Filter，EKF)与反馈控制进行主动对准，在文献[16]中提出了基于EKF的对准控制算法，并且初步验证了算法的有效性.但LED系统通信速率较低，传输距离较近，为了增加传输距离和传输速率，需要使用水下半导体激光器(Laser Diode，LD)无线光通信系统.而目前水下无线光通信的对准解决方案主要集中在LED通信系统的应用中，LD光源发散角和光斑均较小，LD光斑强度分布的高斯标准差更小，接收视场角小，使得接收端的跟踪对准难度增加.

本文主要研究水下无线激光在海水中传输时，探测器在光斑不同位置的接收效率及接收光强分布规律.根据接收端光强分布特性采用EKF算法，用两个连续过程的光强度测量值进行状态估计，得到接收器位置与最大光强处的距离，然后将这个距离信息反馈给接收端的控制器，根据估计的相对值跟踪最大光强.

1 接收功率模型及状态空间模型

1.1 接收功率模型

首先利用文献[17-18]中的蒙特卡洛方法，以高斯光作为光源，基于米氏散射理论，仿真了水下激光通信信道中光子的传输特性.根据接收到的光子权重确定接收位置的光强.通过对通信链路为10 m、15 m、20 m、25 m时接收端光子分布进行仿真分析，得到在不同传输距离下，归一化接收端光强在一维空间与接收机感光面距光斑中心距离r的分布关系，并用高斯函数对其进行拟合，如图 1.表 1为采用的蒙特卡洛数值仿真参数.

图 1. The normalized intensity distribution under different transmission distances

Fig. 1. The normalized intensity distribution under different transmission distances

下载图片 查看所有图片

图1中拟合表达式为

1 $ f\left( r \right) = a \times {{\rm{e}}^{ - {{\left( {\frac{r}{b}} \right)}^2}}} $

式中，a和b为不同距离下的拟合参数，其取值见表 2.

考虑接收机位置对接收光强的影响，结合文献[19-20]中在不考虑接收机位置的情况下(只考虑光斑扩展及衰减)的接收功率模型P，得到接收光强模型为

2 $ {P_{\rm{r}}} = P \times f\left( r \right) $

式中，P_r表示在距离r时接收机的接收光功率.

因此接收机在光斑中心(即r=0)时，具有最大的接收功率.控制接收端跟踪到光斑的中心位置，如图 2所示.

图 2. Alignment diagram of receiver and transmitter of UOWC system

Fig. 2. Alignment diagram of receiver and transmitter of UOWC system

下载图片 查看所有图片

在接收端，当光到达接收器时产生光电流，通过信号处理，最终得到的完整信号模型为

3 $ {V_{\rm{d}}} = C \times {P_{\rm{t}}} \times a_{\rm{r}}^2 \times \exp \left( { - c \times d} \right) \times f\left( r \right)/{\left( {d \times \tan \theta + {a_{\rm{t}}}} \right)^2} $

式中，a_r为接收器的接收孔径半径，a_t为发射器的发射孔径半径，θ为激光光源的光束发散角，c为海水的总衰减系数(包括吸收和散射效应)，d表示通信距离，P_t表示发射端的发射功率值, V_d表示接收到的电压信号，C是比例常数，取决于滤波器和放大器电路相关参数.

1.2 状态空间模型

当接收机中心与光斑中心的距离为r时，控制r来使得接收的性能最优.因此，将状态变量定义为

4 $ \left\{ {\begin{array}{*{20}{l}} {{x_1} = C \times {P_{\rm{t}}} \times a_{\rm{r}}^2 \times \exp \left( { - c \times d} \right)/{{\left( {d \times \tan \theta + {a_{\rm{t}}}} \right)}^2}}\\ {{x_2} = r} \end{array}} \right. $

式中，x₁表示接收器在接收端接收到的最大信号电压，x₂表示接收器与光斑中心的相对位置.在端对端水下激光通信应用场景中，假设发射端固定，接收器相对发射端在r方向上进行缓慢移动，即为动力学中的准静态过程.因此，式(4)中定义的状态的动力学可以在离散时间域中表示为

5 $ \left\{ {\begin{array}{*{20}{l}} {{x_{1,k}} = {x_{1,k - 1}} + {\omega _{1,k - 1}}}\\ {{x_{2,k}} = {x_{2,k - 1}} + {u_k} + {\omega _{2,k - 1}}} \end{array}} \right. $

式中，ω_{1, k－1}和ω_{2, k－1}是系统的过程噪声，假设为独立的高斯白噪声，这些噪声项在某种程度上可理解成x₁和x₂的缓慢变化；x_{1, k}，x_{2, k}表示k时刻噪声干扰下x₁、x₂的离散值；u_k为接收器的相对位移.

将第k个测量值V_{d, k}用状态变量表示，假设包含一个与过程噪声无关的加性高斯白噪声v，式(3)可表示为动力学方程

6 $ {V_{{\rm{d}},k}} = {x_{1,k}}f\left( {{x_{2,k}}} \right) + {v_k} $

根据测量值，目标是估计x_{1, k}和x_{2, k}，通过控制u_k变量，将x₂趋向0，即具有最大接收光强度的位置.

2 估计与控制算法

由于测量模型式(6)是非线性的，采用离散时间EKF^[21]进行估计.根据式(5)，系统状态可写成

7 $ {\mathit{\boldsymbol{x}}_k} = \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{x}}_{k - 1}} + \mathit{\boldsymbol{B}}{u_{k - 1}} + {\mathit{\boldsymbol{\omega }}_{k - 1}} $

式中，$\boldsymbol{A}=\left[ \begin{array}{*{35}{l}} 1 & 0 \\ 0 & 1 \\ \end{array} \right],\boldsymbol{B}=\left[ \begin{array}{*{35}{l}} 0 \\ 1 \\ \end{array} \right],{{\boldsymbol{\omega }}_{k-1}}={{\left[ \begin{array}{*{35}{l}} {{\omega }_{1,k-1}},{{\omega }_{2,k-1}} \\ \end{array} \right]}^{\text{T}}}$

为了保证系统的可观测性，假设系统的线性化输出矩阵为H, 根据可观测性的性质有：H秩必须为2，因此观测时需要至少2次独立的光强度测量^[16].采用扫描技术，接收器通过限定的距离L = L₁, L₂, L₃…L_k移动，获得光强度测量值.根据图 1拟合的接收位置与接收效率的归一化关系，考虑实际的接收端灵敏度以及光斑的几何扩展，预设用于移动的距离L= －3, －6, －9, －12, －9, －6, －3, 0, 3, 6, 9, 12, 9, 6, 3, 0 (mm), 在每次迭代中，从该阵列中顺序选择一个L_k，在L_k和L_k－1处的两个测量形成输出向量y，得到观测方程，即

8 $ {\mathit{\boldsymbol{y}}_k} = \left[ {\begin{array}{*{20}{c}} {{x_{1,k}}f\left( {{x_{2,k}} + {L_k}} \right) + {v_k}}\\ {{x_{1,k - 1}}f\left( {{x_{2,k - 1}} + {L_{k - 1}}} \right) + {v_{k - 1}}} \end{array}} \right] $

通过式(6)和(8)，采用EKF算法估计系统状态，算法控制流程如图 3所示.

图 3. Algorithm control flow chart

Fig. 3. Algorithm control flow chart

下载图片 查看所有图片

在EKF算法中，有三个协方差矩阵，即P、Q和R. P是条件误差协方差矩阵，P_k|k－1表示协方差矩阵的预测. P_k|k－1需要初始化为正定矩阵.最初的x₂的估计值$\overset{\wedge }{\mathop{{{x}_{2}}}}\, $的初值设为0.定义初始估计值$\overset{\wedge }{\mathop{{{x}_{1}}}}\, $为最大强度值的2/3.Q是过程噪声协方差矩阵，为正定矩阵.除此之外，由于假设噪声项是不相关的，因此Q应该是对角矩阵.R是测量噪声协方差矩阵，由于测量的噪声被假定为白噪声，因此R矩阵也是对角矩阵.在初始化这些量之后，算法开始迭代，过程为：

1) 对状态变量进行预测，式(7)可写为

9 $ {{\mathit{\boldsymbol{\hat x}}}_{m,k|k - 1}} = \mathit{\boldsymbol{A}}{{\mathit{\boldsymbol{\hat x}}}_{m,k - 1|k - 1}} + \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{u}}_{k - 1}} $

10 $ {\mathit{\boldsymbol{P}}_{k|k - 1}} = \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{P}}_{k - 1|k - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}} + {\mathit{\boldsymbol{Q}}_k} $

式中，${{\overset{\wedge }{\mathop{\boldsymbol{x}}}\, }_{m, \text{ }k\text{ }|\text{ }k\text{-1}}}$表示第m个状态k时刻的先验状态估计值，是滤波的中间计算结果，即根据上一时刻(k－1时刻)的最优估计预测的k时刻的结果.

2) 估计输出为

11 $ {{\mathit{\boldsymbol{\hat y}}}_k} = \left[ {\begin{array}{*{20}{c}} {{{\hat x}_{1,k|k - 1}}f\left( {{{\hat x}_{2,k|k - 1}} + {L_k}} \right)}\\ {{{\hat x}_{1,k|k - 1}}f\left( {{{\hat x}_{2,k - 1|k - 1}} + {L_{k - 1}}} \right)} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{\hat y}_{1,k}}}\\ {{{\hat y}_{2,k}}} \end{array}} \right] $

把式(9)代入式(11)可得

12 $ {{\mathit{\boldsymbol{\hat y}}}_k} = h\left( {{{\hat x}_{1,k|k - 1}},{{\hat x}_{2,k|k - 1}}} \right) = \left[ {\begin{array}{*{20}{c}} {{{\hat x}_{1,k|k - 1}}f\left( {{{\hat x}_{2,k|k - 1}} + {L_k}} \right)}\\ {{{\hat x}_{1,k|k - 1}}f\left( {{{\hat x}_{2,k - 1|k - 1}} - {u_{k - 1}} + {L_{k - 1}}} \right)} \end{array}} \right] $

则H矩阵为

13 $ {\mathit{\boldsymbol{H}}_k} = \frac{{\partial h\left( {{{\hat x}_{1,k|k - 1}},{{\hat x}_{2,k|k - 1}}} \right)}}{{\partial {{\mathit{\boldsymbol{\hat x}}}_{k|k - 1}}}} $

3) 更新

14 $ {\mathit{\boldsymbol{K}}_k} = {\mathit{\boldsymbol{P}}_{k|k - 1}}\mathit{\boldsymbol{H}}_k^{\rm{T}}{\left( {{\mathit{\boldsymbol{H}}_k}{\mathit{\boldsymbol{P}}_{k|k - 1}}\mathit{\boldsymbol{H}}_k^{\rm{T}} + {\mathit{\boldsymbol{R}}_k}} \right)^{ - 1}} $

15 $ {{\mathit{\boldsymbol{\hat x}}}_{k|k}} = {{\mathit{\boldsymbol{\hat x}}}_{k|k - 1}} + {\mathit{\boldsymbol{K}}_k}\left( {{\mathit{\boldsymbol{y}}_k} - {{\mathit{\boldsymbol{\hat y}}}_k}} \right) $

16 $ {\mathit{\boldsymbol{P}}_{k|k}} = \left( {\mathit{\boldsymbol{I}} - {\mathit{\boldsymbol{K}}_k}{\mathit{\boldsymbol{H}}_k}} \right){\mathit{\boldsymbol{P}}_{k|k - 1}} $

式中，${{\overset{\wedge }{\mathop{\boldsymbol{x}}}\, }_{\text{ }k\text{ }|\text{ }k}}$和P_k|k是更新之后的值.

4) 得到控制变量为

17 $ {u_k} = - \hat G{x_{2,k|k}} $

式中，G为增益系数.将控制变量u_k反馈到接收端将x₂驱动到零，使接收器在短时间内可以接收到最大功率，进而得到最大测量值.

3 实验及仿真结果分析

通过在实验装置上实现与仿真相近的场景来验证接收端光斑光强分布情况.并根据蒙特卡洛仿真的光强分布曲线进行算法仿真，验证算法的有效性.

首先搭建水下光通信实验平台.发射端采用532 nm，发射功率80 mW的激光器，产生高斯光源，如图 4(a)所示.搭建长度为30 m的实验水槽，并装入自来水，水质衰减光在水中传播路径如图 4(b)所示；接收端使用CGU2-500C，CCD相机，利用USB直接与计算机连接，采集光斑图像，可实时观测光斑并拍摄光斑图像，如图 4(c)所示.实验研究在同一水质情况不同距离(15 m，20 m，25 m)下光斑能量分布.

图 4. Spot detection experiment system

Fig. 4. Spot detection experiment system

下载图片 查看所有图片

接收端拍摄的光斑图像如图 5(a)所示.对拍摄到的光斑进行图像处理，沿通过光斑中心的横线做出光强在一维空间的分布，以灰度值的变化表示接收光强的变化，得到像素位置与归一化光强的曲线如图 5(b)所示.根据图 5(b)实验测量的光强分布结果，进行高斯曲线拟合如图 5(c)所示.

图 5. Spot detection experiment system

Fig. 5. Spot detection experiment system

下载图片 查看所有图片

图 5(c)的拟合曲线表达式为

18 $ F\left( x \right) = a \times \exp \left[ { - {{\left( {\frac{{x - b}}{c}} \right)}^2}} \right] $

式中，x表示像素的位置坐标；a、b、c为拟合表达式的参数.用F(x)表归一化后相对接收光强相对于像素位置坐标的函数.

从图 5及式(18)可以看出，在不同距离下光斑能量基本都呈高斯分布，与蒙特卡洛方法仿真得到的规律一致.因此当接收机感光面半径小于光斑半径时，接收机如果不在光斑中心位置接收，那么接收光强将有较大的损失.

因此根据模型式(3)的光强分布进行跟踪算法仿真.系统的状态如式(5)所示，f(r)为传输距离10 m时的函数.扩展卡尔曼滤波初始化参数设置见表 3.

系统的过程噪声ω_{1, k－1}和ω_{2, k－1}的方差分别取0.01和0.001.图 6显示了状态估计$\overset{\wedge }{\mathop{{{x}_{2}}}}\, $与实际状态x₂的变化.

图 6. The convergence of the estimate state and the actual state

Fig. 6. The convergence of the estimate state and the actual state

下载图片 查看所有图片

从图 6中可以看出，估计的状态在大约2 s内收敛到原始状态的邻域并且在整个运行期间趋于稳定.即实际平均值(状态x₂)在大约2~3 s内收敛到零附近，表明接收器在2~3 s内到达光束中心实现对准.

对15 m、20 m、25 m的实验拟合的光强分布模型进行跟踪算法仿真.图 7给出了三种距离下状态x₂的变化，可以看出，在不同距离下系统模型的x₂都可以稳定地收敛到零附近.

图 7. Change of state x₂ at different transmission distances

Fig. 7. Change of state x₂ at different transmission distances

下载图片 查看所有图片

通过计算实际状态(x₂)的标准偏差来评价在不同距离下的对准精度，并计算在此误差下接收光强的效率η，效率计算公式为

19 $ \eta = \frac{{\frac{1}{k}\sum\limits_{k = 1}^n {{x_{1,k}}f\left( {{x_{2,k}}} \right)} }}{{{x_1}f\left( {{x_2} = 0} \right)}} \times 100\% $

图 8给出仿真得到不同传输距离对准误差的变化.表 4列出对准误差、接收效率与传输距离的关系.

图 8. Alignment error at different transmission distances

Fig. 8. Alignment error at different transmission distances

下载图片 查看所有图片

从图 8和表 4可以看出，随着传输距离的增加，对准的平均误差逐渐增大，传输距离小于25 m时，误差小于2 mm，并且接收端处于对准状态时接收效率仍可以达到98%以上.

4 结论

本文首先根据蒙特卡洛仿真结果确定接收端在光斑不同位置的接收效率，得到在不同传输距离时接收效率与接收位置的关系，并通过实验确定水下激光光斑的光强分布特性.根据光强分布特性提出基于EKF的对准控制算法用于水下无线激光通信系统，基于接收端接收光强模型建立状态空间模型，对系统状态进行估计，并结合基本反馈理论用于控制接收器的移动，增加接收端的接收效率.实验与蒙特卡洛仿真表明，高斯激光光束在水下经过较长距离传播后，其接收端光强仍然近似为高斯分布.非线性估计算法仿真表明，基于EKF的估计算法对于水下激光通信系统有较高的精度.并且在长距离的水下激光通信中，接收端对准后的接收效率可达98%，大大提高了水下无线LD光通信接收系统的稳定性以及接收效率.提出的EKF算法以及光学模型可应用于远距离无线激光通信接收端的跟踪对准系统，提高水下无线激光通信系统的实用性.