前言
在波前整形(WFS)领域,虽然LC-SLM能实现聚焦PBR系数 π 4 \frac{\pi}{4} 4π 的相位调制,但其典型值仅为 100 100 100 Hz的低刷新率,导致无论是基于迭代的波前优化还是测TM实现透过复杂介质聚焦均较慢,从而限制了其应用场景。 相比之下,含大量数字微镜阵列的DMD具备高达23KHz的刷新率,其每个微镜都能在“ON”和“OFF”两种状态间切换,能提供聚焦PBR系数 1 2 π \frac{1}{2\pi} 2π1 的二元幅度调制。基于DMD的波前优化算法[1,2] 和传输矩阵™测量用于透过光学散射实现聚焦[3,4,5] 已经有不少研究。这篇博客总结一下基于DMD实现透过多模光纤(MMF)的聚焦。
基于DMD表征MMF的背景研究
同散射介质类似,MMF也是一种复杂介质。光透过MMF传播会发生模式色散因而产生波前畸变,导致其无法直接用于聚焦,或形成object图案用于成像。通常基于全息法让参考光与MMF输出散斑干涉来还原光场相位信息,进而测量MMF的复值TM,但全息法装置复杂,且对系统稳定性依赖高。对于基于DMD表征MMF的非全息法,有:
- 复值TM:基于变分法的相位复原(phase retrieval),仅根据MMF输出光场的强度测量和DMD二值化输入图案估计复值TM[3]。
- 部分TM:包括①对基于Hadamard基的DMD pattern巧妙构造,仅由输出光场强度直接测出复值TM的实部,得到实值TM[6,10];②基于一系列输入掩膜下输出光场强度的变化测量二值TM,用于快速的DMD二元幅度调制[7]。
- DMD微镜开闭的条件概率算法(CPA)[8],用于决定最优DMD pattern实现聚焦。
- 基于DMD的迭代波前优化算法,典型的有遗传算法 (GA)[2]、粒子群优化(PSO)算法、蚁群(AC)算法以及针对波前介质扰动改进后的动态优化算法(DMA)[9]等
针对上述四类测量TM或者迭代优化波前的方法,文献[11]比较了各自的理论聚焦PBR,以及在相同实验装置和条件下产生最优DMD pattern进行聚焦的速度(仿真和实验),同时本文也主要是基于翻译该文献加上个人总结而得来。
测复值TM的理论聚焦PBR
一旦测出复值TM,要在MMF输出光场第 m m m 个mode处聚焦,我们就取TM第 m m m 行元素 t m n ( n = 1 , 2 , . . . N ) t_{mn} (n=1,2,...N) tmn(n=1,2,...N) 二值化处理构成DMD pattern,有两种方式:
- 相位在 ( − π 2 , π 2 ) (-\frac{\pi}{2}, \frac{\pi}{2}) (−2π,2π) 之间(实部>0)的对应1,否则对应0;
- 相位在 ( 0 , π ) (0, \pi) (0,π) 之间(虚部>0)的对应1,否则对应0;
这样DMD输入光场能透过MMF在目标输出位置发生干涉相长形成聚焦。事实上,打开某些微镜对背景光强的比对目标光强的贡献还大,因此,上述两个决定最优DMD二值图案的标准(criteria) 需要修正,当关闭这些不该打开的微镜能进一步提升聚焦PBR,图解如下:
令
N
N
N 为DMD微镜 (input modes) 数,假设
t
m
n
t_{mn}
tmn 服从Rayleigh分布且
∣
t
m
n
∣
2
\left | t_{mn} \right | ^2
∣tmn∣2 遵循指数衰减
∣
t
m
n
∣
2
∼
e
−
1
/
2
σ
2
\left | t_{mn} \right | ^2 \sim e^{-1/2\sigma^2}
∣tmn∣2∼e−1/2σ2 ,其中
2
σ
2
2\sigma^2
2σ2 是
t
m
n
t_{mn}
tmn每个元素的总体平均强度,
φ
\varphi
φ 为绝对相位差异的上确界:
0
⩽
∣
Δ
θ
∣
<
φ
0\leqslant\left|\Delta\theta\right|<\varphi
0⩽∣Δθ∣<φ,则基于DMD二元幅度调制的理论聚焦PBR[12]为:
P
B
R
≈
(
N
−
1
)
s
i
n
2
φ
4
φ
PBR \approx \frac{(N-1)sin^2\varphi}{4\varphi}
PBR≈4φ(N−1)sin2φ
其中绝对相位差异的上确界
φ
=
π
/
2
\varphi=\pi/2
φ=π/2。注意到该理论PBR公式只考虑
t
m
n
t_{mn}
tmn 的相位信息用于决定最优DMD pattern做聚焦,适用于测复值TM;其它测实值TM、条件概率算法CPA以及GA等迭代法能实现的理论聚焦PBR将会与之做对比。
测实值TM的改进算法
测实值TM是由Yongkeun Park在2017年最先提出的[6],这里先简要归纳如下:
利用Hadamard基 h n \mathbf{h_n } hn 作为DMD输入,对应的输出光场为 σ n \mathbf{\sigma_n} σn,输入与输出的线性关系可由复值TM T \mathbf{T} T表示为:
[ σ 1 , σ 2 ⋯ σ N ] = T [ h 1 , h 2 ⋯ h N ] \left [ \mathbf{\sigma_1}, \mathbf{\sigma_2} \cdots \mathbf{\sigma_N} \right ] = \mathbf{T}\left [ \mathbf{h_1}, \mathbf{h_2} \cdots \mathbf{h_N} \right ] [σ1,σ2⋯σN]=T[h1,h2⋯hN]
复值TM实部 R e ( T ) = S H T = S H Re(\mathbf{T}) =\mathbf{S}\mathbf{H}^{T}=\mathbf{S}\mathbf{H} Re(T)=SHT=SH,其中 S = 1 N R e [ σ 1 , σ 2 ⋯ σ n ] \mathbf{S} = \frac{1}{N}Re\left [ \mathbf{\sigma_1}, \mathbf{\sigma_2} \cdots \mathbf{\sigma_n} \right ] S=N1Re[σ1,σ2⋯σn],那么怎么获得实值TM R e ( T ) Re(\mathbf{T}) Re(T)呢?
对于Hadamard矩阵 H = [ h 1 , h 2 ⋯ h n ] \mathbf{H} = \left[ \mathbf{h_1}, \mathbf{h_2} \cdots \mathbf{h_n}\right] H=[h1,h2⋯hn],我们构造如下两个二值化Hadamard基作为DMD pattern:
v n ± = 1 2 ( h 1 ± h n ) , n = 1 , 2 ⋯ N \mathbf{v_{n}^{\pm}} = \frac{1}{2}(\mathbf{h_1} \pm \mathbf{h_n}), n=1,2 \cdots N vn±=21(h1±hn),n=1,2⋯N
对应 v n ± \mathbf{v_{n}^{\pm}} vn± 的输出光强 I n ± \mathbf{I_{n}^{\pm}} In± 为:
I n ± = ∣ T v n ± ∣ 2 = ∣ 1 2 T ( h 1 ± h n ) ∣ 2 = 1 4 ∣ σ 1 ∣ 2 + 1 4 ∣ σ n ∣ 2 ± 1 2 R e ( σ 1 ∗ σ n ) ( 复 数 相 加 的 模 的 平 方 运 算 ) \mathbf{I_{n}^{\pm}} = \left| \mathbf{Tv_{n}^{\pm}}\right|^2 = \left| \frac{1}{2}\mathbf{T(\mathbf{h_1} \pm \mathbf{h_n})}\right|^2 = \frac{1}{4}\left|\mathbf{\sigma_1}\right|^2 + \frac{1}{4}\left|\mathbf{\sigma_n}\right|^2 \pm \frac{1}{2}Re(\mathbf{\sigma^*_1\sigma_n}) (复数相加的模的平方运算) In±=∣∣Tvn±∣∣2=∣∣∣∣21T(h1±hn)∣∣∣∣2=41∣σ1∣2+41∣σn∣2±21Re(σ1∗σn)(复数相加的模的平方运算)
令 σ 1 \mathbf{\sigma_1} σ1 相位为 0 不失一般性,则我们得到输出光场 σ n \mathbf{\sigma_n} σn 的实部为 R e [ σ n ] = I n + − I n − I 1 Re\left[\mathbf{\sigma_n}\right] = \frac{\mathbf{I_{n}^{+} - I_{n}^{-}}}{\sqrt{\mathbf{I_1}}} Re[σn]=I1In+−In− ,其中 I 1 = ∣ T v 1 ∣ 2 = σ 1 2 \mathbf{I_1} = \left| \mathbf{Tv_{1}}\right|^2 = \mathbf{\sigma^2_1} I1=∣Tv1∣2=σ12
这样我们就能通过组合得到二值化的Hadamard基,由输出光场强度直接测出复杂介质的实值TM,要实现透过复杂介质在第 m m m 个mode处聚焦,取实值TM第 m m m 行,按照元素>0的取1否则取0的方法得到最优DMD输入图案。
Tianrui ZHAO等人则进一步提出了实值强度TM ( R V I T M RVITM RVITM),连接二值DMD输入与输出光强[10]。
同Park测实数值TM构造二值化Hadamard基一样,构造 两个大小为 N × N N \times N N×N二值化Hadamard矩阵 H 1 = ( H + 1 ) / 2 , H 2 = ( − H + 1 ) / 2 \mathbf{H_1 = (H + 1) / 2, H_2 =(-H + 1) / 2} H1=(H+1)/2,H2=(−H+1)/2,其中每个列向量作为Hadamard基构成DMD pattern,则二值输入与输出光强的线性关系为:
[ I 1 1 ⋯ I 1 2 N ⋮ ⋱ ⋮ I M 1 ⋯ I M 2 N ] = R V I T M ⋅ [ H 1 , H 2 ] \begin{bmatrix} &I_{1}^{1} &\cdots &I_{1}^{2N} \\ &\vdots &\ddots &\vdots \\ &I_{M}^{1} &\cdots &I_{M}^{2N} \end{bmatrix} = RVITM \cdot \left[\mathbf{H_1, H_2} \right ] ⎣⎢⎡I11⋮IM1⋯⋱⋯I12N⋮IM2N⎦⎥⎤=RVITM⋅[H1,H2]
先对上式都右乘主对角元为2的对角矩阵,再考虑到 I 1 = R V I T M ⋅ h 1 \mathbf{I_1} = RVITM \cdot \mathbf{h_1} I1=RVITM⋅h1, 其中 h 1 \mathbf{h_1} h1 为全1向量,若等式右边矩阵 [ 2 H 1 , 2 H 2 ] \left[\mathbf{2H_1, 2H_2} \right ] [2H1,2H2] 减去1,则有:
[ 2 I 1 1 − I 1 1 ⋯ 2 I 1 2 N − I 1 1 ⋮ ⋱ ⋮ 2 I M 1 − I M 1 ⋯ 2 I M 2 N − I M 1 ] = R V I T M ⋅ [ 2 H 1 − 1 , 2 H 2 − 1 ] = R V I T M ⋅ [ H , − H ] \begin{bmatrix} &2I_{1}^{1}-I_{1}^{1} &\cdots &2I_{1}^{2N}-I_{1}^{1} \\ &\vdots &\ddots &\vdots \\ &2I_{M}^{1}-I_{M}^{1} &\cdots &2I_{M}^{2N}-I_{M}^{1} \end{bmatrix} = RVITM \cdot \left[\mathbf{2H_1-1, 2H_2-1} \right ] = RVITM \cdot \left[\mathbf{H, -H} \right ] ⎣⎢⎡2I11−I11⋮2IM1−IM1⋯⋱⋯2I12N−I11⋮2IM2N−IM1⎦⎥⎤=RVITM⋅[2H1−1,2H2−1]=RVITM⋅[H,−H]
可见等式右侧变成 R V I T M RVITM RVITM 直接与大小为 2 N × N 2N \times N 2N×N 的Hadamard矩阵 [ H , − H ] \left[\mathbf{H, -H} \right ] [H,−H]相乘!因此我们能对二值化Hadamard矩阵 H 1 \mathbf{H_1} H1 和 H 2 \mathbf{H_2} H2 作为输入的MMF输出光强做一些处理,由Hadamard矩阵性质: [ H , − H ] − 1 = 1 2 N ⋅ [ H , − H ] T \left[\mathbf{H, -H} \right ]^{-1} = \frac{1}{2N} \cdot \left[\mathbf{H, -H} \right ]^T [H,−H]−1=2N1⋅[H,−H]T,我们可以直接算出实值强度TM为:
R V I T M = 1 2 N ⋅ [ 2 I 1 1 − I 1 1 ⋯ 2 I 1 2 N − I 1 1 ⋮ ⋱ ⋮ 2 I M 1 − I M 1 ⋯ 2 I M 2 N − I M 1 ] ⋅ [ H , − H ] T RVITM = \frac{1}{2N} \cdot \begin{bmatrix} &2I_{1}^{1}-I_{1}^{1} &\cdots &2I_{1}^{2N}-I_{1}^{1} \\ &\vdots &\ddots &\vdots \\ &2I_{M}^{1}-I_{M}^{1} &\cdots &2I_{M}^{2N}-I_{M}^{1} \end{bmatrix} \cdot \left[\mathbf{H, -H} \right ]^T RVITM=2N1⋅⎣⎢⎡2I11−I11⋮2IM1−IM1⋯⋱⋯2I12N−I11⋮2IM2N−IM1⎦⎥⎤⋅[H,−H]T
进一步,通过对测出的 R V I T M RVITM RVITM求伪逆,能用于从MMF输出散斑中重建输入图像:
I i m a g e = R V I T M + ⋅ I o u t ≈ I i n I_{image} = RVITM^{+} \cdot I_{out} \approx I_{in} Iimage=RVITM+⋅Iout≈Iin
可以发现:实值强度TM ( R V I T M RVITM RVITM)较实值TM的优势是直接连接二值DMD输入与输出光场的强度测量值,而非输出光场的实部,从而能直接从MMF输出散斑测量结果中重建输入图像。
Fig 2. (a) 用RVITM表征MMF图解, (b)基于RVITM求逆从MMF输出散斑图案重建输入图像(KCL), (c) 测RVITM的实验设置
进一步, 实值强度TM ( R V I T M RVITM RVITM) 还能表达成编码相应复值TM的幅度与相位信息的形式[11].
由于输出光场强度与其激发的PA信号信号幅度有 系数为
α
\alpha
α 的正比例关系,因此考虑测量结果为目标处输出光强时,连接第
n
n
n 个DMD微镜处输入光场与
m
m
m 处目标输出光强的
R
V
I
T
M
RVITM
RVITM元素为:
r
v
i
t
m
m
n
=
A
m
n
A
R
c
o
s
(
θ
m
n
−
ϕ
R
)
rvitm_{mn}= A_{mn}A_Rcos(\theta_{mn} - \phi_R)
rvitmmn=AmnARcos(θmn−ϕR)
根据上述推导,当 θ m n \theta_{mn} θmn 在 [ ϕ R − π / 2 , ϕ R + π / 2 ] \left[ \phi_R-\pi/2, \phi_R+\pi/2\right] [ϕR−π/2,ϕR+π/2] 之间,则 r v i t m m n ⩾ 0 rvitm_{mn}\geqslant 0 rvitmmn⩾0,打开相对应的DMD微镜则在目标输出处发生干涉叠加用于聚焦。作者假设对应更大 r v i t m m n rvitm_{mn} rvitmmn值的DMD微镜对焦点处光强比背景光强的贡献更大,因此要优先打开。作者因而按照对应 r v i t m m n rvitm_{mn} rvitmmn 大小对DMD微镜做了降序排列,在仿真和实验上探究只打开前 P % P \% P% 的微镜对聚焦PBR的影响,结果如下:
基于变分法的相位复原估计复值TM
条件概率算法(CPA)
迭代波前优化 (以GA为例)
参考文献
[1] D. Akbulut, T. J. Huisman, E. G. van Putten, W. L. Vos, and A. P. Mosk, “Focusing light through random photonic media by binary amplitude modulation,” Opt. express 19, 4017–4029 (2011).
[2] X. Zhang and P. Kner, “Binary wavefront optimization using a genetic algorithm,” J. Opt. 16, 125704 (2014).
[3] A. Drémeau, A. Liutkus, D. Martina, O. Katz, C. Schülke, F. Krzakala, S. Gigan, and L. Daudet, “Reference-less measurement of the transmission matrix of a highly scattering material using a dmd and phase retrieval techniques,” Opt. express 23, 11898–11911 (2015).
[4] G. Huang, D. Wu, J. Luo, Y. Huang, and Y. Shen, “Retrieving the optical transmission matrix of a multimode fiber using the extended kalman filter,” Opt. Express 28, 9487–9500 (2020).
[5] M. N’Gom, T. B. Norris, E. Michielssen, and R. R. Nadakuditi, “Mode control in a multimode fiber through acquiring its transmission matrix from a reference-less optical system,” Opt. letters 43, 419–422 (2018).
[6] H. Yu, K. Lee, and Y. Park, “Ultrahigh enhancement of light focusing through disordered media controlled by mega-pixel modes,” Opt. express 25, 8036–8047 (2017).
[7] X. Tao, D. Bodington, M. Reinig, and J. Kubby, “High-speed scanning interferometric focusing by fast measurement of binary transmission matrix for channel demixing,” Opt. express 23, 14168–14187 (2015).
[8] T. Zhao, L. Deng, W. Wang, D. S. Elson, and L. Su, “Bayes’ theorem-based binary algorithm for fast reference-less calibration of a multimode fiber,” Opt. express 26, 20368–20378 (2018).
[9] Huanhao Li, etc, “Adaptive optical focusing through perturbed scattering media with a dynamic mutation algorithm,” Photon. Res. 9, 202-212 (2021)
[10] T. Zhao, S. Ourselin, T. Vercauteren, and W. Xia, “Seeing through multimode fibers with real-valued intensity transmission matrices,” Opt. Express 28, 20978–20991 (2020).
[11]. TIANRUI ZHAO, etc. Focusing light through multimode fibres using a digital micromirror device: a comparison study of non-holographic approaches(in press)
[12] D. Wang, E. H. Zhou, J. Brake, H. Ruan, M. Jang, and C. Yang, “Focusing through dynamic tissue with millisecond digital optical phase conjugation,” Optica 2, 728–735 (2015).