CN114692408A - 一种分布式光纤声传感系统的仿真建模方法 - Google Patents

Abstract

一种分布式光纤声传感系统的仿真建模方法，涉及分布式监测技术领域，包括：基于Simulink模块分别构建激光器模型、分束器模型、声光调制器模型、掺铒光纤放大器和衰减器模型、传感光纤模型、光电探测器模型、外界扰动模型和光纤耦合器模型后分别进行封装，根据分布式光纤声传感系统的光路图将对应的功能模块添加到Simulink工作区，利用信号线，将模块的输入信号和其他模块的输出信号名称相同的相连，完成分布式光纤声传感系统的Simulink建模进行仿真运算，为分布式光纤声传感系统的调试、安装、故障诊断以及系统优化提供直观依据，解决了没有分布式光纤声传感系统Simulink建模方法与系统的问题。

Description

一种分布式光纤声传感系统的仿真建模方法

技术领域

本发明涉及分布式监测技术领域，具体涉及一种分布式光纤声传感系统的仿真建模方法。

背景技术

分布式光纤声传感(Distributed fiber Acoustic Sensing，DAS)系统是近十余年来发展最为快速的光纤传感技术之一，该系统基于相位敏感时域反射计(φ-OTDR)原理，由窄线宽激光器向光纤注入窄脉冲光，通过提取光纤内不同时刻的背向瑞利散射光信号之间的相位差，实现光纤上各个位置的震动波场探测。分布式光纤声传感(DAS)系统可以实现现实空间连续的全尺度分布式监测，并且具有监测距离长、测量参数多、灵敏度高、响应速度快、信噪比好、无需外场供电、以及抗电磁干扰等优点，是当前国防、地震、石油、电力等领域迫切需求的技术。

分布式光纤声传感(DAS)系统涉及到的光学器件较多，多为模块化器件，价格昂贵，器件精密度高，如果购买的光学器件参数匹配不合理，会影响系统性能，甚至会烧毁器件，从而增加实验成本和研发周期，因此，有必要开展基于光学器件实际特性的建模仿真研究。此外，基于系统仿真也能够容易的获得光信号经过不同光学器件后的光学现象，这对于分布式光纤声传感(DAS)系统的调试、安装、故障诊断以及系统优化等都具有重要作用。目前尚未见分布式光纤声传感(DAS)系统的Simulink建模方法与系统报道。

发明内容

本发明实施例提供了一种分布式光纤声传感系统的仿真建模方法，基于分布式光纤声传感系统光电器件的数学模型和使用说明书，搭建器件模型，封装后形成Simulink器件模块，并形成器件模块库：使用者从器件模块库中选择合适的器件模块，能够自由的搭建自己需要的分布式光纤声传感系统结构；依据器件模块搭建的器件模块能够仿真器件参数对系统运行的影响，能够验证和优化使用者的设计方案；在器件模块的输出信号部分，容易的连接至To Workspace，获得系统运行后光信号经过不同光学器件后的光学现象，为分布式光纤声传感系统的调试、安装、故障诊断以及系统优化提供直观依据，解决了没有分布式光纤声传感系统Simulink建模方法与系统的问题。

一种分布式光纤声传感系统的仿真建模方法，包括以下步骤：

步骤1，根据激光器数学模型，基于Simulink模块，构建激光器模型，并封装为激光器模块；

步骤2，根据分束器数学模型，基于Simulink模块，构建分束器模型，并封装为分束器模块；

步骤3，根据声光调制器数学模型，基于Simulink模块，构建声光调制器模型，并封装为声光调制器模块；

步骤4，根据掺铒光纤放大器数学模型和衰减器数学模型，基于Simulink模块，构建掺铒光纤放大器和衰减器模型，并封装为掺铒光纤放大器模块；

步骤5，根据传感光纤数学模型，基于Simulink模块，构建传感光纤模型，并封装为传感光纤模块；

步骤6，根据光电探测器数学模型，基于Simulink模块，构建光电探测器模型，并封装为光电探测器模块；

步骤7，根据外界扰动的数学模型，基于Simulink模块，构建外界扰动模型，并封装为外界扰动模块；

步骤8，根据光纤耦合器的数学模型，基于Simulink模块，构建光纤耦合器模型，并封装为光纤耦合器模块；

步骤9，根据分布式光纤声传感系统的光路图，将对应的功能模块添加到Simulink工作区，利用信号线，将模块的输入信号和其他模块的输出信号名称相同的相连，并将模块的输出信号连接到To Workspace，完成分布式光纤声传感系统的Simulink建模；

步骤10，根据分布式光纤声传感系统的参数，计算Simulink系统的仿真时间和采样间隔，设置分布式光纤声传感系统各个功能模块参数，采用固定步长、ode1求解方法，对系统进行仿真运算，从各个模块输出信号的To Workspace获取计算数据，进行分析。

进一步的，所述步骤1中，当考虑激光器的相位噪声、频移特性时，振幅稳定单模激光器的光场为：

式中，E₀是激光器的电场强度，

P₀是激光器的功率，ω＝2πf₀，f₀是激光器中心频率，

是激光器的相位噪声，γ是激光器频移率，其中，

是由自发辐射引起的，可以认为是维纳过程，其增量服从一个均值为0的高斯型正态分布，即：

式中，

是均值为0、方差为2πτΔf的高斯分布，τ是采样间隔，Δf为激光器的线宽。

进一步的，所述步骤2中，激光器的光场E(t)，经过分束其后分成两部分，设分束器的分束比例为g1:g2，则一部分为激光器光场E(t)的g1倍，即为g₁E(t)，另一部分则为激光器光场E(t)的g2倍，即为g₂E(t)。

进一步的，所述步骤3中，激光器输出的光场E(t)经过声光调制器调制后，变成一个宽度为w的脉冲信号，即

其中

为矩形函数，即t≤w当时，

为1，当t＞w时，

为0。

进一步的，所述步骤4中，设掺铒光纤放大器工作在线性放大区域内，其增益为G，则掺铒光纤放大器的输出信号P_o(t)和输入信号P_i(t)之间的关系为：P_o(t)＝P_i(t)10^G/10，设衰减器的调节倍数为At，则经过衰减器调节后的输出光信号为A_tP_o(t)。

