基于FPGA的自适应光学SPGD控制算法的实现方法与流程

本发明属于自适应光学技术领域，涉及到一种无波前探测自适应光学系统，特别是一种基于fpga的自适应光学spgd控制算法的实现方法。

背景技术：

美国陆军研究院的m.a.vorontsov提出随机并行梯度下降(stochasticparallelgradientdescent:spgd)算法，并将此算法应用到自适应光学领域。基于spgd的无波前探测自适应光学系统不需要进行波前测量和相位重构，降低了系统的复杂性。进一步扩展了自适应光学技术的应用领域，尤其是常规自适应光学系统不能应用的领域。

校正速度是自适应光学系校正能力的重要指标之一，spgd控制算法通过迭代的方式寻找最优解。因此，算法的迭代速度对自适应光学系统校正能力起着关键的作用，而算法的实现方式对迭代速度又有着决定性的影响。另外，应用环境的多样性要求算法的实现方式具有一定的通用性。因此spgd算法的实现方式需要同时满足实时性与灵活性两个要求。

自适应光学系统中spgd控制算法的实现方式有两种：一是使用软件编程的方法在通用pc机上实现，该方法具有实现简单、灵活性强等优点，但pc机运行的windows是一款微内核、抢占式多任务、软实时的操作系统，按照线程的优先级高低顺序来进行任务的调度，从而导致整个系统的实时性较弱。二是通过硬件的方式实现。m.a.vorontsov等人使用模拟超大规模集成电路(analogvlsi)实现spgd控制算法，该方法迭代速度快，实时性强，满足自适应光学系统的要求，但是电路设计复杂，芯片成本高，而且缺乏灵活性。

张金宝、张超等人使用fpga实现spgd控制算法。fpga是一类可编程的集成电路，通过硬件描述语言表述其功能，可以将其配置成能够完成特定功能(如spgd控制算法)的硬件电路，具有可重复编程、算法实现灵活等特点，而其专用硬件电路的本质又能够保证实时性的要求。但是，现有的fpga实现spgd控制算法的方式只是针对特定变形镜，应用范围有限。另外，现有技术没有对串行传入fpga内部的图像数据进行矩阵恢复，这将导致部分性能指标函数无法用fpga实现或者实现的方式不精确，降低了fpga实现spgd控制算法方式的通用性。同时，现有技术性能指标函数实现方式不精确，这将影响自适应光学系统的校正能力。

技术实现要素：
：

本发明的所需要解决的技术问题是克服现有技术的不足，提出一种基于fpga的无波前探测自适应光学spgd控制算法的实现方法，其应用范围广适用于任何单元数的自适应光学系统，提高了fpga实现spgd控制算法方式的通用性，保证了无波前探测自适应光学系统的校正能力。

本发明所要解决的技术问题是通过以下的技术方案来实现的。本发明是一种基于fpga的无波前探测自适应光学spgd控制算法的实现方法，其特点是，其步骤如下：

(1)随机扰动扰动模块，借助软件(包括但不限于matlab软件)产生多组伪随机序列，在fpga中将每一组序列乘以扰动电压并添加相应的符号，生成正扰动和负扰动；

(2)选取远场光斑的平均半径mr-meanradius作为性能指标函数,性能指标函数不限于mr；

其中(x0,y0)是光斑的质心；i(x,y)表示对应坐标(x,y)处的像素灰度值；光斑的roi区域是n×n，所以x,y∈[1,n]；

(3)计算变形镜电压控制模块，根据公式u⁽ⁿ⁺¹⁾＝u⁽ⁿ⁾-γδu⁽ⁿ⁾δj⁽ⁿ⁾计算控制电压，性能指标函数mr往极小值方向优化，γ取正值。

本发明所述的基于fpga的无波前探测自适应光学spgd控制算法的实现方法，其进一步优选的技术方案是：

spgd算法中随机扰动电压的特性，对借助matlab软件(不限于matlab)产生多组伪随机序列提出了相应的要求：

a、路数等于变形镜单元数；一个n单元的变形镜需要n路伪随机序列；

b、伪随机序列两两相互独立；相互独立可避免变形镜的各驱动单元间相互耦合；

c、伪随机序列符合伯努利分布，两个样本值出现的概率各为0.5。

