基于多峰高斯拟合的分布式拉曼测温系统设计与应用

张 均，杜 超，于昌智，张 丽，邓 霄,2

（1. 太原理工大学 物理与光电工程学院，太原 030024；2. 太原理工大学 新型传感器与智能控制教育部与山西省重点实验室，太原 030024）

引言

冰凌是冬季高寒地区一种普遍存在的自然现象，常发生在河水从低纬度流向高纬度的河段中[1]。冰凌阻塞河道，上游水位上涨，导致地势低的地区造成淹没，严重时甚至引起大坝决堤，给人民的生命财产安全造成重大损失[2]。水库的调度是防治冰凌灾害的最有效方法之一，但由于冬季气温低下，库区会形成冰盖，冰盖的存在不仅会产生较大的冰压力，对大坝造成损害，而且阻碍了水体和空气的物质能量交换，使水体的溶解氧和水温发生改变，进而对库区生态系统产生不利的影响。

在水库冰情研究中，水库冰盖内部和冰盖下方河水的温度分布对分析冰水热力学模型、库区冰情、冰水温度时空分布提供科学依据[3]。目前，用于冰盖内部的温度分布测量主要采用多个热敏电阻或铂电阻组成的温度传感链[4]。但对于大范围内的温度测量，制作和维护温度传感链相当复杂，所以将其用于冰盖下方河水垂直方向上的温度测量是不现实的。分布式光纤拉曼测温系统的光纤作为传输介质和敏感介质，利用光在光纤中的拉曼散射效应，可适用于长距离大范围的温度测量[5]，也可适用于小尺度的温度测量[6]。

空间分辨率为分布式光纤拉曼测温系统的主要性能指标之一，是指能够准确测量光纤各点温度的最小距离。当温感区域小于空间分辨率时，会使得测量的温度不准确，张磊等提出了一种反卷积校正算法将空间分辨提高了4 倍，修正了空间分辨率不足导致的不准确温度值[7]。Bazzo 等提出总变换反卷积算法，能正确测量温度变化区域为15 cm 的温度值，相比于普遍空间分辨率为1 m 的系统，其空间分辨率提高了约6 倍[8]。石希等对基于热红外遥感影像的河水温度反演方法进行比较，以最佳反演方法获得的河水温度均方根误差为1.0 ℃～1.1 ℃，实现了1 m 的空间分辨率，校正温度的最大误差为2 ℃[9]。Tyler S W 等提出一种分布式拉曼测温系统，该系统可用于测量冰架以及冰架下方800 m 深的海洋垂直方向上的温度分布[10]。Kobs S 等根据分布式拉曼测温系统测得的冰架温度分布，计算了季节性融冰速率[11]。王玎睿等设计了一种具有垂直高分辨率的测温装置，可以精确测量冰盖剖面的温度分布[12]。但是以上方法均不能解决由于空间分辨率不足，导致多个温度值测量不准确的问题；在冬季水库中河道冰盖垂直方向、冰盖下方河水垂直方向上温度分布的测量方法也未提及。因此，使用现有的分布式光纤拉曼测温技术难以同时测量小尺度和大范围环境的温度分布。

本文针对以上问题提出一种基于高斯拟合的多点温度校正方案，可以自动识别多个温度突变点，并对多个温度测量不准确的突变数据进行校正，使系统的温度准确度提高。同时，将分布式拉曼测温系统应用在黄河万家寨水库中测量冰盖内部及冰盖下方河水垂直方向上的温度分布，根据温度分布对万家寨水库中的冰盖冻结过程、冰盖厚度及冰盖下方河水垂直方向上的温度分布进行了分析，证明了系统应用在野外大范围监测河道流水剖面的温度分布的可行性。

1 系统总体设计及温度解调原理

1.1 系统总体设计

分布式拉曼测温系统总体设计如图1 所示。波长为1 550 nm 的光从光源中发出，经波分复用器进入到传感光纤；光在传感光纤中发生拉曼散射，散射回波长为1 450 nm 的拉曼反斯托克斯光和1 660 nm 的拉曼斯托克斯光进入波分复用器；波分复用器将两束光按波长分离出来，传输给光电探测器；光电探测器将两束光的光信号转换为模拟电信号传递给数据采集卡；数据采集卡接收到模拟电信号并将其转换为数字电信号，传递给上位机；上位机对表示斯托克斯光强值和反斯托克斯光强值的数字电信号进行解调，转换为沿光纤的温度分布，并在上位机上显示。

图1 分布式拉曼测温系统总体设计Fig. 1 Overall design of distributed Raman thermometry system

1.2 温度解调原理