进一步的，所述步骤5中给出的传感光纤电场计算式E_cobs(t)数学模型为：

式中，

为背向散射光功率捕获因子，n_f是光纤纤芯的折射率，λ为背向散射光的波长，λ＝c/f₀，c是真空中的光速，f₀是激光器中心频率，r为光纤的模场半径，α为光在光纤的衰减系数，v_g＝c/n_f，N是整个传感光纤的散射体个数，τ_i是时间延迟，l_i是从入射面到第i个散射点的距离，且τ_i＝2l_i/v_g，r_i是第i个散射点的反射率，φ_i是第i个散射点的相位延迟，r_i满足服从瑞利分布，φ_i在(-π,π)内服从均匀分布，

是矩形函数，当

时，

否则矩形函数为0，β为光在光纤中的传播常数，

由于整个光纤被划分为N散射体，设脉冲宽度w下恰好覆盖n_sc个散射体，每个散射体的长度是Δl，则当脉冲上升沿到达第i个散射体时，上述公式修改为传感光纤电场计算式E_cobs(t_i)：

式中，m＝i-n+1，

由此可得背向瑞利散射光的功率P_cobs(t)为：

当光源是非相干光源时，由传感光纤数学模型和背向瑞利散射光功率数学模型并考虑

可得非相干光源光功率P_obs(t)为：

式中，n_sc是脉冲光所覆盖的散射体的个数，

对比时域反射计的光功率P_mobs(t)：

P_mobs(t)＝0.5P₀exp(-αv_gt)Sα_sv_gw

可知，

式中，α_s为瑞利散射系数，w表示光脉冲宽度；

基于Matlab的随机数生成函数rand可以实现瑞利分布和均匀分布，其中，瑞利分布为

r_is为瑞利分布系数，

通过以下方法确定瑞利分布系数r_is：

式中，mean{-2r_islog[1-rand(1,N)]}表示对-2r_islog[1-rand(1,N)]进行平均运算，rand(1,N)是Matlab的随机数生成函数，表示生成了N个(0,1)范围内的随机数，

考虑在坐标l_i处有外界扰动ε(l_i)时，折射率为，

式中标p₁₂光弹常数，取p₁₂＝0.270，t_i＝2l_i/v_g；

将n_fc(l_i)带入传感光纤电场计算式E_cobs(t_i)中的传播常数β，即可实现对外界扰动信号的感知；

根据传感光纤电场计算式E_cobs(t_i)、背向瑞利散射光的功率P_cobs(t)、非相干光源光功率P_obs(t)、时域反射计的光功率P_mobs(t)，基于S-function模块编写Sensingfiber.m文件函数，实现上述公式的计算，Sensingfiber.m模块中输入量分布为系统参数向量fParam，数组长度N1，光纤入射电场信号EAOM和扰动信号Epsilon，输出量分别是背向瑞利散射光的功率P_cobs(t)、时域反射计的光功率P_mobs(t)、非相干光源光功率P_obs(t)、光纤长度离散值l和传感光纤电场计算式E_cobs(t_i)，由于在simulink计算中无法利用历史数据，为此，引入持久变量(即persistent变量)counter，构建复数数组EAOMa，记录历史数据，同时散射点的随机反射率ri和相位延迟phii，也定义为持久变量，在第一次计算中，实现所有散射点的随机反射率和相位延迟的计算；

在S-function模块外部，通过VectorConcatenate将11个光纤参数常量模块转化为系统参数向量fParam，通过名称为N，数值为zeros(1,N)的常数模块产生数组长度N的输入，其中，这11个光纤参数模块分别是：名称为P0、数值为P0的激光光强，名称为L0、数值为L0的光纤长度，名称为lambda、数值为lambda的激光波长，名称为c、数值为c的光速，名称为nf、数值为nf的纤芯折射率，名称为nsc、数值为nsc的脉冲覆盖是单元数，名称为alpha、数值为alpha的光在光纤中的衰减系数，名称为alphas、数值为alphas的瑞利散射系数，名称为w、数值为w的光脉冲宽度，名称为ris、数值为ris的瑞利分布系数和名称为rfiber、数值为rfiber的光纤的模场半径，对上述模型进行封装，形成光纤传感模块，即Fiber模块，将11个光纤参数和数组长度N设置为封面参数。

进一步的，所述步骤6中，光电探测器的作用是将光功率信号P(t)转换为能被数据采集卡采集的电信号VP(t)，不靠滤波特性时，光电探测器的数学模型是：

式中，P_min是最小探测功率，k_p是转换增益，a_p是调整系数，P_s是饱和输入功率，b_p是非线性系数，c_p是衰减系数，当P(t)＝2P_s时，设

于是可得c_p＝0.0465，b_p由P＝a_pP_s处的连续性确定，即b_p＝k_p/c_p，

当考虑总体输出噪声电压时，上式表示为：

VP(t)＝VP'(t)+u_n(t)

式中，u_n(t)＝2u₀·rand-u₀，用于模拟产生±u₀的噪声电压，rand是(0,1)区间的标准均匀分布随机函数；

由光电探测器说明书中的噪声等效功率和带宽的平方根相乘可得最小可测功率P_min，由探测器响应度和夸阻增益相乘可得转换增益k_p，并设a_p＝1.34，再根据说明书中的饱和输入功率，可以实现光电探测器数学模型的计算；

当探测器是双平衡探测器时，有正功率输入P_p(t)和负功率输入P_n(t)，结合光电探测器的数学模型和噪声电压，有双平衡探测器输出数学模型：

VP(t)＝VP_p(t)-VP_n(t)

其中，VP_p(t)＝VP'_p(t)+u_np(t)，VP_n(t)＝VP'_n(t)+u_nn(t)，需要注意的是，当负功率输入P_n(t)为0时，u_nn(t)也为零。

进一步的，所述步骤7中，仅考虑脉冲扰动时，ε(t)＝ε₀δ(t-t₀)，当t＝t₀时，δ(t-t₀)＝1，

