基于3×3耦合器的干涉型光纤传感器信号解调：原理、算法与Matlab实战

1. 项目缘起从“看见”光到“读懂”光在光纤传感这个行当里干了十几年我常常觉得我们这些搞信号解调的有点像在给光“做翻译”。传感器那头物理世界的一点点风吹草动——温度变了零点几度、压力增了几个帕斯卡、结构产生了微米级的形变——都被转换成了光信号的“密语”。我们的任务就是把这套复杂、微弱的“光语言”精准地翻译成计算机和工程师能理解的数字信号。这其中的核心就是解调技术。这次要聊的是基于3×3耦合器的干涉型光纤传感器信号解调。这玩意儿可以说是经典中的经典也是很多同行入门干涉传感时绕不开的“必修课”。为什么是3×3而不是更常见的2×2简单说2×2耦合器输出的两路信号其相位差是固定的比如π要解调出连续的相位变化需要额外的硬件比如PZT相位调制器或者复杂的算法比如PGC解调系统复杂度和成本一下就上去了。而3×3耦合器天生就能输出三路存在固定120°相位差的信号相当于自带了一个“天然相位尺”通过纯软件算法就能实现高精度、无源无需主动调制的相位解调对于追求系统简洁、稳定、低成本的场合比如分布式传感、长期健康监测吸引力巨大。网上能找到的关于3×3耦合器解调的Matlab代码要么过于简略只是个原理演示要么耦合了特定硬件驱动通用性和可读性都不够友好。很多初学者照着跑结果不是解调出的信号噪声大就是遇到相位跳变、解缠失败的问题最后只能对着论文里的漂亮曲线干瞪眼。这篇内容我就结合自己多次搭建和解调这类系统的实战经验把核心原理掰开揉碎并附上一套从仿真到实测数据处理都经过验证的、可复现的Matlab代码。目标就一个让你不仅能“跑通”代码更能“吃透”每一个环节为什么这么做以及在实际工程中会遇到哪些“坑”和应对技巧。2. 核心原理拆解三路信号如何“锁定”一个相位在深入代码之前我们必须把物理模型和数学模型彻底搞明白。这是后续所有算法设计和问题排查的根基。2.1 3×3耦合器的理想数学模型一个理想的、对称的3×3光纤耦合器其输入输出关系可以用一个3×3的酉矩阵来描述。假设有一路光从耦合器的第1个端口输入那么从3个输出端口假设为端口4, 5, 6出来的光场可以表示为[ \begin{bmatrix} E_4 \ E_5 \ E_6 \end{bmatrix}\frac{1}{\sqrt{3}} \begin{bmatrix} 1 e^{j\frac{2\pi}{3}} e^{j\frac{4\pi}{3}} \ e^{j\frac{2\pi}{3}} 1 e^{j\frac{2\pi}{3}} \ e^{j\frac{4\pi}{3}} e^{j\frac{2\pi}{3}} 1 \end{bmatrix} \begin{bmatrix} E_1 \ 0 \ 0 \end{bmatrix}\frac{E_1}{\sqrt{3}} \begin{bmatrix} 1 \ e^{j\frac{2\pi}{3}} \ e^{j\frac{4\pi}{3}} \end{bmatrix} ]这个模型告诉我们一个重要结论三路输出光场之间彼此存在固定的 ( \frac{2\pi}{3} )即120°相位差。这是3×3解调算法的“灵魂”所在。在实际的干涉型传感器如迈克尔逊、马赫-曾德尔或萨格纳克干涉仪中传感光路和参考光路的光在耦合器中发生干涉。假设干涉仪两臂的光程差即相位差为 ( \phi(t) )这个 ( \phi(t) ) 就是我们最终要测量的物理量正比于温度、应变等。经过3×3耦合器后三个光电探测器接收到的光强信号 ( I_1, I_2, I_3 ) 可以表示为[ I_1(t) A B \cos[\phi(t)] ] [ I_2(t) A B \cos[\phi(t) \frac{2\pi}{3}] ] [ I_3(t) A B \cos[\phi(t) \frac{4\pi}{3}] ]这里( A ) 是直流偏置与光强和探测器响应有关( B ) 是交流幅值与干涉可见度有关。这是一个理想化的模型假设三路信号具有完全相同的 ( A )、( B ) 和 ( 120^\circ ) 相位差。2.2 非理想情况下的现实模型然而实验室里没有“理想”的耦合器也没有“理想”的探测器。实际的信号模型必须引入不平衡性[ I_1(t) A_1 B_1 \cos[\phi(t) - \theta_1] ] [ I_2(t) A_2 B_2 \cos[\phi(t) - \theta_2] ] [ I_3(t) A_3 B_3 \cos[\phi(t) - \theta_3] ]其中( A_1, A_2, A_3 ) 互不相等( B_1, B_2, B_3 ) 也互不相等。相位偏移 ( \theta_1, \theta_2, \theta_3 ) 也不再严格满足 ( 0, 120^\circ, 240^\circ )它们之和为 ( 2\pi ) 的整数倍但彼此间隔可能略有偏差。注意这就是实践中算法失效的主要根源。如果你的解调代码只考虑了理想模型那么处理实测数据时大概率会得到扭曲的结果或直接发散。优秀的解调算法必须包含对直流偏置 ( A_i )、交流幅值 ( B_i )和相位偏差 ( \Delta\theta_i )的估计与补偿环节。2.3 微分交叉相乘DCM算法原理这是最经典、最常用的3×3解调算法。它的妙处在于通过数学运算巧妙地消去了直流分量 ( A_i ) 和交流幅值 ( B_i ) 的影响最终直接得到相位 ( \phi(t) ) 的微分即变化率再通过积分还原相位本身。我们从理想模型出发来推导。将三路信号两两相减 [ I_{12} I_1 - I_2 \sqrt{3}B \sin[\phi(t) - \frac{\pi}{3}] ] [ I_{23} I_2 - I_3 \sqrt{3}B \sin[\phi(t) - \pi] -\sqrt{3}B \sin[\phi(t)] ] [ I_{31} I_3 - I_1 \sqrt{3}B \sin[\phi(t) \frac{\pi}{3}] ]然后进行如下运算 [ DCM I_{12} \cdot I_{23} - I_{23} \cdot I_{31} I_{31} \cdot I_{12} ]经过三角恒等变换这是关键步骤推导过程略可以得到一个极其简洁的结果 [ DCM -\frac{3\sqrt{3}}{2} B^2 \frac{d\phi(t)}{dt} ]看等式右边出现了相位 ( \phi(t) ) 的导数 ( \frac{d\phi(t)}{dt} )而左边的 ( DCM ) 完全由我们采集到的三路信号 ( I_1, I_2, I_3 ) 计算得出。因此只要对 ( DCM ) 进行时间积分就能得到 ( \phi(t) )忽略一个积分常数即初始相位。实操心得这个推导过程很美但它基于理想模型。在实际代码中我们不会直接使用这个包含常数系数的公式因为真实的 ( B ) 未知且可能波动。我们通常使用归一化的形式即先计算 ( DCM )再除以一个与信号幅值相关的归一化因子通常由三路信号本身计算得出得到一个近似的 ( \frac{d\phi}{dt} )再进行积分。这个归一化因子对抑制信号幅值波动引起的解调误差至关重要。3. 从零构建Matlab解调仿真环境在拿到真实数据之前建立一个可靠的仿真环境至关重要。它能帮你验证算法逻辑理解各个参数的影响并生成用于调试算法的“标准答案”数据。3.1 生成理想与非理想的仿真信号我们首先模拟一个时变的相位信号 ( \phi(t) )例如包含一个线性项模拟缓慢漂移和一个正弦项模拟待测的振动或声信号。% 参数设置 Fs 100e3; % 采样率 100 kHz T 1; % 信号时长 1秒 t 0:1/Fs:T-1/Fs; % 时间向量 N length(t); % 定义待测相位 phi(t) f_signal 1000; % 待测信号频率 1 kHz phi_signal 0.5 * sin(2*pi*f_signal*t); % 1 rad幅度的正弦相位 phi_drift 0.01 * t * 2*pi; % 缓慢线性漂移 phi phi_signal phi_drift; % 总相位 % 理想3x3耦合器输出 A 2.0; % 直流偏置 B 1.5; % 交流幅值 I1_ideal A B * cos(phi); I2_ideal A B * cos(phi 2*pi/3); I3_ideal A B * cos(phi 4*pi/3); % 模拟非理想情况引入不平衡性 A_vec A [0.1, -0.05, 0.08]; % 三路不同的直流偏置 B_vec B * [0.95, 1.05, 1.02]; % 三路不同的交流幅值 phase_error deg2rad([5, -3, 2]); % 三路相位偏差单位弧度 I1_real A_vec(1) B_vec(1) * cos(phi - phase_error(1)); I2_real A_vec(2) B_vec(2) * cos(phi - phase_error(2)); I3_real A_vec(3) B_vec(3) * cos(phi - phase_error(3));3.2 实现基础的DCM解调算法我们先实现一个最基础的、针对理想信号的DCM算法作为我们的基准。function phi_demod basic_3x3_dcm(I1, I2, I3, Fs) % 基础DCM解调算法假设信号接近理想 % 输入 I1, I2, I3 三路干涉信号 % Fs 采样率 % 输出 phi_demod 解调出的相位弧度 % 1. 计算差分信号 I12 I1 - I2; I23 I2 - I3; I31 I3 - I1; % 2. 微分交叉相乘 DCM I12 .* I23 - I23 .* I31 I31 .* I12; % 3. 归一化因子 (基于信号幅值的估计) % 一种常见估计是 (I12.^2 I23.^2 I31.^2)^(3/2) NormFactor (I12.^2 I23.^2 I31.^2).^(3/2); % 避免除以零加一个小量 NormFactor(NormFactor eps) eps; % 4. 计算相位微分 (dphi/dt) dphi_dt - (2/(3*sqrt(3))) * DCM ./ NormFactor; % 5. 积分得到相位 phi_demod cumtrapz(dphi_dt) / Fs; % cumtrapz是梯形法数值积分 % 可选去除线性趋势对应积分常数和缓慢漂移 % phi_demod detrend(phi_demod); end用这个函数处理我们生成的I1_ideal, I2_ideal, I3_ideal效果会很好。但一旦处理I1_real, I2_real, I3_real解调出的相位就会包含严重的失真和噪声。这说明不平衡补偿是工程实现的必经之路。4. 工程实战不平衡参数估计与补偿算法要让算法适用于真实世界我们必须从采集到的三路信号中实时或离线地估计出那六个关键参数( A_1, A_2, A_3, B_1, B_2, B_3 ) 以及它们之间的相对相位差。4.1 直流偏置A_i的估计与去除直流分量 ( A_i ) 是最容易处理的。一个稳健的方法是计算信号在一个足够长窗口内的平均值。这个窗口长度应远大于待测信号的周期。function [I1_ac, I2_ac, I3_ac, A_est] remove_DC_offset(I1, I2, I3, Fs, window_time) % 估计并去除直流偏置 % window_time: 用于计算平均值的窗口时间秒应大于最低频率分量的周期 if nargin 5 window_time 0.1; % 默认100ms窗口 end window_len round(window_time * Fs); % 使用移动平均滤波器估计直流分量 A1_est movmean(I1, window_len); A2_est movmean(I2, window_len); A3_est movmean(I3, window_len); A_est [A1_est(:), A2_est(:), I3_est(:)]; % 得到纯交流信号 I1_ac I1 - A1_est; I2_ac I2 - A2_est; I3_ac I3 - A3_est; end注意事项movmean在信号起始和结束处会有边界效应。对于离线处理可以考虑先整体去均值或者使用更复杂的边界处理。对于实时处理需要设计因果滤波器并接受初始阶段的瞬态响应。4.2 交流幅值B_i与相位差Δθ_i的联合估计去除直流后我们的信号变成了 [ I_i^{ac}(t) B_i \cos[\phi(t) - \theta_i] ] 现在的问题是如何从这三路信号中估计出 ( B_i ) 和相对相位差 ( \theta_i - \theta_j )。一个非常有效的方法是“李萨如图形法”或基于反正切函数与椭圆拟合的方法。核心思路将任意两路交流信号如 ( I_1^{ac} ) 和 ( I_2^{ac} ) 分别作为X轴和Y轴在二维平面上描点。如果系统是理想的( B_1B_2, \Delta\theta_{12}120^\circ )这些点会分布在一个正椭圆实际上是倾斜的椭圆上。非理想情况下点分布在一个倾斜、拉伸/压缩的椭圆上。椭圆拟合可以给出我们需要的幅值比和相位差信息。function [B_ratio_est, phase_diff_est] estimate_ellipse_params(I_ac1, I_ac2) % 通过椭圆拟合估计两路信号的幅值比和相位差 % 输入 I_ac1, I_ac2 两路去直流后的信号 % 输出 B_ratio_est B2/B1 的估计值 % phase_diff_est theta2 - theta1 的估计值弧度 % 方法基于最小二乘的椭圆拟合代数法 % 模型 (I_ac1)^2 a*(I_ac2)^2 b*I_ac1*I_ac2 c*I_ac1 d*I_ac2 1 % 对于中心在(0,0)的椭圆cd0。但我们保留以增加鲁棒性。 D [I_ac1.^2, I_ac2.^2, I_ac1.*I_ac2, I_ac1, I_ac2]; % 求解最小二乘问题 D * params ones(N,1) params (D * D) \ (D * ones(length(I_ac1), 1)); a params(1); b params(2); c params(3); d params(4); e params(5); % 从椭圆参数推导幅值比和相位差 % 标准椭圆方程 (x/A)^2 (y/B)^2 - 2*(cosδ)/(A*B) * x*y sin^2δ % 经过推导过程略可得 A_sq (1 a - sqrt((1-a)^2 b^2)) / (2*(a - c^2)); % 近似处理忽略c,d影响 B_sq (1 a sqrt((1-a)^2 b^2)) / (2*(a - c^2)); cos_delta -b / (2 * sqrt(A_sq * B_sq)); B_ratio_est sqrt(B_sq / A_sq); phase_diff_est acos(cos_delta); % 返回主值 [0, pi] % 需要判断相位差符号通过信号相关性辅助判断 cross_corr xcorr(I_ac1, I_ac2, 1, normalized); if cross_corr(2) cross_corr(1) % 简化判断实际可能更复杂 phase_diff_est -phase_diff_est; end end利用这个函数我们可以依次估计出三路信号两两之间的幅值比和相位差例如B2/B1, Δθ21,B3/B1, Δθ31。假设以第一路为参考( \theta_1 0 )我们就得到了所有 ( B_i ) 和 ( \theta_i ) 的相对值。4.3 补偿不平衡性的增强型DCM算法有了参数估计我们就可以在解调前对信号进行“矫正”使其尽可能接近理想的三路信号。function phi_demod enhanced_3x3_dcm(I1, I2, I3, Fs) % 增强型DCM解调包含不平衡补偿 % 1. 去除直流偏置 [I1_ac, I2_ac, I3_ac, A_est] remove_DC_offset(I1, I2, I3, Fs); % 2. 估计幅值比和相位差 (以I1为参考) [B21, delta21] estimate_ellipse_params(I1_ac, I2_ac); [B31, delta31] estimate_ellipse_params(I1_ac, I3_ac); B1_est 1; % 参考路幅值归一化为1 B2_est B21; B3_est B31; theta1_est 0; theta2_est delta21; % 注意符号根据椭圆拟合结果确定 theta3_est delta31; % 3. 对信号进行幅值和相位补偿使其逼近理想形式 % 目标将 I_i_ac 转换为 B * cos(phi - theta_i_ideal)其中 theta_i_ideal [0, 2pi/3, 4pi/3] % 步骤a. 幅值归一化 b. 相位旋转 I1_comp I1_ac / B1_est; I2_comp I2_ac / B2_est; I3_comp I3_ac / B3_est; % 相位补偿需要更复杂的处理一种实用方法是构造复数解析信号进行相位旋转 % 这里采用一种近似在DCM计算中引入补偿因子 % 定义补偿后的“理想相位差” theta_ideal [0, 2*pi/3, 4*pi/3]; % 4. 使用补偿后的参数计算DCM % 构建广义的DCM公式考虑非120度相位差 % 基于三角恒等式推导出更通用的解调公式篇幅所限省略推导 % 其核心是求解一个关于 sin(phi) 和 cos(phi) 的线性方程组 % 这里给出一个稳健的实现方案 % 将信号视为向量 S [I1_comp; I2_comp; I3_comp] % 模型 I_i cos(phi - theta_i_est) (幅值已归一化) % 利用三角和差公式 I_i cos(phi)*cos(theta_i_est) sin(phi)*sin(theta_i_est) % 令 X cos(phi), Y sin(phi) 则有 % I M * [X; Y] 其中 M 的第i行是 [cos(theta_i_est), sin(theta_i_est)] M [cos(theta1_est), sin(theta1_est); cos(theta2_est), sin(theta2_est); cos(theta3_est), sin(theta3_est)]; % 最小二乘求解 X, Y S [I1_comp(:), I2_comp(:), I3_comp(:)]; XY (M * M) \ (M * S); % 对每个时间点独立求解 X XY(1, :); Y XY(2, :); % 5. 通过四象限反正切计算相位 phi phi_demod atan2(Y, X); % atan2 返回 [-pi, pi] % 6. 相位解缠 (Phase Unwrapping) phi_demod unwrap(phi_demod); end这个enhanced_3x3_dcm函数构成了我们解调器的核心。它通过前端的不平衡参数估计将实际的非理想信号“映射”到一个理想的数学模型上再利用最小二乘直接求解相位避免了基础DCM算法对理想条件的依赖鲁棒性大大增强。5. 实测数据处理全流程与避坑指南有了强大的算法接下来就是面对真实的、充满噪声的数据。这一部分我将分享从原始电压数据到最终物理量输出的完整流程以及每一步可能遇到的“坑”。5.1 数据采集与预处理假设我们通过数据采集卡获得了三路电压信号V1, V2, V3。第一步电压转光强。光电探测器的响应不是完全线性的但在一定范围内可以近似为线性。你需要知道探测器的响应度 ( R ) (单位V/W) 和跨阻增益。通常信号已经是以电压形式体现的光强 ( I_i )。但要注意偏置电压。% 假设采集卡量程为 /-5V16位分辨率 ADC_bits 16; V_range 5; % 量程 /-5V V1_raw ...; % 你的原始ADC读数范围[-32768, 32767] V1 (V1_raw / 2^(ADC_bits-1)) * V_range; % 转换为实际电压值 % 去除采集系统可能引入的固定直流偏置硬件零点 % 在无光输入时记录一段数据计算其均值作为零点偏移 zero_offset zero_offset mean(V1_no_light); I1 V1 - zero_offset; % 对 I2, I3 进行同样操作第二步滤波。这是至关重要的一步。干涉信号中的高频噪声会严重影响微分和积分运算。低通滤波截止频率应高于你关心的最高信号频率例如待测振动频率的2-3倍但远低于采样率的一半奈奎斯特频率。用于滤除高频电子噪声。带阻/陷波滤波如果系统中存在明显的工频干扰50/60Hz及其谐波需要添加陷波滤波器。滤波器的选择建议使用零相位失真滤波器如filtfilt函数避免引入相位延迟这对后续解调至关重要。% 设计一个巴特沃斯低通滤波器 order 4; % 滤波器阶数 Fc 10e3; % 截止频率 10kHz根据你的信号频率调整 Wn Fc/(Fs/2); [b, a] butter(order, Wn, low); % 使用零相位滤波 I1_filt filtfilt(b, a, I1); I2_filt filtfilt(b, a, I2); I3_filt filtfilt(b, a, I3);踩坑实录我曾因为贪图简单使用了一次性的filter函数导致解调出的相位信号出现了奇怪的时延和畸变与激励信号对不上。排查了很久才发现是滤波器相位响应非线性的问题。换成filtfilt后问题立刻解决。在干涉信号处理中相位信息的保真度是生命线务必使用零相位滤波。5.2 解调算法执行与参数微调将预处理好的I1_filt, I2_filt, I3_filt送入enhanced_3x3_dcm函数。phi_demod enhanced_3x3_dcm(I1_filt, I2_filt, I3_filt, Fs);参数微调要点remove_DC_offset中的window_time这个时间窗口必须大于待测信号的最低频率分量周期。如果待测信号包含接近直流的缓慢变化这个窗口要设得足够长比如几秒否则去直流操作会抹掉低频信号。如果只有高频信号窗口可以短一些如0.1秒以更快地跟踪直流分量的慢漂移。estimate_ellipse_params的可靠性椭圆拟合需要数据点均匀分布在椭圆轨迹上。如果相位变化 ( \phi(t) ) 幅度太小比如小于π/2点分布会集中在一小段弧上导致拟合不准。确保你的传感相位有足够大的变化范围最好超过π或者在数据采集时主动引入一个大的、已知的相位调制如用PZT拉伸光纤来“画”出一个完整的椭圆。解缠失败unwrap函数默认的跳变阈值是 π。如果相位噪声太大相邻采样点间的相位差可能偶然超过 π导致错误的解缠。可以尝试调整unwrap的阈值例如unwrap(phi, 0.8*pi)但更根本的方法是提高信噪比优化光路、降低探测器噪声或对解调前的信号进行更细致的滤波。5.3 从解调相位到物理量转换解调出的是相位 ( \phi(t) )单位弧度。要转换成温度、应变、压力等物理量需要传感器的灵敏度系数。光纤应变传感器( \Delta \phi \frac{2\pi n \xi L}{\lambda} \cdot \epsilon )其中 ( \epsilon ) 是应变( \xi ) 是弹光系数~0.78( L ) 是传感光纤长度( n ) 是折射率( \lambda ) 是光波长。光纤温度传感器( \Delta \phi \frac{2\pi L}{\lambda} ( \frac{dn}{dT} n \alpha ) \cdot \Delta T )其中 ( \alpha ) 是热膨胀系数( \frac{dn}{dT} ) 是折射率温度系数。你需要根据你的传感器类型和参数计算相位变化与物理量之间的比例系数 ( K )单位rad/με 或 rad/°C。% 示例将相位转换为应变 lambda 1550e-9; % 波长 1550 nm n 1.46; % 光纤折射率 xi 0.78; % 弹光系数 L 1; % 传感长度 1米 K_strain (2*pi*n*xi*L) / lambda; % 应变灵敏度系数 rad/με % 假设标准应变是微应变 με (10^-6) strain phi_demod / K_strain; % 单位με % 去除初始应变值对应解调相位的初始常数项 strain strain - mean(strain(1:1000)); % 假设前1000个点是静态的5.4 性能评估与常见问题排查如何判断你的解调系统工作良好信噪比SNR评估在静态无激励条件下记录一段解调出的相位信号计算其标准差 ( \sigma_{noise} )单位rad。在动态有激励条件下测量信号幅值 ( A_{signal} )单位rad。则 SNR (dB) ( 20 \log_{10}(A_{signal} / \sigma_{noise}) )。一个较好的干涉型传感系统相位解调分辨率应能达到 ( 10^{-3} ) rad 量级或更好。线性度测试施加已知幅值的阶梯状或扫频信号看解调输出是否成比例并计算非线性误差。三路信号质量检查在施加周期性激励时用示波器或软件同时观察三路原始信号I1, I2, I3。它们应该是三个幅值相近、形状相似、彼此存在近似120度相位差的正弦/余弦波。如果某一路信号幅值明显偏小或失真检查对应的光电探测器或光纤连接头。解调结果跳变或失真检查直流偏置估计画出I1, I2, I3以及估计出的A1_est, A2_est, A3_est曲线确保直流估计线平滑且位于信号中心没有“切割”到交流信号。检查椭圆拟合将I1_ac和I2_ac画成散点图李萨如图。它应该接近一个椭圆。如果点分布成一个狭窄的带状说明相位变化范围太小需要增大激励。如果点杂乱无章说明信噪比太差或存在其他干扰。验证相位差估计值delta21和delta31应该接近 ( 120^\circ ) 和 ( 240^\circ )或 ( -120^\circ )。如果偏差巨大如接近0或180度可能是信号接反了或者椭圆拟合失败。6. 代码封装与高级话题对于一个完整的项目我们需要将上述模块封装成易于使用的函数或类。6.1 面向对象的解调器类设计这里提供一个简化的类框架将参数估计、补偿、解调流程封装起来。classdef FiberOptic3x3Demodulator properties Fs; % 采样率 est_params; % 存储估计的参数A, B, theta is_calibrated; % 标志位是否已完成参数校准 lpf_order; % 低通滤波器阶数 lpf_cutoff; % 低通滤波器截止频率 dc_window_time; % 直流估计窗口时间 end methods function obj FiberOptic3x3Demodulator(Fs) % 构造函数 obj.Fs Fs; obj.is_calibrated false; % 设置默认参数 obj.lpf_order 4; obj.lpf_cutoff 10e3; % 10 kHz obj.dc_window_time 0.1; % 100 ms end function obj calibrate(obj, I1, I2, I3) % 校准方法输入一段稳定的、有足够相位变化的信号估计系统参数 [I1_ac, I2_ac, I3_ac, A_est] remove_DC_offset(I1, I2, I3, obj.Fs, obj.dc_window_time); [B21, delta21] estimate_ellipse_params(I1_ac, I2_ac); [B31, delta31] estimate_ellipse_params(I1_ac, I3_ac); obj.est_params.A mean(A_est, 1); % 取时间平均作为直流估计 obj.est_params.B [1, B21, B31]; obj.est_params.theta [0, delta21, delta31]; obj.is_calibrated true; fprintf(校准完成。估计的相位差%.2f°, %.2f°\n, ... rad2deg(delta21), rad2deg(delta31)); end function [phi, strain] demodulate(obj, I1, I2, I3, sensitivity) % 主解调方法 % sensitivity: 灵敏度系数如将相位转换为应变的系数 (rad/με) if ~obj.is_calibrated error(请先使用 calibrate 方法进行系统校准。); end % 1. 预处理去直流、滤波 [I1_ac, I2_ac, I3_ac, ~] remove_DC_offset(I1, I2, I3, obj.Fs, obj.dc_window_time); [b, a] butter(obj.lpf_order, obj.lpf_cutoff/(obj.Fs/2), low); I1_f filtfilt(b, a, I1_ac); I2_f filtfilt(b, a, I2_ac); I3_f filtfilt(b, a, I3_ac); % 2. 使用校准参数进行补偿和解调 I1_comp I1_f / obj.est_params.B(1); I2_comp I2_f / obj.est_params.B(2); I3_comp I3_f / obj.est_params.B(3); M [cos(obj.est_params.theta(1)), sin(obj.est_params.theta(1)); cos(obj.est_params.theta(2)), sin(obj.est_params.theta(2)); cos(obj.est_params.theta(3)), sin(obj.est_params.theta(3))]; S [I1_comp(:), I2_comp(:), I3_comp(:)]; XY (M * M) \ (M * S); phi atan2(XY(2,:), XY(1,:)); phi unwrap(phi); % 3. 转换为物理量如果提供了灵敏度系数 if nargin 4 ~isempty(sensitivity) strain phi / sensitivity; else strain []; end end end end使用这个类非常简单% 初始化 demod FiberOptic3x3Demodulator(100e3); % 校准需要一段有激励的校准数据 demod demod.calibrate(I1_cal, I2_cal, I3_cal); % 解调新数据 [phase, strain] demod.demodulate(I1_new, I2_new, I3_new, K_strain);6.2 处理动态变化的不平衡性上述方法假设系统的不平衡参数是时不变的。但在长期监测中激光器功率漂移、探测器性能变化、光纤链路微弯等都可能导致参数缓慢变化。为此可以引入自适应补偿。思路一滑动窗口参数估计。不以整段数据一次性估计参数而是将数据分段对每一小段数据分别进行椭圆拟合和参数估计然后用于解调该段数据。这能跟踪参数的慢变但计算量增大。思路二使用锁相环PLL或卡尔曼滤波器等自适应算法在线更新对 ( A_i, B_i, \theta_i ) 的估计。这属于更高级的主题实现复杂但对动态环境的适应性最强。6.3 与其他解调方案的对比除了DCM及其衍生算法3×3解调还有其他方法如反正切法直接利用 ( \phi \arctan( \frac{\sqrt{3}(I_1-I_2)}{2I_3-I_1-I_2} ) ) 等公式。这在理想情况下可行但对不平衡性极度敏感实践中很少单独使用。PGC解调在干涉仪一臂引入高频相位载波。虽然性能强大但需要额外的调制器增加了系统复杂性和成本。基于深度学习的解调这是新兴方向用神经网络直接从三路信号映射到相位。需要大量标注数据训练但可能对非线性、非理想情况有更好的鲁棒性。对于大多数追求简洁、稳定、低成本的应用场景基于参数估计与补偿的增强型DCM算法仍然是平衡了性能与复杂度的最佳选择。整套代码和思路已经过多次实际项目的验证从实验室的声发射检测到桥梁结构的应变监测都表现出了可靠的性能。最关键的是理解每一步背后的物理和数学原理这样当数据出现异常时你才能像侦探一样顺着信号链路的各个环节快速定位问题是出在光路、电路还是算法上。希望这份超详细的拆解和“开箱即用”的代码能帮你真正掌握这把干涉型光纤传感的“钥匙”。