基于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_1=B_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)的椭圆，c=d=0。但我们保留以增加鲁棒性。 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 )。但要注意偏置电压。

% 假设采集卡量程为 +/-5V，16位分辨率 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);

参数微调要点：

1. remove_DC_offset中的window_time：这个时间窗口必须大于待测信号的最低频率分量周期。如果待测信号包含接近直流的缓慢变化，这个窗口要设得足够长（比如几秒），否则去直流操作会抹掉低频信号。如果只有高频信号，窗口可以短一些（如0.1秒），以更快地跟踪直流分量的慢漂移。
2. estimate_ellipse_params的可靠性：椭圆拟合需要数据点均匀分布在椭圆轨迹上。如果相位变化 ( \phi(t) ) 幅度太小（比如小于π/2），点分布会集中在一小段弧上，导致拟合不准。确保你的传感相位有足够大的变化范围（最好超过π），或者在数据采集时主动引入一个大的、已知的相位调制（如用PZT拉伸光纤）来“画”出一个完整的椭圆。
3. 解缠失败：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 性能评估与常见问题排查

如何判断你的解调系统工作良好？

1. 信噪比（SNR）评估：在静态（无激励）条件下记录一段解调出的相位信号，计算其标准差 ( \sigma_{noise} )（单位：rad）。在动态（有激励）条件下，测量信号幅值 ( A_{signal} )（单位：rad）。则 SNR (dB) = ( 20 \log_{10}(A_{signal} / \sigma_{noise}) )。一个较好的干涉型传感系统，相位解调分辨率应能达到 ( 10^{-3} ) rad 量级或更好。
2. 线性度测试：施加已知幅值的阶梯状或扫频信号，看解调输出是否成比例，并计算非线性误差。
3. 三路信号质量检查：在施加周期性激励时，用示波器或软件同时观察三路原始信号I1, I2, I3。它们应该是三个幅值相近、形状相似、彼此存在近似120度相位差的正弦/余弦波。如果某一路信号幅值明显偏小或失真，检查对应的光电探测器或光纤连接头。
4. 解调结果跳变或失真：
   • 检查直流偏置估计：画出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算法，仍然是平衡了性能与复杂度的最佳选择。

整套代码和思路已经过多次实际项目的验证，从实验室的声发射检测到桥梁结构的应变监测，都表现出了可靠的性能。最关键的是理解每一步背后的物理和数学原理，这样当数据出现异常时，你才能像侦探一样，顺着信号链路的各个环节，快速定位问题是出在光路、电路还是算法上。希望这份超详细的拆解和“开箱即用”的代码，能帮你真正掌握这把干涉型光纤传感的“钥匙”。