通过上述方法，上位机获得两路光强值。因此，需将光强值解调为温度值。在选择解调方法时，考虑到光源的温度会影响入射光的稳定性，从而造成测得的温度分布产生较大误差，为了减小光源所处位置的温度变化对入射光功率的影响，选用双路解调的方式对光强值进行解调[13-14]，进而获得沿光纤的温度分布。其中，双路解调公式如下：

式中：T为沿光纤各个采集点的温度值，单位为K；T0为参考光纤环的温度值，单位为K；k是波尔兹曼常量，单位为J/K；h是普朗克（Planck）常量，单位为J·s；Δv是光纤分子的声子频率，单位为Hz；R为沿光纤采集到的各个点的拉曼比值；R0为参考光纤环内各个点的拉曼比值的均值。

2 温度校正实验及实验结果

2.1 多点温度校正方法

多点温度校正方法可以校正温感区域小于空间分辨率时所导致的不准确温度值，从而提高分布式拉曼测温系统的温度准确度。由双路解调方式获得沿传感光纤的温度分布后，使用极大值和极小值的方式寻找温度分布的所有峰值；剔除掉因噪声的存在使温度波动，被极大值极小值方式误判为由温度变化引起的峰值点，并判定剔除后温度分布的波峰值个数；对温度分布进行平移和翻转处理，并根据波峰值的个数自动调整拟合时所使用的数学方法，获得拟合后的曲线；计算拟合后的曲线所有峰值的半高全宽（full width at half maxima，FWHM），设定半高全宽的阈值范围，判定所有半高全宽所属的阈值范围，使用相应阈值范围内的校正公式对所有不准确的温度值进行校正，获得校正后的温度分布。

1） 高斯拟合

在校正温度之前，需要对温度分布进行高斯拟合，以确定所有温度突变点，并计算对应的半高全宽。当有一个或多个温度段发生变化时，采用极大值、极小值的方式来寻找突变点。但是在实验中，由于噪声存在，温度数据可能存在波动，因此小的波动也会被认为是突变点，所以需要将找寻到的所有极大值极小值与常温数据作差值，差值大于一定阈值的数据点才认为是温度突变点。本文设定的阈值为5 ℃（根据分布式拉曼测温系统精度设定），记录这些温度与常温值差值大于5 ℃点的个数N（即突变点个数），选择对应多峰高斯拟合的函数[15]：

式中：f1(x)为单峰高斯拟合曲线；a1，b1，S1分别为单峰高斯曲线峰值、峰值位置和半宽度信息。

式中：f2(x)为单峰高斯拟合曲线；a1，b1，S1分别为双峰高斯曲线第一个波峰的峰值、峰值位置和半宽度信息；a2，b2，S2分别为双峰高斯曲线第二个波峰的峰值、峰值位置和半宽度信息。

式中：fn(x)为多峰高斯拟合曲线；a1，b1，S1分别为多峰高斯曲线第一个波峰的峰值、峰值位置和半宽度信息；a2，b2，S2分别为多峰高斯曲线第二个波峰的峰值、峰值位置和半宽度信息；an，bn，Sn分别为多峰高斯曲线第n个波峰的峰值、峰值位置和半宽度信息。

（2）式、（3）式和（4）式是分别针对传感光纤温度分布存在1 个突变点，2 个突变点和n个突变点的高斯拟合函数。根据高斯拟合函数，可以计算出每个波峰的半高全宽：

式中：FWHM 是高斯曲线中各个波峰的半高全宽；Sj是高斯曲线中各个波峰的半宽度信息。

2） 校正公式

在校正温度之前，需要确定校正公式，将不准确温度值带入公式中，可获得校正后的温度值。用于确定校正公式的光纤总长为150 m，绕制10 m光纤环放置在恒温环境中用作参考光纤，在其后再分别绕制2 m，1.5 m 光纤环放在恒温槽中，同时将温度精度为0.01 ℃的标准温度计探头放入，准确测量恒温槽的温度值。恒温槽的温度分别设置为：-50 ℃、-40 ℃、-30 ℃、-20 ℃、-10 ℃。测得标准温度与光纤测量温度之间的线性关系如图2所示。

图2 标准温度与光纤测量温度之间的线性关系Fig. 2 Linear relationship between standard temperature and fiber-measured temperature

2 m 和1.5 m 光纤环根据3δ准则设定半高全宽的阈值范围分别为[2 m，3 m]、[1 m，2 m][16]。所以由图2 可得，低温时2 m 和1.5 m 光纤环测得的温度值与标准温度之间的线性关系如下：

式中：Ts为标准温度计测量的温度值；Tm为光纤测得的温度值。

因此，无论2 m 长度的温感区域还是1.5 m 长的温感区域，只要得到光纤测得的偏离真实值的温度值和对应的半高全宽，带入（6）式和（7）式便可以获得校正后准确的温度数据。

2.2 实验过程及结果

