一种基于矩阵运算的光场传输数值仿真方法与流程

1.本发明属于光场传输技术领域，具体涉及一种基于矩阵运算的光场传输数值仿真方法。


背景技术：

2.随着激光技术的不断发展，激光广泛应用于工业、通信、医学和国防等众多领域。激光传输是激光系统应用的重要环节之一，主要研究激光从系统出口发射后在大气或其他介质中传输时所发生的各种物理过程。激光传输问题可以通过物理实验和数值仿真等手段研究激光传输机理，获取规律性认识。(任斌,陈纯毅,杨华民.大气湍流光波传输数值仿真技术研究综述[j].系统仿真学报,2017,29(08):1631-1640.)数值仿真以其易于实现、参数可控、经济高效等鲜明优势，受到广大研究人员的青睐。(孙全,吕品,宁禹,习锋杰,刘文广,许晓军.光学系统仿真软件seelight在自适应光学上的应用[j].光电工程,2018,45(03):130-136.)。在激光传输数值仿真方法中，角谱传播法因其清晰的物理图像而被广泛采用。然而，角谱传播法由于涉及傅里叶变换，不可避免地存在以下局限性：一是有限的离散采样频率会导致在光场突变点或边界处出现信号振荡；二是有限的离散傅里叶变换会出现混淆和拖尾效应；三是角谱传播法对计算网格参数有复杂的要求，从而导致计算不够方便灵活，甚至因网格选取不当导致错误的计算结果；四是角谱传播法要求输入光场和目标光场的网格数目必须一致，这对于只关心斯特列尔比或者桶中功率的仿真场景，计算成本较高(schmidt,jason d.numerical simulation of optical wave propagation with examples in matla[b]，(spie press monograph vol.pm199,2010))。


技术实现要素：

