留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

集成光子储备池的可调谐光滤波器

李文璐 裴丽 白冰 左晓燕 王建帅 郑晶晶 李晶 宁提纲

李文璐, 裴丽, 白冰, 左晓燕, 王建帅, 郑晶晶, 李晶, 宁提纲. 集成光子储备池的可调谐光滤波器[J]. 红外与激光工程, 2023, 52(9): 20220915. doi: 10.3788/IRLA20220915
引用本文: 李文璐, 裴丽, 白冰, 左晓燕, 王建帅, 郑晶晶, 李晶, 宁提纲. 集成光子储备池的可调谐光滤波器[J]. 红外与激光工程, 2023, 52(9): 20220915. doi: 10.3788/IRLA20220915
Li Wenlu, Pei Li, Bai Bing, Zuo Xiaoyan, Wang Jianshuai, Zheng Jingjing, Li Jing, Ning Tigang. Tunable optical filter with integrated photonic reservoir computing[J]. Infrared and Laser Engineering, 2023, 52(9): 20220915. doi: 10.3788/IRLA20220915
Citation: Li Wenlu, Pei Li, Bai Bing, Zuo Xiaoyan, Wang Jianshuai, Zheng Jingjing, Li Jing, Ning Tigang. Tunable optical filter with integrated photonic reservoir computing[J]. Infrared and Laser Engineering, 2023, 52(9): 20220915. doi: 10.3788/IRLA20220915

集成光子储备池的可调谐光滤波器

doi: 10.3788/IRLA20220915
基金项目: 国家自然科学基金项目(62235003)
详细信息
    作者简介:

    李文璐,女,硕士生,主要从事光芯片在通信中应用方面的研究

  • 中图分类号: TN256

Tunable optical filter with integrated photonic reservoir computing