为了验证多点温度校正方法的正确性，将上述实验光纤放入设定值为-40 ℃的恒温槽（BILONW-506S）中。需要指出的是，本文使用的恒温槽温度均匀性可达到±0.1℃，低温时使用无水乙醇作为恒温浴场中的介质，可以稳定-40℃～0℃范围内的温度。为确保恒温槽的温度确实稳定在±0.1℃并测量其准确温度值，使用精度为0.01℃的二等标准铂电阻Fluke 5609 作为标准温度计与光纤一同放入恒温槽中。待标准温度计测得的温度稳定在±0.05℃内时，再采集沿光纤的温度分布，以保证低温情况下的温度稳定性。标准温度计测得恒温槽中的准确温度值为-40.39 ℃，根据解调得到的温度分布，可得2 m光纤环测得对应的温度值为-35.78 ℃，1.5 m 光纤环测得对应的温度值为-20.45 ℃。由上述可得，温感区域为2 m 时，测得的温度偏差为4.61 ℃；温感区域为1.5 m 时，测得的温度偏差为24.55 ℃。因此，需要对2 m 或1.5 m 光纤环测得的温度数据进行校正。

首先，获得沿光纤的温度分布，对其进行求导，求得所有的极小值点，将极小值点对应的极小值与其他未发生突变位置的温度点作差，并求绝对值，记下差值大于5 ℃极小值点的个数为2；其次，将未处理的温度分布进行平移和翻转，由极小值点的个数判定共含有突变点数为2，所以选择双峰高斯拟合公式进行拟合，上述步骤处理前后的效果如图3 所示。

图3 高斯拟合处理前后的温度分布Fig. 3 Temperature distribution before and after Gaussian fitting

最后，拟合得到的2 m 和1.5 m 光纤环对应波峰的半高全宽由（5）式计算约为2.1 m 和1.5 m，分别属于半高全宽阈值范围[2 m，3 m]和[1 m，2 m]。将光纤环对应的波峰值带入（6）式和（7）式，得到2 m和1.5 m 光纤环校正后的温度值分别为-40.70 ℃和-40.89 ℃。校正后的温度值能反映光纤所处环境的温度，温度准确度在1 ℃以内，校正前后的温度分布如图4 所示。

图4 校正前后的温度分布Fig. 4 Temperature distribution before and after correction

3 黄河万家寨现场实验方案及实验结果

3.1 现场实验方案

1） 冰盖剖面温度传感器设计

往年数据显示，黄河万家寨冰盖厚度约为40 cm左右，由于此分布式拉曼测温系统空间分辨率约为2 m 左右，无法直接测量冰盖剖面的温度分布，因此需设计光纤传感器去测量冰盖剖面的温度分布。为了解决传感光纤在布设时易被折断、低温容易对光纤造成破坏的问题，选取带有不锈钢管保护层且可适用于-15℃的铠装光纤作为本系统现场应用时的传感光纤。将直径为3 mm 的多模铠装光纤缠绕在直径为75 mm 的PPR 管上，总共缠绕80 cm，计算得到2 m 空间分辨率对应的垂直分辨率约为2.5 cm[17]，可用于测量冰盖剖面的温度分布。在缠绕好的光纤传感器外部涂覆环氧树脂，避免测量冰盖时静冰压力挤压缠绕光纤使光纤散开。按照上述方式设计的冰盖剖面温度传感器如图5 所示。

图5 冰盖剖面温度传感器Fig. 5 Temperature sensor of ice sheet profile

2） 传感器的布设与数据采集

2022 年1 月11 日至2022 年1 月18 日，在山西省忻州市黄河万家寨水库完成冰盖剖面温度传感器和河水剖面温度传感器的布设。布设过程如下：首先使用电冰钻对水库冰盖进行开孔，孔径为20 cm，大于中小尺度光纤传感器的直径；其次在传感器尾端的接口处固定重物，重物带动水体传感器部分下垂至库底，测量水体的温度分布；将垂直高分辨率传感器固定在冰孔中，测量冰盖冻结过程冰盖内部温度变化情况；最后将分布式光纤拉曼测温系统样机放置在冰盖上方，连接传感器测量温度分布。传感器的布设及数据采集如图6 所示。

图6 传感器的布设与数据采集Fig. 6 Sensor layout and data collection

3.2 现场实验结果

1） 冰盖垂直方向上温度分布现场实验结果