[0003]
本发明的目的在于提供一种基于矩阵运算的光场传输数值仿真方法，解决传统的光场传输数值仿真中的角谱传播法，在有限的离散采样频率会导致在光场突变点或边界处出现信号振荡，且对参与计算的网络参数有较为复杂的要求，从而使得计算不灵活甚至产生错误结果导致计算成本较高的技术问题。
[0004]
为达到上述目的，解决上述技术问题，本发明的技术方案如下：
[0005]
一种基于矩阵运算的光传输数值仿真方法，包括如下步骤：
[0006]
s1：根据任务需求确定光学传输系统的abcd矩阵
[0007]
s2：将目标光场的广义惠更斯-菲涅尔衍射积分转化为离散直角坐标系下的求和形式后，为转化为矩阵乘积形式做准备；
[0008]
目标光场的广义惠更斯-菲涅尔积分如公式(1)所示
[0009][0010]
其中，(xm,yn)为输入光场的笛卡尔直角坐标；m,n∈(1,n1)，(xj,yk)为目标光场的
笛卡尔直角坐标；j,k∈(1,n2)，u为目标光场，u为输入光场，；a、b、d为步骤1确定的光学传输系统abcd矩阵中的参数；
[0011]
通过数值积分将广义惠更斯-菲涅尔衍射积分转化为离散直角坐标系下的求和形式，如公式(2)所示
[0012][0013]
其中w是与数值积分有关的向量，向量w含有n1个元素，δ为输入光场直角坐标系下的网格间隔；
[0014]
数值积分方法选取梯度算法，w＝(0.5,1,1,
…
0.5)
×
δ
[0015]
也可以选择辛普森算法或其它合适的数值积分方法，对于辛普森算法：
[0016]
有
[0017]
对公式(2)进行移项操作，得到如下所示公式(3)
[0018][0019]
由于和有着高度的相似性，因此特别规定：
[0020]
对于(xm,yn)，有x1＝y1,x2＝y2,x3＝y3,x4＝y4，
…
,x
n1
＝y
n1
即输入光场u所在的笛卡尔直角坐标系是网格数量为n1×
n1的矩阵；
[0021]
对于(xj,yk)，有x1＝y1,x2＝y2,x3＝y3,x4＝y4，
…
,x
n2
＝y
n2
；即输出光场u所在的笛卡尔直角坐标系是网格数量为n2×
n2的矩阵；
[0022]
通过此规定，获得到输入光场所在的直角坐标系的横坐标间隔与纵坐标间隔相等，所以数值积分向量w对于输入光场所在的直角坐标系的横纵坐标完全相等；
[0023]
由此公式(3)转化为矩阵乘积的形式如公式(4)所示：
[0024][0025]
其中h为引入的算法矩阵，大小为n2×
n1,h的定义如下：
[0026][0027]
其中xm,yn分别为输入光场的横纵坐标向量，xj,yk分别为目标光场的横坐标向量和纵坐标向量。
[0028]
由此将广义惠更斯-菲涅尔衍射积分转化为离散直角坐标系下的求和形式；
[0029]
s3：将s1中确定好的abcd矩阵参数和与数值积分方法有关的w代入式(4)，求解h；
[0030]
s4：将输入光场u和s3中得到的矩阵h代入式3，计算得到目标光场u如下矩阵乘积的形式，如公式(6)所示：
[0031][0032]
其中，输入光场u为n1×
n1的矩阵，目标光场为u为n2×
n2的矩阵，h为n1×
n2的矩阵。
[0033]
本发明的有益效果为：
[0034]
本发明通过选取合适的数值积分算法，利用广义惠更斯-菲涅尔积分中纵坐标与横坐标的可交换性，将其转化为形式简洁的矩阵乘法运算。该方法避免了傅里叶变换所带来的吉布斯现象、混淆和拖尾效应。此外，该方法可任意选取目标光场的网格数量，具有灵活、经济的特点。
[0035]
本发明相比于传统光场传输仿真方法——傅里叶变换方法(或称角谱传播方法)具有快速、准确、简洁、灵活的优势。具体体现在如下方面：一是避免了傅里叶变换中采样频率有限导致的在光场突变点或边界出现振荡的吉布斯现象；二是避免离散傅里叶变换中频率采样不当导致的混淆和拖尾效应；三是由于数值积分求解方法的操作灵活，从而可简化光传输数值仿真过程中网格参数的选取；四是对于目标光场的计算可局限于某一点或某一区域而不必计算整个目标光场。对于重点关心光场特定物理参量——诸如斯特列尔比或者桶中功率的数值仿真场景，可进一步提升仿真速度。
[0036]
本发明可取代传统角谱传播算法用于快速精确仿真光场在真空、大气或其它介质场中的传输问题。
附图说明
[0037]
图1为本发明技术方案的原理图；
[0038]
图2为矩形光束近场光强分布灰度图；
[0039]
图3为本发明和傅里叶变换法在矩形光束真空传输场景中的数值仿真结果对比；
[0040]
图4为基模高斯光束灰度图；
[0041]
图5为本发明和傅里叶变换法在基模高斯光束真空传输场景中的数值仿真结果对比。
具体实施方式
[0042]
为了更好地理解本发明，本发明的设计原理如图1所示，图1右侧为高斯光束近场光斑图，左侧为高斯光束远场光斑图，从近场到远场的传播距离为300米，中间的式子即为广义惠更斯-菲涅尔衍射积分将其转换为矩阵乘法的详细过程。
[0043]
下面结合附图和具体实例对本发明作进一步地描述。
[0044]
一种基于矩阵运算的光传输数值仿真方法，包括如下步骤：
[0045]
s1：根据任务需求确定光学传输系统的abcd矩阵
[0046]
在几何光学中，光线、光学系统、传输介质有其特定的矩阵表示方法，传输系统变换矩阵可以表示为各个分系统变换矩阵的乘积，该方法属于矩阵光学范畴；
[0047]
s2：将目标光场的广义惠更斯-菲涅尔衍射积分转化为离散直角坐标系下的求和形式后，进一步转化为矩阵乘积的形式；
[0048]
本发明的核心理念就是不用离散傅里叶变换，从而避免傅里叶变换中采样频率有限导致的在光场突变点或边界出现振荡的吉布斯现象；同时也避免离散傅里叶变换中频率
采样不当导致的混淆和拖尾效应；为解决如上所述的技术问题，选择了广义惠更斯-菲涅尔衍射积分，但是如果直接对广义惠更斯-菲涅尔衍射积分并不是妥当的做法也是不现实的，因为计算量太大会耗费巨大的计算成本，所以本发明将设置输入光场以及目标光场的x，y点的坐标一致，然后利用其横纵坐标可互换的性质，将广义惠更斯-菲涅尔积分通过数值积分算法转化为矩阵乘法。
[0049]
因此本发明将目标光场的广义惠更斯-菲涅尔积分通过选取合适的数值积分算法，包括但不限于梯度算法、辛普森算法转化为矩阵乘法；
[0050]
目标光场的广义惠更斯-菲涅尔积分如公式(1)所示
[0051][0052]
其中，(xm,yn)为输入光场的笛卡尔直角坐标；m,n∈(1,n1)，(xj,yk)为目标光场的笛卡尔直角坐标；j,k∈(1,n2)，u为目标光场，u为输入光场，；a、b、c为步骤1确定的光学传输系统abcd矩阵中的参数；
[0053]
通过数值积分将广义惠更斯-菲涅尔衍射积分转化为离散直角坐标系下的求和形式，如公式(2)所示
[0054][0055]
其中w是与数值积分有关的向量，向量w含有n1个元素，δ为输入光场直角坐标系下的网格间隔；
[0056]
对公式(2)进行移相操作，得到如下所示公式(3)
[0057][0058]
这里注意到，和有着高度的相似性，因此特别规定：
[0059]
对于(xm,yn)，有x1＝y1,x2＝y2,x3＝y3,x4＝y4，
…
,x
n1
＝y
n1
即输入光场u所在的笛卡尔直角坐标系是网格数量为n1×
n1的矩阵；
[0060]
对于(xj,yk)，有x1＝y1,x2＝y2,x3＝y3,x4＝y4，
…
,x
n2
＝y
n2
；即输出光场u所在的笛卡尔直角坐标系是网格数量为n2×
n2的矩阵；
[0061]
通过此规定，可以得到输入光场和目标光场所在的直角坐标系的横坐标间隔与纵坐标间隔相等，所以数值积分向量w对于目标光场所在的直角坐标系的横纵坐标完全相等；
[0062]
由此公式(3)转化为矩阵乘积的形式如公式(4)所示：
[0063][0064]
其中h为引入的算法矩阵，大小为n2×
n1,h的定义如下：
[0065][0066]
其中xm,yn分别为输入光场的横纵坐标向量，xj,yk分别为目标光场的横坐标向量和纵坐标向量。
[0067]
由此将目标广场的广义惠更斯-菲涅尔衍射积分转化为离散直角坐标系下的求和
形式；
[0068]
s3：将s1中确定好的abcd矩阵参数(a、b、d)和s2中确定的与数值积分方法有关的w代入式(4)，求解h；
[0069]
例如采用梯度算法展示计算h的过程，w＝(0.5,1,1,
…
0.5)
×
δ，
[0070][0071]
δ代表网格间隔，(x1，x2，
…
，x
n1
)为输入光场的横坐标，(y1，y2，
…
，y
n1
)为输入光场的纵横坐标，(x
’1，x
’2，
…
，x’n2
)为目标光场的横坐标，(y
’1，y
’2，
…
，y’n2
)为目标光场的纵坐标；
[0072]
s4：将输入光场u和s3中得到的矩阵h代入式(4)，计算得到目标光场u如下矩阵乘积的形式，如公式(6)所示：
[0073][0074]
其中，输入光场u为n1×
n1的矩阵，目标光场为u为n2×
n2的矩阵，h为n1×
n2的矩阵。
[0075]
在步骤s2中，数值积分方法选取梯度算法，也可以选择辛普森算法或其它合适的数值积分方法，对于辛普森算法：
[0076]
有
[0077]
实施例1：矩阵光束真空传输
[0078]
利用一种基于矩阵运算的光场传输数值仿真方法进行矩形光束真空传输仿真，包括如下步骤：
[0079]
s1：确定传输系统的abcd矩阵。
[0080]
真空传输菲涅尔衍射的abcd矩阵为
[0081][0082]
s2：将广义惠更斯-菲涅尔积分通过梯度算法转换为矩阵乘法
[0083][0084]
s3：将s1中abcd矩阵参数(a、b、d)和s2中与数值积分方法有关的w代入式2，求解h
n1
×
n2
[0085][0086]
其中为输入光场网格点的横坐标，为输入光场网格点的横坐标，为目标光场网格点的横坐标。
[0087]
在本实例中，n1，n2皆选取为512，δ选取为0.1mm，传输距离z选取为300m。
[0088]
s4：将输入光场u(矩形光束)和步骤三中得到的矩阵h代入式3，计算得到目标光场。本实例矩形光束如图2所示，其中光波长为1μm，光功率为1w，光斑尺寸为8mm。
[0089]
图3为本发明矩阵方法和傅里叶变换方法在矩形光束真空传输场景中的数值仿真结果对比。从中可以看到，矩阵方法相比于傅里叶变换方法，与解析解更加吻合。本实例中，矩阵运算通过blas数值计算工具包实现。对于300m的传输距离，傅里叶变换方法耗时2.239s；而本发明中的矩阵算法，耗时0.343s，快了大约6倍。该实例证明了本发明一种基于矩阵运算的光场传输数值仿真方法具有快速、准确的优势。
[0090]
实施例2：基模高斯光束真空传输
[0091]
利用一种基于矩阵运算的光场传输数值仿真方法进行基模高斯光束真空传输仿真，包括如下步骤：
[0092]
s1：确定传输系统的abcd矩阵。
[0093]
真空传输菲涅尔衍射的abcd矩阵为
[0094][0095]
s2：将广义惠更斯-菲涅尔积分通过梯度算法转换为矩阵乘法
[0096][0097][0098]
w＝(0.5,1,1,
…
0.5)
×
δ
[0099]
s3：将s1中的abcd矩阵参数(a、b、d)和s2中与数值积分方法有关的w代入式2，求解h
n1
×
n2
[0100][0101]
其中为输入光场网格点的横坐标，为输入光场网格点的横坐标，为目标光场网格点的横坐标。
[0102]
在本实例中，n1，n2皆选取为512，δ选取为0.1mm，传输距离z选取为360m。
[0103]
s4：将输入光场u(基模高斯光束)和步骤3中的矩阵h代入式3，计算得到目标光场。本实例基模高斯光束如图3所示，其中光波长为1μm，光功率为1w，光斑尺寸为8mm。
[0104]
图5为本发明矩阵方法和傅里叶变换方法在基模高斯光束真空传输场景中的数值仿真结果对比。从中可以看到，矩阵方法与解析解高度一致，而傅里叶变换方法则出现了吉布斯现象，信号发生了振荡。本实例中，矩阵运算通过blas数值计算工具包实现。对于360m的传输距离，傅里叶变换方法耗时2.199s；而本发明中的矩阵方法耗时0.632s，快了大约3倍。该实例同样证明了本发明一种基于矩阵运算的光场传输数值仿真方法具有快速、准确的优势。