根据上述关系，利用名称为Step1、参数Steptime值为llong*ts、Initialvalue值为eps0、Finalvalue值为0、Sampletime值为ts的阶跃模块和名称为Step2、参数Steptime值为(llong-1)*ts、Initialvalue值为0、Finalvalue值为eps0、Sampletime值为ts的阶跃模块以及乘法模块产生脉冲扰动ε(t)，输出量为EpsilonOut，其中t0＝llong*ts，

对上述模型进行封装，形成外界扰动模块，即Epsilonp模块，将扰动位置点数llong和时间间隔ts设置为封面参数。

进一步的，所述步骤8中，对于输入光E_In1(t)和输入光E_In2(t)经过光纤耦合器后分别再产生两路光：

上述四路光混合后，形成两路拍频信号：

当被光电探测器接受时，转化成光功率信号：

式中，real()函数表示取复数的实部。

进一步的，所述步骤10中，仿真时间和采样间隔计算方法为：

仿真时间为：

T＝2L₀/v_g

散射体单元个数为：

N＝fix(T/2w)，fix()表示对计算结果取整数

采样间隔为：

t_x＝T/(N-1)。

本发明实施例提供的上述技术方案的有益效果至少包括：

本发明基于分布式光纤声传感系统光电器件的数学模型和使用说明书，搭建器件模型，封装后形成Simulink器件模块，并形成器件模块库：使用者从器件模块库中选择合适的器件模块，能够自由的搭建自己需要的分布式光纤声传感系统结构；依据器件模块搭建的器件模块能够仿真器件参数对系统运行的影响，能够验证和优化使用者的设计方案；在器件模块的输出信号部分，容易的连接至To Workspace，获得系统运行后光信号经过不同光学器件后的光学现象，为分布式光纤声传感系统的调试、安装、故障诊断以及系统优化提供直观依据，解决了没有分布式光纤声传感系统Simulink建模方法与系统的问题。

本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在所写的说明书、权利要求书、以及附图中所特别指出的结构来实现和获得。

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。在附图中：

图1为本发明实施例公开的分布式光纤声传感系统的仿真建模方法的流程示意图；

图2为本发明实施例公开的步骤1中激光器模块的结构示意图；

图3为本发明实施例公开的步骤2中分束器模块的结构示意图；

图4为本发明实施例公开的步骤3中声光调制器模块的结构示意图；

图5为本发明实施例公开的步骤4中掺铒光纤放大器模块的结构示意图；

图6为本发明实施例公开的步骤5中传感光纤模块的结构示意图；

图7为本发明实施例公开的步骤6中光电探测器模块的结构示意图；

图8为本发明实施例公开的步骤7中外界扰动模块的结构示意图；

图9为本发明实施例公开的步骤8中光纤耦合器模块的结构示意图；

图10为本发明实施例1中公开的直接探测式分布式光纤声传感系统模型结构示意图；

图11为本发明实施例1中公开的仿真获得光纤被散射光功功率曲线示意图；

图12为本发明实施例1中公开的仿真得到的光电探测器输出电压曲线示意图；

图13为本发明实施例1中公开的加扰动后的光电探测器输出差值电压曲线示意图；

图14为本发明实施例2中公开的外差探测式分布式光纤声传感系统模型的结构示意图；

图15为本发明实施例2中公开的仿真得到的光电探测器输出电压曲线示意图；

图16为本发明实施例2中公开的加扰动后的光电探测器输出差值电压曲线示意图。

具体实施方式

下面将参照附图更详细地描述本公开的示例性实施例。虽然附图中显示了本公开的示例性实施例，然而应当理解，可以以各种形式实现本公开而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更透彻地理解本公开，并且能够将本公开的范围完整的传达给本领域的技术人员。

如图1-9所示，本发明实施例提供一种分布式光纤声传感系统的仿真建模方法，包括以下步骤：

步骤1，根据激光器数学模型，基于Simulink模块，构建激光器模型，并封装为激光器模块；

当考虑激光器的相位噪声、频移特性时，振幅稳定单模激光器的光场为：

式中，E₀是激光器的电场强度，

P₀是激光器的功率，ω＝2πf₀，f₀是激光器中心频率，

是激光器的相位噪声，γ是激光器频移率，其中，

是由自发辐射引起的，可以认为是维纳过程，其增量服从一个均值为0的高斯型正态分布，即：

式中，

是均值为0、方差为2πτΔf的高斯分布，τ是采样间隔，Δf为激光器的线宽；

根据上述公式，基于Simulink模块激光器模型建模方法如下，利用名称为E0、变量值为E0的常数模块产生激光器的电场强度E₀。利用名称为omega、变量值为omega的常数模块，名称为gamma、变量值为gamma的常数模块以及加法模块产生(ω+γ)。时钟模块与名称为t0、变量值为t0的常数模块以及加法模块产生(t+t₀)，进而利用乘法模块产生(ω+γ)(t+t₀)。利用名称为fs、变量值为fs的常数模块、白噪声模块、乘法模块和名称为phi0、变量值为phi0的常数模以及块积分模块产生包含激光器的相位噪声

其中白噪声模块中的Noisepower参数值为ts*(2*pi*ts*fw)，Sampletime参数为ts，seed参数为[23341]。(ω+γ)(t+t₀)与

通过加法模块产生

激光器的电场强度E₀作为复数的模，激光器的(ω+γ)(t+t₀)作为幅角，基于Magnitude-AngletoComplex模块产生激光器的光场

输出接口为光场输出Ecout和相位输出phi。对上述模型进行封装，形成激光器模块，即Laser模块，将电场强度E0、激光器中心频率omega、频移率gamma、采样间隔ts、线宽fw、开始时间t0和初始相位phi0设置为封面参数。

步骤2，根据分束器数学模型，基于Simulink模块，构建分束器模型，并封装为分束器模块；

激光器的光场E(t)，经过分束其后分成两部分，设分束器的分束比例为g1:g2，则一部分为激光器光场E(t)的g1倍，即为g₁E(t)，另一部分则为激光器光场E(t)的g2倍，即为g₂E(t)；

