CN116402696A - 基于非局部哈尔变换图像去噪光频域反射分布式传感方法 - Google Patents

Abstract

本发明涉及一种基于非局部哈尔变换图像去噪光频域反射分布式传感方法，包括下列步骤：在可调谐激光器的触发下，在高温炉升温过程中利用OFDR系统连续采集在不同温度下的瑞利散射信号，取相邻两组信号分别进行处理，得到局部参考瑞利散射光谱和局部测量瑞利散射光谱；去除直流分量；分别进行互相关计算，并对其结果进行归一化处理，每一组归一化互相关结果都对应于光纤上的一个位置，得到原始二维互相关强度图y；将与参考行欧氏距离最小的行组成矩阵；执行可分的二维哈尔变换，得到变换系数；双硬阈值收缩；获得初步去噪图；执行维纳滤波；得到对应于温度变化T₀～T_n的光谱频移分布。

Description

基于非局部哈尔变换图像去噪光频域反射分布式传感方法

技术领域

本发明涉及光纤分布式传感技术领域，尤其涉及一种利用非局部哈尔变换图像去噪光频域反射分布式传感方法。

背景技术

光纤因其具有抗电磁干扰、耐腐蚀、耐高低温等特点而得以用于恶劣环境中进行传感测量。在钢铁工业、固体氧化物燃料电池、核电站等领域，分布式光纤传感器可以通过嵌入、插入或直接接触的方式部署在测试结构上，实现准确、快速的分布式高温测量，光纤分布式高温传感展现出十分广阔的应用前景。在现有的光纤温度传感技术中，基于拉曼散射和布里渊散射的光纤分布式温度传感技术通常采用光时域反射测量方法，因此不适合用于要求高空间分辨率的场合；光纤布拉格光栅也被用于温度传感，并且能够达到厘米量级的空间分辨率，但是由于光纤光栅制作成本较高，并且在高于700℃时，会出现光栅擦除现象，限制了其在更高温度场景下的应用。基于瑞利散射的光频域反射(Optical FrequencyDomain Reflectometry,OFDR)方法是另一种光纤分布式高温传感方案。Patrick Bulot^[1]等人报道了利用ZrO₂掺杂光纤实现了800℃的测量，测试距离为40cm，空间分辨率为1cm；A.K.Sang^[2]等人利用镀金光纤进行分布式温度测量，最高可测量温度为850℃，测试距离为1m，空间分辨率为厘米量级。然而，以上基于瑞利散射的OFDR高温传感方法均没有采用普通光纤，制作成本较高，并且无法实现高温度分辨率、长距离以及更高温度下的温度传感。

光纤中的瑞利散射来源于光纤中折射率的随机波动，当光纤某处有温度变化时，便会引起该位置瑞利散射光的变化，进而获得瑞利散射光谱(Rayleigh backscatteringspectra,RBS)的波长漂移，利用温度变化前后的瑞利散射信号进行互相关运算得到的频谱漂移从而实现温度的测量。然而在1000℃附近时，时间间隔长的两组瑞利散射信号的相关性会急剧下降，造成互相关峰值位置的偏移，从而解调出错误的温度值，这在Chen^[3]所作的报道中也有所提及。

近年来，图像去噪方法已被广泛应用于OFDR分布式传感中。由于OFDR系统所测量的多为信息具有较高的相似性和冗余度，因此可以使用图像处理的方法来对二维图像进行降噪。例如小波去噪、高斯滤波^[4]、非局部均值^[5]、中值滤波^[6]、形状自适应主成分分析BM3D^[7]都已被用于OFDR分布式应变传感中，以提高传感性能和测量精度。在这几种方法中，性能最好的是形状自适应主成分分析BM3D，但其由于引入了大量复杂的小波变换和离散余弦变换，导致运算时间较长，限制了其在实际场景中的使用。另外，利用图像去噪来提高OFDR分布式高温传感的性能还未曾报道。

非局部哈尔变换(Non-Local Haar Transform,NLH)图像去噪方法是由YingkunHou^[8]等人提出的，通过将图像的块级非局部自相似性(Non-local selfsimilarity,NSS)提升至像素级NSS，更大程度利用了NSS先验，对最相似的像素群组执行简单的哈尔变换即可实现有效的图像去噪，并且能够较好地保留图像细节。整个算法不需要任何矩阵乘法运算，与形状自适应主成分分析BM3D方法相比大大降低了运算复杂度。

参考文献：

[1]Bulot,P.,Bernard,R.,Cieslikiewicz-Bouet,M.,Laffont,G.,and Douay,M.“Performance study of a zirconia-doped fiber for distributed temperaturesensing by OFDR at 800℃.”Sensors,21(11),article no.3788.

[2]A.K.Sang,D.K.Gifford,B.D.Dickerson,B.F.Fielder,and M.E.Froggatt,“One centimeter spatial resolution temperature measurements in a nuclearreactor using Rayleigh scatter in optical fiber,”Proc.SPIE 6619,66193D(2007).

[3]Chen,C.；Chen,L.；Bao,X.,“Distributed temperature profile inhydrogen flame measured by telecom fiber and its durability under flame byOFDR.”Opt.Express 2022,30,19390.

[4]S.Qu et al.,"High Spatial Resolution Investigation of OFDR Basedon Image Denoising Methods,"IEEE Sensors Journal,vol.21,no.17,pp.18871-18876,1Sept.1,2021.

[5]S.Zhao,et al.,“Accuracy improvement in OFDR based distributedsensing system by image processing,”Opt.,Lasers Eng.,vol.124,p.105824,Jan.2020.

[6]Q.Wang et al.,"Improving OFDR Distributed Fiber Sensing by FibersWith Enhanced Rayleigh Backscattering and Image Processing,"IEEE SensorsJournal,vol.22,no.19,pp.18471-18478,1Oct.1,2022.

[7]M.Pan et al.,"Long Distance Distributed Strain Sensing in OFDR byBM3D-SAPCA Image Denoising,"Journal ofLightwave Technology,2022.

