江嘉乔，鲁冰新，肖天亮，翟 锦

（北京航空航天大学化学学院,仿生智能科学与技术教育部重点实验室,北京 100191）

纳米离子通道具有离子选择性、离子整流性和离子门控性3个重要的特征.其中，离子整流性（即离子在纳米通道内定向运输）是指在一定的外加电势下，离子优先从通道的一侧向另一侧进行输运，反过来几乎不能通过，表现为类似二极管的整流性能.当纳米离子通道用于构造具有整流性的器件时，可在nA级别调控电流，此时离子通道体系在相反的外加电势下可以测得非对称的离子电流，因此，也被称为纳米流体二极管（Nanofluidic diode）［1～3］.基于连续电荷的假设，纳米离子通道研究有许多理论模型.其中，Siwy等［4～7］对锥形纳米通道电势分布进行了研究并提出静电场的棘轮模型；Woermann等［8］则对锥形单纳米通道的离子整流行为提出了一种半定量的理论模型，提出了拟合实验数据电流的公式；Cervera等［9］提出了基于泊松-能斯特-普朗克（Poisson-Nernst-Planck，PNP）方程组的模型，该模型定量描述锥形通道中离子电流与电压的关系.利用计算机在离子通道二维模型中求解PNP连续方程组最早由Daiguji 等［10］提出，它可以给出模型中电势、离子浓度等物理量的分布.Cheng 等［11］对具有纳米流体二极管性质的多种离子通道实验体系和理论做了介绍和总结.

在描述离子通道性质的PNP方程模型中，泊松（Poisson）方程是描述静电场分布的通用方程.而在恒温条件下，对于纳米离子通道，使用能斯特-普朗克（Nernst-Planck）方程进行描述溶液中载流子的定向运动，则是基于电化学液相传质理论中的扩散-电迁移模型［12，13］.根据经典的半导体载流子输运理论，在恒温条件下，半导体导电是由于电子、空穴两种载流子通过浓度梯度形成的扩散作用和电场下的漂移作用产生定向运动形成电流，该模型称为扩散-漂移（Diffusion-drift）模型［14，15］.比较两种模型，半导体的扩散电流与溶液的扩散电流之间、半导体的漂移电流和溶液的电迁移电流之间有相同的描述形式，即各自的载流子输运都可以通过能斯特-普朗克方程描述.

PN结二极管是一种典型的半导体［14～16］.对于理想二极管模型，通过二极管的电流对其两端电压呈指数函数，当反向偏压逐渐增大且未击穿时，对应反向电流将趋于饱和电流（Saturated current）.仿真电路模拟器（Simulation Program with Integrated Circuit Emphasis，SPICE）模型是一种通用的电路模拟模型，它可以用于求解包含二极管的电路，且已被用于模拟神经元离子通道的等效电路［17，18］.在实际的固体半导体器件中，一般都存在偏离理想二极管模型的情况，使得在反向偏压下存在一定的额外电流［15］.Antognetti等［19］提出在二极管电路的SPICE模型中，对于偏离理想模型的电流，可以由包含反向饱和电流（Reverse saturated current）的项描述.

在离子通道整流性测定实验中，通常将具有离子通道的薄膜夹在两个装有相同溶液的电解池之间组成实验装置，并在一定的外加电势范围（如-2～+2 V）内测量I-V特性曲线，该曲线与半导体二极管电流-电压（I-V）特性曲线有相同的非对称和通过原点的性质.因此，本文比较固体PN结二极管半导体电流和锥形离子通道实验装置离子电流的产生原理，提出可以参考PN结二极管半导体载流子输运的理论，使用SPICE模型建立的等效二极管电路来拟合离子通道整流性的I-V曲线.离子通道实验测量I-V曲线来表征整流性时，在实验装置施加外加电势范围内通常不会出现较大反向偏压下的击穿效应和较大正向偏压下的大电流饱和效应.在实验的反向偏压范围内，此时离子通道对离子输运的阻力较大，因此测得较小的离子电流.可以注意到随着反向偏压增大，离子电流线性增大，而不是类似于理想PN结二极管的I-V曲线趋近于下界的特性.这种I-V曲线的特性既呈现在聚合物锥孔离子通道，也呈现在多种类型材料复合而成的异质结离子通道，且可以存在于低整流比和高整流比的体系［20～27］.

