企业官网建设流程全解析

做干涉型光纤传感解调的人，多半都经历过这种尴尬：光强抖一下、调制深度偏一点，辛辛苦苦写好的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做过一组比较典型的验证，参数如下：

这里故意让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), btype='low') 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, w=50): return np.convolve(x, np.ones(w)/w, mode='same') 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, initial=0) # 高通滤除积分漂移 b_high, a_high = butter(2, 100/(fs/2), btype='high') 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绝对值得认真试一次。