根据上述关系，利用名称为Gain1、变量值为g1的比例模块产生g1倍光场的分束，输出量为SPbOut，利用名称为Gain2、变量值为g2的比例模块产生g2倍光场的分束，输出量为SPsOut。对上述模型进行封装，形成分束器模块，即Beam Splitter模块，将分束比例g1、分束比例g2设置为封面参数。

步骤3，根据声光调制器数学模型，基于Simulink模块，构建声光调制器模型，并封装为声光调制器模块；

激光器输出的光场E(t)经过声光调制器调制后，变成一个宽度为w的脉冲信号，即

其中

为矩形函数，即t≤w当时，

为1，当t＞w时，

为0；

根据上述关系，利用名称为Step、参数Step time值为nsc*ts、Initial value值为1、Final value值为0、Sample time值为ts的阶跃模块和名称为Product的乘法模块产生声光调制器调制后的输出

输出量为EAOMOut。对上述模型进行封装，形成声光调制器模块，即AOM模块，将每个脉冲覆盖的光纤单元个数nsc和时间间隔ts设置为封面参数。

步骤4，根据掺铒光纤放大器数学模型和衰减器数学模型，基于Simulink模块，构建掺铒光纤放大器和衰减器模型，并封装为掺铒光纤放大器模块；

设掺铒光纤放大器工作在线性放大区域内，其增益为G，则掺铒光纤放大器的输出信号P_o(t)和输入信号P_i(t)之间的关系为：P_o(t)＝P_i(t)10^G/10，设衰减器的调节倍数为At，则经过衰减器调节后的输出光信号为A_tP_o(t)；

根据上述关系，利用名称为Constant、参数为G的常数模块，名称为Gain1、参数为1/10的放大模块与名称为MathFuntion、函数为10u的函数模块产生10^G/10；再通过与名称为Product的乘法模块与输入信号EAOMIn相乘形成产生P_i(t)10^G/10；进而再经过名称为Gain2、参数为At的放大模块产生衰减器的调节。输出信号为EDFAOut；

对上述模型进行封装，形成脉冲掺铒光纤放大器模块，即EDFA模块，将掺铒光纤放大器增益G和衰减器的调节倍数At为封面参数。

步骤5，根据传感光纤数学模型，基于Simulink模块，构建传感光纤模型，并封装为传感光纤模块；

给出的传感光纤电场计算式E_cobs(t)数学模型为：

式中，

为背向散射光功率捕获因子，n_f是光纤纤芯的折射率，λ为背向散射光的波长，λ＝c/f₀，c是真空中的光速，f₀是激光器中心频率，r为光纤的模场半径，α为光在光纤的衰减系数，v_g＝c/n_f，N是整个传感光纤的散射体个数，τ_i是时间延迟，l_i是从入射面到第i个散射点的距离，且τ_i＝2l_i/v_g，r_i是第i个散射点的反射率，φ_i是第i个散射点的相位延迟，r_i满足服从瑞利分布，φ_i在(-π,π)内服从均匀分布，

是矩形函数，当

时，

否则矩形函数为0，β为光在光纤中的传播常数，

由于整个光纤被划分为N散射体，设脉冲宽度w下恰好覆盖n_sc个散射体，每个散射体的长度是Δl，则当脉冲上升沿到达第i个散射体时，上述公式修改为传感光纤电场计算式E_cobs(t_i)：

式中，m＝i-n+1，

由此可得背向瑞利散射光的功率P_cobs(t)为：

当光源是非相干光源时，由传感光纤数学模型和背向瑞利散射光功率数学模型并考虑

可得非相干光源光功率P_obs(t)为：

式中，n_sc是脉冲光所覆盖的散射体的个数，

对比时域反射计的光功率P_mobs(t)：

P_mobs(t)＝0.5P₀exp(-αv_gt)Sα_sv_gw

可知，

式中，α_s为瑞利散射系数，w表示光脉冲宽度；

本方案提出的传感光纤电场计算式E_cobs(t)数学模型克服了常用背向瑞利散射一维脉冲响应模型中不考虑背向散射光功率捕获因子S、光的传播相位β以及当光源是非相干光时光功率也不能退化为时域反射计通用公式P_mobs(t)的不足；

基于Matlab的随机数生成函数rand可以实现瑞利分布和均匀分布，其中，瑞利分布为

r_is为瑞利分布系数，

通过以下方法确定瑞利分布系数r_is：

式中，mean{-2r_islog[1-rand(1,N)]}表示对-2r_islog[1-rand(1,N)]进行平均运算，rand(1,N)是Matlab的随机数生成函数，表示生成了N个(0,1)范围内的随机数，

考虑在坐标l_i处有外界扰动ε(l_i)时，折射率为，

式中标p₁₂光弹常数，取p₁₂＝0.270，t_i＝2l_i/v_g；

将n_fc(l_i)带入传感光纤电场计算式E_cobs(t_i)中的传播常数β，即可实现对外界扰动信号的感知；

根据传感光纤电场计算式E_cobs(t_i)、背向瑞利散射光的功率P_cobs(t)、非相干光源光功率P_obs(t)、时域反射计的光功率P_mobs(t)，基于S-function模块编写Sensingfiber.m文件函数，实现上述公式的计算，Sensingfiber.m模块中输入量分布为系统参数向量fParam，数组长度N1，光纤入射电场信号EAOM和扰动信号Epsilon，输出量分别是背向瑞利散射光的功率P_cobs(t)、时域反射计的光功率P_mobs(t)、非相干光源光功率P_obs(t)、光纤长度离散值l和传感光纤电场计算式E_cobs(t_i)，由于在simulink计算中无法利用历史数据，为此，引入持久变量(即persistent变量)counter，构建复数数组EAOMa，记录历史数据，同时散射点的随机反射率ri和相位延迟phii，也定义为持久变量，在第一次计算中，实现所有散射点的随机反射率和相位延迟的计算；

