做干涉型光纤传感解调的人多半都经历过这种尴尬光强抖一下、调制深度偏一点辛辛苦苦写好的PGC解调算法立刻给你“脸色看”输出不是带毛刺就是谐波失真。今天要聊的PGC-SDD-DSM算法全称是单路径微分相除和微分自乘相减Phase Generated Carrier - Single-path Differential Division and Differential Self-Multiply Subtraction这类改进型解调方案的核心目的就是治“参数漂移”的毛病。我最早接触它是在一套光纤水听器项目上当时被光强扰动坑得够呛换成这个思路之后问题基本一次解决。这篇文章不绕弯子直接讲它的原理推导、实现步骤和我在调试中攒下的经验适合正在做PGC解调算法仿真、FPGA或DSP实时解调的工程师以及研究干涉型光纤传感方向的研究生参考。1. PGC解调为什么难从干涉信号到两路正交分量1.1 干涉仪输出一个被载波隐藏的相位先回到最基础的物理模型。干涉型光纤传感器比如水听器、地震计、光纤麦克风的输出光强可以写成I(t) A B·cos(C·cos(ω0·t) φ(t))这里的A是直流光强项B是干涉条纹的对比度也就是有效光强幅度C是相位调制深度由PZT调制器或者激光器波长调制引入ω0是载波角频率φ(t)才是我们真正想测的待测相位信号。PGC的基本思想是把φ(t)“搬”到高频载波的边带上然后用混频和低通滤波把它解调出来。对上述式子做贝塞尔函数展开会得到一串以ω0为基频的谐波组合。其中比较重要的是载波基频cos(ω0·t)这一项的幅度正比于J1(C)·sin(φ)载波二倍频cos(2ω0·t)这一项的幅度正比于J2(C)·cos(φ)。这就是整个PGC解调能成立的数学基础——待测相位信息同时出现在了基频和二倍频两个“通道”里而且这两个通道天然正交。为什么需要载波呢因为干涉仪的输出光强本身对φ(t)是余弦关系直接解调会遇到灵敏度为零的点也就是所谓“余弦响应退化”。加了高频载波之后φ(t)被调制到载波边带上解调时无论φ怎么变化都能通过两个正交分量恢复出来。这个思路听起来很清晰但真正实现时B、C这些非理想参数会让人非常头疼。1.2 混频和低通如何拿到sin项和cos项把I(t)分别乘以cos(ω0·t)和cos(2ω0·t)再用低通滤波器滤掉高频分量会得到两路基带信号。忽略系数中不影响原理的常数因子它们可以写成X K1·sin(φ(t))Y K2·cos(φ(t))其中K1正比于B·J1(C)K2正比于B·J2(C)。这里的J1、J2分别是第一类一阶、二阶贝塞尔函数。要注意K1和K2一般不相等除非C恰好取一个特殊值。这就是PGC解调的“分水岭”你已经拿到了两个携带相位信息的分量但它们不仅包含共同的B还分别掺入了不同的J1(C)和J2(C)。后续算法的任务就是从这个X和Y里干净地恢复出φ(t)同时尽量不让B和C的波动影响结果。1.3 DCM和Arctan的“老毛病”传统方案里最有名的是PGC-DCM微分交叉相乘和PGC-Arctan。先说DCM它的核心操作是S X·dY/dt - Y·dX/dt把X K1·sinφ和Y K2·cosφ代进去会得到S K1·K2·φ(t)看起来很好但问题在于K1·K2里面含有B²。一旦光强抖动比如光纤弯曲损耗变化、光源功率漂移B一变S的增益就跟着变。解调输出的幅度误差是实打实的不是滤个波就能消掉的。而且DCM需要同时做两路微分和三路乘法在FPGA上占的资源不算少。PGC-Arctan的思路则是直接用两个低通分量做反正切φ(t) arctan(X / Y)这里有两个隐患。第一如果直接用X/Y隐含假设了K1 K2但实际中C值稍微偏离2.63左右J1和J2就不相等反正切出来就是带谐波失真的相位。即使提前标定好比值C漂移之后误差又回来了。第二反正切函数在相位跨越±π/2时会有跳变后续需要额外的解包裹逻辑信噪比差的时候解包裹很容易出错一个跳变就是一大段毛刺。所以做PGC解调的人一直想要一种方案不显式依赖B对C漂移不敏感不要在相位边界处搞跳变。PGC-SDD-DSM正是往这个方向走的。2. 微分自乘相减把光强和调制深度“吃”掉2.1 先平方再微分两条新支路这个算法里最妙的一步不是直接去解φ而是先构造两个“派生信号”把B和C的影响单独拎出来。具体做法是把刚才低通滤波得到的X、Y分别平方然后再对时间求导得到两条新支路d(X²)/dt 2·K1²·sin(φ)·cos(φ)·φ(t)d(Y²)/dt -2·K2²·sin(φ)·cos(φ)·φ(t)注意这里用到了链式法则所以结果里多了一个φ(t)这一点在后面的“单路径微分相除”环节会非常有用。先记下这两个式子然后做加减法A_sum d(X²)/dt d(Y²)/dt 2·(K1² - K2²)·sin(φ)·cos(φ)·φ(t)A_diff d(X²)/dt - d(Y²)/dt 2·(K1² K2²)·sin(φ)·cos(φ)·φ(t)这两个式子放在一起会发现问题已经被“约”掉了一大批东西。B呢藏在K1、K2里但没有单独出现。φ本身呢变成了sin·cos·φ的组合。只要sin(φ)·cos(φ)不为零φ(t)也不为零两个式子里的动态项就可以被约掉。2.2 相除算出比值反推调制深度把A_sum和A_diff相除得到R A_sum / A_diff (K1² - K2²) / (K1² K2²)这一步很有意思分母里的sin·cos·φ被约掉了分子里剩下的只是一个与K1、K2比值有关的数。换句话说R不包含待测相位信息也不包含光强幅度B它只跟C值下J1(C)和J2(C)的比例有关。反解一下设k² K1² / K2²那么k² (1 R) / (1 - R)这里的k就是J1(C)/J2(C)的绝对值严格来说还带有符号后面会讲。于是就有了一个不需要预先标定C、不需要知道光强B就能实时估计出J1/J2比值的方法。这在工程上意义非常大你不再需要相信“C值出厂设好了就不会变”因为温度、应力、驱动电压漂移都会让C偏掉而算法可以自己在线感知这个偏差。需要注意实际计算R时分子分母都可能因为sin(φ)·cos(φ)接近零而变得很小直接除法会产生极大的毛刺。所以工程实现里一定要对R或者k²做平滑滤波最好加上状态保护比如在能量过低时保持上一拍的k值等有效信号来了再更新。2.3 这一步到底消掉了什么把这一阶段做的事情梳理一下输入是两路低通后的准正交信号X、Y输出是一个比值系数k。这个k的用处就是替代传统Arctan算法里“查表估计C值”的角色而且是实时更新的不需要人去现场重新标定。我举个实际例子。某次实验里光纤跳线被人无意中弯了一下光强B直接掉了30%。用传统DCM解调输出幅度立刻缩水波形幅度变得忽大忽小用带固定补偿系数的Arctan解调谐波失真明显增加。但用SDD-DSM的框架只要低通后的X、Y本身没有饱和、没有被噪声淹没R的估计几乎不受B变化影响后续解调出来的相位幅度基本稳得住。这就是“把光强和调制深度吃掉”的含义。当然这个环节不是没有代价。它需要两次平方、两次微分、一次除法计算量比DCM多了一大截在DSP上写可能还好在FPGA上就要认真规划流水线和资源。另外如果sin(φ)·cos(φ)长期接近零比如待测信号极小、相位在某个固定值附近基本不动那么R的估计会非常不稳定。实际使用时建议对k做分时段平滑别让它被几个噪声尖峰带跑。3. 单路径微分相除一路微分、一路相除积分还原相位3.1 一路信号的微分除以另一路拿到k之后接下来要恢复φ(t)。很多经典算法到这里会选择做反正切但SDD-DSM用的是“微分相除”。取X的导数dX/dt K1·cos(φ)·φ(t)然后把这个导数和Y相除V (dX/dt) / Y (K1·cos(φ)·φ(t)) / (K2·cos(φ)) k·φ(t)你看cos(φ)也约掉了只剩一个与相位导数成正比的量V前面的系数正好是上一步估计出来的k。这个操作只需要对一路信号做微分另一路直接用来做除法所以叫“单路径”微分相除相比DCM的两路同时微分再交叉相乘结构和计算路径都简单一些。如果Y的表达式中cos(φ)恰好过零除法会出问题这个我后面专门讲。但在正常工作点附近这个除法输出就是一个和φ成线性关系的干净信号。3.2 为什么能直接积分还原既然拿到了φ(t)剩下的就顺理成章了φ(t) V / k对时间做一次积分就能恢复出φ(t)φ(t) ∫(V/k) dt积分之后通常会带一个直流漂移分量这其实是积分器的初始常数和低通残余直流造成的。实际解调系统里待测信号一般是交流信号振动、声压变化或者我们关心的是相对变化量所以积分后加一个高通滤波器把直流和超低频漂移去掉即可。这个过程与传统Arctan相比最大的优势是整个恢复过程是连续的、单值的没有±π/2边界不需要解包裹。用Arctan方案时相位一旦接近±π/2反正切函数就会跳到另一个分支解包裹逻辑要小心翼翼而SDD-DSM这种“先求导再积分”的思路天然绕开了相位跳变问题。和DCM相比它的输出又经过k修正对C漂移不那么敏感。3.3 过零“陷阱”和符号问题工程上永远绕不开边界情况。第一个问题是Y K2·cos(φ)过零时V (dX/dt)/Y会变成无穷大。实际信号里φ通常是一个围绕某点波动的量比如φ(t) A·sin(ωt)那么cos(φ)确实会周期性地接近零。处理办法有几个在分母里加一个很小的常数ε牺牲一点精度避免除零。对V做滑窗中值滤波把尖峰去掉。在|Y|低于某个门限时改用另一路的组合比如W (dY/dt)/X -φ/k两条路径互为补充总有一路在正常工作区附近。这就是实践里常说的“双路径互相备份”思想。第二个问题是符号问题。前面提到k² K1²/K2²开方后默认k为正但J1(C)和J2(C)并不总是同号。当C小范围变化时比如在2到3.5之间两者基本都是正数问题不大但如果C漂移过大跨过了某个贝塞尔函数零点J1可能变号k的符号就会反解调结果会整体反相相位看起来像翻了个个儿。这是所有依赖贝塞尔函数比值算法都躲不过的坑不是SDD-DSM特有。工程上最简单的做法是把调制深度C控制在2到3.5附近不要让它跑到贝塞尔函数零点附近。4. 仿真与实现从公式到可跑通的代码4.1 仿真参数怎么定纸上谈兵没用直接看仿真里怎么布参数。我用Python做过一组比较典型的验证参数如下参数取值说明载波频率 f010 kHzPGC载波频率决定边带位置采样率 fs1 MHz过采样方便低通滤波和微分待测信号频率 fsig1 kHz典型声学/振动信号频段待测信号幅度 A_phi1.2 rad相位幅度超过1 rad才能体现非线性问题调制深度 C1.0 ~ 4.0测试中做扫掠观察C漂移影响光强B1.0 0.3·sin(2π·50t)加入50 Hz光强扰动模拟电源/光纤损耗波动噪声高斯白噪声SNR30 dB模拟探测器噪声这里故意让C从1到4变化光强加了50 Hz扰动就是要看算法在“不理想环境”下还能不能抗住。实际项目中C往往会因为PZT驱动电压漂移而变动光强扰动更是家常便饭所以仿真时别把参数调得太“干净”。4.2 混频低通与微分器的工程实现代码之前先说说工程实现上的几个关键点这些直接决定算法能不能跑稳。混频这一步要生成参考信号cos(ω0·t)和cos(2ω0·t)。注意参考信号和载波之间最好保持相位锁定如果有相位偏差解调出来的X、Y会串扰。工程上用锁相环或直接与调制信号源同步产生参考信号比较稳妥。低通滤波器的设计要重点注意群延迟。X和Y两路必须用完全相同的滤波器保证延迟一致否则后面做除法和比值时会引入固定的相位误差。我在项目里用的是FIR等纹波滤波器阶数取128通带截止频率设为待测信号最大频率的1.2倍左右阻带衰减做到60 dB以上。实测下来FIR的线性相位特性对这类算法特别友好不像IIR滤波器那样群延迟随频率抖。微分器看似简单实际是噪声放大器。理想微分器的幅频响应随频率线性上升高频噪声会被成倍放大。工程上建议用中心差分或者设计一个带限微分器。中心差分本质上也是一个高通特性所以微分之前最好再对X、Y做一次轻滤波。我一般会在微分环节后面加一个10阶左右的FIR低通专门压高频噪声。4.3 Python代码示例下面是一个从生成干涉信号到解调相位的Python示例核心步骤都做了注释。这个代码不是完整工程版本但拿来跑通算法流程、理解数据流足够了。import numpy as np from scipy.signal import butter, lfilter, filtfilt fs 1_000_000 t np.arange(0, 0.05, 1/fs) f0 10_000 f_sig 1_000 C 2.63 B 1.0 0.3*np.sin(2*np.pi*50*t) phi 1.2*np.sin(2*np.pi*f_sig*t) # 干涉信号 I B * (1 np.cos(C*np.cos(2*np.pi*f0*t) phi)) # 混频 mix1 I * np.cos(2*np.pi*f0*t) mix2 I * np.cos(2*np.pi*2*f0*t) # 低通滤波两路使用相同滤波器 b, a butter(4, 20000/(fs/2), btypelow) X filtfilt(b, a, mix1) Y filtfilt(b, a, mix2) # 消除常数增益差异简化处理 X X / np.std(X) Y Y / np.std(Y) # 微分 dX np.gradient(X, t) dY np.gradient(Y, t) # 微分自乘相减估计 k d_X2 np.gradient(X**2, t) d_Y2 np.gradient(Y**2, t) A_sum d_X2 d_Y2 A_diff d_X2 - d_Y2 # 移动平均平滑防止尖峰 def moving_average(x, w50): return np.convolve(x, np.ones(w)/w, modesame) R moving_average(A_sum / (A_diff 1e-12)) R np.clip(R, -0.99, 0.99) k2 (1 R) / (1 - R) k np.sqrt(np.clip(k2, 0.01, None)) # 单路径微分相除 V dX / (Y 1e-12) # 还原相位导数并积分 phi_dot V / k phi_est np.cumtrapz(phi_dot, t, initial0) # 高通滤除积分漂移 b_high, a_high butter(2, 100/(fs/2), btypehigh) phi_est filtfilt(b_high, a_high, phi_est)这段代码里我故意把低通滤波改用了filtfilt保证零相位失真。实际实时系统没法用零相位滤波要改用因果FIR并忍受一截固定延迟。如果你在MATLAB里写思路完全一样把np.gradient换成diff加适当对齐就行。4.4 结果怎么看、如何评估仿真好不好不能只看波形像不像。我一般看三个指标相位误差把解调出来的phi_est和原始phi做差去掉直流后看标准差。理想情况下应该小于0.01 rad。总谐波失真THD对解调输出做FFT看二次、三次谐波相对于基波的幅度。传统Arctan在C偏离2.63时THD会明显上涨而SDD-DSM在C从2.2到3.2变化时THD应该保持很低。抗光强扰动能力解调输出幅度是否随B的50 Hz波动而变化。好的算法输出应该基本是一条干净的等幅包络。我实测下来SDD-DSM在光强扰动抑制上比DCM强很多THD表现也优于Arctan。代价是计算量上来了在普通PC上仿真毫无压力但要移植到FPGA时就得好好设计流水线了。5. 常见问题与调试技巧实录5.1 分母太小尖峰爆炸这是单路径微分相除最典型的坑。Y接近零时V dX/Y会产生巨大的尖峰积分之后就是一大段错误波形。我的解决办法是组合拳第一除法之前给分母加一个跟信号幅度相关的自适应ε第二对V用滑窗中值滤波窗口取5到7个采样点就行第三如果信号本身信噪比好还可以加一个|Y|门限低于门限时保持上一拍的V输出。这三个手段叠加尖峰基本能压到可以接受的程度。5.2 微分噪声放大微分运算是PGC这类“微分-积分”结构绕不开的痛点。采样率越高微分噪声越容易失控。我的经验是在微分之前先把X、Y做一次轻度的低通平滑不要把高频噪声喂给微分器。另外中心差分本身只有两阶精度如果噪声还是大可以用Savitzky-Golay滤波器来做数值微分平滑和微分同时完成效果会好很多。不过SG滤波器的窗口长度不能太长否则会把信号本身的快速变化也磨平了。5.3 积分漂移phi ∫(V/k)dt之后直流漂移几乎是必然的。原因有三低通滤波后的X、Y本身可能残留直流偏置除法输出V微小的直流分量会被积分累积数值积分的初始常数也没法预知。最简单有效的处理是加一个高通滤波器把几十赫兹以下的分量切掉。要注意高通截止频率不能高于待测信号最低频率否则信号本身会被衰减得很难看。5.4 C漂移过大导致符号反相前面提到了k的符号问题。如果你把C从2.63一路调到4.0附近J1(C)会过零点J1/J2变成负数。此时k的估算如果只取正根解调结果的相位会整体反相听起来像“声音倒放”一样振动方向完全反了。排查时如果发现解调结果在某个调制电压之后突然反相不要怀疑滤波器和微分器先去看贝塞尔函数值。实用建议是加一个最小调制深度监测用RMS值跟踪载波谐波幅度如果发现C漂移过大提醒现场调驱动电压或者用符号校验逻辑自动翻转。5.5 调试中的经验顺序最后分享一点调试顺序能省不少时间。第一步先用干净信号跑通流程不掺噪声、不加光强扰动确认公式和代码方向没错。第二步加入光强扰动观察解调输出幅度有没有明显起伏验证DSM部分的功效。第三步扫C值看THD和符号确认k估计是否跟得上。第四步加入噪声和实际信号特征慢慢拧参数。别一上来就加噪声否则哪里出了问题你都分不清是公式错还是滤波没调好。这套算法用下来我最直观的感受是它把PGC解调里最烦的两个“环境变量”——光强B和调制深度C——变成了可在线估计、可补偿的量而不是让系统去“赌”它们不变。再加上它天然绕开了反正切的相位跳变问题做实时解调时省掉了解包裹的很多麻烦。如果你手头的系统光路环境不稳定、现场又不好频繁标定PGSDD-DSM绝对值得认真试一次。