[8]Y.Hou et al.,"NLH:A Blind Pixel-Level Non-Local Method for Real-World Image Denoising,"IEEE Transactions on Image Processing,vol.29.

[9]I.Daubechies and W.Sweldens,“Factoring wavelet transforms intolifting steps,”J.FourierAnal.Appl.,vol.4,no.3,pp.247–269,1998.

[10]W.Sweldens,“The lifting scheme:A custom-design construction ofbiorthogonal wavelets,”Appl.Comput.Harmon.Anal.,vol.3,no.2,pp.186–200,Apr.1996.

发明内容

本发明提供了利用NLH图像去噪光频域反射分布式高温测量方法，技术方案如下：

一种基于非局部哈尔变换图像去噪光频域反射分布式传感方法，包括下列步骤：

第一步，在可调谐激光器的触发下，在高温炉升温过程中利用OFDR系统连续采集在不同温度T₀,T₁,…,T_n下的瑞利散射信号，其中n为总测量组数，温度步长为Δt＞0，即Δt＝T₁-T₀＝…＝T_n-T_n-1，取相邻两组信号分别进行处理，共n-1组，温度较低组为参考信号，温度较高组为测量信号；对于取T₀和T₁两个温度下的瑞利散射信号，T₁＞T₀，T₀温度下的瑞利散射信号即为参考信号S_r，T₁温度下的瑞利散射信号即为测量信号S_m，以下第二步到第十二步，是对T₀和T₁两个温度下的瑞利散射信号的处理，其余n-2组信号均做相同处理；

第二步，分别对T₀温度下的参考信号S_r和T₁温度下的测量信号S_m进行快速傅里叶变换，得到T₀和T₁两个温度下对应光纤各个位置的距离域信号R_r和R_m；对R_r和R_m分别进行取窗补零操作得到R_ri和R_mi，其中取窗点数为N，每个窗的补零点数为M，i∈[1,D]，D为段数，M和D之间满足M＝有效点数/D；补零过程将每一段R_ri和R_mi填充成了具有(M+N)个点的数据段，每一段R_ri和R_mi都对应着光纤上相应位置在T₀和T₁两个温度下的瑞利散射信息；随后利用快速逆傅里叶变换将每一段R_ri和R_mi反变换回波长域，得到局部参考瑞利散射光谱S_ri和局部测量瑞利散射光谱S_mi；取S_ri和S_mi的幅值，并去除直流分量得到S_ri'和S_mi'；

第三步，将得到的D组S_ri'和S_mi'分别进行互相关计算，并对其结果进行归一化处理，每一组归一化互相关结果都对应于光纤上的一个位置；将D组归一化互相关结果沿光纤对应位置排布，即可得到沿光纤各个位置的互相关三维图，再将其投影到二维平面上，即可得到原始二维互相关强度图y；

第四步，在第三步中得到的原始二维互相关强度图y上按一个指定的步长N_step提取大小为N₁×N₁的图像块作为参考块，然后在以该参考块为中心的一个大小为N_s×N_s的邻域内进行基于欧氏距离的块匹配获得包括参考块本身的数量为N₂的相似图像块，将N₂个大小为N₁×N₁的相似图像块进行列扫描转换成N₂个大小为N₁ ²×1的矢量，将N₂个大小为N₁ ²×1的矢量拼接成一个大小为

的矩阵Y_l，在Y_l上将每一行作为参考行与其余的所有行计算欧氏距离ED，将与参考行欧氏距离ED最小的包括参考行在内的N₃行组成一个大小为N₃×N₂的矩阵Y_l'；

第五步，将第四步获得的大小为N₃×N₂的矩阵Y_l'执行可分的二维哈尔变换，即纵向和横向各自执行哈尔变换，在执行变换过程中应用提升哈尔变换，即只对行或列之间执行加权的平均，得到变换系数；

第六步，双硬阈值收缩：在第五步对大小为N₃×N₂的矩阵Y_l'中欧氏距离最小的两行应用提升哈尔变换的结果相当于对这两行的加权求差，这些加权求差后的结果系数基本全部来自于噪声，对这样的系数采用全部设置为零的方法；对其余哈尔变换的系数采用系数硬阈值进行系数收缩，系数收缩硬阈值参数为Thr；

第七步，对第六步执行完双硬阈值收缩后的哈尔变换系数执行逆哈尔变换，获得去噪后的二维矩阵，再将此二维矩阵扫描回原来的图像块，最终将这些图像块加权放回到图像的原始位置获得初步去噪图像；为了最好保留图像细节，采用了迭代策略，迭代次数为K，经过K次迭代得到的初步去噪图像记为y_K；用于第k次迭代的含噪图像，通过第k-1次迭代得到的去噪图像y_k-1按照y_k＝λy_k-1+(1-λ)y计算得到，其中k＝1,…,K,λ为正则化参数；

第八步，将第七步中经过K次迭代的初步去噪图像y_K和第三步中得到的原始含噪图像y同步执行块匹配和行匹配操作；按一个指定的步长N_step在y_K中提取大小为N₁×N₁的图像块作为参考块，然后在以该参考块为中心的一个大小为N_s×N_s的邻域内进行基于欧氏距离的块匹配获得包括参考块在内的数量为N₂的相似图像块，将N₂个大小为N₁×N₁的相似图像块进行列扫描转换成N₂个大小为N₁ ²×1的矢量，将N₂个大小为N₁ ²×1的矢量拼接成一个大小为

的矩阵Y_K，在Y_K上将每一行作为参考行与其余的所有行计算欧氏距离ED，将与包括参考行在内的参考行欧式距离最小的N₃行组成一个大小为N₃×N₂的矩阵Y′_K；原始含噪图像y进行同步的同参数的操作得到一个大小为N₃×N₂的矩阵Y_l'；经过块匹配和行匹配后获得两个大小为N₃×N₂的矩阵Y′_K和Y_l'；

第九步，将第八步获得的两个大小为N₃×N₂的矩阵Y′_K和Y_l'分别执行同样的可分的二维哈尔变换，即纵向和横向各自执行哈尔变换，在执行变换过程中应用了提升哈尔变换，即只对行或列之间执行加权的平均；其中Y′_K的变换系数为θ₁，Y_l'的变换系数为θ₂；