Funds: National Natural Science Foundation of China (62235003)
  • 摘要: 为了适应滤波、波分复用等不同的应用场景,光子滤波器需要具备可调谐以及滤波形状可变的功能。提出一种集成光子储备池结构的新型可调谐光滤波器。该结构由输入层、储层、读出层三部分组成,输入层由2×1多模干涉仪组成,储层由定向耦合器以及波导组成,按螺旋拓扑规律相互连接组成梅花型网格,匹配粒子群寻优算法进行全光域训练,利用读出层的热光调制器和相移器对光信号的幅度和相位进行调整,实现无限冲激响应型与有限冲激响应型可调谐光滤波器。以无限冲激响应型光滤波器仿真结果为例,分析其自由光谱范围与储层中波导长度的关系,以及定向耦合器分光比对透射端透过率的影响。通过设置储层中波导的相位在0~3/2π范围内变化,实现了光滤波器中心波长在一个自由光谱范围(1.18 nm)内的连续可调。该滤波器采用集成光子储备池与粒子群算法结合方案设计,无需考虑储层光网络传输路径的权重配置,仅对读出层权重进行训练,大大简化了网络训练过程。同时,该滤波器具有尺寸小、功耗低、灵活性高、算法可控滤波参数的优势,广泛应用于集成微波光子学和光通信领域中。
  • 图  1  光子储备池结构。(a)整体结构;(b)输入层结构

    Figure  1.  Structure of photonic reservoir computing. (a) Overall structure; (b) Input layer

    图  2  光子滤波器系统实验装置图

    Figure  2.  Diagram of photonic filter system experimental setup

    图  3  光子储备池中器件模型示意图。(a)波导; (b)定向耦合器;(c)光调制器

    Figure  3.  Schematic diagram of device model in photonic reservoir computing. (a) Waveguide; (b) Directional coupler; (c) Optical modulation

    图  4  储层输出的8路信号光谱图。(a)第1路;(b)第2路;(c)第3路;(d)第4路;(e)第5路;(f)第6路;(g)第7路;(h)第8路

    Figure  4.  8 channel signal spectra of reservoir output. (a) Rout 1; (b) Rout 2; (c) Rout 3; (d) Rout 4; (e) Rout 5; (f) Rout 6; (g) Rout 7; (h) Rout 8

    图  5  IIR型光滤波器仿真训练效果。(a) IIR滤波谱拟合效果图;(b)Fitness曲线

    Figure  5.  Simulation training effect of IIR optical filter. (a) IIR filter spectrum fitting rendering; (b) Fitness curve

    图  6  FIR型光滤波器仿真训练效果。(a) FIR滤波谱拟合效果图;(b) Fitness曲线

    Figure  6.  Simulation training effect of FIR optical filter. (a) FIR filter spectrum fitting rendering; (b) Fitness curve

    图  7  储层中波导在不同长度下的IIR滤波谱。(a)波导长度为400 μm ;(b)波导长度为500 μm ;(c)波导长度为600 μm

    Figure  7.  IIR filter spectrum of waveguide at different lengths in reservoir. (a) Waveguide length of 400 μm; (b) The waveguide length of 500 μm; (c) The waveguide length of 600 μm

    图  8  FSR随波导长度变化

    Figure  8.  FSR changing with waveguide length

    图  9  不同定向耦合器分光比滤波器输出功率随波长变化

    Figure  9.  Output power of filter changing with the wavelength in different directional coupler radio

    图  10  不同相移量滤波器输出功率随波长变化

    Figure  10.  Output power of filter changing with the wavelength in different phase shifts

    表  1  IIR型光滤波器对应α值的设定参数

    Table  1.   Setting parameters of the corresponding α of IIR optical filter

    |ylabel|max|ylabel|minabx
    10.054.510.95
    下载: 导出CSV

    表  2  IIR型光滤波器8路信号权重训练值

    Table  2.   IIR optical filter 8 channel signal weight training value

    IIR1st2nd3rd4th5th6th7th8th
    WAmp0.4960.0150.8470.3580.2400.6490.8650.641
    φw4.7170.5935.011.1764.0720.1732.1801.935
    下载: 导出CSV

    表  3  FIR型光滤波器对应α值的设定参数

    Table  3.   Setting parameters of the corresponding α of FIR optical filter

    |y’label|max|y’label|mina’b’x’
    0.810240.4
    下载: 导出CSV

    表  4  FIR型光滤波器8路信号权重训练值

    Table  4.   FIR optical filter 8 channel signal weight training value

    FIR1st2nd3rd4th5th6th7th8th
    WAmp0.9660.4520.9980.9280.6850.9910.9970.766
    φw3.2582.8613.0975.1893.0800.2596.2825.246
    下载: 导出CSV
  • [1] Madsen C K, Zhao J H. Optical Filter Design And Analysis: a Signal Processing Approach[M]. New York: Wily, 1999.
    [2] 钟昌锦, 刘斌, 付益, 等. 密集波分复用(DWDM)技术实现原理分析[J]. 广东通信技术, 2021, 41(4): 66-68. doi:  10.3969/j.issn.1006-6403.2021.04.015
    [3] 顾佩芸, 张娟, 周益. 基于各向异性介质的多通道平顶偏振滤波器的设计[J]. 光电子激光, 2012, 23(6): 1043-1050.
    [4] Butt M A, Khonina S N, Kazansiy N L. A compact design of a modified Bragg grating filter based on a metal-insulator-metal waveguide for filtering and temperature sensing applications [J]. Optik, 2022, 251: 168466. doi:  10.1016/j.ijleo.2021.168466
    [5] 廖莎莎, 包航, 冯玉婷, 等. 基于级联啁啾亚波长光栅辅助反向耦合器的超宽带可调滤波器[J]. 光学学报, 2022, 42(14): 1405003. doi:  10.3788/AOS202242.1405003

    Liao Shasha, Bao Hang, Feng Yuting, et al. Ultra-Broadband tunable filter based on cascaded chirped subwavelength grating assisted contra-directional coupler [J]. Acta Optica Sinica, 2022, 42(14): 1405003. (in Chinese) doi:  10.3788/AOS202242.1405003
    [6] 杨双月, 谷一英, 胡晶晶, 等. 基于亚波长超窄带滤波器设计[J]. 红外与激光工程, 2020, 49(S1): 20200134.

    Yang Shuangyue, Gu Yiying, Hu Jingjing, et al. Design of sub-wavelength ultra-narrow band filter [J]. Infrared and Laser Engineering, 2020, 49(S1): 20200134. (in Chinese)
    [7] 陈彧芳, 沈骁, 周权, 等. 基于磷酸盐玻璃微球腔的全光调谐光纤滤波器[J]. 中国激光, 2021, 48(1): 0106003. doi:  10.3788/CJL202148.0106003

    Chen Yufang, Shen Xiao, Zhou Quan, et al. All-optical tunable fiber filter based on phosphate glass microspheres [J]. Chinese Journal of Lasers, 2021, 48(1): 0106003. (in Chinese) doi:  10.3788/CJL202148.0106003
    [8] 唐钊, 张钧翔, 付士杰, 等. 基于MMI滤波器的可调谐连续光全光纤OPO[J]. 红外与激光工程, 2019, 48(5): 520002-0520002(6). doi:  10.3788/IRLA201948.0520002

    Tang Zhao, Zhang Junxiang, Fu Shijie, et al. Tunable CW all-fiber optical parametric oscillator based on the multimode interference filter [J]. Infrared and Laser Engineering, 2019, 48(5): 0520002. (in Chinese) doi:  10.3788/IRLA201948.0520002
    [9] Mokhtari M. Tunable optical filter design with ring resonator based Sagnac loop [J]. Optik, 2021, 242: 167068. doi:  10.1016/j.ijleo.2021.167068
    [10] 崔文翔, 周雪芳, 胡淼, 等. 双Sagnac环滤波器的传输特性分析和实验研究[J]. 中国激光, 2022, 49(4): 0406006. doi:  10.3788/CJL202249.0406006

    Cui Wenxiang, Zhou Xuefang, Hu Miao, et al. Analysis and experimental study on transmission characteristics of double Sagnac loop filter [J]. Chinese Journal of Lasers, 2022, 49(4): 0406006. (in Chinese) doi:  10.3788/CJL202249.0406006
    [11] Deng Lin, Li Dezhao, Liu Zilong, et al. Tunable optical filter using second-order micro-ring resonator [J]. Chinese Physics B, 2017, 26(2): 024209. doi:  10.1088/1674-1056/26/2/024209
    [12] Wang Haoyan, Dai Jincheng, Jia Hao, et al. Polarization-independent tunable optical filter with variable bandwidth based on silicon-on-insulator waveguides [J]. Nanophotonics, 2018, 7(8): 1469-1477. doi:  10.1515/nanoph-2018-0058
    [13] 刘大建, 赵伟科, 张龙, 等. 高性能无源硅光波导器件: 发展与挑战[J]. 光学学报, 2022, 42(17): 1713001. doi:  10.3788/AOS202242.1713001

    Liu Dajian, Zhao Weike, Zhang Long, et al. High-performance passive silicon photonic waveguide devices: progress and challenges [J]. Acta Optica Sinica, 2022, 42(17): 1713001. (in Chinese) doi:  10.3788/AOS202242.1713001
    [14] Zhou Hailong, Zhao Yuhe, Wang Xu, et al. Self-configuring and reconfigurable silicon photonic signal processor [J]. ACS Photonics, 2020, 7(3): 792-799. doi:  10.1021/acsphotonics.9b01673
    [15] Fandiño J S, Muñoz P, Doménech D, et al. A monolithic integrated photonic microwave filter [J]. Nature Photonics, 2017, 11(2): 124-129.
    [16] Li Jiachen, Yang Sigang, Chen Hongwei, et al. Reconfigurable rectangular filter with continuously tunable bandwidth and wavelength [J]. IEEE Photonics Journal, 2020, 12(4): 1-9.
    [17] Katumba A. Energy-efficient photonic neuromorphic computing for telecommunication applications[D]. Chicago: Ghent University, 2019.
  • [1] 刘云哲, 董岩, 王伟, 宋建林.  光电跟踪系统的摩擦模型辨识与补偿策略研究 . 红外与激光工程, 2023, 52(11): 20230151-1-20230151-9. doi: 10.3788/IRLA20230151
    [2] 丁国建, 王晓晖, 冯琦, 于萍, 贾海强, 陈弘, 汪洋.  高效铌酸锂薄膜波导模斑转换器设计 . 红外与激光工程, 2023, 52(9): 20220897-1-20220897-9. doi: 10.3788/IRLA20220897
    [3] 刘鹏飞, 任麟昊, 闻浩, 施雷, 张新亮.  集成电光频率梳研究进展(特邀) . 红外与激光工程, 2022, 51(5): 20220381-1-20220381-18. doi: 10.3788/IRLA20220381
    [4] 陈沁, 南向红, 梁文跃, 郑麒麟, 孙志伟, 文龙.  片上集成光学传感检测技术的研究进展(特邀) . 红外与激光工程, 2022, 51(1): 20210671-1-20210671-18. doi: 10.3788/IRLA20210671
    [5] 王希, 刘英杰, 张子萌, 王嘉宁, 姚勇, 宋清海, 徐科.  2 μm波段片上光子集成器件的研究进展(特邀) . 红外与激光工程, 2022, 51(3): 20220087-1-20220087-12. doi: 10.3788/IRLA20220087
    [6] 罗强, 薄方, 孔勇发, 张国权, 许京军.  铌酸锂薄膜微腔激光器研究进展(特邀) . 红外与激光工程, 2021, 50(11): 20210546-1-20210546-13. doi: 10.3788/IRLA20210546
    [7] 吕桓林, 梁宇鑫, 韩秀友, 谷一英, 武震林, 赵明山.  基于狭缝波导的聚合物基微环折射率传感器研究 . 红外与激光工程, 2020, 49(1): 0118001-0118001(6). doi: 10.3788/IRLA202049.0118001
    [8] 丁君珂, 陈浩, 蒋建光, 孟浩然, 刘欣悦, 郝寅雷.  集成光学移相器波长相关性的比较研究 . 红外与激光工程, 2019, 48(5): 520001-0520001(5). doi: 10.3788/IRLA201948.0520001
    [9] 郝寅雷, 丁君珂, 陈浩, 蒋建光, 孟浩然, 刘欣悦.  集成光学移相干涉仪的研制与性能表征 . 红外与激光工程, 2019, 48(4): 420001-0420001(5). doi: 10.3788/IRLA201948.0420001
    [10] 刘鑫, 孔梅, 徐亚萌, 王雪萍.  微环谐振器中各参数对光速控制输出脉冲畸变的影响仿真分析 . 红外与激光工程, 2019, 48(9): 918002-0918002(6). doi: 10.3788/IRLA201948.0918002
    [11] 管磊, 王卓然, 袁国慧, 陈昱任, 董礼, 彭真明.  微环差分光子生物传感器的传感性能 . 红外与激光工程, 2018, 47(2): 222002-0222002(6). doi: 10.3788/IRLA201847.0222002
    [12] 谭巧, 徐启峰, 黄奕钒, 项宇锴.  一种基于径向偏振解调的线性光学电流传感器 . 红外与激光工程, 2018, 47(2): 222003-0222003(6). doi: 10.3788/IRLA201847.0222003
    [13] 杨旭, 李亚明, 郭肃丽, 李晶, 刘旭东.  拉曼增益对回音壁模式光学微腔的全光调制 . 红外与激光工程, 2017, 46(11): 1122003-1122003(5). doi: 10.3788/IRLA201746.1122003
    [14] 张尧, 王宏力, 陆敬辉, 何贻洋, 姜伟.  基于粒子群算法的星敏感器光学误差标定方法 . 红外与激光工程, 2017, 46(10): 1017002-1017002(8). doi: 10.3788/IRLA201770.1017002
    [15] 陈明, 赵永乐, 牛奔, 宋华.  基于铌酸锂光子线波长分裂器的研究 . 红外与激光工程, 2016, 45(6): 620003-0620003(4). doi: 10.3788/IRLA201645.0620003
    [16] 薛彬, 邾继贵, 郑迎亚.  工作空间测量定位系统最佳测量点的确定方法 . 红外与激光工程, 2015, 44(4): 1218-1222.
    [17] 秦华, 冯东太, 刘波, 吴国栋.  用粒子群算法设计非球面准直透镜的方法 . 红外与激光工程, 2015, 44(6): 1811-1817.
    [18] 许兆美, 刘永志, 杨刚, 王庆安.  粒子群优化BP神经网络的激光铣削质量预测模型 . 红外与激光工程, 2013, 42(9): 2370-2374.
    [19] 秦华, 王立刚, 张静华, 类成新, 韩克祯.  用粒子群算法设计的5镜系统与ZEBASE中5镜系统比较 . 红外与激光工程, 2013, 42(10): 2724-2731.
    [20] 高旭, 万秋华, 杨守旺, 陈伟, 赵长海.  提高光电轴角编码器细分精度的改进粒子群算法 . 红外与激光工程, 2013, 42(6): 1508-1513.
  • 加载中