对于锥形离子通道体系，基于PNP 方程组的模型可以较好地模拟实验测得的I-V特性曲线［20，21，28～31］.因此，本文以PNP模型模拟一定几何外型、表面电荷密度均匀的锥形离子通道在多种条件下的I-V特性曲线，并通过PNP模型的模拟结果分析离子通道体系根据能斯特-普朗克方程描述的离子电流的组成.以此为依据，进一步在SPICE模型中将离子通道装置作为电路中的一个二极管器件模拟I-V特性曲线，并将该结果同PNP 模型的结果进行比较，阐述离子通道体系类似PN 结二极管的特性.离子通道SPICE电路模型的建立可用于研究纳米流体二极管作为一个器件在电路中的应用.

1 实验部分

1.1 PN结二极管半导体I-V特性曲线SPICE模型

固体半导体载流子输运的扩散-漂移模型，对于电子n和空穴p各自的扩散流密度Ji（mol·m-2·s-1）由下式描述：

电子n和空穴p各自在电场E作用下形成的漂移流密度Ji由下式描述：

式中：Di（m2/s）为扩散系数；Ni（mol/m3）为电子/空穴浓度；ei为电子/空穴电荷数；μi（m2·V-1·s-1）为电子/空穴迁移率；i表示电子/空穴种类；φ（V）表示电势.自由电子贡献电流与通量方向相反，自由空穴贡献电流与通量方向相同.

在PN结二极管中，无外加电势时，PN结界面附近会自发形成电子-空穴平衡.对于理想的固体PN结半导体，当处于外加电势作用下，界面处的电子-空穴平衡被打破而产生扩散电流.在一定的外加电势范围内，不考虑因外加电势过大导致反向击穿电压效应或正向注入大电流饱和效应的情况，理论推导的通过理想PN结二极管器件电流I（A）表达式如下：

式（3）中，q=1代表电子n和空穴p具有基本电荷量；对于理想模型，参数N=1；kB为玻尔兹曼常数；T（K）为温度，此处固定为室温293.15 K；Va（V）为PN结两侧的外加电势；IS（A）为饱和电流.式（4）中，Dn和Dp分别为电子n和空穴p的扩散系数；Ln（m）和Lp（m）分别为电子n和空穴p的扩散长度；和分别为P型半导体平衡状态下的电子n的浓度和N型半导体平衡状态下的空穴p的浓度.对于一定温度下特定的PN结材料，它们属于固有属性，决定了饱和电流密度j0（A/m2），式（4）中括号内两项分别描述电子n向P型区的扩散和空穴p向N型区的扩散行为.Sarea（m2）为具体器件的横截面积［14，16］.

在SPICE模型语言中，二极管的串联电阻可以作为二极管器件的属性描述，因此仅需在电路模型中定义一个直流电源器件、一个二极管器件以模拟I-V曲线.对于理想PN结二极管器件，不考虑反向击穿电压（模型符号为BV）以及正向注入大电流饱和（模型符号为IKF），需要定义IS和定值串联电阻［RS（Ω）］两个参数，使用式（3）对通过理想PN结二极管器件的电流I进行计算［17，19］.Scheme 1给出了理想PN 结二极管器件的工作原理及求解所得的I-V特性曲线，可见，当反向偏压逐渐增大且未击穿时，通过PN结的反向电流将趋于IS.

Scheme 1 Schematic diagram of interface electric field and charge carrier transportation of ideal PN junction on negative bias(A) and positive bias(B),and corresponding I-V characteristic curve(C)

在实际的固体半导体器件中，通常存在偏离理想二极管模型的情况，因而，在反向偏压下存在一定的额外电流［15］.此外，如果存在一定的外界作用使得较大数量的少数载流子（P区中的电子n和N区中的空穴p）并非由于PN结界面平衡作用生成，且可以到达PN结的界面，被PN结的强电场拉向界面的另一侧从而导向外电路，则此时在反向偏压下也可以形成较大的电流.两个常用器件的例子可以说明这种情况：（1）PNP型三极管集电区部分的PN结，如Scheme 2（A）所示，对于PNP型三极管而言，中间的N型基区很薄，在工作状态时，发射区PN结处于正向偏压下，生成大量载流子空穴pE向集电区发射并能避免在N型基区复合，而集电区PN结处于反向偏压下，其界面的强电场可以将载流子空穴pE快速拉向另一侧P区并导出形成大电流，可见此时集电区PN结区域有大量载流子由外部引入，并在电场施加的漂移作用下运动；（2）PN结光生伏特效应，如Scheme 2（B）所示，PN结界面吸收光照受到激发而产生电子-空穴对并产生电势，电子-空穴对被PN结界面强电场分离形成光电流，并在电势作用下通过漂移作用分别导向外电路，其光生电子nE、光生空穴pE的运动方向相对于PN 结界面电场的方向与Scheme 1中理想PN结反向偏压下的情况相同［14，16］.