第十步，按公式(1)执行维纳滤波：

第十一步，为了提高去噪性能，将Y′_K的变换系数θ₁与维纳滤波得到的系数

再执行一次第十步所示维纳滤波，即执行2次迭代；将执行完2次迭代得到的系数执行逆哈尔变换，获取去噪后的二维矩阵，再将此二维矩阵扫描回原来的图像块，最终将这些图像块加权放回到图像中的原始位置获得最终去噪后的图像y_denoised；

第十二步，将去噪后的二维互相关强度图y_denoised转换成三维，即可得到去噪后的对应于温度变化范围T₀～T₁光纤各个位置的互相关三维图；再通过光纤每一个位置的互相关结果寻找其峰值位置，即可得到温度变化范围T₀～T₁光谱频移分布；

第十三步，n-1组相邻温度值的瑞利散射信号均执行完第二步～第十二步后，得到n-1个光谱频移分布，分别对应于n-1温度变化区间T₀～T₁,T₁～T₂,…,T_n-1～T_n,再将n-1个光谱频移分布相加，即可得到对应于温度变化T₀～T_n的光谱频移分布。

本发明提供的技术方案的有益效果是：

1、利用细径光纤实现了测试距离100m，空间分辨率2cm的分布式高温传感；

2、实现了温度分辨率1℃，温度测量范围950℃～1050℃的分布式高温测量，经过NLH图像去噪算法处理后得到了R-square值为0.9999的高线性度OFDR系统；

3、克服了光纤用于1000℃测量时，由于温度间隔大导致信号相似度差，从而导致无法测量的问题。

附图说明

图1为基于光频域反射的分布式高温传感装置；

图2为三维互相关分布利用NLH去噪前后对比图及通过寻找峰值得到的光谱频移结果利用NLH去噪前后对比图；

图3为1000℃～1010℃，温度间隔1℃利用NLH去噪前后累加结果对比；

图4为空间分辨率2cm，温度测量范围950℃～1050℃利用NLH去噪前后对比及NLH去噪后线性拟合图；

图5为基于NLH图像去噪OFDR实现分布式高温传感方法流程图；

附图中，各标号所代表的部件列表如下：

1：可调谐激光器； 2：第一平衡探测器；

3：20:80保偏耦合器； 4：1:99保偏光分束器；

5：50:50耦合器； 6：延迟光纤；

7：第一法拉第转镜； 8：第二法拉第转镜；

9：第一环形器； 10：第二环形器；

11：参考臂； 12：测试臂；

13：光混合器； 14：高温炉；

15：模数采集装置； 16：第二平衡探测器；

17：第三平衡探测器； 18：传感光纤；

19：计算机； 20：USB控制线；

21：主干涉仪； 22：基于附加干涉仪的时钟触发装置

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面对本发明实施方式做进一步详细描述。

在本专利中，我们在高温炉的升温过程中，温度间隔为1℃，测量范围从950℃到1050℃，利用OFDR系统采集细径光纤在这101个温度值下的背向瑞利散射信号。在此过程中，由于温度高，时间间隔久的两组背向瑞利散射信号的相似度会变低，因此高温炉升温要尽量快，我们取1℃为温度间隔，时间间隔约4s，既可以保证相邻两个温度值下的两组背向瑞利散射信号较高的相似度，又可以确保在每一个温度值有足够的时间采集信号。我们将相邻两个温度值下(即间隔1℃)的背向瑞利散射信号归为1组(共100组)进行后续的处理，这里仅以950℃和951℃为例进行说明，其余相邻两个温度值下的信号均以相同方法进行处理。我们将温度值较低即950℃的那组信号记为参考信号，将温度值较高即951℃的那组信号记为测量信号，将这两组信号通过FFT、取窗补零、IFFT、去直流和互相关处理得到光纤局部的互相关分布，将得到的各个窗位置处的局部互相关分布沿光纤按顺序排列得到三维互相关分布，将得到的三维互相关分布投影到二维平面转化为二维图像，随后再利用NLH图像去噪算法对该二维图像进行处理，从而有效地去除随机噪声，再将此二维图像转换成三维图像，通过寻找各个窗的峰值位置，即可得到950℃～951℃的光谱频移分布。在将100组相邻温度值的背向瑞利散射信号都按上述方法处理后，将950℃～951℃的光谱频移分布和951℃～952℃的光谱频移分布相加，即可得到950℃～952℃的光谱频移分布；以此类推，将100组相邻温度值的背向瑞利散射信号处理完后得到的光谱频移分布相加，即可得到950℃～1050℃的光谱频移分布。基于该方法，我们利用细径光纤实现了950℃～1050℃的温度测量范围，空间分辨率为2cm，温度分辨率为1℃，经拟合得到的光谱频移值与温度变化呈线性关系，R-square为0.9999。

实施例一：

本发明采用的OFDR装置，与丁振扬课题组之前申请的专利所采用的OFDR装置基本一致，但为了保持本案的易于理解性，下面对采用的装置做介绍。

本实例包括基于光频域反射系统分布式光纤传感装置

光频域反射的分布式光纤传感装置包括：可调谐激光器1、1:99保偏光分束器4、计算机19、USB控制线20、模数采集装置15、基于附加干涉仪的时钟触发装置22、主干涉仪21。

其中，基于辅助干涉仪的时钟触发装置22包括：第一平衡探测器2、50:50耦合器5、延迟光纤6、第一法拉第转镜7、第二法拉第转镜8和第一环形器9。基于附加干涉仪的时钟触发装置22用于实现等光频间隔采样，其目的是补偿光源的非线性扫频。

其中，主干涉仪21包括：20:80保偏耦合器3、第二环形器10、光混合器13、第二平衡探测器16、第三平衡探测器17、参考臂11、测试臂12、高温炉14和传感光纤18，其为细径光纤。主干涉仪21是光频域反射的分布式光纤传感装置的核心，其为改进型马赫泽德干涉仪。