图(10) / 表(4)
计量
  • 文章访问数:  79
  • HTML全文浏览量:  28
  • PDF下载量:  28
  • 被引次数: 0
出版历程
  • 收稿日期:  2023-01-03
  • 修回日期:  2023-04-06
  • 刊出日期:  2023-09-28

集成光子储备池的可调谐光滤波器

doi: 10.3788/IRLA20220915
    作者简介:

    李文璐,女,硕士生,主要从事光芯片在通信中应用方面的研究

基金项目:  国家自然科学基金项目(62235003)
  • 中图分类号: TN256

摘要: 为了适应滤波、波分复用等不同的应用场景,光子滤波器需要具备可调谐以及滤波形状可变的功能。提出一种集成光子储备池结构的新型可调谐光滤波器。该结构由输入层、储层、读出层三部分组成,输入层由2×1多模干涉仪组成,储层由定向耦合器以及波导组成,按螺旋拓扑规律相互连接组成梅花型网格,匹配粒子群寻优算法进行全光域训练,利用读出层的热光调制器和相移器对光信号的幅度和相位进行调整,实现无限冲激响应型与有限冲激响应型可调谐光滤波器。以无限冲激响应型光滤波器仿真结果为例,分析其自由光谱范围与储层中波导长度的关系,以及定向耦合器分光比对透射端透过率的影响。通过设置储层中波导的相位在0~3/2π范围内变化,实现了光滤波器中心波长在一个自由光谱范围(1.18 nm)内的连续可调。该滤波器采用集成光子储备池与粒子群算法结合方案设计,无需考虑储层光网络传输路径的权重配置,仅对读出层权重进行训练,大大简化了网络训练过程。同时,该滤波器具有尺寸小、功耗低、灵活性高、算法可控滤波参数的优势,广泛应用于集成微波光子学和光通信领域中。