Scheme 2 Scheme of electron and hole transportation in a working PNP transistor(A),scheme of electron and hole transportation in a working photovoltaic PN junction(B)

对于PN 结二极管偏离理想模型的电流，在SPICE 电路模型中，Antognetti 等［19］提出的包含反向饱和电流ISR（A）项，由下式描述：

式中：q=1；参数NR=2；T固定为室温293.15 K；VJ（V）为PN 结内建电场的电势；M为PN 结的梯度系数，M＞1/2表示PN结为突变结，M≤1/3表示PN结为缓变结；反向饱和电流ISR参数由具体器件决定.

式（5）属于构造公式，并不明确地对应实际的物理意义.但它需要给出VJ，M和ISR，同样属于器件的固有属性.该式表明，在较大的反向偏压下PN结界面的电场明显增大，由ISR描述的项会被PN结界面电场快速分离，避免了电子-空穴复合效应的发生，通过二极管器件的电流会逐渐增大［17，19］.

在一定的外加电势范围内，不考虑较大反向偏压下的击穿效应以及较大正向偏压下的大电流饱和效应，由式（3）和式（5）两部分共同构成恒温条件下SPICE 模型描述的二极管器件I-V曲线的电流［JD（A）］，如下式所示：

SPICE模型可由COMSOL Multiphysics软件使用“电路”模块依据电路原理图建立和求解，其作为通用模型也可以由其它软件求解.

1.2 锥形离子通道体系PNP模型

对于描述离子通道性质的PNP方程模型，以常用的纳米离子通道在KCl浓度为1～100 mmol/L的中性水溶液体系为例，此时，H+/OH-的浓度可以忽略不计，电解质溶液中的离子电流大小主要由K+和Cl-两种载流子决定，能斯特-普朗克方程描述的通量密度J由扩散通量密度和电迁移通量密度两部分组成，如下所示：

式中：ui（m2·V-1·s-1）为离子淌度；ci（mol/m3）为离子浓度；zi为离子电荷数；i表示离子种类；F为法拉第常数；R为理想气体常数；T固定为室温293.15 K.将离子通量换算为离子电流时，正离子贡献电流与通量方向相同，负离子贡献电流与通量方向相反.

对于锥孔离子通道在10 mmol/L 中性KCl 水溶液体系的模型，假设锥孔通道表面电荷密度均匀分布，且不与K+或Cl-发生化学反应.讨论该体系时，改变离子通道表面电荷密度，而不改变锥孔通道的几何外形.此时，离子通道体系的载流子为K+和Cl-，其分布主要受到通道表面电荷和外加电势的影响.

锥形离子通道PNP模型使用COMSOL Multiphysics软件建立.如图1所示，其几何外型为一个锥孔通道的二维模型，参考常用于制备锥孔的聚对苯二甲酸乙二醇酯（PET）和聚酰亚胺（PI）薄膜，包括一个长度10 μm、左侧小孔端为10 nm、右侧大孔端为400 nm的锥孔［20～22］.

锥形通道两侧为两个长500 nm、宽1000 nm的缓冲池.二维模型的面外厚度设置为1 m，以此换算模型的输出电流.其中，模型左侧边界定义为工作电极电势（φL）及其所在电解池体相溶液浓度（cL），右侧边界定义为对电极接地电势（φR）及其所在电解池体相溶液浓度（cR）；代表离子通道表面的边界具有均匀的表面电荷密度（σ）.

在COMSOL Multiphysics软件中求解时，使用软件的“静电”和“稀物质传递”模块.离子通道在KCl中性水溶液体系中，数值模拟使用PNP方程组与时间无关的稳态形式，如下所示：