USB控制线20的输出端与可调谐激光器1的输入端相连；USB控制线20的输入端与计算机19的输出端相连；可调谐激光器1与1:99保偏光分束器4的a端口相连；1:99保偏光分束器4的b端口即1％分光口与第一环形器9的a端口相连；1:99保偏光分束器4的c端口即99％分光口与20:80保偏耦合器3的a端口相连；第一环形器9的b端口与50:50耦合器5的a端口相连；第一环形器9的c端口与第一平衡探测器2的输入端相连；50:50耦合器5的b端口与第一平衡探测器2的输入端相连；50:50耦合器5的c端口通过延迟光纤6与第一法拉第转镜7相连；50:50耦合器5的d端口与第二法拉第转镜8相连；第一平衡探测器2的输出端与模数采集装置15的输入端相连；20:80保偏耦合器3的c端口即20％分光口通过参考臂11与光混合器13的输入端a相连；20:80保偏耦合器3的d端口即80％分光口通过测试臂12与第二环形器10的a端口相连；第二环形器10的c端口与传感光纤18相连；第二环形器10的b端口与光混合器13的输入端b相连；光混合器13的输出端c和输出端d与第二平衡探测器16的两个输入端相连；光混合器13的输出端e和输出端f与第三平衡探测器17的两个输入端相连；第二平衡探测器16的输出端与模数采集装置15的输入端相连；第三平衡探测器17的输出端与模数采集装置15的输入端相连；模数采集装置15的输出端与计算机19的输入端相连。

装置工作时，计算机19通过USB控制线20控制可调谐激光器1控制调谐速度、中心波长、调谐启动等；可调谐激光器1的出射光由1:99保偏光分束器4的a端口进入，从1:99保偏光分束器4的b端口即1％分光口进入第一环形器9的a端口，经第一环形器9的b端口进入50:50耦合器5的a端口，50％从耦合器5的c端口出射，50％从耦合器5的d端口出射，来自50:50耦合器5的c端口的出射光经延迟光纤6被第一法拉第旋转镜7反射，返回至50:50耦合器5的c端口，来自50:50耦合器5的d端口的出射光被第二法拉第旋转镜8反射，返回至50:50耦合器5的d端口，两束光在50:50耦合器5中发生干涉，从50:50耦合器5的b端口输出；该输出光进入第一平衡探测器2，第一平衡探测器2将探测到的光信号转换为干涉拍频信号传输至模数采集装置15，作为模数采集装置15的外部时钟信号。

从1:99光分束器4的c端口即99％分光口进入20:80保偏耦合器3的a端口；经过20:80保偏耦合器3从c端口即20％分光口进入参考臂11，从d端口即80％分光口进入测试臂12，光从第二环形器10的a端口进入，从第二环形器10的c端口进入传感光纤18，传感光纤18的背向散射光从第二环行器10端口c端口进入，从第二环行器10的b端口输出；参考臂11中输出的参考光进入光混合器13的a端口与第二环行器10端口b端口进入光混合器13的b端口的测试光进行合束，形成拍频干涉；光混合器13通过c端口、d端口输入到第二平衡探测器16，第二平衡探测器16将探测到的光信号转换为干涉拍频信号传输至模数采集装置15；光混合器13通过e端口、f端口输入到第三平衡探测器17，第二平衡探测器17将探测到的光信号转换为干涉拍频信号传输至模数采集装置15；模数采集装置15在附加干涉仪的时钟触发装置22形成的外部时钟信号作用下将采集到的模拟电信号传输至计算机19。

USB控制线20用于计算机19通过其控制可调谐激光器1。

可调谐激光器1为光频域反射系统提供光源，其光频能够进行线性扫描。

第一环形器9防止附加干涉仪中50:50耦合器5的b端口反射光进入激光器。

50:50耦合器5用于光干涉。

延迟光纤6用于实现非等臂的拍频干涉，能够根据拍频和延迟光纤长度得到光频。

高温炉14用于使传感光纤18产生温度的变化。

第一法拉第转镜7和第二法拉第转镜8用于为干涉仪提供反射，且能够消除干涉仪的偏振衰落现象。

光混合器13完成对信号进行偏振分束，使参考光和测试光在偏振分束时两个正交方向上光强基本一致，消除偏振衰落噪声的影响，实现参考光和测试光的合束，形成拍频干涉。

计算机19：对模数采集装置15采集的干涉信号进行数据处理，实现基于光频域反射利用细径光纤分布式测量高温的光纤传感。

本发明的利用NLH图像去噪光频域反射分布式高温测量方法，具体步骤如下：

第一步，在可调谐激光器的触发下，在高温炉升温过程中利用OFDR系统连续采集在不同温度T₀,T₁,…,T_n下的瑞利散射信号，其中n为总测量组数，温度步长为Δt＞0，即Δt＝T₁-T₀＝…＝T_n-T_n-1，取相邻两组信号进行处理(共n-1组)，温度较低组为参考信号，温度较高组为测量信号，如取T₀和T₁两个温度下的瑞利散射信号，T₁＞T₀，T₀温度下的瑞利散射信号即为参考信号S_r，T₁温度下的瑞利散射信号即为测量信号S_m。以下步骤中仅以T₀和T₁两个温度下的瑞利散射信号为例，其余n-2组信号均做相同处理。

第二步，分别对T₀温度下的参考信号S_r和T₁温度下的测量信号S_m进行快速傅里叶变换，得到T₀和T₁两个温度下对应光纤各个位置的距离域信号R_r和R_m。对R_r和R_m分别进行取窗补零操作得到R_ri和R_mi，其中取窗点数为N，每个窗的补零点数为M，i∈[1,D]，D为段数，N和D之间满足M＝有效点数/D。补零过程将每一段R_ri和R_mi填充成了具有(M+N)个点的数据段，每一段R_ri和R_mi都对应着光纤上相应位置在T₀和T₁两个温度下的瑞利散射信息。随后利用快速逆傅里叶变换将每一段R_ri和R_mi反变换回波长域，得到局部参考瑞利散射光谱S_ri和局部测量瑞利散射光谱S_mi。取S_ri和S_mi的幅值，并去除直流分量得到S_ri'和S_mi'。