在S-function模块外部，通过VectorConcatenate将11个光纤参数常量模块转化为系统参数向量fParam，通过名称为N，数值为zeros(1,N)的常数模块产生数组长度N的输入，其中，这11个光纤参数模块分别是：名称为P0、数值为P0的激光光强，名称为L0、数值为L0的光纤长度，名称为lambda、数值为lambda的激光波长，名称为c、数值为c的光速，名称为nf、数值为nf的纤芯折射率，名称为nsc、数值为nsc的脉冲覆盖是单元数，名称为alpha、数值为alpha的光在光纤中的衰减系数，名称为alphas、数值为alphas的瑞利散射系数，名称为w、数值为w的光脉冲宽度，名称为ris、数值为ris的瑞利分布系数和名称为rfiber、数值为rfiber的光纤的模场半径，对上述模型进行封装，形成光纤传感模块，即Fiber模块，将11个光纤参数和数组长度N设置为封面参数；

步骤6，根据光电探测器数学模型，基于Simulink模块，构建光电探测器模型，并封装为光电探测器模块；

光电探测器的作用是将光功率信号P(t)转换为能被数据采集卡采集的电信号VP(t)，不靠滤波特性时，光电探测器的数学模型是：

式中，P_min是最小探测功率，k_p是转换增益，a_p是调整系数，P_s是饱和输入功率，b_p是非线性系数，c_p是衰减系数，当P(t)＝2P_s时，设

于是可得c_p＝0.0465，b_p由P＝a_pP_s处的连续性确定，即b_p＝k_p/c_p，

当考虑总体输出噪声电压时，上式表示为：

VP(t)＝VP'(t)+u_n(t)

式中，u_n(t)＝2u₀·rand-u₀，用于模拟产生±u₀的噪声电压，rand是(0,1)区间的标准均匀分布随机函数；

由光电探测器说明书中的噪声等效功率和带宽的平方根相乘可得最小可测功率P_min，由探测器响应度和夸阻增益相乘可得转换增益k_p，并设a_p＝1.34，再根据说明书中的饱和输入功率，可以实现光电探测器数学模型的计算；

当探测器是双平衡探测器时，有正功率输入P_p(t)和负功率输入P_n(t)，结合光电探测器的数学模型和噪声电压，有双平衡探测器输出数学模型：

VP(t)＝VP_p(t)-VP_n(t)

其中，VP_p(t)＝VP'_p(t)+u_np(t)，VP_n(t)＝VP'_n(t)+u_nn(t)，需要注意的是，当负功率输入P_n(t)为0时，u_nn(t)也为零；

根据光电探测器的数学模型，基于S-function模块编写Photoelectricdetector.m文件函数，实现上述公式的计算。S-function模块的输入分别是光电探测器参数向量DParam、正功率输入Pp和负功率输入Pn，输出则是差值功率Vpcobs。在S-function模块外部通过，通过VectorConcatenate将4个光电探测器参数常量模块转化为光电探测器参数向量DParam，这四个常数模块分别是名称为Pmin、数值为Pmin的最小探测功率，名称为Ps、数值为Ps的饱和功率，名称为ap、数值为ap的调整系数和名称为kp、数值为kp的转换增益；

对上述模型进行封装，形成光电探测器模块，即Photoelectric detector模块，将4个光电探测器参数设置为封面参数。

步骤7，根据外界扰动的数学模型，基于Simulink模块，构建外界扰动模型，并封装为外界扰动模块；

仅考虑脉冲扰动时，ε(t)＝ε₀δ(t-t₀)，当t＝t₀时，δ(t-t₀)＝1，

根据上述关系，利用名称为Step1、参数Steptime值为llong*ts、Initialvalue值为eps0、Finalvalue值为0、Sampletime值为ts的阶跃模块和名称为Step2、参数Steptime值为(llong-1)*ts、Initialvalue值为0、Finalvalue值为eps0、Sampletime值为ts的阶跃模块以及乘法模块产生脉冲扰动ε(t)，输出量为EpsilonOut，其中t0＝llong*ts，

对上述模型进行封装，形成外界扰动模块，即Epsilonp模块，将扰动位置点数llong和时间间隔ts设置为封面参数；

步骤8，根据光纤耦合器的数学模型，基于Simulink模块，构建光纤耦合器模型，并封装为光纤耦合器模块；

对于输入光E_In1(t)和输入光E_In2(t)经过光纤耦合器后分别再产生两路光：

上述四路光混合后，形成两路拍频信号：

当被光电探测器接受时，转化成光功率信号：

式中，real()函数表示取复数的实部；

根据上述关系，利用名称为Constant1、数值为0的常数模块为实部，名称为Constant2、数值为1的常数模块为虚部以及复数产生函数Real-Imag to Complex1生成复数

进而通过乘法模块Product1与输入光E_In1(t)相乘，产生E_In11(t)；利用名称为Constant2、数值为0的常数模块为实部，名称为Constant3、数值为1的常数模块为虚部以及复数产生函数Real-Imag to Complex2生成复数

进而通过乘法模块Product2与输入光E_In2(t)相乘，产生E_In21(t)；输入光E_In1(t)和输入光E_In21(t)经过加法模块Add1相加后产生E₁(t)，然后分成两路，一路经过复数转化实部虚部模块Comlex to Real-Imag1，产生实部信号和虚部信号，其中虚部信号经过数值为-1的比例模块Gain1后，再通过实部虚部转化复数模块Real-Imag to Comlex3产生输入光E₁(t)的共轭复数

再通过乘法模块Product3与E₁(t)直接相乘，形成

再经过复数转化实部虚部模块Comlex toReal-Imag3取实部，产生输出光功率信号P_Out1(t)；输入光E_In11(t)和输入光E_In2(t)经过加法模块Add2相加后生产E₂(t)，然后分成两路，一路经过复数转化实部虚部模块Comlex toReal-Imag2产生实部信号和虚部信号，其中虚部信号经过数值为-1的比例模块Gain2后，再通过实部虚部转化复数模块Real-Imagto Comlex4输入光E₂(t)的共轭复数

再通过乘法模块Product4与E₂(t)直接相乘，形成

再经过复数转化实部虚部模块Comlex to Real-Imag4取实部，产生输出光功率信号P_Out2(t)。

对上述模型进行封装，形成光纤耦合器模块，即Coupler模块。