式中：K+和Cl-的离子扩散系数分别为1.96×10-9和2.03×10-9m2/s；ε为水溶液介电常数，相对介电常数为80；T固定为室温293.15 K.

模型左侧（图1红色）边界条件代表工作电极和左侧电解池KCl体相溶液浓度，电势定义为仪器的外加电势φA（V），电解池的体相溶液浓度为c0（mol/m3），边界条件由如下两式描述：

而右侧（图1蓝色）边界代表接地电极的电势和右侧电解池体相浓度，边界条件如下两式描述：

其中，工作电极电势φA在-2～+2 V 范围之间，体相溶液浓度c0根据具体模型分别设定为1，10，100 mmol/L.

模型的输出电流I（A）为以上左侧或右侧边界各离子通量积分并换算后的总和，如下式所示：

式中：l（m2）为考虑面外厚度的模型左侧边界.

给定代表离子通道各表面的模型边界（图1黄色）具有均匀的表面电荷密度σ0（C/m2），该离子通道表面对溶液体系产生的电势梯度符合静电学高斯（Gauss）定律，由如下两式描述：

式中：σ（C/m2）为离子通道的表面电荷密度；n为该处通道表面的单位法向量.其余模型边界（图1 黑色）表面电荷密度和法向通量均为零.

2 结果与讨论

2.1 10 mmol/L KCl溶液锥形离子通道体系PNP模型对I-V曲线的模拟

设置电解池KCl溶液体相浓度为10 mmol/L、锥形通道表面电荷密度为-5 mC/m2，软件对PNP模型求解可得该离子通道体系在外加电势-2～+2 V范围内的I-V特性曲线，另外，模拟结果还给出了由能斯特-普朗克方程得到的扩散电流和电迁移电流的组成，如图2（A）所示.在±2 V外加电势下离子通道体系轴线上的电势和离子浓度分布如图2（B）所示.

Fig.2 Simulated I-V curve and its composition of diffusion and electrophoresis currents(A),electric voltage and ions concentrations distributions on the axis of channel system under ±2 V applied voltages of the 10 mmol/L bulk concentration,-5 mC/m2 surface charge density conical ion channel PNP model(B)

由PNP 模型I-V曲线模拟结果可知，-2 V 时总离子电流为-0.849 μA，其中，扩散电流为0.004 μA，电迁移电流为-0.853 μA；2 V时总离子电流为8.405 μA，其中，扩散电流为-0.007 μA，电迁移电流为8.412 μA.可见，在-2～+2 V外加电势范围内，扩散作用贡献的离子电流十分微小，而电迁移作用贡献的离子电流占据主导地位.

如果考察轴线上电势和离子浓度的分布，在-2 V 外加电势下通道内小孔端附近正负离子逐渐排空，而电势逐渐升高；在2 V下通道内小孔端附近两种离子都迅速富集，且此处电势快速下降，表明此时电场强度较大.

Nernst-Planck方程扩散电流大小由两种离子浓度变化之差决定，电迁移电流大小由两种离子浓度之和与电场强度的乘积决定.由于K+和Cl-浓度在-2 V或2 V下的分布没有显著差别，因此扩散电流较小；而±2 V下小孔端附近K+和Cl-浓度之和与电场强度的乘积描述的电迁移电流可以描述相反偏压下的整流性.

对比PN结二极管半导体中始终存在快速的电子-空穴平衡反应，电子在P区或空穴在N区运动时容易受到复合作用的影响.在10 mmol/L KCl 中性水溶液中则几乎不存在因化学平衡改变而产生的电流，且离子通道表面电荷吸附相反离子的电容效应而产生的电流也可以忽略不计，因此扩散电流的贡献十分微小；另一方面，离子通道两侧的电解池是浓度稳定的载流子的外部来源，在外加电势下，它们的电迁移运动成为了离子电流的主要来源.这一行为十分类似偏离理想模型的PN结中外部来源的载流子形成漂移电流的行为.

进一步改变该PNP 模型的表面电荷密度（-10，-5，-2 和-1 mC/m2）时的I-V曲线如图3（A）所示，相应±2 V 外加电势下整流比为14.1，9.90，4.69 和2.60.对于这一小孔端直径10 nm 的锥形通道在10 mmol/L KCl中性水溶液的体系而言，较大的离子通道表面电荷密度可以产生较大的整流比.同时，这些条件下电迁移作用贡献的离子电流占据主导地位（图S1，见本文支持信息）.