为了分析冰盖在测量过程中每天的变化情况，分别选取正在冻结（1 月14 日和1 月15 日）和冻结完成（1 月16 日和1 月17 日）的4 天温度数据作分析，温度数据如图7 所示。其中，h为冰盖剖面温度传感器深度，h=0 处为冰盖剖面温度传感器缠绕光纤第1 个点的位置，h=80 cm 处为最后1 个点的位置，新开孔的冰面不再与原始冰面齐平，记录第1 个点的位置与新冰盖上表面垂直方向的距离为7.47 cm。从图7（a）中可以看出，h=-43.16 cm 为1 月14 日每组数据的拐点所在位置；从图7（b）中可以看出，h=-47.31 cm 为1 月15 日每组数据的拐点所在位置；从图7（c）和图7（d）中可以看出，h=49.80 cm为1 月16 日和1 月17 日每组数据的拐点所在位置，拐点位置即为冰盖与河水的分界面。1 月13 日至1 月14 日与1 月14 至1 月15 日拐点位置变化较大，即冰盖厚度变化较大；1 月15 日至1 月16 日拐点位置变化较小，即冰盖厚度变化较小；1 月16 日至1 月17 日拐点位置基本不变，即冰盖厚度未变化。从图7 中每一天的温度数据也可以发现，冰盖内部温度数据随时间变化，这是因为冰盖内部温度与太阳辐射和空气温度有关。

图7 正在冻结和冻结完成的温度数据Fig. 7 Freezing and frozen temperature data

为了分析冰盖的厚度变化过程，将1 月11 日至1 月18 日的温度数据放在一起，如图8（a）所示，从图中数据计算得1 月11 日至1 月18 日冰盖厚度分别约为0，16.60 cm，29.05 cm，35.69 cm，39.84 cm，42.33 cm，42.33 cm，42.33 cm。综上所述，冰盖剖面温度传感器所在位置处在测量期间新生长的冰盖厚度约为42.33 cm，经5 个夜晚，新冰盖的下界面长到与原始冰盖下界面齐平的位置，冰盖8 天的生长趋势如图8（b）所示。

图8 冰盖垂直方向上的温度分布Fig. 8 Temperature distribution in vertical direction of ice sheet

为了验证上述冰盖剖面温度传感器测得的厚度是否正确，使用冰尺测得原始冰盖厚度约为48.4 cm，如图9（a）所示；又测量原始冰盖上界面与新冰盖上界面垂直方向上的距离约为4.4 cm，如图9（b）所示。计算得新生长的冰盖厚度约为44 cm，与冰盖剖面温度传感器测得的新生长冰盖厚度相差1.67 cm。

图9 人工测量的新冰盖厚度Fig. 9 New ice sheet thickness by manual measurement

2） 河水垂直方向上温度分布现场实验结果

为了分析河水垂直方向上温度分布的变化情况，分别选取1 月14 日空气与河水温度数据与8 天河水垂直方向上的温度数据作分析，如图10所示。其中，H为河水剖面温度传感器深度，记录H=0 处为空气与河水的界面位置，H=36 m 处为河底最后一个点的位置。从图10（a）中可以看出，H=0处为1 月14 日每组数据的拐点所在位置，其上下温度分布差异较大，说明H=0 处是空气与河水的分界点，这与上述记录的位置数据一致。从1 月14 日的温度数据还可以看出，在测量时间段内湖底温度数据比较稳定。2022 年1 月11 日至2022年1 月18 日8 天的温度数据如图10（b）所示。从图中可以看出，空气与河水的温度界限较明显，如黑色虚线所示。在黄河万家寨水库中，由于冰盖的存在阻碍了空气与和水的热交换，所以冰盖下的河水温度在8 天内变化不大，水温基本约在0.32～0.90 ℃内波动。

图10 河水垂直方向上的温度分布Fig. 10 Temperature distribution in vertical direction of river water

4 结论

针对温感区域小于空间分辨率引起的温度测量不准确问题，提出了多点校正方法，设计了温度校正实验，对温感区域为2 m 和1.5 m 的情况进行温度校正，校正后的温度测量准确度在1 ℃之内，证明了多点校正方法的正确性。将分布式拉曼测温系统应用在万家寨水库中冰盖及冰盖下方河水垂直方向上的温度测量，分析了冰盖的生长过程。数据显示，白天冰盖基本不生长，新冰盖经6 个夜晚的生长后与原始冰盖底部齐平，生长速度与夜晚空气温度有关。根据温度分布，计算出新生长的冰盖厚度约为42.33 cm，冰尺测量出新生长的冰盖厚度约为44 cm，光纤传感器测量冰厚的误差约为1.67 cm，证明了光纤拉曼测温系统用于冰盖厚度测量的可行性。将分布式拉曼测温系统用于测量河水垂直方向上的温度分布，发现由于冰盖的存在基本阻断了空气与河水之间的热交换，所以测量期间万家寨水库的温度基本稳定，水温约在0.32℃～0.90 ℃内波动。综上所述，分布式拉曼测温系统虽温度准确度为1 ℃，但可用于野外大范围内监测河道的温度分布。