步骤9，根据分布式光纤声传感系统的光路图，将对应的功能模块添加到Simulink工作区，利用信号线，将模块的输入信号和其他模块的输出信号名称相同的相连，并将模块的输出信号连接到ToWorkspace，完成分布式光纤声传感系统的Simulink建模；

步骤10，根据分布式光纤声传感系统的参数，计算Simulink系统的仿真时间和采样间隔，设置分布式光纤声传感系统各个功能模块参数，采用固定步长、ode1求解方法，对系统进行仿真运算，从各个模块输出信号的To Workspace获取计算数据，进行分析；

仿真时间和采样间隔计算方法为：

仿真时间为：

T＝2L₀/v_g

散射体单元个数为：

N＝fix(T/2w)，fix()表示对计算结果取整数

采样间隔为：

t_x＝T/(N-1)。

实施例1

如图10-13所示，直接探测式分布式光纤声传感系统模型与仿真方法：

1.搭建直接探测式分布式光纤声传感系统模型

将需要的功能模块从模型库中拖入Simulink模型工作区，用信号线连接各个模块，搭建直接探测式分布式光纤声传感系统模型，如图10所示；

2.参数计算

系统参数为：L0＝10000；lambda＝1550*1e-9；c＝299792458；nf＝1.4685；ndu＝20；则vg＝c/nf；

计算可得：

T＝9.7968e-05，N＝8163，ts＝1.2003e-08；

采用固定步长、ode1求解方法，其中，仿真时间为T，固定步长长度为ts；

3.模块参数设置

其它还有参数为：P0＝20/1000；alpha＝0.0495/1000；alphas＝1.020/1000；ris＝6.1465e-04；rfiber＝10.5*10^-6/2；G＝13；At＝0.2；eps0＝1*10^-6；Ps＝150；ap＝1.34；kp＝28.5。在激光器模块中：t0＝0；phi0＝0。

其中手动开关模块Manual Switch模块，选择Epsilon0模块，数值为0，表示没有外界扰动输入。

系统仿真完成后，以l.data为横坐标，Pcobs.data、Pmobs.data和Pobs.data为纵坐标获得光功率曲线如图11所示；

以l.data为横坐标，VPcobs.data为纵坐标获得光电探测器输出电压曲线如图12所示。同时设VPcobs1＝VPcobs；

4.改变参数做差值仿真

在激光器模块中：t0＝9.7968e-05；phi0＝0.7225。

其中手动开关模块Manual Switch模块，选择Epsilon1模块，表示有外界扰动。其中，扰动模块即Epsilonp的参数为llong＝820,ts＝9.7968e-05；lshort＝810。

仿真完成后，设VPcobs2＝VPcobs。以l.data为横坐标，VPcobs2.data-VPcobs1.data为纵坐标获得加扰动后的光电探测器输出差值电压曲线如图13所示。可见在820*ts*vg＝1000m处有一个明显的扰动。

实施例2

如图14-16所示，外差探测式分布式光纤声传感(DAS)系统模型与仿真方法，搭建的外差探测式分布式光纤声传感(DAS)系统模型如图14所示；

分束器模块，即Beam Splider，参数g1＝0.9,g2＝0.1。以l.data为横坐标，VPcobs.data为纵坐标获得光电探测器输出电压曲线如图15所示；

扰动模块(即Epsilonp)的参数为llong＝504,ts＝9.7968e-05；lshort＝500。以l.data为横坐标，VPcobs2.data-VPcobs1.data为纵坐标获得加扰动后的光电探测器输出差值电压曲线如图16所示。可见在504*ts*vg＝610m处有一个明显的扰动。

本发明公开的一种分布式光纤声传感系统的仿真建模方法，基于分布式光纤声传感系统光电器件的数学模型和使用说明书，搭建器件模型，封装后形成Simulink器件模块，并形成器件模块库：使用者从器件模块库中选择合适的器件模块，能够自由的搭建自己需要的分布式光纤声传感系统结构；依据器件模块搭建的器件模块能够仿真器件参数对系统运行的影响，能够验证和优化使用者的设计方案；在器件模块的输出信号部分，容易的连接至To Workspace，获得系统运行后光信号经过不同光学器件后的光学现象，为分布式光纤声传感系统的调试、安装、故障诊断以及系统优化提供直观依据，解决了没有分布式光纤声传感系统Simulink建模方法与系统的问题。

应该明白，公开的过程中的步骤的特定顺序或层次是示例性方法的实例。基于设计偏好，应该理解，过程中的步骤的特定顺序或层次可以在不脱离本公开的保护范围的情况下得到重新安排。所附的方法权利要求以示例性的顺序给出了各种步骤的要素，并且不是要限于所述的特定顺序或层次。

在上述的详细描述中，各种特征一起组合在单个的实施方案中，以简化本公开。不应该将这种公开方法解释为反映了这样的意图，即，所要求保护的主题的实施方案需要清楚地在每个权利要求中所陈述的特征更多的特征。相反，如所附的权利要求书所反映的那样，本发明处于比所公开的单个实施方案的全部特征少的状态。因此，所附的权利要求书特此清楚地被并入详细描述中，其中每项权利要求独自作为本发明单独的优选实施方案。

本领域技术人员还应当理解，结合本文的实施例描述的各种说明性的逻辑框、模块、电路和算法步骤均可以实现成电子硬件、计算机软件或其组合。为了清楚地说明硬件和软件之间的可交换性，上面对各种说明性的部件、框、模块、电路和步骤均围绕其功能进行了一般地描述。至于这种功能是实现成硬件还是实现成软件，取决于特定的应用和对整个系统所施加的设计约束条件。熟练的技术人员可以针对每个特定应用，以变通的方式实现所描述的功能，但是，这种实现决策不应解释为背离本公开的保护范围。

结合本文的实施例所描述的方法或者算法的步骤可直接体现为硬件、由处理器执行的软件模块或其组合。软件模块可以位于RAM存储器、闪存、ROM存储器、EPROM存储器、EEPROM存储器、寄存器、硬盘、移动磁盘、CD-ROM或者本领域熟知的任何其它形式的存储介质中。一种示例性的存储介质连接至处理器，从而使处理器能够从该存储介质读取信息，且可向该存储介质写入信息。当然，存储介质也可以是处理器的组成部分。处理器和存储介质可以位于ASIC中。该ASIC可以位于用户终端中。当然，处理器和存储介质也可以作为分立组件存在于用户终端中。