Fig.3 Simulated I-V curves(A) and corresponding electric voltage distributions on the axis of channel system(B) under 0 V applied voltages of the 10 mmol/L bulk concentration conical ion channel PNP model with different surface charge densities

2.2 10 mmol/L KCl溶液锥形离子通道体系SPICE模型对I-V曲线的模拟

由以上讨论可知，在本文离子通道模型中，与理想PN结相似的因载流子平衡被打破而产生的扩散电流十分微小，电迁移作用贡献的离子电流占据主导地位.因此，为了简化需要拟合的参数，在应用二极管器件SPICE模型构造本锥形离子通道I-V曲线测定的电路模型［图4（A）］时，考虑设定饱和电流IS为零，而仅使用含反向饱和电流ISR的项进行描述.

本文SPICE电路模型由电压源V1和二极管D1两个器件组成.其中，电压源V1表示离子通道装置的外加电势，范围设定为-2～+2 V；二极管D1表示离子通道体系，其中包括一个恒定的串联电阻RS.模型的电路示意图如图4（A）所示.此时，二极管D1器件的电流特性仅包含反向饱和电流ISR的项，如下式所示：

根据此假设，在COMSOL Multiphysics软件中求解SPICE模型，对锥形离子通道I-V曲线进行模拟，已知参数正电荷载流子K+和负电荷载流子Cl-的电荷数q=1，NR=2，玻尔兹曼常数kB，恒定温度（T=293.15 K）条件.另外，还需要给出PN结内建电场电势VJ，PN结梯度系数M，反向饱和电流ISR以及串联电阻RS.

对于内建电势，PNP 模型给出0 V外加电势下轴线上的电势分布，如图3（B）所示，在离子通道内小孔端附近电势会快速下降，这个电势差可以用于估计二极管模型的内建电势VJ.对于串联电阻RS，它在二极管器件SPICE 模型中为一个定值，考虑到在2 V 外加电势下电流较大，此时PN 结阻值较小，缓冲池阻值有重要作用，因此可以粗略估计RS≤2 V/I2V，通过拟合PNP模型在2 V外加电势时计算得到的电流确定串联电阻RS.当处于反向偏压范围内，需要通过设置反向饱和电流ISR和PN 结梯度系数M求解SPICE模型的I-V曲线，从而对PNP模型的计算结果进行拟合.

Fig.4 Scheme of SPICE model for the conical ion channel(A),simulated I-V curves of the 10 mmol/L bulk concentration conical ion channel SPICE model with different surface charge densities(B)

对于表面电荷密度为-5 mC/m2的锥形通道，其串联电阻RS通过PNP模型在2 V外加电势下总离子电流8.405 μA估计在小于238 kΩ的范围附近为215 kΩ，内建电势VJ通过图3（B）得到为7.4 mV.而后指定ISR为0.054 μA、M为0.5进行求解得到SPICE模型的I-V曲线.本SPICE模型的导出网表见本文支持信息.当表面电荷密度较大时，0 V外加电势下离子通道轴线上电势变化更陡峭，因此，考虑较大的PN结梯度系数M进行拟合；相应表面电荷密度较小时，降低PN结梯度系数M.同样方法通过SPICE模型拟合出的锥形通道表面电荷密度为-10，-5，-2和-1 mC/m2时的I-V曲线［图4（B）］.

体相浓度10 mmol/L下不同表面电荷密度的PNP模型计算得到的±2，±1，±0.5 V下的离子电流，以及对应SPICE模型求解使用的参数和计算得到的±2，±1，±0.5 V下的电流列于表1中.可见，特别是在反向偏压下离子输运阻力较大时，由于离子通道内的电阻远大于电解池的电阻，仅考虑反向饱和电流项的SPICE二极管模型对该10 mmol/L浓度下锥形离子通道PNP模型的I-V曲线有很好的拟合结果.在正向偏压下，由于实际实验中缓冲池尺寸远大于PNP模型的设置，但又非SPICE模型中定值串联电阻的情况，I-V曲线介于两种模型之间［20，21］.

Table 1 Parameters for SPICE model and simulated results of PNP and SPICE models of the 10 mmol/L bulk concentration conical ion channel

