下载:
-
可调谐半导体激光吸收光谱(TDLAS)技术具有高灵敏度、非入侵、高测量精度和工业普适性等优点, 可实现对燃烧场和推进流场温度、组分浓度的实时测量[1,2]. 为了准确获得非均匀分布流场温度和浓度的空间分布, 将TDLAS技术与计算机层析成像(CT)技术相结合, 即激光吸收光谱层析成像(LAST)技术[3-9]. LAST对于研究燃烧过程和火焰结构等诊断工具具有十分重要的意义, 是流场测量的理想方法之一. LAST基本原理是将感兴趣区域(ROI)离散为一定数量的网格, 假设其物理参数不变, 然后通过测量到的投影数据, 建立参数方程求解逆问题, 从而获得网格中的参数[10-12]. 在经典的吸收层析成像测量中, 积分吸光度(投影数据)是长度、浓度和压力的线性函数, 这满足重建算法的条件. 因此, 准确的积分吸光度测量结果是LAST重建的一个重要因素[13].
LAST重建通常采用直接吸收光谱法(DAS)和波长调制光谱(WMS)两种方案获得积分吸光度[14]. DAS获得积分吸光度的一个关键步骤是确定基线, 从透射强度中计算出入射强度. 然而, 在恶劣的环境中, 入射激光强度受振动和窗口污损而发生变化, 这在计算积分吸光度时引入了误差, 进而这些误差被引入层析成像的重建中. 即使对重建算法和模型进行优化, 也很难消除积分吸光度中的误差. 相反WMS可以达到非常高的灵敏度, 对背景噪声有很强的鲁棒性, 被广泛应用在恶劣的环境中. Rieker等[9,15]开发了一种免标定的WMS方法, 并采用一次谐波对二次谐波进行归一化(WMS-2f/1f )的方式消除光强波动的影响. Sun等[16]在此基础上提出了一种基于WMS-2f/1f线型拟合的测量方法(WMS-Fit), 极大地拓展了免标定波长调制光谱方法的适用范围, 不再受光学深度、调制指数和激光波长的限制. 然而WMS-Fit方法需要大量的迭代计算, 每一步迭代都需要对信号滤波和拟合处理, 计算效率较低. 针对复杂燃烧场的重建过程, 过长的计算积分吸光度时间, 无法保证重建场的实时性, 只能通过离线采集, 后期再完成计算, 这无法为复杂燃烧场提供一些实时参数. 此外, WMS-Fit算法对拟合的参数初始值很敏感, 不恰当的初始值可能会增加计算时间, 甚至使解落入局部最优. 在之前的工作中, 我们提出了一种基于非零偶次谐波的免标定波长调制测量算法[17], 只在低频实验信号、均匀温度状态下进行了视线测量, 然而, 并没有对非均匀复杂燃烧环境场重建进一步说明.
本文重点针对非均匀复杂燃烧环境场重建, 呈现了一种快速免标定波长调制光谱法获取积分吸光度的方法(R-WMS-IA), 并且进一步将该算法应用于非均匀复杂燃烧场温度和浓度重建工作. 所采用的方法只涉及简单的代数运算来计算积分吸光度, 从而消除了传统拟合算法中耗时的线型拟合过程. 通过仿真和实验室丁烷喷灯复杂燃烧火焰对该方法的重建质量和重建效率进行验证.
-
采用基于偶次谐波的波长调制方法实现积分吸光度的快速计算, 该方法的更多细节可见文献[17], 本节只做简要介绍. 当频率为ν (单位: cm–1)的一束激光在长度为L (单位: cm)的路径上穿透吸收的气体时, 其强度的一部分被吸收. 根据Beer-Lambert定律, 与频率有关的吸光度可定义为
$\begin{split}\;& \alpha (\nu) = - \ln \left[ {{{{I_{\text{t}}}(\nu)}/ {{I_0}(\nu)}}} \right] \\ = \;&\int_0^L {P\left( l \right)X\left( l \right){S_\nu}\left[ {T(l)} \right]\phi (\nu,l){\text{d}}l} = A\phi (\nu) ,\end{split} $ 式中I0(ν)和It(ν)分别是入射和透射的激光强度, l是沿着路径的位置, P(l) (单位: atm, 1 atm = 1.013×105 Pa)是总压力, X(l)表示气体摩尔分数, T(l ) (单位: K)是气体温度, S(T )(单位: cm–2·atm–1)是随温度变化的线强, φ(ν) (单位: cm–1)是线型函数, 为洛伦兹(Lorenz)和高斯(Gauss)线型函数的卷积形式, 采用Voigt线型函数描述. 其具体数学表达式为[18-20]
$\begin{split} \varphi (\nu) =\;& \frac{{{c_{\text{L}}}}}{{\text{π }}}\frac{{{\text{δ}}{{\nu}_{\text{V}}}}}{{{{\left( {{\nu} - {{\nu}_0}} \right)}^2} + {\text{δ}}{\nu}_{\text{V}}^{\text{2}}}} + \frac{{{c_{\text{G}}}}}{{{\text{δ}}{{\nu}_{\text{V}}}}}{\left( {\frac{{\ln 2}}{{\text{π }}}} \right)^{1/2}} \\ &\times \exp \left[ {\frac{{ - \ln 2 \cdot {{\left( {{\nu} - {{\nu}_0}} \right)}^2}}}{{{\text{δ}}{\nu}_{\text{V}}^{\text{2}}}}} \right] \text{, } \end{split}$ 式中, ν0为谱线中心频率,
${\text{δ}}{{\nu}_{\text{V}}}$ 为Voigt半高全宽,${\text{δ}}{{\nu}_{\text{L}}}$ 和${\text{δ}}{{\nu}_{\text{G}}}$ 分别为Lorenz和Gauss的半高全宽, cL和cG分别对应其权重系数. 相关的表达式为[21]$ {\text{δ}}{{\nu}_{\text{V}}} = 0.5346{\text{δ}}{{\nu}_{\text{L}}} + \sqrt {0.2166{\text{δ}} {\nu}_{\text{L}}^{\text{2}} + {\text{δ}}{\nu}_{\text{G}}^2} , $ $\begin{split} {c_{\text{L}}} =\;& 0.6818817 + 0.6129331d \\ &- 0.1838439{d^2} - 0.1156844{d^3} ,\end{split} $ $ \begin{split}{c_{\text{G}}} = \;&0.3246017 - 0.6182531d \\ & + 0.1768139{d^2} + 0.1210944{d^3} \text{, } \end{split}$ 其中, d是线型参数,
$ d = ({\text{δ}}{{\nu}_{\text{L}}} - {\text{δ}}{{\nu}_{\text{G}}})/({\text{δ}}{{\nu}_{\text{L}}} + {\text{δ}}{{\nu}_{\text{G}}}) . $ 根据文献[17]理论, 光谱吸光度的谐波中心峰高为
$ \begin{split}&{h}_{n} = \\ & \left\{ \begin{array}{ll} \dfrac{4A}{\text{πδ}{{\nu}}_{\text{V}}{c}_{\text{L}}{p}_{n}(m)} + \dfrac{4A}{\text{πδ}{{\nu}}_{\text{V}}{c}_{\text{G}}{q}_{n}(m)}, & n为偶数,\\ 0, & n为奇数,\end{array} \right.\end{split}$ 式中m为调制系数, In/2(z) = in/2Jn/2(iz), Jn/2为第一类n/2阶贝塞尔函数, pn(m)和 qn(m)分别为Lorenz和Gauss对应谐波中心峰值的权重系数, 具体公式为
$ {p_n}(m) = \frac{{{\varepsilon _n}}}{{2{m^n}}}\frac{{{{(\sqrt {1 + {m^2}} - 1)}^n}}}{{\sqrt {1 + {m^2}} }} \text{, } $ $ {q_n}(m) = \sqrt {{\text{π}}\ln 2} \exp (-{m^2}\ln 2/2){I_{{n/2}}}({m^2}\ln 2/2). $ (7)式可以写成矩阵形式:
$ {\boldsymbol{H}} = {k_p}{\boldsymbol{P}}(m) + {k_q}{\boldsymbol{Q}}(m) \text{, } $ 其中,
${\boldsymbol{H}}\in {\mathbb{R}} ^N$ 是各非零偶次谐波中心峰高的向量值, 即H = (h2, h4, ···, h2N)T; P(m)$\in {\mathbb{R}} ^N $ 和 Q(m)$\in {\mathbb{R}} ^N $ 是与调制系数m相关的基向量, 分别定义为P(m) = [p2(m), p4(m), ···, p2N(m)]T ,
Q(m) = [q2(m), q4(m), ···, q2N(m)]T.
待求的向量系数kp和kq定义为kp = 4AcL/(πλ) 和 kq = 4AcG/(πλ). H是矩阵M = [P(m), Q(m)]列空间中的一个元素. 因此向量H与其在矩阵M列空间中的投影向量H'严格相等. 向量H与投影向量H'满足等式:
$ {\boldsymbol{H}} - {{\boldsymbol{H}}'} = \big[ {{\boldsymbol{E}} - {\boldsymbol{M}}{{({{\boldsymbol{M}}^{\text{T}}}{\boldsymbol{M}})}^{ - 1}}{{\boldsymbol{M}}^{\text{T}}}} \big]{\boldsymbol{H}}={\boldsymbol{0}} \text{, } $ 其中E是N维单位矩阵. 调制系数m*可通过求解非线性最优化问题来求解, 具体为
$ {m^*} = \mathop {{\text{arg min}}}\limits_m ||({\boldsymbol{E}} - {\boldsymbol{M}}{({{\boldsymbol{M}}^{\text{T}}}{\boldsymbol{M}})^{ - 1}}{{\boldsymbol{M}}^{\text{T}}}){\boldsymbol{H}}|{|_2} . $ 根据(10)式, kp/kq与线型参数d的关系可以改成:
$ \frac{{{k_p}}}{{{k_q}}} = \frac{{{c_{\text{L}}}}}{{{c_{\text{G}}}}} =\frac{{0.6818817 + 0.6129331d - 0.1838439{d^2} - 0.1156844{d^3}}}{{0.3246017 - 0.6182531d+ 0.1768139{d^2}+ 0.1210944{d^3}}} . $ 根据获得的调制系数m*, 利用最小二乘法的原理求解(10)式中的系数kp和kq :
$ \left[ {\begin{array}{*{20}{c}} {{k_p}} \\ {{k_q}} \end{array}} \right]={({{\boldsymbol{M}}^{*{\text{T}}}}{{\boldsymbol{M}}^*})^{ - 1}}{{\boldsymbol{M}}^{{\text{*T}}}}{\boldsymbol{H}} \text{, } $ 式中矩阵M*通过计算得到调制系数m*代入M = [P(m), Q(m)]获得. 在获得调制系数m*后,
$ {\text{δ}}{{\nu}_{\text{V}}} $ 为Voigt半高全宽, 可通过下式得到:$ {\text{δ}}{{\nu}_{\text{V}}} = 2a/{m^*} , $ 其中a为调制深度. 根据最小二乘法的原理计算得到的kp和kq, 积分吸光度可通过(16)式获得:
$ A=\text{πδ}{{\nu}}_{\text{V}}{k}_{p}/(4{c}_{\text{L}})\text{ }或\text{ }A=\text{πδ}{{\nu}}_{\text{V}}{k}_{q}/(4{c}_{\text{G}}) . $ -
由于线性函数可以归一化为
$\displaystyle \int\nolimits_{ - \infty }^{ + \infty } {\phi \left( {\nu} \right){\text{d}}\nu} \equiv 1$ , 则路径积分吸光度A为$ \begin{split}A =\;& \int_{ - \infty }^\infty {\alpha \left( {\nu} \right){\text{d}}\nu} = \int_0^L {P\left( l \right)X\left( l \right)S\left[ {T\left( l \right)} \right]{\text{d}}l} \\ =\;& \int_0^L {u\left( l \right){\text{d}}l} ,\end{split} $ 式中u为路径上局部吸光度. 对于吸收光谱层析成像过程中, (17)式可以描述为
$ {\boldsymbol{Lu}} = {\boldsymbol{A}}, $ 式中A
$\in {\mathbb{R}} $ M×1表示从M个视线测量得到的路径积分吸光度向量. L$\in {\mathbb{R}} $ M×N是灵敏度矩阵, 其元素li, j代表第i束激光通过第j个像素的激光路径段的长度. i (i = 1, 2, ···, M)和j (j = 1, 2, ···, N)分别是激光束和像素的索引. u$\in {\mathbb{R}} $ N×1是A的局部吸光度向量, 其元素uj = Pj Xj S(Tj).当通过测量获得目标气体在两个或多个波数的局部吸光度u, 从中可以提取温度和气体浓度的火焰分布参数. 第j个网格的温度可以从两条预选的吸收谱线(ν1和ν2)下重建的局部吸光度uν1, j和uν2
, j的比率计算出来, 具体计算温度公式为 $ \begin{split}{T_j} =\;& \frac{{hc}}{k}\left( {{E''_2} - {E''_1}} \right) \bigg[ \ln\frac{{{u_{{\nu}1,j}}}}{{{u_{{\nu}2,j}}}} \\ &\times \ln\frac{{S{}_{{\nu}2}\left( {{T_0}} \right)}}{{S{}_{{\nu}1}\left( {{T_0}} \right)}} + \frac{{hc}}{k}\left( {\frac{{{E''_{{\nu}2}} - {E''_{{\nu}1}}}}{{{T_0}}}} \right) \bigg]^{-1}, \end{split} $ 式中, ν1和ν2代表两条吸收谱线频率. 当温度值获得后, 气体浓度Xj表达式为
$ {X_j} = \frac{{{u_j}}}{{{P_j}{l_j}S\left( {{T_j}} \right)}}\text{, } $ 式中, 局部积分吸光度u可以通过迭代算法求解逆问题来获得. 由于代数迭代重建算法(ART)收敛速度较慢, 且松弛因子μ一般为常数, 不能反映当前网格信息, 本文采用修正ART(MAART)算法[13,22], 将吸收强度引入到迭代关系中, 对松弛因子μ进行了修正, 具体为
$ {u^{(k + 1,i)}} = {u^{(k,i)}} - \mu \big[ {{L_i}{u^{(k,i)}} - {A_i}} \;\big] \text{, } $ $ \mu = \beta \times \frac{{u_j^{(k)}{L_{ij}}}}{{\displaystyle\sum\nolimits_{j = 1}^N {u_j^{(k)}{L_{ij}}} }} \text{, } $ 式中k表示当前迭代次数, Ai表示第i条激光束获得的积分吸光度. 松弛因子μ用来控制收敛率, 通常范围为(0, 2). 该算法(22)式中的β通过最小二乘法寻优获得, 最优值为0.01. 图1为基于激光吸收光谱快速工业级层析成像的整个算法流程. 针对复杂燃烧场, 首先将燃烧场离散化, 通过光路切换获得每一路激光束的入射光强, 在经过燃烧场则获得每一路的透射光强, 通过FFT光谱分析法获得每一路的非零偶次谐波中心峰值. 通过偶次谐波算法计算出积分吸光度, 再通过求解逆问题获得局部吸光度, 从而反演出燃烧场的温度和浓度分布.
图 1 在线成像算法的基本框架
Figure 1. Basic framework of on-line imaging algorithms.
-
水蒸气是碳氢化合物燃烧的主要产物, 有很强的近红外吸收光谱, 因此选择水蒸气作为目标吸收气体. 本文选择7185.60 cm–1和7444.36 cm–1两条谱线, 通过双线法来推断气体温度. 根据HITRAN数据库[23], 谱线对参数见表1, 列出了T0 = 296 K时的线强度(S(T0))、自展宽(ξself)和空气(ξair)的压力增宽引起的半宽、低态的能量(E'')和温度相关增宽系数(nair).
Line index wavenumber/cm–1 S(T0)/(cm–2·atm–1) ξself/(cm–1·atm–1) ξair/(cm–1·atm–1) E''/cm–1 nair 1 7185.596 0.00490 0.371 0.0342 1045.0583 0.63 7185.597 0.0147 0.195 0.0413 1045.058 0.65 2 7444.351 0.000541 0.366 0.0199 1774.750 0.44 7444.368 0.000154 0.250 0.0188 1806.670 0.41 7444.371 0.000462 0.194 0.0153 1806.669 0.41 表 1 7185.60 cm–1和7444.36 cm–1中心谱线处的参数
Table 1. Parameters of the selected transitions at around 7185.60 cm–1 and 7444.36 cm–1.
3个非均匀温度和H2O浓度分布区域如图2所示, 其中3个非均匀区域分别是由1个、2个和3个高斯曲线拟合所得. 在碳氢化合物火焰中, H2O浓度通常与温度有很好的相关性. 因此, 每个分布区域中的H2O浓度分布与温度分布相似. 为了验证所提方法的有效性, 本文采用三维高斯分布来模拟非均匀的复杂燃烧场分布. 重建范围为200 mm × 200 mm, 网格划分为120×120, 网格尺寸为1.67 mm × 1.67 mm. 具体的场分布可以用(23)式和(24)式来计算:
图 2 非均匀温度和H2O浓度分布区域 (a), (b) 分别表示模型1温度和浓度分布; (c), (d) 分别表示模型2温度和浓度分布; (e), (f) 分别表示模型3温度和浓度分布
Figure 2. Non-uniform temperature and H2O concentration distribution regions: (a), (b) Indicate the temperature, concentration distribution of model 1; (c), (d) indicate the temperature, concentration distribution of model 2; (e), (f) indicate the temperature, concentration distribution of model 3, respectively.
$\begin{split} T =\;& 400 + \sum\limits_{k = 1}^K 1200 \cdot {\eta ^k} \\ &\times\exp \left[ { - \frac{{{{\left( {x - x_{\text{c}}^k} \right)}^2} + {{\left( {y - y_{\text{c}}^k} \right)}^2}}}{{{\sigma ^2}}}} \right] \text{, }\end{split} $ $ \begin{split}X =\;& 0.015 + \sum\limits_{k = 1}^K 0.04 \cdot {\eta ^k} \\ &\times\exp \left[ { - \frac{{{{\left( {x - x_{\text{c}}^k} \right)}^2} + {{\left( {y - y_{\text{c}}^k} \right)}^2}}}{{{\sigma ^2}}}} \right] \text{, } \end{split}$ 其中x和y分别表示重建区域的横坐标和纵坐标. (
$ x_{\text{c}}^k $ ,$ y_{\text{c}}^k $ ) 是第k个高斯峰的中心, η是高斯分布峰值参数, σ是标准差. 3个分布区域的详细参数见表2.分布模型 η ($ x_{\text{c}}^k $, $ y_{\text{c}}^k $)/cm σ/cm 1 0.40 (12, 8) 4 2 0.40
0.35(6, 14)
(14, 6)4 3 0.40
0.35
0.25(6, 14)
(14, 14)
(9, 6)4 表 2 三种不同分布模型的详细参数
Table 2. Detailed parameters of the three different distribution models.
整套测量系统的光路布置如图3所示. 结合整套测量系统传感器的体积和布置位置, 同时满足激光束尽可能多角度布置[24,25]的原则, 文中100路激光束, 分别布置在待测区域四边, 单边25路, 交叉入射到对边2个积分球中, 从而实现了最优光路布置方案.
图 3 (a)测量系统中光路布置的位置关系; (b)测量区域中光路布置
Figure 3. (a) Position relationship of the optical path arrangement in the measurement system; (b) optical paths arrangement in the measurement area.
为了验证该方法的有效性, 用于评估重建质量的表达式为
$ e = {{{{\left\| {{{\vec u}^{\,\text{p}}} - {{\vec u}^{\,\text{r}}}} \right\|}_2}} / {{{\left\| {{{\vec u}^{\,\text{p}}}} \right\|}_2}}} \text{, } $ 式中,
${\vec u^{\,\text{p}}}$ 和${\vec u^{\,\text{r}}}$ 分别表示设定值和重建值, ‖·‖2是向量的二范数. -
仿真中, 将R-WMS-IA方法和文献[16]所提算法WMS-Fit与理论值进行了对比分析, 结果如图4所示, 三种不同分布模型在谱线7185.60 cm–1计算获得的积分吸光度. 图4各图上部分为积分吸光度计算结果, 下部分为两种方法与理论值的相对偏差, 其中激光束编号与图3相对应. 由图4不难看出, R-WMS-IA方法获得积分吸光度相对偏差比WMS-Fit方案大, 但均控制在3%以内. 该结果表明R-WMS-IA方法在计算非均匀燃烧场的可行性. 谱线7444.36 cm–1计算结果与谱线7185.60 cm–1相似.
图 4 三种不同分布模型通过R-WMS-IA和WMS-Fit算法获得积分吸光度的计算结果及相对误差 (a) 分布模型1; (b) 分布模型2; (c) 分布模型3
Figure 4. Calculation results and relative error of integrated absorbance by R-WMS-IA and WMS-Fit algorithms for three different distribution models: (a) Distribution model 1; (b) distribution model 2; (c) distribution model 3.
三种不同分布模型重建结果分别见图5—图7. 在三种不同的分布模型中, R-WMS-IA和WMS-Fit算法获得的积分吸光度用于温度、浓度重建结果接近于理论仿真结果. 分布模型1重建偏差范围在5.00%—5.58%之间(图5), 分布模型2重建偏差范围在5.74%—6.36%之间(图6), 分布模型3重建偏差范围在4.40%—4.87%之间(图7). 这些偏差主要由求解逆问题的重建算法引起, 由于积分吸光度偏差导致重建偏差很小. 温度是由两条谱线的积分吸光度的比率获得. 然而, 在浓度重建中不能消除积分吸光度偏差, 因为浓度测量依赖于温度重建和一个过渡的积分吸光度. 因此, 浓度重建误差要比温度重建大.
图 5 分布模型1中不同方法的重建结果 (a) e = 5.34%; (b) e = 5.38%; (c) e = 5.53%; (d) e = 5.42%; (e) e = 5.58%; (f) e = 5.00%
Figure 5. Reconstruction results of different methods in distribution model 1: (a) e = 5.34%; (b) e = 5.38%; (c) e = 5.53%; (d) e = 5.42%; (e) e = 5.58%; (f) e = 5.00%.
图 6 分布模型2中不同方法的重建结果 (a) e = 5.74%; (b) e = 6.20%; (c) e = 5.91%; (d) e = 5.91%; (e) e = 6.33%; (f) e = 6.36%
Figure 6. Reconstruction results of different methods in distribution model 2: (a) e = 5.74%; (b) e = 6.20%; (c) e = 5.91%; (d) e = 5.91%; (e) e = 6.33%; (f) e = 6.36%.
图 7 分布模型3中不同方法的重建结果 (a) e = 4.46%; (b) e = 4.70%; (c) e = 4.61%; (d) e = 4.74%; (e) e = 4.87%; (f) e = 4.40%
Figure 7. Reconstruction results of different methods in distribution model 3: (a) e = 4.46%; (b) e = 4.70%; (c) e = 4.61%; (d) e = 4.74%; (e) e = 4.87%; (f) e = 4.40%.
-
3.2节分析了所提算法R-WMS-IA获得积分吸光度的能力及成像可行性. 本节分析所提算法获得积分吸光度的时间, 如图8所示. R-WMS-IA和WMS-Fit计算积分吸光度的方法在重建误差方面的偏差相对较小, 但计算积分吸光度所用时间的差异特别大. 所涉及的计算机CPU为Intel Core i7-9700 @ 3.00 GHz. 对于100路激光束的非均匀场, R-WMS-IA方法的平均计算时间为0.14 s. 然而, WMS-Fit方法的平均计算时间为19.80 s. 换言之, 所提算法的计算效率至少提高了2个数量级.
图 8 仿真分析通过两种方法计算积分吸光度时间
Figure 8. Simulation analysis calculates the integral absorbance time by two methods.
-
实验测试系统如图9所示, 整个系统主要包括激光源、丁烷喷灯(BRS-75)、K-型热电偶(HH-K-24-SLE)、光开关(定制)、信号传感器和主机. 重建范围为200 mm×200 mm, 网格划分为120×120, 网格尺寸为1.67 mm×1.67 mm.
图 9 实验测试系统
Figure 9. Experimental test system.
数据采集设备(NI-DAQ, PCI6115)同时具备信号发生功能, 产生1 kHz正弦扫描电压信号与100 kHz正弦调制的叠加电压信号. 两个分布式反馈(DFB)二极管的激光器7185.60 cm–1(NTT, NLK1E5GAAA)和7444.36 cm–1(NTT, NLK1B5EAAA)在H2O分子吸收的中心谱线附近工作, 以时分复用(TDM)模式运行, 通过控制器 (LDC501, SRS)对其出光波长进行调制. 两束激光合并后由一个2×2的单模光纤耦合器分开, 一条接入一个自由光谱范围(FSR)为0.01 cm–1的标准具(PMMZI-1400—300 MHz), 用于波长监测; 另一条则接入光开关, 通过光开关切换光路, 将每一路输出激光束入射到对应的积分球. 光开关可以实现激光束1路输入100路输出. 通过内部切换器将输入激光束切换到不同的输出激光束, 切换时间为10 ms. 将100路激光束中的每25路安装在待测区的一侧, 交叉入射到对边2个积分球中, 实现多角度光路布置.
为了实现激光束多自由度的调节, 设计了如图9所示的小型化激光发射模块. 整个传感器安装有100条激光束, 完成一幅图像的扫描需要1 s. 每束激光的信号接收器(IS-Detector)主要由一个积分球(Thorlabs, IS200)和一个光电探测器(Thorlabs, PDA20CS-EC)组成. 激光束入射到直径为12 mm的积分球, 这意味着激光束的接收面被扩大为一个直径为12 mm的圆面积, 同时进入积分球的激光被多次漫反射, 因此在端口处实现了光强的均化, 使得探测器接收的信号更稳定. 为了检测各个角度的激光束, 传感器配备了8套信号接收器, 分别布置在待测试区域的四边, 具体布置如图9所示. DAQ可实现同步发生和采集信号. 吸收信号采集的数字化采样率为107 sample/s. 所提信号接收系统不仅可以大大降低信号接收和采集的成本, 而且还可以提高空间分辨率.
-
为了进一步验证该方法的有效性, 通过丁烷喷灯燃烧实验重建燃烧场的温度分布图像[26,27]. 丁烷喷灯燃烧器实验平台3D建模和局部视图如图10(a), (b)所示, 实验燃烧火焰和局部视图如图10(c), (d)所示. 在与吸收光谱测量相同截面均匀布置9个K-型热电偶, 布置间距为20 mm, 以便与所提方法重建结果进行温度比较, 从而可以对测量精度进行评估.
图 10 丁烷喷灯燃烧器的实验平台 (a)实验平台3D建模; (b)建模局部视图; (c)实验照片; (d)局部火焰视图
Figure 10. Experimental platform of butane burner: (a) 3D modeling structure of experimental platform; (b) partial view of modeling; (c) experimental photo; (d) partial view of flame.
当燃烧火焰稳定后, 通过火焰的激光束被光电探测器接收, 其电压信号用来计算积分吸光度进而重建温度场. 以测量系统中第50路激光束为例, 探测器收到透射光强的电压信号和通过数字锁相技术提取的谐波信号, 如图11所示. 以7185.60 cm–1和7444.36 cm–1为中心的光谱范围的投射光强曲线分别为绿色和紫色. 由于入射激光的波长被100 kHz的正弦信号调制, 所以以7185.60 cm–1和7444.36 cm–1为中心的吸收光谱解调的偶次谐波受到的干扰很小. 通过获得的偶次谐波峰值高度进一步计算积分吸光度.
图 11 第50路透射光强的电压信号和以7185.60 cm–1(绿色)、7444.36 cm–1(紫色)光谱吸收为中心的2、4、6次谐波信号
Figure 11. Voltage signal of the 50th transmitted light intensity and 2nd, 4th and 6th harmonic signals centered on 7185.60 cm–1 (green) and 7444.36 cm–1 (violet) spectral absorption.
激光经过丁烷喷灯燃烧场, 通过R-WMS-IA和WMS-Fit算法获得的积分吸光度计算结果, 如图12所示. 由于实验过程中受到噪声的影响, 两种算法计算的积分吸光度较仿真结果偏差增大,但整体偏差控制在5.5%以内.
图 12 通过R-WMS-IA和WMS-Fit算法获得积分吸光度的计算结果 (a) 谱线中心为7185.60 cm–1积分吸光度; (b) 谱线中心为7444.36 cm–1积分吸光度
Figure 12. Results of the integrated absorbance calculations were obtained by R-WMS-IA and WMS-Fit algorithms: (a) Spectral line centered at 7185.60 cm–1 integrated absorbance; (b) spectral line centered at 7444.36 cm–1 integrated absorbance.
通过两条谱线积分吸光度求解逆问题, 获得温度场重建结果, 具体如图13所示. 热电偶与光电探测器安装在同一平面上, 方便与所提算法重建结果进行温度比较, 从而可以对测量精度进行评估. 两种算法重建温度值与热电偶对应截面测量结果如图14所示, 两种算法的温度重建结果基本相同, 最大偏差仅为0.94%, 与热电偶测量值相比最大偏差为3.5%, 说明了所提算法的有效性.
图 13 分别通过两种算法获得积分吸光度的温度场重建结果 (a) R-WMS-IA方法; (b) WMS-Fit算法
Figure 13. The temperature field reconstruction results of the integrated absorbance were obtained by two algorithms, respectively: (a) R-WMS-IA algorithm; (b) WMS-Fit algorithm.
图 14 热电偶与两种算法的温度测量值
Figure 14. Temperature curves measured by the thermocouple and the two algorithms respectively.
-
与仿真分析重建时间一致, 采用相同的计算平台, 本节通过实验分析所提算法获得积分吸光度的时间, 如图15所示. 对于100路激光束的非均匀场, R-WMS-IA方法的平均计算时间为0.153 s. 然而, WMS-Fit方法的平均计算时间为21.1 s. R-WMS-IA方法的计算效率与仿真分析结果一致, 与传统WMS-Fit方法相比至少提高了2个数量级.
图 15 实验分析通过两种方法计算积分吸光度时间
Figure 15. Experiment analysis calculates the integral absorbance time by two methods.
-
本文采用了一种快速免标定波长调制光谱法获取积分吸光度的方法, 并且将该算法应用于非均匀复杂燃烧场. 通过数值模拟和丁烷喷灯燃烧火焰的实验对所提算法的有效性进行了验证.
仿真结果表明, 三种不同分布模型中R-WMS-IA方法获得积分吸光度与理论值相比的相对偏差高于WMS-Fit方法, 但均控制在3%以内. 通过两种方法获得的积分吸光度进行温度、浓度重建. 分布模型1重建偏差范围在5.00%—5.58%之间, 分布模型2重建偏差范围在5.74%—6.36%之间, 分布模型3重建偏差范围在4.40%—4.87%之间. 这些偏差主要是由求解逆问题的重建算法引起, 由于积分吸光度偏差导致重建偏差很小. 这表明R-WMS-IA方法在计算非均匀燃烧场的可行性. 但在计算100路激光束积分吸光度所用的时间时, R-WMS-IA方法的平均计算时间为0.14 s; 而WMS-Fit方法的平均计算时间为19.8 s. 换言之, 所提算法的计算效率至少提高了2个数量级.
实验结果表明, 与传统的WMS-Fit方法相比, 两种方法的重建分布基本一致, 最大测量偏差仅为0.94%, 与热电偶测量值相比最大偏差为3.5%, 同样表明了该方法的有效性. 在重建精度相当的前提下, R-WMS-IA方法消除了耗时的拟合程序, 大大减少了计算积分吸光度的时间, 不需要先采集好数据, 再离线计算积分吸光度. 与仿真结果相同, 计算效率至少提高了2个数量级, 为快速实现工业级重建燃烧场的温度浓度分布提供了有效的方法.
-
传统的免标定波长调制光谱方法一般需要结合光谱数据库和激光调制参数进行复杂的吸收光谱模拟, 对先验光谱参数的准确度和硬件参数提出了很高的要求, 同时不合适的初值会增加计算时间, 甚至会导致陷入局部最优解. 为提高计算效率, 本文引入一种快速免标定波长调制光谱技术获取积分吸光度. 该方法对光谱数据库的依赖性低, 计算效率高, 同时解决了传统方法在高温高压下由于吸收谱线展宽变大而导致的谐波信号不完整问题. 进一步将该方法应用于非均匀复杂燃烧场层析成像, 并结合所提成像系统实现了快速在线重建温度、浓度分布. 通过数值模拟和丁烷喷灯燃烧火焰的实验验证该方法获得积分吸光度的准确性和计算效率. 结果表明, 与传统的波长调制方法相比重建分布基本一致, 最大测量相对偏差仅为0.94%, 与热电偶测量值相比最大相对偏差为3.5%, 验证了该方法的准确性. 在重建精度相当的前提下, 分析两种方法获得积分吸光度的计算效率. 所引方法和传统方法平均每路计算时间分别为0.15 s和21.10 s. 所引方法的计算效率比传统方法至少提高了2个数量级, 为实现在线重建燃烧场的温度、浓度分布提供了快速可靠的研究方法和技术手段.Conventional calibration-free wavelength modulation spectroscopy generally requires complex absorption spectrum simulations in combination with spectral databases and laser modulation parameters, placing high demands on the accuracy of a priori spectral parameters and hardware parameters. Meanwhile, inappropriate initial values can increase the computation time and even lead to local optimal solutions. In order to improve the computational efficiency, a rapid calibration-free wavelength modulation spectroscopy to obtain the integrated absorbance is presented in this work. First, this method is computationally efficient, requiring only algebraic calculations by using the 2nd, 4th, and 6th harmonic center peak height parameters to obtain the integrated absorbance, eliminating the need for computationally intensive harmonic fitting calculations. Secondly, this method has low dependence on the spectral database, requiring only line intensity and low-state energy level spectral parameters. Finally, this method is highly adaptable and does not require scanning the complete absorption spectral line shape, which solves the problem of incomplete harmonic signals caused by the conventional method at high temperature and high pressure due to the broadening of the absorption spectral line. This method has previously been used only for line-of-sight measurements at low-frequency experimental signals in stable environments, and for calculating the integrated absorbance at average temperature, concentration and pressure states. In this work, the method is applied to non-uniform complex combustion field tomography and combined with the proposed tomographic system to achieve online reconstructing temperature and concentration distributions. The accuracy and computational efficiency of the method in obtaining the integrated absorbance are verified by numerical simulations and experiments on the butane burner flame. The results show that the presented method is consistent with the reconstructed distribution compared with the conventional wavelength modulation method, with a maximum relative deviation of only 0.94% from the measurement and 3.5% from the thermocouple measurement, verifying the accuracy of the method. The computational efficiencies of the two methods for obtaining the integrated absorbance are analyzed. The average calculation time per path is 0.15 s for the present method and 21.10 s for the conventional method. The calculation efficiency of the present method is at least two orders of magnitude higher than that of the conventional method, which provides a fast and reliable research method and technical means to realize the industrial-grade online reconstruction of temperature and concentration distribution of combustion fields.
-
Keywords:
- online tomography /
- laser absorption spectroscopy /
- wavelength modulation spectroscopy /
- complex combustion field
[1] 黄安, 许振宇, 夏晖晖, 姚路, 阮俊, 胡佳屹, 臧益鹏, 阚瑞峰 2021 光谱学与光谱分析 41 1144
Google Scholar
Hang A, Xu Z Y, Xia H H, Yao L, Ruan J, Hu J Y, Zang Y P, Kan R F 2021 Spectrosc. Spect. Anal. 41 1144
Google Scholar
[2] Wang Y, Zhou B, Liu C 2021 IEEE Photonics Technol. Lett. 33 1487
Google Scholar
[3] Liu C, Xu L, Li F, Cao Z, Tsekenis S A, McCann H 2015 Appl. Phys. B 120 407
Google Scholar
[4] 宋俊玲, 洪延姬, 王广宇, 潘虎 2012 物理学报 61 124
Google Scholar
Song J L, Hong Y J, Wang G Y, Pan H 2012 Acta Phys. Sin. 61 124
Google Scholar
[5] Liu C, Xu L 2019 Appl. Spectrosc. Rev. 54 1
Google Scholar
[6] 臧益鹏, 许振宇, 黄安, 艾苏曼, 夏晖晖, 阚瑞峰 2021 物理学报 70 134205
Google Scholar
Zang Y P, Xu Z Y, Hang A, Ai S M, Xia H H, Kan R F 2021 Acta Phys. Sin. 70 134205
Google Scholar
[7] Liu C, Cao Z, Li F, Lin Y, Xu L 2017 Meas. Sci. Technol. 28 054002
Google Scholar
[8] 屈东胜, 洪延姬, 朱晓辉 2021 光谱学与光谱分析 41 1072
Qu D S, Hong Y J, Zhu X H 2021 Spectrosc. Spect. Anal. 41 1072
[9] Rieker G B, Jeffries J B, Hanson R K 2009 Appl. Opt. 48 5546
Google Scholar
[10] Grauer S J, Emmert J, Sanders S T, Wagner S, Daun K J 2019 Meas. Sci. Technol. 30 105401
Google Scholar
[11] Zhang R, Si J, Enemali G, Bao Y, Liu C 2022 IEEE Sens. J. 22 12728
Google Scholar
[12] Shui C, Wang Y, Cai W, Zhou B 2021 Opt. Express 29 20889
Google Scholar
[13] Song J, Xin M, Rao W, Hong Y, Feng G 2021 Appl. Opt. 60 5056
Google Scholar
[14] Peng D, Jin Y, Zhai C, Yang J 2018 Spectrosc. Lett. 51 7
Google Scholar
[15] Rieker G B, Li H, Liu X, Jeffries J B, Hanson R K, Allen M G, Wehe S D, Mulhall P A, Kindle H S 2007 Meas. Sci. Technol. 18 1195
Google Scholar
[16] Sun K, Chao X, Sur R, Goldenstein C S, Jeffries J B, Hanson R K 2013 Meas. Sci. Technol. 24 125203
Google Scholar
[17] Wang Y, Zhou B, Liu C 2021 Opt. Express 29 26618
Google Scholar
[18] 张书锋, 蓝丽娟, 丁艳军, 贾军伟, 彭志敏 2015 物理学报 64 053301
Google Scholar
Zhang S F, Lan L J, Ding Y J, Jia J W, Peng Z M 2015 Acta Phys. Sin. 64 053301
Google Scholar
[19] Wang Y, Zhou B, Wang B, Zhao R, Liu Q, Dai M 2022 Mathematics 10 308
Google Scholar
[20] Wang Y, Zhou B, Zhao R, Wang B, Liu Q, Dai M 2022 Mathematics 10 210
Google Scholar
[21] Liu Y, Lin J, Huang G, Guo Y, Duan C 2001 J. Opt. Soc. Am. B 18 666
Google Scholar
[22] Li N, Weng C 2011 Chin. Opt. Lett. 9 061201
Google Scholar
[23] Gordon I E, Rothman L S, Hill C, et al. 2017 J. Quant. Spectrosc. Radiat. Transf. 203 3
Google Scholar
[24] Terzija N, Davidson J L, Garcia-Stewart C A, Wright P, Ozanyan K B, Pegrum S, Litt T J, McCann H 2008 Meas. Sci. Technol. 19 094007
Google Scholar
[25] Grauer S J, Rice K M, Donbar J M, Bisek N J, France J J, Ochs B A, Steinberg A M 2022 AIAA J. 60 1
Google Scholar
[26] Ma L, Lau L Y, Ren W 2017 Appl. Phys. B 123 83
Google Scholar
[27] Huang A, Cao Z, Zhao W, Zhang H, Xu L 2020 IEEE Trans. Instrum. Meas. 69 9087
Google Scholar
-
图 2 非均匀温度和H2O浓度分布区域 (a), (b) 分别表示模型1温度和浓度分布; (c), (d) 分别表示模型2温度和浓度分布; (e), (f) 分别表示模型3温度和浓度分布
Fig. 2. Non-uniform temperature and H2O concentration distribution regions: (a), (b) Indicate the temperature, concentration distribution of model 1; (c), (d) indicate the temperature, concentration distribution of model 2; (e), (f) indicate the temperature, concentration distribution of model 3, respectively.
图 3 (a)测量系统中光路布置的位置关系; (b)测量区域中光路布置
Fig. 3. (a) Position relationship of the optical path arrangement in the measurement system; (b) optical paths arrangement in the measurement area.
图 4 三种不同分布模型通过R-WMS-IA和WMS-Fit算法获得积分吸光度的计算结果及相对误差 (a) 分布模型1; (b) 分布模型2; (c) 分布模型3
Fig. 4. Calculation results and relative error of integrated absorbance by R-WMS-IA and WMS-Fit algorithms for three different distribution models: (a) Distribution model 1; (b) distribution model 2; (c) distribution model 3.
图 5 分布模型1中不同方法的重建结果 (a) e = 5.34%; (b) e = 5.38%; (c) e = 5.53%; (d) e = 5.42%; (e) e = 5.58%; (f) e = 5.00%
Fig. 5. Reconstruction results of different methods in distribution model 1: (a) e = 5.34%; (b) e = 5.38%; (c) e = 5.53%; (d) e = 5.42%; (e) e = 5.58%; (f) e = 5.00%.
图 6 分布模型2中不同方法的重建结果 (a) e = 5.74%; (b) e = 6.20%; (c) e = 5.91%; (d) e = 5.91%; (e) e = 6.33%; (f) e = 6.36%
Fig. 6. Reconstruction results of different methods in distribution model 2: (a) e = 5.74%; (b) e = 6.20%; (c) e = 5.91%; (d) e = 5.91%; (e) e = 6.33%; (f) e = 6.36%.
图 7 分布模型3中不同方法的重建结果 (a) e = 4.46%; (b) e = 4.70%; (c) e = 4.61%; (d) e = 4.74%; (e) e = 4.87%; (f) e = 4.40%
Fig. 7. Reconstruction results of different methods in distribution model 3: (a) e = 4.46%; (b) e = 4.70%; (c) e = 4.61%; (d) e = 4.74%; (e) e = 4.87%; (f) e = 4.40%.
图 8 仿真分析通过两种方法计算积分吸光度时间
Fig. 8. Simulation analysis calculates the integral absorbance time by two methods.
图 10 丁烷喷灯燃烧器的实验平台 (a)实验平台3D建模; (b)建模局部视图; (c)实验照片; (d)局部火焰视图
Fig. 10. Experimental platform of butane burner: (a) 3D modeling structure of experimental platform; (b) partial view of modeling; (c) experimental photo; (d) partial view of flame.
图 11 第50路透射光强的电压信号和以7185.60 cm–1(绿色)、7444.36 cm–1(紫色)光谱吸收为中心的2、4、6次谐波信号
Fig. 11. Voltage signal of the 50th transmitted light intensity and 2nd, 4th and 6th harmonic signals centered on 7185.60 cm–1 (green) and 7444.36 cm–1 (violet) spectral absorption.
图 12 通过R-WMS-IA和WMS-Fit算法获得积分吸光度的计算结果 (a) 谱线中心为7185.60 cm–1积分吸光度; (b) 谱线中心为7444.36 cm–1积分吸光度
Fig. 12. Results of the integrated absorbance calculations were obtained by R-WMS-IA and WMS-Fit algorithms: (a) Spectral line centered at 7185.60 cm–1 integrated absorbance; (b) spectral line centered at 7444.36 cm–1 integrated absorbance.
图 13 分别通过两种算法获得积分吸光度的温度场重建结果 (a) R-WMS-IA方法; (b) WMS-Fit算法
Fig. 13. The temperature field reconstruction results of the integrated absorbance were obtained by two algorithms, respectively: (a) R-WMS-IA algorithm; (b) WMS-Fit algorithm.
图 14 热电偶与两种算法的温度测量值
Fig. 14. Temperature curves measured by the thermocouple and the two algorithms respectively.
图 15 实验分析通过两种方法计算积分吸光度时间
Fig. 15. Experiment analysis calculates the integral absorbance time by two methods.
表 1 7185.60 cm–1和7444.36 cm–1中心谱线处的参数
Table 1. Parameters of the selected transitions at around 7185.60 cm–1 and 7444.36 cm–1.
Line index wavenumber/cm–1 S(T0)/(cm–2·atm–1) ξself/(cm–1·atm–1) ξair/(cm–1·atm–1) E''/cm–1 nair 1 7185.596 0.00490 0.371 0.0342 1045.0583 0.63 7185.597 0.0147 0.195 0.0413 1045.058 0.65 2 7444.351 0.000541 0.366 0.0199 1774.750 0.44 7444.368 0.000154 0.250 0.0188 1806.670 0.41 7444.371 0.000462 0.194 0.0153 1806.669 0.41 
下载: 导出CSV
表 2 三种不同分布模型的详细参数
Table 2. Detailed parameters of the three different distribution models.
分布模型 η ($ x_{\text{c}}^k $, $ y_{\text{c}}^k $)/cm σ/cm 1 0.40 (12, 8) 4 2 0.40
0.35(6, 14)
(14, 6)4 3 0.40
0.35
0.25(6, 14)
(14, 14)
(9, 6)4
下载: 导出CSV
-
[1] 黄安, 许振宇, 夏晖晖, 姚路, 阮俊, 胡佳屹, 臧益鹏, 阚瑞峰 2021 光谱学与光谱分析 41 1144
Google Scholar
Hang A, Xu Z Y, Xia H H, Yao L, Ruan J, Hu J Y, Zang Y P, Kan R F 2021 Spectrosc. Spect. Anal. 41 1144
Google Scholar
[2] Wang Y, Zhou B, Liu C 2021 IEEE Photonics Technol. Lett. 33 1487
Google Scholar
[3] Liu C, Xu L, Li F, Cao Z, Tsekenis S A, McCann H 2015 Appl. Phys. B 120 407
Google Scholar
[4] 宋俊玲, 洪延姬, 王广宇, 潘虎 2012 物理学报 61 124
Google Scholar
Song J L, Hong Y J, Wang G Y, Pan H 2012 Acta Phys. Sin. 61 124
Google Scholar
[5] Liu C, Xu L 2019 Appl. Spectrosc. Rev. 54 1
Google Scholar
[6] 臧益鹏, 许振宇, 黄安, 艾苏曼, 夏晖晖, 阚瑞峰 2021 物理学报 70 134205
Google Scholar
Zang Y P, Xu Z Y, Hang A, Ai S M, Xia H H, Kan R F 2021 Acta Phys. Sin. 70 134205
Google Scholar
[7] Liu C, Cao Z, Li F, Lin Y, Xu L 2017 Meas. Sci. Technol. 28 054002
Google Scholar
[8] 屈东胜, 洪延姬, 朱晓辉 2021 光谱学与光谱分析 41 1072
Qu D S, Hong Y J, Zhu X H 2021 Spectrosc. Spect. Anal. 41 1072
[9] Rieker G B, Jeffries J B, Hanson R K 2009 Appl. Opt. 48 5546
Google Scholar
[10] Grauer S J, Emmert J, Sanders S T, Wagner S, Daun K J 2019 Meas. Sci. Technol. 30 105401
Google Scholar
[11] Zhang R, Si J, Enemali G, Bao Y, Liu C 2022 IEEE Sens. J. 22 12728
Google Scholar
[12] Shui C, Wang Y, Cai W, Zhou B 2021 Opt. Express 29 20889
Google Scholar
[13] Song J, Xin M, Rao W, Hong Y, Feng G 2021 Appl. Opt. 60 5056
Google Scholar
[14] Peng D, Jin Y, Zhai C, Yang J 2018 Spectrosc. Lett. 51 7
Google Scholar
[15] Rieker G B, Li H, Liu X, Jeffries J B, Hanson R K, Allen M G, Wehe S D, Mulhall P A, Kindle H S 2007 Meas. Sci. Technol. 18 1195
Google Scholar
[16] Sun K, Chao X, Sur R, Goldenstein C S, Jeffries J B, Hanson R K 2013 Meas. Sci. Technol. 24 125203
Google Scholar
[17] Wang Y, Zhou B, Liu C 2021 Opt. Express 29 26618
Google Scholar
[18] 张书锋, 蓝丽娟, 丁艳军, 贾军伟, 彭志敏 2015 物理学报 64 053301
Google Scholar
Zhang S F, Lan L J, Ding Y J, Jia J W, Peng Z M 2015 Acta Phys. Sin. 64 053301
Google Scholar
[19] Wang Y, Zhou B, Wang B, Zhao R, Liu Q, Dai M 2022 Mathematics 10 308
Google Scholar
[20] Wang Y, Zhou B, Zhao R, Wang B, Liu Q, Dai M 2022 Mathematics 10 210
Google Scholar
[21] Liu Y, Lin J, Huang G, Guo Y, Duan C 2001 J. Opt. Soc. Am. B 18 666
Google Scholar
[22] Li N, Weng C 2011 Chin. Opt. Lett. 9 061201
Google Scholar
[23] Gordon I E, Rothman L S, Hill C, et al. 2017 J. Quant. Spectrosc. Radiat. Transf. 203 3
Google Scholar
[24] Terzija N, Davidson J L, Garcia-Stewart C A, Wright P, Ozanyan K B, Pegrum S, Litt T J, McCann H 2008 Meas. Sci. Technol. 19 094007
Google Scholar
[25] Grauer S J, Rice K M, Donbar J M, Bisek N J, France J J, Ochs B A, Steinberg A M 2022 AIAA J. 60 1
Google Scholar
[26] Ma L, Lau L Y, Ren W 2017 Appl. Phys. B 123 83
Google Scholar
[27] Huang A, Cao Z, Zhao W, Zhang H, Xu L 2020 IEEE Trans. Instrum. Meas. 69 9087
Google Scholar
计量
- 文章访问数: 891
- PDF下载量: 46
- 被引次数: 0