对于软件实现，本申请中描述的技术可用执行本申请所述功能的模块(例如，过程、函数等)来实现。这些软件代码可以存储在存储器单元并由处理器执行。存储器单元可以实现在处理器内，也可以实现在处理器外，在后一种情况下，它经由各种手段以通信方式耦合到处理器，这些都是本领域中所公知的。

上文的描述包括一个或多个实施例的举例。当然，为了描述上述实施例而描述部件或方法的所有可能的结合是不可能的，但是本领域普通技术人员应该认识到，各个实施例可以做进一步的组合和排列。因此，本文中描述的实施例旨在涵盖落入所附权利要求书的保护范围内的所有这样的改变、修改和变型。此外，就说明书或权利要求书中使用的术语“包含”，该词的涵盖方式类似于术语“包括”，就如同“包括，”在权利要求中用作衔接词所解释的那样。此外，使用在权利要求书的说明书中的任何一个术语“或者”是要表示“非排它性的或者”。

Claims (10)

1.一种分布式光纤声传感系统的仿真建模方法，其特征在于，包括以下步骤：

步骤1，根据激光器数学模型，基于Simulink模块，构建激光器模型，并封装为激光器模块；

步骤2，根据分束器数学模型，基于Simulink模块，构建分束器模型，并封装为分束器模块；

步骤3，根据声光调制器数学模型，基于Simulink模块，构建声光调制器模型，并封装为声光调制器模块；

步骤4，根据掺铒光纤放大器数学模型和衰减器数学模型，基于Simulink模块，构建掺铒光纤放大器和衰减器模型，并封装为掺铒光纤放大器模块；

步骤5，根据传感光纤数学模型，基于Simulink模块，构建传感光纤模型，并封装为传感光纤模块；

步骤6，根据光电探测器数学模型，基于Simulink模块，构建光电探测器模型，并封装为光电探测器模块；

步骤7，根据外界扰动的数学模型，基于Simulink模块，构建外界扰动模型，并封装为外界扰动模块；

步骤8，根据光纤耦合器的数学模型，基于Simulink模块，构建光纤耦合器模型，并封装为光纤耦合器模块；

步骤9，根据分布式光纤声传感系统的光路图，将对应的功能模块添加到Simulink工作区，利用信号线，将模块的输入信号和其他模块的输出信号名称相同的相连，并将模块的输出信号连接到To Workspace，完成分布式光纤声传感系统的Simulink建模；

步骤10，根据分布式光纤声传感系统的参数，计算Simulink系统的仿真时间和采样间隔，设置分布式光纤声传感系统各个功能模块参数，采用固定步长、ode1求解方法，对系统进行仿真运算，从各个模块输出信号的To Workspace获取计算数据，进行分析。

2.如权利要求1所述的一种分布式光纤声传感系统的仿真建模方法，其特征在于，所述步骤1中，当考虑激光器的相位噪声、频移特性时，振幅稳定单模激光器的光场为：

式中，E₀是激光器的电场强度，

P₀是激光器的功率，ω＝2πf₀，f₀是激光器中心频率，

是激光器的相位噪声，γ是激光器频移率，其中，

是由自发辐射引起的，可以认为是维纳过程，其增量服从一个均值为0的高斯型正态分布，即：

式中，

是均值为0、方差为2πτΔf的高斯分布，τ是采样间隔，Δf为激光器的线宽。

3.如权利要求1所述的一种分布式光纤声传感系统的仿真建模方法，其特征在于，所述步骤2中，激光器的光场E(t)，经过分束其后分成两部分，设分束器的分束比例为g1:g2，则一部分为激光器光场E(t)的g1倍，即为g₁E(t)，另一部分则为激光器光场E(t)的g2倍，即为g₂E(t)。

4.如权利要求1所述的一种分布式光纤声传感系统的仿真建模方法，其特征在于，所述步骤3中，激光器输出的光场E(t)经过声光调制器调制后，变成一个宽度为w的脉冲信号，即

其中

为矩形函数，即t≤w当时，

为1，当t＞w时，

为0。

5.如权利要求1所述的一种分布式光纤声传感系统的仿真建模方法，其特征在于，所述步骤4中，设掺铒光纤放大器工作在线性放大区域内，其增益为G，则掺铒光纤放大器的输出信号P_o(t)和输入信号P_i(t)之间的关系为：P_o(t)＝P_i(t)10^G/10，设衰减器的调节倍数为At，则经过衰减器调节后的输出光信号为A_tP_o(t)。

6.如权利要求1所述的一种分布式光纤声传感系统的仿真建模方法，其特征在于，所述步骤5中给出的传感光纤电场计算式E_cobs(t)数学模型为：

式中，

为背向散射光功率捕获因子，n_f是光纤纤芯的折射率，λ为背向散射光的波长，λ＝c/f₀，c是真空中的光速，f₀是激光器中心频率，r为光纤的模场半径，α为光在光纤的衰减系数，v_g＝c/n_f，N是整个传感光纤的散射体个数，τ_i是时间延迟，l_i是从入射面到第i个散射点的距离，且τ_i＝2l_i/v_g，r_i是第i个散射点的反射率，φ_i是第i个散射点的相位延迟，r_i满足服从瑞利分布，φ_i在(-π,π)内服从均匀分布，

是矩形函数，当

时，

否则矩形函数为0，β为光在光纤中的传播常数，

由于整个光纤被划分为N散射体，设脉冲宽度w下恰好覆盖n_sc个散射体，每个散射体的长度是Δl，则当脉冲上升沿到达第i个散射体时，上述公式修改为传感光纤电场计算式E_cobs(t_i)：

式中，m＝i-n+1，

由此可得背向瑞利散射光的功率P_cobs(t)为：

当光源是非相干光源时，由传感光纤数学模型和背向瑞利散射光功率数学模型并考虑

可得非相干光源光功率P_obs(t)为：