第三步，将得到的D组S_ri'和S_mi'分别进行互相关计算，并对其结果进行归一化处理，每一组归一化互相关结果都对应于光纤上的一个位置。将D组归一化互相关结果沿光纤对应位置排布，即可得到沿光纤各个位置的互相关三维图，再将其投影到二维平面上，即可得到原始二维互相关强度图y。

第四步，在第三步中得到的原始二维互相关强度图y上按一个指定的步长N_step提取大小为N₁×N₁的图像块作为参考块，然后在以该参考块为中心的一个大小为N_s×N_s的邻域内进行基于欧氏距离的块匹配获得数量为N₂的相似图像块(包括参考块本身)，将N₂个大小为N₁×N₁的相似图像块进行列扫描转换成N₂个大小为N₁ ²×1的矢量，将N₂个大小为N₁ ²×1的矢量拼接成一个大小为

的矩阵Y_l，在Y_l上将每一行作为参考行与其余的所有行计算欧氏距离ED，将与参考行欧氏距离ED最小的N₃行(包括参考行本身)组成一个大小为N₃×N₂的矩阵Y_l'。

第五步，将第四步获得的大小为N₃×N₂的矩阵Y_l'执行可分的二维哈尔变换，即纵向和横向各自执行哈尔变换，在执行变换过程中应用了提升哈尔变换，即只对行或列之间执行加权的平均，得到变换系数。

第六步，双硬阈值收缩。在第五步对大小为N₃×N₂的矩阵Y_l'中欧氏距离最小的两行应用提升哈尔变换的结果相当于对这两行的加权求差，这些求差后的结果系数基本全部来自于噪声，对这样的系数采用了全部设置为零的方法；对其余哈尔变换的系数采用系数硬阈值进行系数收缩，系数收缩硬阈值参数为Thr。

第七步，对第六步执行完双硬阈值收缩后的哈尔变换系数执行逆哈尔变换，获得去噪后的二维矩阵，再将此二维矩阵扫描回原来的图像块，最终将这些图像块加权放回到图像的原始位置获得初步去噪图像。由于每次参与去噪的最相似像素数目是有限的，为了最好保留图像细节，Thr不能选得很大，仅执行一次上述步骤很难将图像噪声全部去除，因此采用了迭代策略，迭代次数为K，经过K次迭代得到的初步去噪图像记为y_K。用于第k次迭代的含噪图像，通过第k-1次迭代得到的去噪图像y_k-1按照y_k＝λy_k-1+(1-λ)y计算得到，其中k＝1,…,K,λ为正则化参数。

第八步，将第七步中经过K次迭代的初步去噪图像y_K和第三步中得到的原始含噪图像y同步执行块匹配和行匹配操作。按一个指定的步长N_step在y_K中提取大小为N₁×N₁的图像块作为参考块，然后在以该参考块为中心的一个大小为N_s×N_s的邻域内进行基于欧氏距离的块匹配获得数量为N₂的相似图像块(包括参考块本身)，将N₂个大小为N₁×N₁的相似图像块进行列扫描转换成N₂个大小为N₁ ²×1的矢量，将N₂个大小为N₁ ²×1的矢量拼接成一个大小为

的矩阵Y_K，在Y_K上将每一行作为参考行与其余的所有行计算欧氏距离ED，将与参考行欧式距离最小的N₃行(包括参考行本身)组成一个大小为N₃×N₂的矩阵Y′_K。原始含噪图像y进行同步的同参数的操作得到一个大小为N₃×N₂的矩阵Y_l′。经过块匹配和行匹配后获得两个大小为N₃×N₂的矩阵Y′_K和Y_l′。

第九步，将第八步获得的两个大小为N₃×N₂的矩阵Y′_K和Y_l′分别执行同样的可分的二维哈尔变换，即纵向和横向各自执行哈尔变换，在执行变换过程中应用了提升哈尔变换，即只对行或列之间执行加权的平均。其中Y′_K的变换系数为θ₁，Y_l′的变换系数为θ₂。

第十步，按公式(1)执行维纳滤波。

第十一步，为了提高去噪性能，将Y′_K的变换系数θ₁与维纳滤波得到的系数

再执行一次第十步所示维纳滤波，即执行2次迭代。将执行完2次迭代得到的系数执行逆哈尔变换，获取去噪后的二维矩阵，再将此二维矩阵扫描回原来的图像块，最终将这些图像块加权放回到图像中的原始位置获得最终去噪后的图像y_denoised。

第十二步，将去噪后的二维互相关强度图y_denoised转换成三维，即可得到去噪后的对应于温度变化范围T₀～T₁光纤各个位置的互相关三维图。再通过光纤每一个位置的互相关结果寻找其峰值位置，即可得到温度变化范围T₀～T₁光谱频移分布。

第十三步，n-1＝100组相邻温度值的瑞利散射信号均执行完第二步～第十二步后，即可得到n-1＝100个光谱频移分布，分别对应于n-1＝100个温度变化区间T₀～T₁即950℃～951℃，T₁～T₂即951℃～952℃，...，T_n-1～T_n即1049℃～1050℃，再将n-1＝100个光谱频移分布相加，即可得到对应于温度变化区间T₀～T_n即950℃～1050℃的光谱频移分布。

实施例二：

本发明提供了利用NLH图像去噪光频域反射分布式高温测量方法，技术方案如下：