2.3 1 mmol/L和100 mmol/L KCl溶液锥形离子通道体系I-V曲线模拟结果比较

使用相同方法，对低浓度1 mmol/L KCl溶液体系的PNP模型进行模拟计算.离子通道几何模型不变，设定其表面电荷密度为-1，-0.5，-0.2和-0.1 mC/m2时的I-V曲线如图5（A）所示，相应±2 V下整流比为13.7，9.12，4.22，2.32，此时，较小的表面电荷密度可以得到相近的整流效应.模拟结果也表明本组模型电迁移作用贡献的离子电流占据主导地位（图S2，见本文支持信息）.

Fig.5 Simulated I-V curves of PNP model(A),corresponding electric voltage distribution on the axis of channel system(B) under 0 V applied voltages,simulated I-V curves of SPICE model(C) of the 1 mmol/L bulk concentration conical ion channel with different surface charge densities,simulated I-V curves of PNP model(D),corresponding electric voltage distributions on the axis of channel system(E) under 0 V applied voltages,simulated I-V curves of SPICE model(F) of the 100 mmol/L bulk concentration conical ion channel with different surface charge densities

同样，应用SPICE 模型求解I-V曲线，需要的内建电势由PNP 模型在0 V 外加电势下给出，如图5（B）所示，得到的I-V曲线如图5（C）所示.体相浓度1 mmol/L下不同表面电荷密度的PNP模型计算得到的±2，±1，±0.5 V 下的离子电流，以及对应SPICE 模型求解使用的参数和计算得到的±2，±1，±0.5 V下的电流列于表2.可见，SPICE模型同样对这一条件下的PNP模型有较好的拟合结果.

Table 2 Parameters for SPICE model and simulated results of PNP and SPICE models of the 1 mmol/L bulk concentration conical ion channel

对高浓度100 mmol/L KCl溶液体系的PNP模型进行模拟计算.离子通道几何模型不变，设定其表面电荷密度为-10，-5，-2，-1 mC/m2时的I-V曲线如图5（D）所示，相应±2 V 下整流比为2.56，1.67，1.23，1.11，此时，锥孔端10 nm的直径相对100 mmol/L溶液体相浓度而言较大，因此，这一离子通道体系难以达到较高的整流比.而模拟结果同样表明本组模型电迁移作用贡献的离子电流占据主导地位（图S3，见本文支持信息）.

用相同方法应用SPICE 模型模拟I-V曲线，图5（E）给出PNP 模型模拟的0 V 外加电势下的内建电势，得到的I-V曲线如图5（F）所示.

体相浓度100 mmol/L 下不同表面电荷密度的PNP 模型计算得到的±2，±1，±0.5 V下的离子电流，以及对应SPICE模型求解使用的参数和计算得到的±2，±1，±0.5 V下的电流列于表3.可见，在低整流比下SPICE模型同样对PNP模型有较好的拟合结果.

3 结论

考察了常用离子通道体系I-V曲线偏离理想PN 结二极管的特性，通过PNP 模型模拟典型的锥形离子通道实验装置，发现离子电流中电迁移电流占据主导地位.通过SPICE 模型将该离子通道作为电路中的一个器件进行研究时，比较固体PN结二极管器件的电流成因，提出其在反向偏压下仍然可以由电解池引入离子电流，因此，使用包含反向饱和电流项的二极管I-V曲线模型进行拟合.通过对KCl体相浓度为1，10和100 mmol/L且具有不同表面电荷密度的锥形离子通道体系进行I-V特性曲线数值模拟可知，SPICE模型的模拟结果与PNP模型结果比较，特别是在反向偏压下离子通道对离子输运阻力占主要作用时的I-V曲线有很好的拟合.本文初步给出了一定条件范围内的锥形离子通道I-V特性曲线的SPICE模型，后续的研究可以引入本文省略的参数.如比较理想PN结二极管的电流成因，引入饱和电流项可以讨论化学平衡存在明显作用的体系，而反向击穿电压和正向注入大电流饱和的引入可以更精细地刻画不同种类离子通道体系的I-V曲线，此外，电路中的电容项可用于模拟I-t曲线.使用SPICE通用模型描述可以进一步探索纳米流体二极管与其它元件相连接成电路时在不同条件下的工作状态.

支持信息见http://www.cjcu.jlu.edu.cn/CN/10.7503/cjcu20230384.