spgd控制算法通过对多路的控制参数加入随机并行扰动，使用性能指标测量值的变化量与控制参数的变化量进行控制参数的搜索。

本发明所述的基于fpga的无波前探测自适应光学spgd控制方法，该方法采用以下步骤实现spgd控制算法：

(1)借助软件生成随机数序列，并存入fpga的开发板；

(2)fpga读取随机序列并生成扰动电压{δuⁿ},(j＝1,...,n,|δu⁽ⁿ⁾|＝0.2,sgn(δu⁽ⁿ⁾)＝±1)；

(3)fpga向变形镜施加正扰动u⁽ⁿ⁾+δu⁽ⁿ⁾，u⁽ⁿ⁾是第n-1迭代计算出的控制电压，采集图像数据，计算性能指标函数

(4)fpga向变形镜施加负扰动u⁽ⁿ⁾-δu⁽ⁿ⁾，采集图像数据，计算性能指标函数

(5)求出步骤(3)、(4)性能指标函数的变化量根据变化量计算梯度γδu⁽ⁿ⁾δj⁽ⁿ⁾，根据公式u⁽ⁿ⁺¹⁾＝u⁽ⁿ⁾-γδu⁽ⁿ⁾δj⁽ⁿ⁾完成变形镜控制电压的计算；

重复步骤(3)～(5)，进行多次迭代，直到寻找到变形镜的最佳控制电压，使得初始像差(ψ(r))与变形镜的校正相位(m(r))的和趋近于0。

与现有技术相比，本发明具有以下有益效果：

1.本发明利用fpga实现spgd控制算法，实现方便且具有一定的灵活性，又能够保证实算法实时性的要求。

2.本发明利用图像数据存储在fpga内部rom核的地址位，设定判断条件，对图像数据进行矩阵恢复，常用的性能指标函数都能用fpga准确实现，提高fpga实现spgd控制算法方式的通用性，也不会因为性能指标函数的计算影响基于fpga的spgd控制算法的无波前探测自适应光学系统的校正能力。

3.本发明在矩阵恢复完成后，用两个计数器和一些辅助信号对图像数据进行精确定位，在fpga内部规定图像数据与矩阵坐标之间的关系，按照质心公式完成质心计算过程的编译，准确的计算出图像矩阵的质心。本发明用fpga实现mr计算模块的方式精确性更高，基于fpga的spgd控制算法的无波前探测自适应光学系统的校正能力更强。

4.通过理论分析，可以知道增益系数和扰动幅度都是正常数。本发明在预处理部分将两者相乘，结果为一个正常数，在fpga内部直接定义这个常数，减少了乘法器的使用，节约fpga的内部资源，同时减少了算法的复杂性，一定程度上提高了算法的计算速度。

附图说明

图1为生成随机扰动的时序图；

图2为计算性能指标函数mr的流程图；

图3为性能指标函数mr分子分步计算的时序图；

图4为计算变形镜控制电压的流程图；

图5为计算变形镜控制电压部分时序图；

图6为基于fpga的无波前探测自适应光学spgd控制算法的实验系统的结构框图。

具体实施方式

以下参照附图，进一步描述本发明的具体技术方案，以便本领域的技术人员进一步理解本发明。

以下实施方式以97单元变形镜为例：

一生成随机扰动电压模块：

图1为本发明随机扰动模块的时序图，clk为系统时钟，具体实施时按如下步骤进行。

1.预处理部分，借助matlab软件但不限于matlab，利用rand函数生成若干组97路伪随机序列(由0,1组合而成)。将生成的若干组97路伪随机序列存储备用。

2.以其中一组伪随机序列为例计算正负扰动，将这组伪随机序列读取并存储到fpga的rom核中。设置计数器cnt，并将cnt作为地址位读取保存在rom里面的伪随机序列(cnt等于0时，对应的是97路伪随机序列的第1个值，依次类推)。

3.计数器cnt的数值与97路伪随机序列的对应关系如图1所示，q表示从rom里面读取的随机序列，cnt等于0时，对应的q会在cnt等于2的时候输出(从rom里面读取数据有一定的延迟)。图中两条虚线之间的q就是从rom里面读取的97路伪随机序列。

4.q等于0时，对应的扰动为1，q等于1时对应的扰动为-1，具体结果如图1的信号d所示，共97个1，-1的随机组合。由于要先判断q的值，所以对应的d值会延迟，借助一定的辅助信号，确定两个点划线(标注为开始生成扰动和完成97路扰动)之间的d为所需要的97路随机扰动，扰动幅度为1。