第一步，在可调谐激光器的触发下，在高温炉升温过程中利用OFDR系统连续采集在不同温度T₀，T₁，…，T_n下的瑞利散射信号，其中n为总测量组数，温度步长为Δt＞0，即Δt＝T₁-T₀＝…＝T_n-T_n-1，取相邻两组信号进行处理(共n-1组)，温度较低组为参考信号，温度较高组为测量信号，如取T₀和T₁两个温度下的瑞利散射信号，T₁＞T₀，T₀温度下的瑞利散射信号即为参考信号S_r，T₁温度下的瑞利散射信号即为测量信号S_m。在本专利中，T₀＝950℃，T₁＝951℃，T_n-1＝1049℃，T_n＝1050℃，n＝101，Δt＝1℃，T₀,T₁,…,T_n-1,T_n的选择需要考虑温度测量范围以及测量精度，Δt的选择需要考虑相邻两组信号能够保持较高的相似性以及高温炉的温度调节精度，n与温度测量范围和温度测量间隔有关。以下步骤中仅以T₀＝950℃和T₁＝951℃两个温度下的瑞利散射信号为例，其余n-2组信号均做相同处理。

第二步，分别对T₀＝950℃温度下的参考信号S_r和T₁＝951℃温度下的测量信号S_m进行快速傅里叶变换，得到T₀＝950℃和T₁＝951℃两个温度下对应光纤各个位置的距离域信号R_r和R_m。对R_r和R_m分别进行取窗补零操作得到R_ri和R_mi，其中取窗点数为N，每个窗的补零点数为M，i∈[1,D]，D为段数，N和D之间满足M＝有效点数/D。在本专利中，N＝200，M＝5050，N和M的选择需要考虑OFDR系统的空间分辨率和频谱精细度，有效点数＝1050000，D＝5250,补零过程将每一段R_ri和R_mi填充成了具有M+N＝5250个点的数据段，每一段R_ri和R_mi都对应着光纤上相应位置在T₀＝950℃和T₁＝951℃两个温度下的瑞利散射信息。随后利用快速逆傅里叶变换将每一段R_ri和R_mi反变换回波长域，得到局部参考瑞利散射光谱S_ri和局部测量瑞利散射光谱S_mi。取S_ri和S_mi的幅值，并去除直流分量得到S_ri'和S_mi'。

第三步，将得到的D＝5250组S_ri'和S_mi'分别进行互相关计算，并对其结果进行归一化处理，每一组归一化互相关结果都对应于光纤上的一个位置。将D＝5250组归一化互相关结果沿光纤对应位置排布，即可得到沿光纤各个位置的互相关三维图，再将其投影到二维平面上，即可得到原始二维互相关强度图y。

第四步，在第三步中得到的原始二维互相关强度图y上按一个指定的步长N_step提取大小为N₁×N₁的图像块作为参考块，然后在以该参考块为中心的一个大小为N_s×N_s的邻域内进行基于欧氏距离的块匹配获得数量为N₂的相似图像块(包括参考块本身)，将N₂个大小为N₁×N₁的相似图像块进行列扫描转换成N₂个大小为N₁ ²×1的矢量，将N₂个大小为N₁ ²×1的矢量拼接成一个大小为

的矩阵Y_l，在Y_l上将每一行作为参考行与其余的所有行计算欧氏距离ED，将与参考行欧氏距离ED最小的N₃行(包括参考行本身)组成一个大小为N₃×N₂的矩阵Y_l'。在本专利中，N_step＝6，N₁＝9，N_s＝43，N₂＝16，N₃＝4，N_step、N₁、N_s、N₂、N₃的选择需要考虑原始含噪图像的噪声水平。

第五步，将第四步获得的大小为N₃×N₂即4×16的矩阵Y_l'执行可分的二维哈尔变换，即纵向和横向各自执行哈尔变换，在执行变换过程中应用了提升哈尔变换^[9][10]，即只对行或列之间执行加权的平均，得到变换系数。

第六步，双硬阈值收缩。在第五步对大小为N₃×N₂即4×16的矩阵Y_l'中欧氏距离最小的两行应用提升哈尔变换的结果相当于对这两行的加权求差，这些求差后的结果系数基本全部来自于噪声，对这样的系数采用了全部设置为零的方法；对其余哈尔变换的系数采用系数硬阈值进行系数收缩，系数收缩硬阈值参数为Thr。在本专利中，Thr＝1.65，Thr的选择需要考虑原始含噪图像的噪声水平。

第七步，对第六步执行完双硬阈值收缩后的哈尔变换系数执行逆哈尔变换，获得去噪后的二维矩阵，再将此二维矩阵扫描回原来的图像块，最终将这些图像块加权放回到图像的原始位置获得初步去噪图像。由于每次参与去噪的最相似像素数目是有限的，为了最好保留图像细节，Thr不能选得很大，仅执行一次上述步骤很难将图像噪声全部去除，因此采用了迭代策略，迭代次数为K，经过K次迭代得到的初步去噪图像记为y_K。用于第k次迭代的含噪图像，通过第k-1次迭代得到的去噪图像y_k-1按照y_k＝λy_k-1+(1-λ)y计算得到，其中k＝1,…,K,λ为正则化参数。在本专利中，K＝3，λ＝0.6，K和λ的选择需要考虑原始含噪图像的噪声水平。

第八步，将第七步中经过K＝3次迭代的初步去噪图像y₃和第三步中得到的原始含噪图像y同步执行块匹配和行匹配操作。按一个指定的步长N_step＝6在y₃中提取大小为N₁×N₁即9×9的图像块作为参考块，然后在以该参考块为中心的一个大小为N_s×N_s即43×43的邻域内进行基于欧氏距离的块匹配获得数量为N₂＝16的相似图像块(包括参考块本身)，将N₂＝16个大小为N₁×N₁即9×9的相似图像块进行列扫描转换成N₂＝16个大小为N₁ ²×1即81×1的矢量，将N₂＝16个大小为N₁ ²×1即81×1的矢量拼接成一个大小为

即81×16的矩阵Y₃，在Y₃上将每一行作为参考行与其余的所有行计算欧氏距离ED，将与参考行欧式距离最小的N₃＝4行(包括参考行本身)组成一个大小为N₃×N₂即4×16的矩阵Y₃′。原始含噪图像y进行同步的同参数的操作得到一个大小为N₃×N₂即4×16的矩阵Y_l′。经过块匹配和行匹配后获得两个大小为N₃×N₂即4×16的矩阵Y₃′和Y_l′。

