简介面向光纤传感与分布式测量领域的科研人员和工程师该压缩包提供了一套基于OFDR光学频率域反射技术的MATLAB仿真与数据处理示例可帮助理解高分辨率分布式光纤传感的实现思路。包内共5个源码文件均为m脚本总大小仅2KB但已覆盖核心算法环节啁啾脉冲生成与回波解调、温度变化引起的频移分析、空间分辨率计算以及测量距离与精度处理基本串联起OFDR系统从信号激励到参数提取的完整流程。目前已有1877人学习反映出该主题在电力电缆热监测、桥梁健康监测、油气管道泄漏检测等应用场景中的关注热度。通过学习这些示例代码读者能够快速掌握OFDR的关键数学模型与编程方法也可为后续基于MATLAB或LabVIEW开展二次开发、算法验证和教学演示提供直接参考。该资源适合作为课程设计、科研预研或工程应用的入门资料。1. OFDR 光纤传感的数据链路为什么绕不开 Matlab 和 LabVIEW一套 OFDR 系统跑起来之后采集卡灌进来的原始数据是几 GB 的干涉波形。要把这些波形变成沿光纤分布的应变和温度曲线中间隔着一条从光域到数字域的换算链路等光频重采样、FFT、相位解调、标定系数修正每一步都有参数要调、有坑要躲。实际项目里最常见的分工方式是 LabVIEW 负责采集控制、触发同步和数据落盘Matlab 负责离线重采样、频谱分析和相位解缠。两者接力不是习惯问题而是 LabVIEW 的硬件定时和 DMA 缓冲管理比 Matlab 脚本可控得多Matlab 的矩阵运算和绘图生态又比 LabVIEW 的 G 语言高效得多。这篇把这条链路上关键的原理、参数和排错方法摊开适合正在搭 OFDR 系统、或者准备接手光纤传感数据处理的人参考。2. OFDR 测量原理与信号链路的三个关键环节2.1 拍频位置映射扫频范围决定分辨率扫频速率决定测量速度OFDR 的基本结构并不复杂。可调谐激光器输出线性扫频光经耦合器分成参考臂和测量臂测量臂的后向瑞利散射光回到探测器与参考光干涉。由于散射点位置不同参考光与散射光之间存在时延 τ而激光器扫频使得不同时延对应不同瞬时频率差探测器上就出现了拍频信号。拍频频率与散射点位置的关系是 OFDR 一切数据处理的基础fb γ × τ τ 2nz / c其中 γ 是激光器扫频速率Hz/sn 是光纤有效折射率z 是散射点距参考位置的距离。频率轴可以直接映射到距离轴。距离分辨率由扫频范围 ΔF 决定Δz c / (2n × ΔF)这个公式和 OTDR 完全不同。OTDR 分辨率取决于脉冲宽度增大测量距离就要降低分辨率OFDR 的分辨率只跟扫频范围挂钩扫频范围做到几十纳米时分辨率能到微米级这是 OFDR 在分布式光纤传感里的核心价值。系统设计时几个关键参数会互相牵制扫频速率越高测量速度越快但同一位置的拍频频率也越高对采集卡采样率和模拟带宽的要求随之上升扫频范围越大分辨率越好但激光器调谐线性度和波长重复性压力也越大。典型配置如下表。参数典型值对测量的影响扫频范围 ΔF100 GHz ~ 20 THz决定距离分辨率范围越大分辨率越高扫频速率 γ1 ~ 100 nm/s决定单次扫频耗时和拍频上限辅助干涉仪延迟 L_MZI50 ~ 200 m决定重采样频率间隔需与扫频速度匹配激光器线宽1 ~ 100 kHz决定相干长度影响最大测量距离采集卡采样率 fs5 倍以上拍频上限采样不足直接产生频谱混叠计算出拍频上限后采样率宁高勿低。实际调试中我一般先按理论值的 5 倍设置采样率再观察距离域曲线远端噪声基底是否随幅度抬升以此判断是否存在混叠。2.2 辅助干涉仪与等光频重采样抑制扫频非线性的标准做法半导体可调谐激光器的扫频过程并不是理想的线性调谐瞬时扫频速率 γ 会随波长和时间漂移。如果直接按等时间间隔采样做 FFT相同的散射点位置会因为 γ 的变化产生不同的拍频频率距离域峰值会展宽严重时几十米外的峰直接看不见。行业里的标准做法是用辅助干涉仪做等光频重采样。辅助干涉仪本质上是一个马赫-曾德尔干涉仪其中一臂有几十到两百米的延迟光纤。激光器扫频时辅助干涉仪输出一个正弦干涉信号其过零点对应光频变化量恒定的时刻。光频间隔与辅助干涉仪延迟的关系为Δf_mzi c / (n × L_MZI)假设 L_MZI 为 100 mn 取 1.468则每个过零对应的光频变化约为 2.04 MHz。以这些过零时刻为采样点对主干涉仪信号重新插值就得到了光频域等间隔采样的序列。这个序列再去做 FFT扫频非线性带来的展宽基本被消除。重采样后的等效采样率不再是采集卡的物理采样率而是由辅助干涉仪决定的“光频采样率”。计算距离轴时不能直接用物理采样率换算否则距离标定会整体偏移。这一点是很多刚接触 OFDR 数据处理的人最容易忽略的。2.3 从复数散射幅度到相位信号OFDR 解调的基本操作对象重采样和 FFT 之后得到的是复数频谱。每个距离点上的幅度代表该处瑞利散射的强度相位则携带了光纤轴向应变和温度的信息。这里的关键认知是OFDR 分布式传感解调的不是某个单点反射峰的相位而是整段光纤的瑞利散射图案。两次扫描之间如果光纤某一段发生应变该段散射信号的相位会整体平移相位变化量与该段的应变近似成正比。解调思路就是沿距离方向取一个滑动窗口对参考信号和当前信号的复数频谱做互相关用互相关峰值的相位差作为该位置的测量值。这个思路在数据实现上要处理两个问题一是相干衰落导致的相位不确定性二是窗口长度与空间分辨率的取舍。窗口越长互相关信噪比越好但空间分辨率越差。工程上一般先根据需要的空间分辨率确定窗口长度再做参数扫描找到相位噪声的拐点。3. 用 Matlab 把 OFDR 原始数据变成距离域散射曲线核心算法与参数3.1 数据预处理去直流、归一化、裁剪扫频首尾区段采集卡记录的数据通常包含两到三个通道主干涉仪信号、辅助干涉仪 MZI 信号有时还有激光器功率监测通道。预处理的第一步是确认每个通道的物理含义尤其是 MZI 通道直流偏置是否稳定。主干涉仪信号的直流成分来自参考光的恒定功率直接影响后续过零点检测和 FFT 的零频分量。处理时先减去均值再除以标准差做归一化。归一化不是必须的但当激光器功率在扫频过程中有缓慢漂移时归一化能减少能量泄漏带来的旁瓣抬升。激光器扫频的起始段和结尾段通常存在调谐非线性剧烈变化即使做重采样也难以完全修正。我的做法是在重采样前先裁剪掉首尾各 2% 到 5% 的数据点虽然损失一点测量距离但能显著改善距离域曲线的基底平坦度。裁剪比例需要根据激光器实际调谐曲线调整没有统一标准。function [sig_c, mzi_c] preprocess(sig_main, sig_mzi, fs) % 预处理去直流、归一化、裁剪首尾 % fs 为采集卡采样率裁剪比例为每端 3% n length(sig_main); cut round(n * 0.03); sig_main sig_main(cut1:end-cut); sig_mzi sig_mzi(cut1:end-cut); % 去直流并归一化幅度 sig_main sig_main - mean(sig_main); sig_main sig_main / std(sig_main); sig_mzi sig_mzi - mean(sig_mzi); sig_mzi sig_mzi / std(sig_mzi); sig_c sig_main; mzi_c sig_mzi; end裁剪比例放在函数开头集中管理方便对不同激光器做标定时快速调整。去直流的均值计算要用整段数据不要用滑动平均否则会引入低频畸变。3.2 基于 MZI 过零点的等光频重采样实现等光频重采样的核心是找到 MZI 信号的过零点时刻然后在过零点时刻对主干涉仪信号做插值。Matlab 里可以用diff(sign())找符号变化位置但直接取整数索引会引入最多一个采样周期的抖动因此需要用线性插值把过零点时刻细化到子采样精度。function [sig_res, t_zero] resample_by_mzi(sig_main, sig_mzi, fs) % 等光频重采样 % 返回重采样后的主干涉仪信号和过零点时刻序列 s sign(sig_mzi); idx find(diff(s) ~ 0); % 符号变化位置 t_zero zeros(length(idx), 1); for k 1:length(idx) i0 idx(k); t0 (i0-1)/fs; t1 i0/fs; y0 sig_mzi(i0); y1 sig_mzi(i01); t_zero(k) t0 - y0*(t1-t0)/(y1-y0); % 线性求根 end % 去掉首尾不完整的 MZI 周期 t_zero t_zero(2:end-1); t_orig (0:length(sig_main)-1)/fs; sig_res interp1(t_orig, sig_main, t_zero, linear); end这里interp1使用线性插值已经足够因为重采样后的等效光频采样率远高于信号带宽高阶插值带来的计算开销大于收益。过零点数量就是重采样后的点数如果激光器扫频范围是 20 THz、MZI 光频间隔是 2 MHz理论上得到约 1000 万个点。过零点检测对 MZI 信号质量很敏感。如果 MZI 信号幅值波动太大建议先做一次带通滤波中心频率设为 MZI 信号的拍频频率。滤波器的带宽要足够宽否则过零点位置会被滤波器的群延迟扭曲。3.3 FFT 加窗与距离轴标定窗函数参数和坐标换算重采样完成后对信号做 FFT 得到距离域的复数散射曲线。距离轴不能直接按物理采样率换算而是要用有效扫频范围和重采样点数来算N length(sig_res); % 使用凯泽窗控制频谱泄漏 beta 6; win kaiser(N, beta); sig_w sig_res .* win; S fft(sig_w); S S(1:floor(N/2)); % 距离轴标定总扫频范围对应总距离范围 df_total 2e13; % 有效扫频范围Hz需按激光器实际值修改 d_max c / (2 * n * df_total / N); % 单点距离间隔 z_axis (0:length(S)-1) * (c / (2 * n * df_total));凯泽窗的 β 值控制主瓣宽度和旁瓣抑制的折中。β 取 0 相当于矩形窗旁瓣高β 取 8 以上旁瓣极低但主瓣明显展宽。对 OFDR 距离域谱来说β6 左右通常是比较好的起点既能把距离域峰值的旁瓣压下去又不会让空间分辨率恶化太多。窗函数旁瓣抑制主瓣宽度适用场景Hamming-43 dB适中常规散射曲线显示Kaiser β6-50 dB 左右略宽弱散射点检测需要压低旁瓣Kaiser β10-70 dB 左右明显展宽强反射峰附近的弱信号分析距离轴的标定误差主要来自有效扫频范围 df_total 的取值。如果激光器实际扫频范围与设定值有偏差距离域峰的位置会整体偏移。标定方法是用一段已知长度的光纤端面反射峰去反推扫频范围这个校准系数在每次实验前都应该重新确认。3.4 相位解调与应变换算滑动互相关例程得到复数频谱后两次扫描之间的相位差通过滑动互相关来计算。直接对每个距离点做angle(S2) - angle(S1)的问题是散射信号幅度接近零的位置相位噪声极大。滑动窗口互相关可以抑制这种噪声窗口内所有散射点的贡献被平均。function strain_curve demod_strain(S_ref, S_cur, win_len, lambda, n, L_gauge) % 滑动互相关解调相位差换算应变 % win_len 窗口长度点数lambda 中心波长mL_gauge 窗口对应光纤长度m N length(S_ref); half floor(win_len/2); phase_diff zeros(N, 1); for z half1 : N-half seg_ref S_ref(z-half : zhalf); seg_cur S_cur(z-half : zhalf); % 互相关并取相位 cross sum(seg_ref .* conj(seg_cur)); phase_diff(z) angle(cross); end % 沿距离方向解缠避免 ±π 跳变 phase_diff unwrap(phase_diff); % 应变换算Δφ (4πnL/λ) × ε × 0.78 strain phase_diff * lambda / (4 * pi * n * L_gauge * 0.78); end0.78 是标准单模光纤有效光弹系数的近似值物理上接近 (1 - n²/2 × (p12 - ν(p11p12))) 的数值结果。不同纤芯掺杂的光纤会有差异精确标定需要用已知应变的试验段反推实际系数。窗口长度 win_len 的选择直接决定应变曲线的空间分辨率和噪声水平。窗口越短空间细节越好但相位噪声越大窗口越长曲线越平滑但会抹掉短距离内的应变集中区域。实用做法是先按设计空间分辨率换算窗口点数再做一倍和两倍窗口的对比看应变曲线波峰位置是否变化。4. 用 LabVIEW 搭 OFDR 实时采集链路别在数据落盘上翻车4.1 采集卡触发与同步MZI 信号该走硬件触发还是软件触发OFDR 采集的同步核心是让采集窗口与激光器扫频起始点对齐。一般做法是把激光器的触发输出接到采集卡的 Start Trigger上升沿到来时采集开始单次扫频记录整段原始波形。MZI 信号的处理有两种路线。一种是把 MZI 信号作为普通模拟通道与主干涉仪信号同步采样后续在 Matlab 里做软件重采样。这种做法的好处是硬件简单且保留原始信息重采样参数随时可调坏处是数据量大且对 MZI 通道和主通道的相位一致性要求高两条模拟通道的延迟差异必须校准。另一种路线是 MZI 信号经过比较器整形为 TTL接入采集卡的 PFI 引脚作为外部采样时钟。每个 MZI 过零点触发一次采样硬件直接输出等光频间隔的数据无需软件重采样。这个方案对采集卡要求高需要支持外部时钟或重触发模式而且一旦 MZI 信号幅值波动导致比较器误触发整帧数据作废。对实验室环境来说我更推荐软件重采样路线灵活性和排查便利性都更好。4.2 生产者-消费者架构正常扫频速度下不要直接在采集回调里写文件OFDR 一次扫频产生的数据量非常可观。假设采样率 250 MS/s、16 bit、双通道单次扫频 0.1 秒数据量就有 100 MB 左右。在 LabVIEW 里用 DAQmx 读回数据后直接写 TDMS 文件虽然单次操作耗时不大但界面刷新、波形显示和文件写入混在同一个循环里时任何磁盘抖动都会造成下一次读取超时。常见做法是生产者-消费者架构采集循环只负责把数据块推入队列消费者循环从队列取数据执行 TDMS 写入和 UI 更新。队列深度按 10 帧数据估算内存占用在 GB 级别时要用 DVR 或者直接分配预缓冲区避免队列拷贝导致内存翻倍。// LabVIEW 框图中的关键节点文字描述 // 1. DAQmx Create Virtual Channel (AI, 双通道) // 2. DAQmx Timing (Sample Clock, 连续采样模式, 采样率fs) // 3. DAQmx Start Task // 4. While 循环: // DAQmx Read (N 通道 N 采样, 超时5000ms) // → 生产者队列 Enqueue 元素 // 5. 消费者 While 循环: // Dequeue 元素 → TDMS 文件写入 → 前面板波形更新采集循环的超时设置很关键。超时太短会在磁盘忙时误报错误太长则可能积累延迟导致缓冲溢出。5000 ms 的读取超时配合队列消费者通常够用如果频繁超时优先排查磁盘写入速度而不是调整超时。4.3 与 Matlab 的数据衔接MATLAB Script 节点与 TDMS 属性传递LabVIEW 写 TDMS 文件时把每次扫频的系统参数同时写入属性。这样 Matlab 读取时不需要手动维护一份参数表避免参数错位的问题。% 在 LabVIEW 的 MATLAB Script 节点中调用 filePath D:\ofdr_data\scan_0042.tdms; info TDMS_getStructInfo(filePath); % 读取 TDMS 属性 fs info.Property(SampleRate); df_total info.Property(SweepRangeHz); n_fiber info.Property(RefractiveIndex); data TDMS_readChannel(filePath, AI0); % 读主干涉仪通道TDMS 读取函数在 File Exchange 上有多实现版本选一个带属性读取功能的即可。注意区分 TDMS 文件中的通道名和属性名通道名在写入时要固定属性名建议统一大写驼峰风格减少跨平台大小写问题。LabVIEW 调用 Matlab 有三种常见方式MATLAB Script 节点、MathScript 节点、编译 DLL。MATLAB Script 节点要求本机装有完整的 Matlab 许可证启动速度慢但调试直观适合开发阶段。MathScript 节点语法兼容性有限很多向量化函数不支持不推荐用于 OFDR 这类重算法场景。编译成 DLL 部署最干净不依赖 Matlab 图形界面但参数接口改一次就要重编译一次。调用方式依赖优点缺点MATLAB Script 节点本机完整 Matlab调试方便支持全部函数启动慢许可证绑定MathScript 节点LabVIEW 自带无额外依赖语法兼容差不支持部分高级函数MATLAB Compiler SDK DLLMCR 运行时发布方便无界面算法修改需重新编译5. 数据格式与离线重处理LabVIEW 记录、Matlab 分析的典型分工5.1 为什么实时采集段交给 LabVIEW而不是直接在 Matlab 里循环Matlab 的采集工具箱也能驱动部分采集卡但它的内存管理和实时调度不如 LabVIEW DAQmx 稳定。说得直白一点Matlab 脚本在采集循环里如果遇到一次绘图刷新或者 JIT 编译停顿缓冲溢出就来了。LabVIEW 的 DAQmx 驱动在内核层面做 DMA 环形缓冲管理连续高速采集时丢帧概率远低于脚本环境。另一个原因是工程现场通常不需要在采集的同时做完整解调。实时显示距离域散射曲线可以用简单的 FFT 预览但完整的滑动互相关解调计算量大放在线路上做会造成采集循环阻塞。最佳实践是采集与存储跑在 LabVIEW 里解调算法在 Matlab 里离线重处理。5.2 数据格式选择TDMS 还是自定义二进制TDMS 是 NI 的二进制流式格式兼具随机访问和属性存储能力Matlab 读取也有现成工具。对 OFDR 这种单次扫频百 MB 级的数据量TDMS 的属性机制特别方便每次扫频的激光器参数、时间戳、实验备注都跟着文件走不会丢上下文。自定义二进制格式的优势是体积更小、读取速度更快但需要自己管理数据字典。如果团队里同时有 LabVIEW 和 Matlab 两套环境我建议直接用 TDMS省下格式维护的成本。这里列一个 TDMS 写入时的属性设置清单照着填可以规避大部分对不上的问题属性名内容示例SampleRate采集卡物理采样率250000000SweepRangeHz有效扫频范围用于距离标定2e13RefractiveIndex光纤有效折射率1.468MZIDelay_m辅助干涉仪延迟光纤长度100ScanIndex扫频序号42Timestamp采集时刻2025-06-01 10:00:005.3 字节序与路径乱码跨环境数据处理的隐性坑LabVIEW 和 Matlab 都跑在 x86 平台上时二进制数据的字节序默认小端一般不会出问题。但如果以后迁移到嵌入式平台或使用国产化采集卡大端小端不一致会导致数据完全错误。自定义二进制格式时在文件头写入一个魔数和字节序标志是基本操作不要相信环境默认。中文路径是另一个高频坑。LabVIEW 安装在带中文的路径下或者 TDMS 文件存放在中文目录中Matlab 读取时可能出现无法找到文件或者音符乱码的情况。尤其是 Matlab 脚本文件本身如果用了中文注释在某些系统区域设置下会显示乱码把脚本统一存成 UTF-8并在 LabVIEW 的 MATLAB Script 节点中避免写入非 ASCII 字符的参数名。给一个用 Matlab 自定义二进制读取的例子这个格式与 TDMS 可以并存% 读取自定义二进制 OFDR 数据帧 fid fopen(scan_0042.bin, rb); magic fread(fid, 4, uint8); % 验证文件头 assert(isequal(magic, [79 70 68 82]), 文件头不匹配); byteOrder fread(fid, 1, uint8); % 1小端, 0大端 machineFormat ieee-le; % 或 ieee-be if byteOrder 0, machineFormat ieee-be; end nPoints fread(fid, 1, uint64, 0, machineFormat); data fread(fid, nPoints, int16, 0, machineFormat); fclose(fid);直接读 TDMS 固然简单但理解底层二进制布局对排查问题和制定备份策略有实际帮助。文件头做魔数校验能在分卷拷贝或磁盘损坏时第一时间发现数据异常。6. OFDR 排错进阶相干衰落、偏振衰落与 PZT 标定验证6.1 相干衰落导致相位跳变不是算法错是信号本身劣化距离域幅度谱上某些位置会出现幅度接近零的暗点这些暗点的相位完全随机滑动互相关解调在这里会产生突然的相位跳变。这种现象叫相干衰落根源是多点散射光随机干涉叠加幅度出现瑞利分布的低值。处理手段不是单纯加大窗口而是用移频方法在参考臂串入声光移频器让拍频信号载波频率偏移每次扫频的散射图案发生变化对多次解调结果做平均能有效降低暗点概率。6.2 偏振衰落的表现与验证手段散射光的偏振态沿光纤随机旋转在某一段可能与参考光偏振正交导致干涉对比度消失。表现是距离域曲线中一段持续的幅度凹陷且改变扰动光纤的弯曲状态时凹陷位置会移动。验证手段很简单用手轻轻扰动怀疑段的光纤如果幅度谱发生明显变化就说明偏振衰落占主导。系统层面解决需要偏振分集接收在接收端分三路正交偏振态独立解调再合成。6.3 用 PZT 标定整条链路的完整验证流程标定是最容易被跳过的环节但 OFDR 信号链上任何一环节的系数错误都会让应变数据整体偏差。简单的做法是把 1 到 2 米光纤紧密绕制在 PZT 压电陶瓷管上PZT 通以已知频率和幅值的正弦电压产生周期性应变。用解调出的相位信号与驱动信号对比就能检查整条链路是否有信号丢失或相位跳变。% PZT 标定后的结果验证对比驱动信号与解调相位 t (0:length(phase_diff)-1) / frame_rate; ref_signal A * sin(2 * pi * f_pzt * t); % 驱动参考 % 在 PZT 段取平均减小随机噪声 seg_phase mean(phase_diff(z_start:z_end)); plot(t, seg_phase, t, ref_signal); % 计算相关系数低于 0.9 说明链路存在明显噪声或非线性 r corrcoef(seg_phase, ref_signal);相关系数低于 0.9 时优先检查是否存在明显的相干衰落暗点落在 PZT 段内。如果暗点无法避开把 PZT 绕制的光纤段移到距离域谱上幅度较高的位置或延长绕制长度以覆盖更多独立散射点通常就能改善。PZT 标定这一步跑通了后续换激光器、换采集卡再也不用从头摸参数。本文还有配套的精品资源点击获取