English Abstract

    • 可调谐光滤波器,通过对中心滤波波长的实时调谐,实现对所需光波长的选择性通过,逐渐发展成为密集波分复用网络信号控制和处理的关键技术[1]。随着波分复用系统信道数的增加、信道间隔向更窄带宽发展[2],对可调谐光滤波器提出了更大调谐范围、更低损耗且小型化的要求,基于介质膜滤光片[3]、光纤光栅[4-6]结构的传统滤波器已不能满足通信系统高速发展的需求。光纤滤波器[7-8]偏振不敏感、插入损耗小、结构简单,且易于与光纤通信系统结合实现全光纤可调谐波长滤波,被广泛应用于光通信领域。近年来,随着光子集成技术的发展,采用集成方式成为未来实现高速率、大容量以及低损耗光网络的重要手段。光纤滤波器难以满足集成化发展的趋势,而集成化的光滤波器具有损耗低、尺寸小等优势成为研究热点。目前国内外对集成可调谐光滤波器的研究主要基于以下几种类型:Sagnac环干涉结构[9-10]、微环谐振器(Micro-ring Resonator, MRR)[11-13]、马赫-曾德尔干涉仪(Mach-Zehnder Interferometer, MZI)[14]、环辅助马赫-曾德尔干涉仪(Ring-assisted Mach-Zehnder Interferometer, RAMZI)[15-16]等。其中基于Sagnac环具有结构简单、输入偏振无关的优势,但调谐范围有限;基于MRR级联结构的集成光滤波器结构简单紧凑,具有良好的波长选择性、带外信号抑制比高、功耗低等优势,但制造MRR的谐振腔尺寸精确度要求较高,并且滤波参数精确度控制的实现较为困难;华中科技大学研究团队提出一种由MZI级联构成的自配置集成可调谐光滤波器,将人工智能技术和可编程光信号处理芯片相结合,实现具有学习能力的自配置可编程光信号处理器,但采用的MZI器件对波长不敏感,调谐灵敏度比较低,且训练过程较为复杂;提出了基于硅基RAMZI型滤波器结构的光滤波器,通过调节两个微环的谐振波长,可以实现波长和带宽的同时可调,并且具有较高的消光比,但滤波稳定性较差,且仅能实现单一滤波功能。虽然基于上述器件的光滤波器已做验证,但采用新型神经网络,例如光子储备池的全光滤波器还处于空白,仍需探索研究。

      文中为了验证利用光子储备池神经网络与粒子群寻优算法(Particle Swarm Optimization, PSO)结合方案实现可调谐光滤波器的可行性,对光子储备池的全光滤波器的实现进行前沿探索,提出一种基于集成光子储备池的可调谐光滤波器,利用散射矩阵对光信号的处理过程进行分析,得到光子储备池端口间的映射关系; PSO算法对读出层中各路光信号相位和幅度参数进行训练,实现滤波形状可选;改变储层中波导的相位参数,实现中心滤波波长可调。采用光子储备池的优势在于并行处理信号和波长敏感特性,支持光信号的前馈传输和反馈配置,利用粒子群算法进行训练,实现定制化滤波功能,对滤波参数精确控制。且仅对网络部分连接权重训练,减少训练过程计算负担,通过算法训练进行权重寻优,实现给定的多种滤波波形。

    • 图1(a)所示,光子储备池整体结构由输入层、储层和读出层(全光训练读出)三部分组成。输入层结构如图1(b)所示,由耦合比为0.5的3个2×1多模干涉仪(Multi-mode Interference, MMI)组成,光信号经输入层均分4路进入储层;储层由12个定向耦合器(Directional Coupler, DC)作为节点,两节点间由波导(Waveguide, WG)相连,采用螺旋拓扑规律组成梅花型网格结构,具有波长敏感特性,可提高滤波器的调谐灵敏度。输入层输出的4路光信号由储层中心4个节点进入,由于中心节点的4个端口均与波导相连,需在中心节点处额外引入4个DC作为输入节点,输入节点的4个端口分别与输入信号、终端、中心节点端口以及波导相连,其中输入节点分光比设置为0.9,偶数路输出节点分光比设置为0.45,其余节点分光比设置为0.5,储层中波导长度均设置为500 μm,光信号经输入节点进入储层,按拓扑规律沿不同路径传输,并在节点处发生干涉并重新组合,重组后的光信号由储层外围的8个节点输出进入读出层;读出层由8个光调制器(Optical Modulator , OM)以及7个2×1 MMI组成,信号进入读出层的OM后,经过对每路光信号的幅度以及相位进行加权调节,调节后的光信号由2×1 MMI拟合成一路光信号,经光放大器进行功率放大后输出。其中在储层波导上添加微型加热电极,如图1(a)中红色部分,基于硅材料的热光效应,改变波导的相对折射率引起信号的相位变化;在读出层光调制器上下臂上同样添加微型加热电极改变其相位值,实现读出权重的赋值。

      图  1  光子储备池结构。(a)整体结构;(b)输入层结构

      Figure 1.  Structure of photonic reservoir computing. (a) Overall structure; (b) Input layer

      基于集成光子储备池的光子滤波器系统实验装置如图2所示,激光器输出功率为1 mW、中心波长为1550 nm的光信号,经单模光纤(Single Mode Fiber, SMF)和偏振控制器(Polarization Controller, PC)对光信号的偏振态进行控制,采用垂直耦合方式,将输入光通过如图1(b)中输出层的TE光栅耦合器进入光子储备池芯片中,经过光子储备池对信号处理后,通过光栅耦合到输出光纤,最后由光谱分析仪(Optical Spectrum Analyzer, OSA)对光滤波器的FSR、中心波长以及滤波波形进一步分析,芯片中的相移器均由电压源阵列(Voltage Source Array, VSA)进行驱动。

      图  2  光子滤波器系统实验装置图

      Figure 2.  Diagram of photonic filter system experimental setup

    • 引入散射矩阵,对光子储备池各个端口之间的映射关系进行描述[117]。如图1所示,该结构分为输入层、储层、读出层三部分,不考虑端口处的反射现象,每一部分均可用唯一对应的散射矩阵描述。

      1)输入层:输入层由2×1 MMI器件连接而成,MMI是一种无源器件,主要起到功率分配作用,即可作为合束器又可作为分束器,其输入输出端口均为单模波导,具有结构紧凑、易于制作、损耗小、制作容差性好、偏振相关性小等优点。同时可以通过级联多个MMI实现光功率的按比例分配。在输入层中通过级联3个1×2 MMI将一路光信号的功率均分四路进入储层。

      采用的2×1 MMI的耦合比为0.5,1个2×1 MMI具有3个端口,其端口间的映射关系表示为:

      $$ \left( {\begin{array}{*{20}{c}} {{E_{M,t1}}} \\ {{E_{M,t2}}} \\ {{E_{M,t3}}} \end{array}} \right) = \sqrt {{{10}^{0.1{L_M}}}} \left( {\begin{array}{*{20}{c}} 0&{\dfrac{1}{2}}&{\dfrac{1}{2}} \\ {\dfrac{1}{2}}&0&0 \\ {\dfrac{1}{2}}&0&0 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{E_{M,i1}}} \\ {{E_{M,i2}}} \\ {{E_{M,i3}}} \end{array}} \right) $$ (1)

      式中:${E_{M,t1}}$、${E_{M,t2}}$、${E_{M,t3}}$为2×1 MMI端口处的透射光场;${E_{M,i1}}$、${E_{M,i2}}$、${E_{M,i3}}$为2×1 MMI端口处的输入光场;${L_M}$为2×1 MMI器件的损耗,单位为dB。

      根据公式(1)可得到2×1 MMI作为分束器时,其两个输出端口与输入端口的映射关系为:${E_{{\text{out}}1}} = {E_{{\rm{out}}2}} = \dfrac{1}{2}\sqrt {{{10}^{0.1{L_M}}}} {E_{{\rm{in}}}}$,作为分束器时,其输入输出端口的映射关系为:${E_{{\rm{out}}}} = \dfrac{1}{2}\sqrt {{{10}^{0.1{L_M}}}} ({E_{{\rm{in}}1}} + {E_{{\rm{in}}2}})$。

      信号经输入层一分四路,对输入层的1个输入端口、4个输出端口之间的映射关系描述为:

      $$ \left( {\begin{array}{*{20}{c}} {{E_{t1}}} \\ \vdots \\ {{E_{t4}}} \end{array}} \right) = {W_{{\rm{in}}}}{E_{{\rm{in}}}} $$ (2)

      式中:${E_{t1}}$、${E_{t2}}$、${E_{t3}}$、${E_{t4}}$为输入层输出端口处的透射光场;${E_{{\rm{in}}}}$为光子储备池输入端口处的输入光场;${W_{{\rm{in}}}}$为光子储备池输入权重,其值随机产生并固定不变,无需进行训练。根据公式(1),结合输入层中2×1 MMI的连接方式,由传递矩阵法计算${W_{{\rm{in}}}}$得:

      $$ {W_{{\rm{in}}}} = {10^{0.1{L_M}}}\left( {\begin{array}{*{20}{c}} {{{\text{1}} \mathord{\left/ {\vphantom {{\text{1}} {\text{4}}}} \right. } {\text{4}}}} \\ {{{\text{1}} \mathord{\left/ {\vphantom {{\text{1}} {\text{4}}}} \right. } {\text{4}}}} \\ {{{\text{1}} \mathord{\left/ {\vphantom {{\text{1}} {\text{4}}}} \right. } {\text{4}}}} \\ {{{\text{1}} \mathord{\left/ {\vphantom {{\text{1}} {\text{4}}}} \right. } {\text{4}}}} \end{array}} \right) $$ (3)

      2)储层:输入层输出的四路信号进入储层,储层由WG、DC器件连接而成,如图3(a)所示,WG器件具有2个端口,端口间的映射关系表示为:

      图  3  光子储备池中器件模型示意图。(a)波导; (b)定向耦合器;(c)光调制器

      Figure 3.  Schematic diagram of device model in photonic reservoir computing. (a) Waveguide; (b) Directional coupler; (c) Optical modulation

      $$ \left( {\begin{array}{*{20}{c}} {{E_{W,t1}}} \\ {{E_{W,t2}}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} 0&{t\exp (j\phi )} \\ {t\exp (j\phi )}&0 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{E_{W,i1}}} \\ {{E_{W,i2}}} \end{array}} \right) $$ (4)
      $$ t = 1{0^{\tfrac{{ - {L_W} * L * 1{0^{ - 4}}}}{{20}}}} $$ (5)
      $$ \phi = {\phi _2}{{ - }}{\phi _1} = 2\pi L{n_{{\rm{eff}}}}/\lambda $$ (6)

      式中:${E_{W,t1}}$、${E_{W,t2}}$为波导端口处的透射光场;${E_{W,i1}}$、${E_{W,i2}}$为波导端口处输入光场;$t$为波导的透射率;${L_W}$为波导的损耗;$L$为波导长度;$\phi $为波导的相位变化;${\phi _1}$和${\phi _2}$分别为波导两端的相位值;$\lambda $为信号的中心波长;${n_{{\rm{eff}}}}$为波导在$\lambda $处的有效折射率。

      根据公式(4)可以看出信号经过波导后,信号的幅度变化与透射率t有关,而t的值由波导的损耗以及波导长度决定;信号的相位变化与传输信号的中心波长、波导长度以及波导的有效折射率有关,因而通过热调改变硅波导的折射率可对光信号的相位进行调制。

      图3(b)所示,DC器件具有4个端口,各个端口间的映射关系表示为:

      $$ \left( {\begin{array}{*{20}{c}} {{E_{D,t1}}} \\ {{E_{D,t2}}} \\ {{E_{D,t3}}} \\ {{E_{D,t4}}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} 0&{\sqrt {1 - \kappa } }&{\sqrt \kappa \exp \left( {j\dfrac{\pi }{2}} \right)}&0 \\ {\sqrt {1 - \kappa } }&0&0&{\sqrt \kappa \exp \left( {j\dfrac{\pi }{2}} \right)} \\ {\sqrt \kappa \exp\left( {j\dfrac{\pi }{2}} \right)}&0&0&{\sqrt {1 - \kappa } } \\ 0&{\sqrt \kappa \exp \left( {j\dfrac{\pi }{2}} \right)}&{\sqrt {1 - \kappa } }&0 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{E_{D,i1}}} \\ {{E_{D,i2}}} \\ {{E_{D,i3}}} \\ {{E_{D,i4}}} \end{array}} \right) $$ (7)

      式中:${E_{D,t1}}$、${E_{D,t2}}$、${E_{D,t3}}$、${E_{D,t4}}$为DC端口处的透射光场;${E_{D,i1}}$、${E_{D,i2}}$、${E_{D,i3}}$、${E_{D,i4}}$为DC端口处输入光场;$\kappa $为DC耦合系数。

      光信号经储层中心4节点进入储层,按螺旋拓扑路径传输后,经外围的8节点输出。利用散射矩阵,对储层4个输入端口、8个输出端口间的频域关系描述为:

      $$ \left( {\begin{array}{*{20}{c}} {{E_{t5}}} \\ \vdots \\ {{E_{t12}}} \end{array}} \right) = {W_{{\rm{res}}}}\left( {\begin{array}{*{20}{c}} {{E_{i1}}} \\ \vdots \\ {{E_{i4}}} \end{array}} \right) $$ (8)

      式中:${E_{t5}} \cdots {E_{t12}}$为储层输出端口处的透射光场;${E_{i1}} \cdots {E_{i4}}$为储层输入端口处的输入光场;即输入层输出端口处的透射光场;${W_{{\rm{res}}}}$矩阵为光子储备池连接权重,其矩阵大小为(8×4),一经产生就固定不变,同样不需要进行训练。根据信号在储层的传输路径以及公式(4)、(7),由传递矩阵法进行计算,由于计算量非常庞大,${W_{{\rm{res}}}}$具体数值计算需由计算机辅助完成。

      3)读出层:读出层由OM及2×1 MMI组成,其中2×1 MMI作为合束器使用,将7个2×1 MMI进行连接,将输出的8路信号拟合成一路。如图3(c)所示,OM器件具有2个端口,端口间的映射关系表示为:

      $$ \left( {\begin{array}{*{20}{c}} {{E_{O,t1}}} \\ {{E_{O,t2}}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} 0&{\cos\dfrac{{({\varphi _a} - {\varphi _b})}}{2}\exp \left(j\dfrac{{({\varphi _a} + {\varphi _b})}}{2}\right)} \\ {\cos\dfrac{{({\varphi _a} - {\varphi _b})}}{2}\exp \left(j\dfrac{{({\varphi _a} + {\varphi _b})}}{2}\right)}&0 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{E_{O,i1}}} \\ {{E_{O,i2}}} \end{array}} \right) $$ (9)

      式中:${E_{O,t1}}$、${E_{O,t2}}$为OM端口处的透射光场;${E_{O,i1}}$、${E_{O,i2}}$为OM端口处输入光场;${\varphi _a}$、${\varphi _b}$分别为OM上、下干涉臂的相位值。

      根据公式(9)可以看出信号经过光调制器后,信号的幅度变化与调制器上下臂相位差值有关,信号的相位变化与调制器上下臂相位和有关,因而通过热调改变调制器上下臂中硅波导的折射率可实现光信号的相位和幅度的调整,从而将训练权重赋值到储层输出的各路光信号。

      由储层输出的8路信号进入读出层后,通过OM对每路信号分别赋予合适的权重,经2×1 MMI拟合成1路输出,对读出层输入输出端口间的映射关系描述为:

      $$ {E_{{\rm{out}}}} = {W_{{\rm{out}}}}\left( {\begin{array}{*{20}{c}} {{E_{i5}}} \\ \vdots \\ {{E_{i12}}} \end{array}} \right) $$ (10)

      式中:${E_{{\rm{out}}}}$为光子储备池输出端口处的光场;${E_{i5}} \cdots {E_{i12}}$为读出层输入端口的光场,即储层输出端口的透射光场;${W_{{\rm{out}}}}$为光子储备池的读出权重,其值需要对储层输出信号进行训练得到。根据公式(1)、(9)以及器件在读出层的连接方式,计算${W_{{\rm{out}}}}$得:

      $$ {W_{{\rm{out}}}} = \frac{1}{8}{\left( {{{10}^{0.1{L_M}}}} \right)^{\tfrac{3}{2}}}{\left( {\begin{array}{*{20}{c}} {{W_{{\rm{Amp1}}}}\exp \left( {j{\varphi _{w1}}} \right)} \\ \vdots \\ {{W_{{\rm{Amp8}}}}\exp \left( {j{\varphi _{w8}}} \right)} \end{array}} \right)^{\rm{T}}} $$ (11)
      $$ {W_{{\rm{Amp}}i}} = \cos \frac{{({\varphi _{{\rm{a}}i}} - {\varphi _{{\rm{b}}i}})}}{2} $$ (12)
      $$ {\varphi _{{\rm{w}}i}} = \frac{{({\varphi _{{\rm{a}}i}} + {\varphi _{{\rm{b}}i}})}}{2} $$ (13)

      式中:${W_{{\rm{Amp}}i}}$、${\varphi _{{\rm{w}}i}}$分别为读出层中第$i$路光信号赋值的幅度权重及相位权重;${\varphi _{{\rm{a}}i}}$、${\varphi _{{\rm{b}}i}}$分别为第$i$路OM上、下干涉臂的相位值。利用训练得到${W_{{\rm{Amp}}i}}$及${\varphi _{{\rm{w}}i}}$值,通过热调对OM器件上、下臂的${\varphi _{{\rm{a}}i}}$、${\varphi _{{\rm{b}}i}}$参数进行调整。

      结合公式(2)、(8)、(10),整个光子储备池输入端口与输出端口之间传输响应关系表示为:

      $$ {E_{{\rm{out}}}} = {W_{{\rm{out}}}}{W_{{\rm{res}}}}{W_{{\rm{in}}}}{E_{{\rm{in}}}} $$ (14)

      由上式可知,通过对光子储备池内读出权重、连接权重、输入权重的确定,得到整个结构输入端与输出端的频域响应。

    • 利用PSO优化算法对读出权重的幅度${W_{{\rm{Amp}}i}}$、相位${\varphi _{{\rm{w}}i}}$共16个参数进行寻优,其中${W_{{\rm{Amp}}i}}$的参数寻优范围为(0~1),${\varphi _{{\rm{w}}i}}$参数寻优范围为(0~2π),设置初始种群数为50,惯性权重为0.8,加速常数均为2,空间维数为16,通过不断更新粒子速度和位置寻找损失函数的最小值,得到最优的权重参数值。利用损失函数定量判断训练效果的优劣,当损失函数值越小,训练得到的权重值越接近理想值,将损失函数定义为:

      $$ Fitness = \frac{1}{n}\min \sum {{{\left[ {\alpha (\left| {{y_{{\rm{pred}}}}} \right| - \left| {{y_{{\rm{label}}}}} \right|)} \right]}^2}} $$ (15)

      式中:$n$为总的采样点数目;$\left| {{y_{{\rm{pred}}}}} \right|$为训练波形幅值;$\left| {{y_{{\rm{label}}}}} \right|$为理想波形幅值;$\alpha $为调节系数,将理想波形幅度值与训练波形幅度值之间的误差进行放大或缩小,最终拟合出理想波形效果。其$\alpha $值描述为:

      $$ \alpha = \left\{ \begin{gathered} a,{\left| {{y_{{\rm{label}}}}} \right|_{\max }} \geqslant \left| {{y_{{\rm{label}}}}} \right| \geqslant x \\ b,x > \left| {{y_{{\rm{label}}}}} \right| \geqslant {\left| {{y_{{\rm{label}}}}} \right|_{\min }} \\ \end{gathered} \right. $$ (16)

      式中:${\left| {{y_{{\rm{label}}}}} \right|_{\max }}$为理想幅值的最大值;${\left| {{y_{{\rm{label}}}}} \right|_{\min }}$为理想幅值的最小值;$a$与$b$为$\alpha $参数在$\left| {{y_{{\rm{label}}}}} \right|$不同范围内的设定值;$x$为设定在${\left| {{y_{{\rm{label}}}}} \right|_{\max }}$与${\left| {{y_{{\rm{label}}}}} \right|_{\min }}$之间的任意实数值。在拟合不同的波形时,无需改变光子储备池的结构以及器件参数值,只需调整$a$、$b$、$x$的取值,通过对权重进行训练优化,得到所需要的滤波效果,从而实现基于光子储备池的光滤波器可重构特性。

      与其他光波导网络的算法训练过程相比,仅对光子储备池的读出权重${W_{{\rm{out}}}}$进行训练,无需训练其输入权重及储层连接权重,大大减少了训练过程中训练负担。其次,采用粒子群算法训练时,可通过当前搜索的最优值寻找全局最优,这种算法具有训练规则简单、易实现、精度高、收敛快的优势。

    • 根据图1搭建仿真模型,考虑储层中波导的实验制造误差,在波导上添加−π~π范围内的随机相位,对储层输出的8路信号进行仿真。

      图4所示,每路信号光谱的FSR相同,均为1.18 nm,其中第1路信号与第2路信号光谱形状相同,但第2路信号光功率高于第1路信号。由图1可知,信号在储层传输过程中,第1路输出信号的传输路径与第2路输出信号路径相似,差别在于两路输出节点处连接的波导路径,由于波导固有损耗,使得两路光谱所对应的光功率不同,但光谱形状不变。由图4(c)~(h)可知,第3、4路、第5、6路、第7、8路均有相同的规律。

      图  4  储层输出的8路信号光谱图。(a)第1路;(b)第2路;(c)第3路;(d)第4路;(e)第5路;(f)第6路;(g)第7路;(h)第8路

      Figure 4.  8 channel signal spectra of reservoir output. (a) Rout 1; (b) Rout 2; (c) Rout 3; (d) Rout 4; (e) Rout 5; (f) Rout 6; (g) Rout 7; (h) Rout 8

      为解决目前滤波器所实现的滤波波形单一问题,文中提出的基于光子储备池的光滤波器利用粒子群算法对读出层训练,寻找多组读出权重参数,实现两种滤波波形,此节以经典的IIR型及FIR型两种不同的滤波功能实现为例进行说明。对图4中的8路信号进行线性拟合,通过对读出层权重训练以及$\alpha $值的相关参数调整,实现IIR型和FIR型光滤波器。如图5图6所示,分别为IIR型光滤波器与FIR型光滤波器的仿真结果。

      图  5  IIR型光滤波器仿真训练效果。(a) IIR滤波谱拟合效果图;(b)Fitness曲线

      Figure 5.  Simulation training effect of IIR optical filter. (a) IIR filter spectrum fitting rendering; (b) Fitness curve

      图  6  FIR型光滤波器仿真训练效果。(a) FIR滤波谱拟合效果图;(b) Fitness曲线

      Figure 6.  Simulation training effect of FIR optical filter. (a) FIR filter spectrum fitting rendering; (b) Fitness curve

      对IIR型光滤波器进行分析,由公式(16)知,$\alpha $值由$a$、$b$、$x$ 3个参数确定,通过多次调整$a$、$b$、$x$的值并进行仿真,找到最小Fitness值所对应的$\alpha $值以及权重值,如表1表2所示。

      表 1  IIR型光滤波器对应α值的设定参数

      Table 1.  Setting parameters of the corresponding α of IIR optical filter

      |ylabel|max|ylabel|minabx
      10.054.510.95

      图5(a)所示,理想波形与训练波形基本完全贴合,在波谷处有少许波动,但仍能实现IIR型滤波器的滤波功能,实现的IIR光滤波器的FSR为1.18 nm,3 dB带宽为0.18 nm,消光比为30 dB左右。图5(b)为误差值随迭代次数的变化,随着迭代次数的增加,误差值呈阶梯状下降趋势。在训练初始阶段,由于粒子本身具有较强的扩展搜索空间能力,其寻找最优解的速度比较快,随着迭代次数增加,误差值快速下降;在迭代过程中,粒子群算法容易陷入局部最优空间,随着迭代次数增加,Fitness值下降缓慢,甚至保持不变,当粒子跳出局部最优时,Fitness值进一步下降;在迭代结束阶段,Fitness值在0.022附近基本保持平稳。由此可以看出,PSO算法的参数搜索速度快,但容易陷入局部最优,需要对加速常数值以及惯性权重进行合理配置,平衡粒子的全局和局部搜索能力。

      对FIR型光滤波器进行权重拟合,理想波形的幅度范围在(0~0.81)之间,通过多次仿真对比,找到最小Fitness值所对应的$\alpha $相关参数以及权重值,如表3表4所示。

      表 2  IIR型光滤波器8路信号权重训练值

      Table 2.  IIR optical filter 8 channel signal weight training value

      IIR1st2nd3rd4th5th6th7th8th
      WAmp0.4960.0150.8470.3580.2400.6490.8650.641
      φw4.7170.5935.011.1764.0720.1732.1801.935

      表 3  FIR型光滤波器对应α值的设定参数

      Table 3.  Setting parameters of the corresponding α of FIR optical filter

      |y’label|max|y’label|mina’b’x’
      0.810240.4

      表 4  FIR型光滤波器8路信号权重训练值

      Table 4.  FIR optical filter 8 channel signal weight training value

      FIR1st2nd3rd4th5th6th7th8th
      WAmp0.9660.4520.9980.9280.6850.9910.9970.766
      φw3.2582.8613.0975.1893.0800.2596.2825.246

      图6所示,为训练得到的FIR型滤波谱以及Fitness曲线。由FIR光谱图看出,训练得到的光谱波形与理想波形相比,在波形顶部两者相差较大,其余部分贴合较好。实现的FIR光滤波器的FSR为1.18 nm,3 dB带宽为0.77 nm,消光比为50 dB左右。顶部毛刺较多的原因在于光信号在储层中既有前向传输又有反馈配置,不同路径的信号在DC节点处发生多次干涉并重组。

      图4所示,由储层输出的每路光信号均具有不同程度的毛刺,仅通过读出层的光调制器对8路信号进行权重赋值并线性拟合,其训练波形的波动很难消除。由Fitness曲线图看出,其下降趋势与IIR型滤波器的Fitness曲线相近,误差值随迭代次数的增加呈阶梯状趋势下降,Fitness值在0.11附近基本保持平稳。对比图5(a)与图6(a),可以看出,FIR型训练波形顶部波动比IIR型训练波形底部波动大,原因在于波形顶部相较于底部所对应的光功率较大,权重值对于波形顶部影响较大。

      通过对IIR型及FIR型光滤波器仿真结果分析可知,在给定理想目标波形情况下,利用粒子群算法对储层的输出信号进行权重训练,使训练得到的滤波波形不断逼近理想波形,进而获得最优权重系数,当滤波波形无限接近时,所实现的光滤波器相关滤波参数便与期望实现的滤波器参数基本一致,从而所设计的该光滤波器的相关滤波参数精确度较高。

    • 以IIR型的训练波形为例,通过对储层中WG、DC器件参数调整,分析WG长度、DC分光比对滤波谱参数的影响。

      图7所示,给出在不同波导长度下IIR型光滤波谱的变化情况,从光谱图看出,随着波导长度的变化,滤波谱形状保持不变,滤波谱的FSR随着波导长度的增加而减小。对图7的仿真数据进一步分析,得到IIR型光滤波谱FSR随波导长度的变化。如图8所示,波导长度与FSR成负相关,随着波导长度的增加,相应的FSR变小。对数据点进行多项式拟合,得到FSR随波导长度变化曲线,可以看出,FSR与波导长度并不符合严格的线性关系。由此可见,波导长度是影响FSR的主要因素,根据实际应用要求,通过对储层中波导长度的调整实现对FSR的灵活控制。

      图  7  储层中波导在不同长度下的IIR滤波谱。(a)波导长度为400 μm ;(b)波导长度为500 μm ;(c)波导长度为600 μm

      Figure 7.  IIR filter spectrum of waveguide at different lengths in reservoir. (a) Waveguide length of 400 μm; (b) The waveguide length of 500 μm; (c) The waveguide length of 600 μm

      图  8  FSR随波导长度变化

      Figure 8.  FSR changing with waveguide length

      保持输入节点DC分光比与储层奇数路输出节点DC分光比不变,对储层偶数路输出节点DC的分光比进行调整。由图9可以看出,DC分光比由0.45变化至0.75过程中,光谱形状保持不变,但对于波形峰值功率影响较大,光功率由1 mW降低至0.67 mW。由此可见,根据实际需要,通过对储层偶数路输出节点分光比进行设置,实现不同滤波强度的调节。

      图  9  不同定向耦合器分光比滤波器输出功率随波长变化

      Figure 9.  Output power of filter changing with the wavelength in different directional coupler radio

    • 在光滤波器滤波波长可调谐仿真中,通过对光滤波器的透射谱分析,得到波导上的相位对滤波波长偏移的影响。

      图10给出在不同相移量下光滤波器透射端光功率随波长的变化。从光谱图中显示,添加在波导上相位变化(0~3π/2)时,波形形状保持不变、自由光谱范围保持为1.18 nm,但滤波中心波长实现了在一个自由光谱范围内连续可调谐。与目前基于Sagnac环结构实现的光滤波器中心波长调谐范围在0.5 nm内相比,该滤波器的中心波长的调谐范围提高两倍。因此,利用热调方法对储层波导上的相位进行控制,可以实现光滤波器中心滤波波长的可调谐。

      图  10  不同相移量滤波器输出功率随波长变化

      Figure 10.  Output power of filter changing with the wavelength in different phase shifts

    • 文中构建了一种基于12节点梅花型储备池的新型可调谐光滤波器,利用散射矩阵分析储备池输入输出端口之间的映射关系,通过PSO算法对读出层的各路光信号幅度和相位参数进行权重训练并拟合,实现IIR型与FIR型滤波器功能。通过对储层波导长度进行调整,实现FSR的灵活控制;调整储层中波导上相位分布(0~3π/2),在滤波波形不变的情况下实现光滤波器滤波波长在一个自由光谱范围内的连续可调。文中通过理论与仿真对该光滤波器进行可行性验证,并在仿真过程中充分考虑实际实验中的数据偏差。与MRR型光滤波器相比,文中设计的基于片上集成光子储备池的光滤波器由波导进行连接而成,具有制作难度小,更易于落地生产的优势。光子储备池相比于MZI型网络,在实现不同滤波功能时,光信号在储层中的传输路径保持不变,避免了更改路径配置增加的系统额外功耗,下一步也将对读出层中MZI光调制器的原理和技术进一步改进,找到克服相关技术的方法。由于目前光子储备池芯片正处于流片阶段,且周期较长,在未来1~2年内,笔者将进行光滤波器芯片的实验测试验证。

参考文献 (17)

目录

    /

    返回文章
    返回