第九步，将第八步获得的两个大小为N₃×N₂即4×16的矩阵Y′_K和Y_l′分别执行同样的可分的二维哈尔变换，即纵向和横向各自执行哈尔变换，在执行变换过程中应用了提升哈尔变换，即只对行或列之间执行加权的平均。其中Y′_K的变换系数为θ₁，Y_l′的变换系数为θ₂。

第十步，按公式(1)执行维纳滤波。

在本专利中，σ＝25，σ的选择需要考虑图像的去噪效果和细节的保留程度。

第十一步，为了提高去噪性能，将Y′₃的变换系数θ₁与维纳滤波得到的系数

再执行一次第十步所示维纳滤波，即执行2次迭代。将执行完2次迭代得到的系数执行逆哈尔变换，获取去噪后的二维矩阵，再将此二维矩阵扫描回原来的图像块，最终将这些图像块加权放回到图像中的原始位置获得最终去噪后的图像y_denoised。

第十二步，将去噪后的二维互相关强度图y_denoised转换成三维，即可得到去噪后的对应于温度变化范围T₀～T₁即950℃～951℃光纤各个位置的互相关三维图。再通过光纤每一个位置的互相关结果寻找其峰值位置，即可得到温度变化范围T₀～T₁即950℃～951℃的光谱频移分布。

第十三步，n-1＝100组相邻温度值的瑞利散射信号均执行完第二步～第十二步后，即可得到n-1＝100个光谱频移分布，分别对应于n-1＝100个温度变化区间T₀～T₁即950℃～951℃，T₁～T₂即951℃～952℃，...，T_n-1～T_n即1049℃～1050℃，再将n-1＝100个光谱频移分布相加，即可得到对应于温度变化区间T₀～T_n即950℃～1050℃的光谱频移分布。

第十三步，n-1＝100组相邻温度值的瑞利散射信号均执行完第二步～第十二步后，即可得到n-1＝100个光谱频移分布，分别对应于n-1＝100个温度变化区间T₀～T₁即950℃～951℃，T₁～T₂即951℃～952℃，…，T_n-1～T_n即1049℃～1050℃，再将n-1＝100个光谱频移分布相加，即可得到对应于温度变化区间T₀～T_n即950℃～1050℃的光谱频移分布。

下面结合具体的试验对NLH对OFDR分布式高温测量系统的去噪效果进行可行性验证，详见下文描述：

实施例三：

本发明实施例验证采用的传感光纤18为细径光纤，传感光纤总长度约为102.5m。在实验中使用OFDR系统的调谐范围为8.33nm，扫频速率为100nm/s,采样点数为5M，附加干涉仪光纤长度为500m。取窗点数N为200，对应着2cm的空间分辨率，补零点数M为5050。

图2a是间隔1℃原始互相关三维图像，图2b是通过寻找峰值位置得到的光谱频移结果。可以看出，由于受到噪声的影响，尽管光纤各个位置所处环境温度变化相同，但是光谱频移大小出现了不规则的抖动。图2c是经过NLH图像去噪之后得到的互相关三维图像，图2d是经过NLH图像去噪之后得到的光谱频移结果。可以看出，经过NLH图像去噪之后，光纤上温度变化分布变得平稳，这说明NLH图像去噪提高了OFDR测温的性能。

为了展示该方法对1000℃的分布式测量效果，如前所述，我们将1001℃～1002℃两组瑞利散射信号得到的去噪后的光谱频移分布加到1000℃～1001℃两组瑞利散射信号得到的去噪后的光谱频移分布，即可得到1000℃～1002℃两组瑞利散射信号得到的去噪后的光谱频移分布。以此类推，即可得到1000℃～1010℃，温度分辨率1℃，空间分辨率2cm的累加结果。图3a为1000℃～1010℃，温度间隔1℃，空间分辨率2cm的原始累加结果，图3b为1000℃～1010℃，温度间隔1℃，空间分辨率2cm的NLH去噪后累加结果。

实施例四：

为了验证利用NLH图像去噪后OFDR系统得到的光谱移动大小与温度变化之间的线性关系，我们通过累加的方式得到了950℃～960℃(间隔10℃)、950℃～970℃(间隔20℃)、950℃～980℃(间隔30℃)、950℃～990℃(间隔40℃)、950℃～1000℃(间隔50℃)、950℃～1010℃(间隔60℃)、950℃～1020℃(间隔70℃)、950℃～1030℃(间隔80℃)、950℃～1040℃(间隔90℃)、950℃～1050℃(间隔100℃)的光谱移动结果。利用NLH去噪前的温度分布如图4a所示，去噪后的温度分布如图4b所示，由图4b的结果进行拟合得到图4c的线性拟合结果。从图4c中可以看出，去噪后系统的线性程度很高，R-square值达到0.9999，光谱频移与温度变化呈线性变化关系。

综上所述，我们NLH图像去噪算法对分布式高温测量的二维互相关分布图进行去噪，实现了测量100m细径光纤末端，空间分辨率为2cm，温度分辨率为1℃且温度测量范围为950℃～1050℃的OFDR系统。去噪后的OFDR系统呈现出很高的线性度，R-square值达到0.9999。

本发明实施例对各器件的型号除做特殊说明的以外，其他器件的型号不做限制，只要能完成上述功能的器件均可。

本领域技术人员可以理解附图只是一个优选实施例的示意图，上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种基于非局部哈尔变换图像去噪光频域反射分布式传感方法，包括下列步骤：

第一步，在可调谐激光器的触发下，在高温炉升温过程中利用OFDR系统连续采集在不同温度T₀,T₁,…,T_n下的瑞利散射信号，其中n为总测量组数，温度步长为Δt＞0，即Δt＝T₁-T₀＝…＝T_n-T_n-1，取相邻两组信号分别进行处理，共n-1组，温度较低组为参考信号，温度较高组为测量信号；对于取T₀和T₁两个温度下的瑞利散射信号，T₁＞T₀，T₀温度下的瑞利散射信号即为参考信号S_r，T₁温度下的瑞利散射信号即为测量信号S_m，以下第二步到第十二步，是对T₀和T₁两个温度下的瑞利散射信号的处理，其余n-2组信号均做相同处理；