式中，n_sc是脉冲光所覆盖的散射体的个数，

对比时域反射计的光功率P_mobs(t)：

P_mobs(t)＝0.5P₀exp(-αv_gt)Sα_sv_gw

可知，

式中，α_s为瑞利散射系数，w表示光脉冲宽度；

基于Matlab的随机数生成函数rand可以实现瑞利分布和均匀分布，其中，瑞利分布为

r_is为瑞利分布系数，

通过以下方法确定瑞利分布系数r_is：

式中，mean{-2r_islog[1-rand(1,N)]}表示对-2r_islog[1-rand(1,N)]进行平均运算，rand(1,N)是Matlab的随机数生成函数，表示生成了N个(0,1)范围内的随机数，

考虑在坐标l_i处有外界扰动ε(l_i)时，折射率为，

式中标p₁₂光弹常数，取p₁₂＝0.270，t_i＝2l_i/v_g；

将n_fc(l_i)带入传感光纤电场计算式E_cobs(t_i)中的传播常数β，即可实现对外界扰动信号的感知；

根据传感光纤电场计算式E_cobs(t_i)、背向瑞利散射光的功率P_cobs(t)、非相干光源光功率P_obs(t)、时域反射计的光功率P_mobs(t)，基于S-function模块编写Sensingfiber.m文件函数，实现上述公式的计算，Sensingfiber.m模块中输入量分布为系统参数向量fParam，数组长度N1，光纤入射电场信号EAOM和扰动信号Epsilon，输出量分别是背向瑞利散射光的功率P_cobs(t)、时域反射计的光功率P_mobs(t)、非相干光源光功率P_obs(t)、光纤长度离散值l和传感光纤电场计算式E_cobs(t_i)，由于在simulink计算中无法利用历史数据，为此，引入持久变量(即persistent变量)counter，构建复数数组EAOMa，记录历史数据，同时散射点的随机反射率ri和相位延迟phii，也定义为持久变量，在第一次计算中，实现所有散射点的随机反射率和相位延迟的计算；

在S-function模块外部，通过VectorConcatenate将11个光纤参数常量模块转化为系统参数向量fParam，通过名称为N，数值为zeros(1,N)的常数模块产生数组长度N的输入，其中，这11个光纤参数模块分别是：名称为P0、数值为P0的激光光强，名称为L0、数值为L0的光纤长度，名称为lambda、数值为lambda的激光波长，名称为c、数值为c的光速，名称为nf、数值为nf的纤芯折射率，名称为nsc、数值为nsc的脉冲覆盖是单元数，名称为alpha、数值为alpha的光在光纤中的衰减系数，名称为alphas、数值为alphas的瑞利散射系数，名称为w、数值为w的光脉冲宽度，名称为ris、数值为ris的瑞利分布系数和名称为rfiber、数值为rfiber的光纤的模场半径，对上述模型进行封装，形成光纤传感模块，即Fiber模块，将11个光纤参数和数组长度N设置为封面参数。

7.如权利要求1所述的一种分布式光纤声传感系统的仿真建模方法，其特征在于，所述步骤6中，光电探测器的作用是将光功率信号P(t)转换为能被数据采集卡采集的电信号VP(t)，不靠滤波特性时，光电探测器的数学模型是：

式中，P_min是最小探测功率，k_p是转换增益，a_p是调整系数，P_s是饱和输入功率，b_p是非线性系数，c_p是衰减系数，当P(t)＝2P_s时，设

于是可得c_p＝0.0465，b_p由P＝a_pP_s处的连续性确定，即b_p＝k_p/c_p，

当考虑总体输出噪声电压时，上式表示为：

VP(t)＝VP'(t)+u_n(t)

式中，u_n(t)＝2u₀·rand-u₀，用于模拟产生±u₀的噪声电压，rand是(0,1)区间的标准均匀分布随机函数；

由光电探测器说明书中的噪声等效功率和带宽的平方根相乘可得最小可测功率P_min，由探测器响应度和夸阻增益相乘可得转换增益k_p，并设a_p＝1.34，再根据说明书中的饱和输入功率，可以实现光电探测器数学模型的计算；

当探测器是双平衡探测器时，有正功率输入P_p(t)和负功率输入P_n(t)，结合光电探测器的数学模型和噪声电压，有双平衡探测器输出数学模型：

VP(t)＝VP_p(t)-VP_n(t)

其中，VP_p(t)＝VP'_p(t)+u_np(t)，VP_n(t)＝VP'_n(t)+u_nn(t)，需要注意的是，当负功率输入P_n(t)为0时，u_nn(t)也为零。

8.如权利要求1所述的一种分布式光纤声传感系统的仿真建模方法，其特征在于，所述步骤7中，仅考虑脉冲扰动时，ε(t)＝ε₀δ(t-t₀)，当t＝t₀时，δ(t-t₀)＝1，

根据上述关系，利用名称为Step1、参数Steptime值为llong*ts、Initialvalue值为eps0、Finalvalue值为0、Sampletime值为ts的阶跃模块和名称为Step2、参数Steptime值为(llong-1)*ts、Initialvalue值为0、Finalvalue值为eps0、Sampletime值为ts的阶跃模块以及乘法模块产生脉冲扰动ε(t)，输出量为EpsilonOut，其中t0＝llong*ts，

对上述模型进行封装，形成外界扰动模块，即Epsilonp模块，将扰动位置点数llong和时间间隔ts设置为封面参数。

9.如权利要求1所述的一种分布式光纤声传感系统的仿真建模方法，其特征在于，所述步骤8中，对于输入光E_In1(t)和输入光E_In2(t)经过光纤耦合器后分别再产生两路光：

上述四路光混合后，形成两路拍频信号：

当被光电探测器接受时，转化成光功率信号：

式中，real()函数表示取复数的实部。

10.如权利要求1所述的一种分布式光纤声传感系统的仿真建模方法，其特征在于，所述步骤10中，仿真时间和采样间隔计算方法为：

仿真时间为：

T＝2L₀/v_g

散射体单元个数为：

N＝fix(T/2w)，fix()表示对计算结果取整数

采样间隔为：

t_x＝T/(N-1)。

Images