5.在步骤4的基础上，将d乘以200，得到图1中的dv1(正扰动)。200就是扰动幅度，扰动幅度定义为后10位是小数位的二进制定点小数，所以200对应的实际扰动幅度为0.1953(近似为0.2)。

6.在dv1的基础上每个扰动加上负号得到图1中的数据dv2(负扰动)。

二性能指标的计算：

图2为性能指标函数mr的计算流程图。具体实施时，按如下步骤进行：

1.采集图像数据(光斑图像的灰度值，roi区域为n×n)，由于图像数据是串行传输到fpga内部的，所以计算之前需要进行矩阵恢复。以地址位作为判断条件将图像数据恢复成n×n的矩阵形式。

2.计算质心，质心的公式是根据公式质心的计算按以下步骤进行计算：

(1)按照公式(2)(3)计算坐标值与之对应的灰度值的乘积：

式中向量x1表示的是横坐标1与矩阵第1行的所有元素(图像灰度值)的乘积，依次类推，向量x2～xn表示的都是横坐标与其对应的所有元素的乘积。向量y1表示纵坐标1与矩阵第1列的所有元素的乘积，依次类推，向量y2～yn表示的都是纵坐标与其对应的所有元素的乘积。

(2)将向量x1～xn内的所有元素相加得到公式∑x,yxi(x,y)的计算结果为sun_xi，将向量y1～yn内的所有元素相加得到公式∑x,yyi(x,y)的计算结果为sun_yi。

(3)将矩阵所有元素(图像灰度值)相加得到公式∑x,yi(x,y)的计算结果为sun_q。

(4)利用除法ip核进行两次除法运算，被除数为sun_xi，除数为sun_q得到质心坐标x0，被除数为sun_yi，除数为sun_q得到质心坐标y0。

3.计算mr的分子式中i(x,y)表示的坐标为(x,y)的灰度值。性能指标函数mr分子分步计算的流程图如图2所示，部分时序图如图3所示，具体计算过程如下：

(1)根据时钟clk生成另一个时钟sclk，sclk是时钟clk的n分频，两者的关系如图3所示，sclk时钟的1个周期包含n个clk时钟周期。

(2)以clk为时序，设计y_cnt计数器，计数器从0开始计数，clk的上升沿触发，计数器加1，当计数到n时，y_cnt从0开始重新计数。以sclk为时序，设计x_cnt计数器，计数器从0开始计数，当计数到n时，停止计数。两个计数的部分计数过程如图3的x_cnt和y_cnt所示。两个计数器主要是为了对矩阵元素(图像灰度值)进行定位，x_cnt定位横坐标，y_cnt定位纵坐标，矩阵的横坐标x等于x_cnt的数值加上1，矩阵的,纵坐标y等于y_cnt的数值加上1，如x_cnt等于0，y_cnt等于1表示矩阵第1行第2列的元素。

(3)计算(x-x0)²-(y-y0)²，先计算x-x0，y-y0，然后按照公式(x-x0)×(x-x0)-(y-y0)×(y-y0)计算，运算结果会延迟1个时钟周期输出，如x等于1，y等于3(对应的位置为x_cnt等于0，y_cnt等于2)，计算的出的结果会在x_cnt等于0，y_cnt等于3定位的位置输出。矩阵第1行所有位置的计算结果在a和b两条虚线间输出。

(4)计算将(3)计算出的结果输入开方ip核，延迟一个时钟周期输出运算结果sqr_q。综合(3)和(4)计算过程，时钟共延迟了2个时钟周期。如x等于1，y等于1(对应的位置为x_cnt等于0，y_cnt等于0)计算的出的结果会在x_cnt等于0，y_cnt等于2的位置输出。矩阵第1行所有坐标的计算结果在c和d两条点划线间输出。

(5)按照时钟clk的时序，设置计数器address，开始计数的使能信号和y_cnt计数器保持一致，address计数器部分时序图如图3所示。address计数器的功能是从rom核中读取图像灰度值。从地址位address的发出，到输出数据有两个时钟周期的延迟，address等0的图像灰度值，会在address等于2的位置被读出。由于address计数器和y_cnt计数器的关系，address等0的位置定位是x_cnt等于0，y_cnt等于0，输出的灰度值定位在x_cnt等于0，y_cnt等于2。具体对应关系如图3所示。