第二步，分别对T₀温度下的参考信号S_r和T₁温度下的测量信号S_m进行快速傅里叶变换，得到T₀和T₁两个温度下对应光纤各个位置的距离域信号R_r和R_m；对R_r和R_m分别进行取窗补零操作得到R_ri和R_mi，其中取窗点数为N，每个窗的补零点数为M，i∈[1,D]，D为段数，M和D之间满足M＝有效点数/D；补零过程将每一段R_ri和R_mi填充成了具有(M+N)个点的数据段，每一段R_ri和R_mi都对应着光纤上相应位置在T₀和T₁两个温度下的瑞利散射信息；随后利用快速逆傅里叶变换将每一段R_ri和R_mi反变换回波长域，得到局部参考瑞利散射光谱S_ri和局部测量瑞利散射光谱S_mi；取S_ri和S_mi的幅值，并去除直流分量得到S_ri′和S_mi′；

第三步，将得到的D组S_ri′和S_mi′分别进行互相关计算，并对其结果进行归一化处理，每一组归一化互相关结果都对应于光纤上的一个位置；将D组归一化互相关结果沿光纤对应位置排布，即可得到沿光纤各个位置的互相关三维图，再将其投影到二维平面上，即可得到原始二维互相关强度图y；

第四步，在第三步中得到的原始二维互相关强度图y上按一个指定的步长N_step提取大小为N₁×N₁的图像块作为参考块，然后在以该参考块为中心的一个大小为N_s×N_s的邻域内进行基于欧氏距离的块匹配获得包括参考块本身的数量为N₂的相似图像块，将N₂个大小为N₁×N₁的相似图像块进行列扫描转换成N₂个大小为N₁ ²×1的矢量，将N₂个大小为N₁ ²×1的矢量拼接成一个大小为N₁ ²×N₂的矩阵Y_l，在Y_l上将每一行作为参考行与其余的所有行计算欧氏距离ED，将与参考行欧氏距离ED最小的包括参考行在内的N₃行组成一个大小为N₃×N₂的矩阵Y_l′；

第五步，将第四步获得的大小为N₃×N₂的矩阵Y_l′执行可分的二维哈尔变换，即纵向和横向各自执行哈尔变换，在执行变换过程中应用提升哈尔变换，即只对行或列之间执行加权的平均，得到变换系数；

第六步，双硬阈值收缩：在第五步对大小为N₃×N₂的矩阵Y_l′中欧氏距离最小的两行应用提升哈尔变换的结果相当于对这两行的加权求差，这些加权求差后的结果系数基本全部来自于噪声，对这样的系数采用全部设置为零的方法；对其余哈尔变换的系数采用系数硬阈值进行系数收缩，系数收缩硬阈值参数为Thr；

第七步，对第六步执行完双硬阈值收缩后的哈尔变换系数执行逆哈尔变换，获得去噪后的二维矩阵，再将此二维矩阵扫描回原来的图像块，最终将这些图像块加权放回到图像的原始位置获得初步去噪图像；为了最好保留图像细节，采用了迭代策略，迭代次数为K，经过K次迭代得到的初步去噪图像记为y_K；用于第k次迭代的含噪图像，通过第k-1次迭代得到的去噪图像y_k-1按照y_k＝λy_k-1+(1-λ)y计算得到，其中k＝1,…,K,λ为正则化参数；

第八步，将第七步中经过K次迭代的初步去噪图像y_K和第三步中得到的原始含噪图像y同步执行块匹配和行匹配操作；按一个指定的步长N_step在y_K中提取大小为N₁×N₁的图像块作为参考块，然后在以该参考块为中心的一个大小为N_s×N_s的邻域内进行基于欧氏距离的块匹配获得包括参考块在内的数量为N₂的相似图像块，将N₂个大小为N₁×N₁的相似图像块进行列扫描转换成N₂个大小为N₁ ²×1的矢量，将N₂个大小为N₁ ²×1的矢量拼接成一个大小为N₁ ²×N₂的矩阵Y_K，在Y_K上将每一行作为参考行与其余的所有行计算欧氏距离ED，将与包括参考行在内的参考行欧式距离最小的N₃行组成一个大小为N₃×N₂的矩阵Y′_K；原始含噪图像y进行同步的同参数的操作得到一个大小为N₃×N₂的矩阵Y_l′；经过块匹配和行匹配后获得两个大小为N₃×N₂的矩阵Y′_K和Y_l′；

第九步，将第八步获得的两个大小为N₃×N₂的矩阵Y′_K和Y_l′分别执行同样的可分的二维哈尔变换，即纵向和横向各自执行哈尔变换，在执行变换过程中应用了提升哈尔变换，即只对行或列之间执行加权的平均；其中Y′_K的变换系数为θ₁，Y_l′的变换系数为θ₂；

第十步，按公式(1)执行维纳滤波：

第十一步，为了提高去噪性能，将Y′_K的变换系数θ₁与维纳滤波得到的系数

再执行一次第十步所示维纳滤波，即执行2次迭代；将执行完2次迭代得到的系数执行逆哈尔变换，获取去噪后的二维矩阵，再将此二维矩阵扫描回原来的图像块，最终将这些图像块加权放回到图像中的原始位置获得最终去噪后的图像y_denoised；

第十二步，将去噪后的二维互相关强度图y_denoised转换成三维，即可得到去噪后的对应于温度变化范围T₀～T₁光纤各个位置的互相关三维图；再通过光纤每一个位置的互相关结果寻找其峰值位置，即可得到温度变化范围T₀～T₁光谱频移分布；

第十三步，n-1组相邻温度值的瑞利散射信号均执行完第二步～第十二步后，得到n-1个光谱频移分布，分别对应于n-1温度变化区间T₀～T₁,T₁～T₂,…,T_n-1～T_n,再将n-1个光谱频移分布相加，即可得到对应于温度变化T₀～T_n的光谱频移分布。

Images