(6)从图3可以公式看出的计算结果和坐标为(x,y)的灰度值是两个完全同步的信号，可以直接将两个结果相乘，完成公式的计算过程。得到的结果记为向量sqr_x1。

(7)按照上述的计算过程可以算出第2行至第n行的所有的结果(sqr_x2～sqr_xn)，然后所有向量里面的元素结果相加得出的结果是mr的分子。

4.将步骤(3)中的(7)计算的结果与矩阵所有元素(图像灰度值)的和相除，计算得到mr。

三变形镜控制电压的计算：

变形镜控制电压的计算公式为:u⁽ⁿ⁺¹⁾＝uⁿ-γδuⁿδj⁽ⁿ⁾，式中γδuⁿδj⁽ⁿ⁾为梯度估计，γ为增益系数，δuⁿ为随机扰动，δj⁽ⁿ⁾为性能指标的变化量性能指标函数mr是往极小值方向寻找最优解，所以当δuⁿ与δj⁽ⁿ⁾的符号一致时，减梯度估计结果，当δuⁿ与δj⁽ⁿ⁾的符号相反时，加梯度估计结果。本发明按照公式(4)实现变形镜控制电压的计算。

u⁽ⁿ⁺¹⁾＝uⁿ-γδuⁿδj⁽ⁿ⁾(4)

式中：γ为增益系数取正值。控制电压的计算模块如图4所示，具体实施步骤如下：

1.正扰动施加至变形镜后，采集图像数据，利用性能指标函数计算模块计算性能指标函数mr，记为mr1。负扰动施加至变形镜后，采集图像数据，利用性能指标函数计算模块计算性能指标函数mr2。

2.计算性能指标函数的变化量δmr＝mr1-mr2。

3.计算梯度。

4.按照公式(4)，计算控制电压，为了实现的方便，将公式转化为u⁽ⁿ⁺¹⁾＝uⁿ-asign(δuⁿ)δj⁽ⁿ⁾,其中a表示的是增益系数和随机扰动幅度(正值)的乘积，sign(δuⁿ)等于正负1，图1随机扰动的时序图中的d表示的就是不同时序sign(δuⁿ)的值。

5.计算完成控制电压后，加上第n+1次的随机扰动电压，然后再将控制电压输出至变形镜，采集图像数据。减去第n+1次的随机扰动电压，然后再将控制电压输出至变形镜，再次采集图像数据。

图5给出了控制电压计算过程的时序图，mr_done信号高电平，表示性能指标函数模块计算完成，mr_done信号触发后，时钟clk的下1个上升沿，输出mr1mr2，同步计算出mr1和mr2的差值δmr。sju_en使能信号是触发随机扰动模块，输出随机扰动模块里面的d信号。为确定d信号输出所在的位置，设计一个0～98的计数器cnt。从图5可以看出信号d从cnt等于3的时候输出，cnt记数完成后，延时1个时钟周期输出完毕，共有97个d输出。预处理部分完成随机扰动幅度和增益系数的乘积a。t_d信号(梯度)是δmr，d，a三者乘积输出的结果。同步完成当前电压减去梯度的计算过程，得到u⁽ⁿ⁺¹⁾，并保存。sju_en同时也是下一次随机扰动模块的触发信号，在cnt等于3的位置随机扰动dv开始输出，cnt记数完成后，延时1个时钟周期输出完毕。将控制电压u⁽ⁿ⁺¹⁾加上随机扰动输出给变形镜和减去随机扰动，得到该次迭代所要输出的控制电压。

四基于fpga的无波前探测自适应光学spgd控制算法的实验系统：

系统的结构框图如图6所示，成像设备(如各种接口的相机、带有coms和ccd感光芯片的图像采集系统)采集光斑图像数据，采集的图像数据通过数据线(如usb、光纤)传输并保存到fpga开发板中。用fpga实现的spgd控制算法对图像数据进行处理，计算出波前校正器的控制电压，将计算出的控制电压通过usb或光纤等数据线(取决于波前校正器接受数据的接口)传输到波前校正器(如各种单元数的变形镜、空间光调制器)，达到控制波前校正器对光斑图像进行校正的目的。