基于黎曼函数的宽带布里渊光谱信号解析算法

针对布里渊激光雷达领域现在普遍使用的光源过于昂贵、笨重的问题，对可促成使用廉价紧凑光源的宽带布里渊信号的解谱算法进行了研究.在宽带布里渊光谱信号的处理过程中，首先构造了黎曼函数的一个变体，再利用该函数在布里渊频移处取极小值性质，从泵浦光-布里渊光的1：1混合光谱信号中复原出了泵浦光的频谱并得到相对应的布里渊频移.对测试数据集的实验表明，在没有噪声的理想情况下算法准确率接近100%；而非理想状态下只有当噪声小于-30 dB，且泵浦光和布里渊光幅值之比小于1.05时，算法才能得到准确结果.

Abstract

As the seed-injected lasers widely used in Brillouin Lidar area are too expensive and huge, the algorithm for a broadband Brillouin signal should be developed, which could accelerate the application of compact and economic diode lasers. A variation of Thomae's function was found during the signal processing of broadband Brillouin lidar, and the property of this function-reaching minimum value at the Brillouin frequency shift was used to recover the original signal spectrum and the corresponding frequency shift from a 1:1 superposition of pump light and Brillouin light spectra. Experiments using test data show nearly 100% accuracies for ideal cases; but for non-ideal cases, only when the noise level is less than -30 dB and the ratio of pump light to Brillouin light is less than 1.05 can this algorithm obtain accurate result.

0 引言

近十年来，利用布里渊散射对大气、海洋等进行探测成为一种新兴的遥感方式 [1] .布里渊雷达通常包括窄线宽激光器、准直聚焦系统、高分辨率光谱分辨装置.其测量过程和普通激光雷达大致相同：光源发射探测脉冲，经过准直聚焦系统照射介质，然后接收介质反射(散射)的回波进行探测.布里渊雷达直接获得的信号是后向布里渊散射光.海水的布里渊信号经过光谱解析后，从布里渊频移可以得到盐度、密度、声速等信息 [2 - 4] ，从布里渊线宽可以得到粘度信息 [5] .大气的布里渊信号解析后还能用于测量一些痕量气体 [6] .

布里渊频移的数值很小，一般气体的频移是几百MHz [7] ，对液体、固体是几个GHz [8 - 10] ，因此对光源和探测器都提出了很高的要求.在探测器方面，获得布里渊光谱图需要有很高的光谱分辨精度，目前主要的方法有法布里-珀罗(Fabry-Pero, F-P)标准具+增强电荷耦合器(Intensified Charge Coupled Device, ICCD) [11] 、扫描F-P干涉仪+光电倍增管 [12] 、碘分子窄带滤波器+边缘探测技术 [13 - 14] 、利用高频微波源产生并测量拍频 [15] 、基于塞曼效应进行调谐的超窄带可调谐光学滤波器(Excited State Faraday Anomalous Dispersion Optical Filter, ESFADOFs) [16] 、布里渊光时域分析(Brillouin Optical Time Domain Analysis，BOTDA) [17 - 19] .BOTDA常用于分布式光纤传感，需要预先铺设固定光纤，故通常只能用于固定结构的桥梁、石油设备、水利工程构筑物的检测，而不适用于大规模普查性质的遥感检测.其他方法中，扫描F-P干涉仪和可调谐窄带滤波器采集光谱需要时间，无法做到实时检测；碘分子吸收池吸收线位置过于固定，对泵浦光种子的稳频特性提出了非常高的要求，给实际应用带来不便.综上，F-P标准具结合ICCD摄谱在布里渊雷达遥感中是最有优势的 [20] .除了探测器，在光源方面布里渊雷达通常要求泵浦光是窄线宽的.这是因为窄线宽的激发光源有助于在布里渊频移特别微小时分别不同的谱峰.

近年来，已经有报道用宽带飞秒激光作为激发源获得布里渊光谱，已经成功用宽带激光测量了声速、热效应和材料的纳米级缺陷 [21 - 23] .这些工作均使用钛宝石飞秒激光器，利用近年流行的泵浦-探测技术(pump-probe)获得后向布里渊散射光谱.通过改变泵浦-探测的时间延迟，可以发现探测光的反射率随度时间延迟存在周期为数皮秒的振荡，对该振荡曲线进行傅里叶变换，得到的就是布里渊光谱.由于泵浦-探测技术要求有皮秒分辨的光脉冲延迟器，这项技术是无法应用到纳秒或准连续激光脉冲上的.综上，现有报道的布里渊光谱技术即便不用窄线宽激光器，也需要体积庞大、价格高昂的钛宝石飞秒激光.如果一个普通的宽带光源(例如二极管激光、二极管泵浦固体激光)也可以用在布里渊雷达中，那么整个系统会变得更加紧凑、经济.因此，解析宽带激光产生的布里渊光谱很具吸引力.

本文利用数学方法解决通过F-P标准具和ICCD所获得的宽带布里渊光谱的解析问题.在解析过程中，首先用光谱信号构造出黎曼函数(Riemann Function，国外通常称Thomae's Function)的一个变体，再借助黎曼函数的某些特殊性质来解谱.黎曼函数长期以来被很少有实际的物理应用.自2011年起，对黎曼函数应用的研究取得突破，该函数先后在高通量基因测序 [24] 、人体T细胞免疫修复的组织 [25] 、三维显示技术的摩尔纹 [26] 、量子多体 [27] 、俄罗斯联邦选举作弊 [28] 中出现.尽管黎曼函数在自然科学领域的出现已经引起了很多物理学家的注意，但是关于这方面的报道仍较少，且不够系统，已经报道的案例大都是关于离散系统的.数学研究表明，虽然这些离散系统的情况各不一样，但是关于其中的黎曼函数的应用存在一个一般性的概率论解释：它可以被认为是一个稀有事件几率的分布函数 [29] .和已有报道不同，本文试图把黎曼函数应用到一个连续统中，这个系统和概率、统计完全无关.

1 基本原理

布里渊雷达有两种，一种是收集自发布里渊散射；另一种是收集受激布里渊散射(Stimulated Brillouin Scattering, SBS).若检测的是自发布里渊散射，其光谱图中会同时出现瑞利散射、布里渊散射的Stokes线和反Stokes线，一共三个谱峰，这三个峰之间存在确定的强度关系.而当泵浦光强度超过阈值时，则发生受激布里渊散射.此时反Stokes线很弱，乃至消失；只剩下泵浦光的界面反射光和SBS的Stokes线，且这种方法信噪比更高 [20] .本研究只针对受激布里渊散射的情形. 图 1 是两张受激布里渊散射的光谱图.这种光谱是由F-P标准具拍摄得到的.标准具有很高的分辨率，但是自由光谱程(Free Spectral Range, FSR)很窄.为了方便起见，本文兼用FSR作为光谱图 x 轴的单位.由于标准具产生的光谱是以FSR为周期的周期函数，图中只需要画出一个区间[0, 1]即可.

图 1. 窄带和宽带激光泵浦所得的典型布里渊光谱及相应的干涉图样

Fig. 1. Typical images of Brillouin spectra that pumped with narrow linewidth laser and broadband laser, and the corresponding fringe patterns

下载图片查看所有图片

由图可知，在 图 1(a) 中泵浦光频率、布里渊光频率和相应的频移很容易分辨，但是 图 1(b) 中的曲线看起来却像是随机曲线.这是因为图中的泵浦光和布里渊光都包含很多个纵模.利用F-P标准具进行分光时，由于其自由光谱程不足，多个纵模必然会重叠到一个自由光谱程内，就形成了一个叠加的宽峰.

在宽峰的内部，泵浦光和SBS的纵模存在一一对应的关系，每一个泵浦光纵模激发一个SBS的Stokes峰.实际上，这些SBS的峰不但在频率上是和泵浦光纵模一一对应，在强度上也是成线性关系的.文献[30]中推导了经典框架下受激布里渊的小信号增益表达式

式中， E 2 是泵浦光的幅值， γ 是压电效应系数， β 和 β s 都是唯象的衰减常数， η 是介质弹性振动的弛豫系数， ε 是介电常数， ρ 是介质密度， v s 是声速.

可以发现，随着在泵浦光电场强度的增长可以分多个阶段：

1) 泵浦光小于阈值时， g \approx0.

2) 泵浦光大于阈值，但仍然满足小信号近似，则后向SBS光呈现非线性放大.增益系数为$ \sqrt {\frac{{{\gamma ^2}{k_1}{k_s}{{\left| {{E_2}} \right|}^2}}}{{8\rho \upsilon _s^2\varepsilon }}} \propto \left| {{E_2}} \right|$，即增益g正比于泵浦光电场强度.而SBS光和增益系数为指数函数关系

3) 泵浦光远大于阈值时，无法满足小信号近似，出现增益饱和，此时SBS光强和泵浦光光强近似呈现线性关系.此时，对于固定的介质和实验装置可以得出一个最佳的SBS能量转化率.

4) 当泵浦光强进一步增长时，SBS由于SRS、激光诱导击穿、克尔透镜等过程的竞争被削弱，斜率效率下降.

以上四个阶段可以参考文献[31, 32]中的仿真和实验结果.文献[32]报道了水中宽带布里渊散射(Wideband Stimulated Brillouin Scattering, WSBS)和前向喇曼散射(Forward Stimulated Raman Scattering)的实验结果.从该文献中的实验结果可知：水的WSBS阈值不高，非线性增长阶段极短；在达到阈值之后较宽的能量范围(可达阈值能量的几十倍)，WSBS的能量都和泵浦光能量成正比.宽带布里渊激光雷达正是工作在这个线性阶段，因此每一个纵模，SBS光的能量都和泵浦光正比；故SBS光谱函数曲线的整体形状和泵浦光相同，仅仅幅值有差异.在实际应用中，通常SBS反射率比界面反射率要大很多，故通过一片可调石英片增强泵浦光的反射，可以把SBS和泵浦光的平均幅值调节到几乎相同.这样，SBS的光谱可以看作泵浦光光谱平移得到的；若泵浦光光谱为 f (x)，则SBS光谱为 f (x + b)，其中 b 是布里渊频移(Brillouin Frequency Shift, BFS).

通过这种方法得到的光谱曲线是泵浦光光谱和SBS光谱之和，这一点和窄带SBS光谱是相同的.但是，由于两个光谱都是宽光谱，布里渊频移本身无法直接从中读出来.本文研究的目标是，从这个看起来像是随机信号的曲线中识别出布里渊频移.为了让表述更加清晰，令泵浦光光谱为 f (x)，SBS光谱为 f (x + b)，标准具采集到的混合光谱为 g (x).研究目标可以用数学语言写为

问题：从一个已知函数

求解函数 f (x) (周期= 1)和相应的平移量 b .

从信息论的角度来看，这个问题需要求解的信息量大于已知条件的信息量： f (x)和 g (x)中的信息量是相同的，解中的平移量 b 是一个多出来的额外信息，因此这个问题是不可能解决的.然而，需要注意的是，这个问题并不是一个纯粹的数学问题.如果多考虑一些 f (x)的物理性质，求解这个额外信息 b 是可能的.

由于 f (x)是一个周期函数，式(3)可以展开为傅里叶级数

式中， F n 是 f (x)的傅里叶系数. g (x)作为周期函数也可以这样展开，因此可以建立关于这两个傅里叶系数 G n 和 F n 之间的关系

为了防止0作为分母出现，引入一个小的正数 ε .

式中右边部分，只有位移 b 是未知的.不妨把 b 看作是一个尝试值 b trial ，把 g (x)中的真实的位移记作 b real .这样，还未解出的傅里叶系数 F n 就可以看作是以尝试值 b trial 为宗量的函数.为了简化式(6)，将其重写为

其中定义

式(7)意味着，对于每一个尝试值 b trial 都会有一套对应的傅里叶系数{F n}.如果把 b trial 遍历所有可能的数值，会得到一个包含多套傅里叶系数{F n}的大集合.在这个集合中，只有最接近真实情况的一套系数应该被识别并挑选出来.在这一步，问题已经转化为：寻找一个有效的方法来评价{F n}是否合理.

为了找到这个方法，仔细分析式(7).考虑到 ε \to0 + ，如果有exp(－i2π nb trial)=－1，X n (b trial)的值会发散到1/ ε \to\infty.这个发散条件exp(－i2π nb trial)=－1可以改写为

如果把n从－N/2遍历到N/2，符合条件n=kp/2并且k∈odd的n每隔p才会出现一次.利用“地板函数” $\left\lfloor x \right\rfloor $表示小于x的最大整数，可以把符合条件的n的数量表达为$ 2 \times \left\lfloor {(N - p)/2p + 1} \right\rfloor $.另一方面，当ε→0+时，只有发散的部分需要在求和时考虑，其他都可以作为小量忽略不计了.利用这两个结论，可以导出和式

其极限为

对于${b_{{\text{trial}}}} \in \mathbb{R} $，式(11)的完整表达形式为

把式(12)和黎曼函数的原始形式

进行比较，可以发现它们非常相似.仅有的区别是，当分母 p 是奇数，黎曼函数得到的是1/ p ，而式(12)得到的是0.这两个函数的曲线如 图 2 ，为了清晰显示，图中用圆点标注分母较小的有理数.注意到，分母是奇数的点，例如1/3、2/3、1/5、2/5，在左边的子图里都没有出现.很显然，两个曲线有类似分形的结构.

图 2. 式(12)的函数曲线和黎曼函数的曲线

Fig. 2. Curves of Eq. (12) and Thomae's function

下载图片查看所有图片

现在已经揭示了\sum X n (b trial)和黎曼函数拥有相似的形式.利用这个结论来考虑对式(7)中的 n 求和.如果 G n 恒不为0，这个求和\sum G n X n (b trial)就会在所有{b trial | b trial = q / p, p \in odd}上发散. b trial 只是一个尝试值，真正的位移是 b real .利用 b real 可以重新展开式(5)

式中，没有解出的“真实” F n 记作 F n, real ，来区分真实值 F n, real 和从尝试值 b trial 解出的 F n (b trial).将式(14)代入\sum| F n (b trial)|= \sum| G n X n (b trial)|，可以得到

显然，如果选择 b real 作为尝试值，那么使得 X n (b trial)发散的 n 也会使得(1+e －i2π nb real)变成0.这是因为两者都要求 n 满足exp(－i2π nb)=－1.因此，整个求和式(15)可以视为 b trial 的函数，并且该函数除了在 b trial = b real 时理论上为0，其他地方都随着 ε \to0 +, N\to\infty趋向于发散.这里的论述中，已经考虑到计算机在数值遍历的过程中只能取有理数.虽然 X n (b trial)在无理点不发散，但实际上计算机是遍历不到的.

利用这样一个特殊的性质，可以找到一种可能途径来解决本节开头提出的问题.首先，计算g(x)的傅里叶系数Gn；然后用式(7)来计算$\sum\limits_{\begin{array}{*{20}{l}} {n{\text{ }} = {\text{ }} - N/2} \end{array}}^{\begin{array}{*{20}{l}} {N/2} \end{array}} {\left| {{F_n}({b_{{\text{trial}}}}){\text{ }}} \right|{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} } $，把它看作一个关于btrial的函数；画出这个函数在$ {b_{{\text{trial}}}} \in \mathbb{Q}$上的所有点，连成曲线并寻找最小值.寻找最小值的工作可以限定在一个合理的范围进行.如，对海水可以在7~9GHz的范围内寻找极值点，以减少计算量.

2 实验研究

2.1 测试数据集

为了测试算法的准确度，用一系列测试数据进行验证.测试数据是由真实的宽带激光光谱中生成.采用的激光器是Nd:YAG调Q激光(光谱物理公司Quanta-Ray 290).实验中使用的波长为532 nm，单脉冲平均能量为38 mJ，脉冲半高宽为12 ns.在无种子注入的工作模式下，激光器的标称线宽为1 cm －1 .激光束散射到一个固体标准具(CVI ET-25.4-4.00-UV)上，然后用一个镜头系统(Takumar 1:4/150)聚焦到一个CCD相机(PointGrey GRAS-20S4M-C, 1/1.8”)上.整个实验装置如 图 3 .由于数据传输速率的限制，在1 600 pixel\times1 200 pixel的分辨率下，CCD相机的采集速率最高只能到7.5 Hz，而激光器的重复频率为30 Hz.因此，用分频器把相机的触发频率调整到6 Hz，用Stanford DG-645时间延迟器来同步激光脉冲和CCD快门的触发时间.标准具标称的透过波前误差是 λ /10，其中 λ =632.8 nm，标称楔角为 α < 1″.标准具基片采用紫外石英制作，折射率为1.461，双面覆盖99.5%高反射率介质膜，中心波长532 nm，理论上的自由光谱程是25.67 GHz.以上所述的光谱数据采集方式，也是实际的布里渊激光雷达所常用的.

图 3. 实验装置

Fig. 3. Experimental configuration

下载图片查看所有图片

CCD相机拍摄到的是二维衍射图样，即一系列同心圆环.用数据折叠(data folding)和除去噪声的方法处理 [33 - 34] ，可以获得高分辨率一维光谱曲线.因为激光器不是在单纵模状态运行，每一发脉冲在各个纵模上的能量分布都是不相同的，所以拍摄的光谱看起来都像是低频的随机信号.为了构造具有各种不同频移的布里渊-泵浦光混合光谱，把光谱先做平移(离散化处理之后实际上是圆移(circular shift))，然后再和自身混合.该处理过程的流程如 图 4 . 图 4(a) 中样本1是以宽带激光光源照射固体标准具拍摄的宽带激光源的光谱， 图 4(b) 中的样本2是样本1的圆移.样本1和样本2之和绘制在 图 4(c) 的图上. 图 4(c) 中的曲线来自于光谱的平移、叠加，可以用作测试数据.不断重复这个流程可以产生一系列具有不同 b real 的测试数据.在本次研究中，对200个不同的激光脉冲进行摄谱，并分别加上8个不同的频移： b real =0.28, 0.29, 0.30, 0.31, 0.32, 0.33, 0.34, 0.35，这里频移的单位为FSR.这样，一共是8\times200=1 600个测试数据.

图 4. 测试数据准备流程

Fig. 4. The synthesis of test data

下载图片查看所有图片

2.2 方法与结果

对于测试数据中的每一个光谱g(x)，计算相应的一整套傅里叶系数Gn.利用式(7)来计算$ \sum {\left| {{F_n}({b_{{\text{trial}}}})} \right|} $

式中，N和ε的数值分别选择为104和0.1，尝试值btrial是一个从0.0到0.999 9的数列，步长为0.000 1.图 5为一个$ \sum {\left| {{F_n}({b_{{\text{trial}}}})} \right|} $的典型曲线，所用的测试数据位移为breal=0.35.其中分母较小的有理点用小圆圈做了标记.breal的位置在曲线上用交叉线来标记，这个位置附近的曲线细节呈现在右上角的插图中.显然，在区间[0.1, 0.5]内，曲线在0.35处达到最低点，这刚好精确等于真实的位移breal.其他测试数据得到的∑|Fn|曲线也表现出相似的特征.

图 5. 典型的∑|F_n|关于b_trial的函数曲线

Fig. 5. Typical ∑|F_n|to b_trial curve

下载图片查看所有图片

对所有的测试数据求出了\sum| F n |曲线在[0.1, 0.5]区间的最低点.该最低点的横坐标 b min 就是算法识别出的SBS频移.若 b min 和真实频移 b real 之差小于0.000 4个FSR(约0.01 GHz)，可认为结果准确.故“识别准确率”可定义为：归类为“正确”(| b min － b real | < 0.0004)的样本占样本总数的比例.

对每个不同的 b real 均使用200个样本进行统计，统计结果显示：在没有添加噪声的测试数据中，除了 b real =0.29和 b real =0.31识别的准确率是99.0%，其他情况的识别准确率都是100.0%.对于标准具来说，该判据的精度已经达到或超过标准具加工误差带来的精度限制.

3 讨论

3.1 额外信息的来源

从信息论的角度来理解， g (x)和 f (x)所含有的信息量是相同的，因此不可能从混合信号 g (x)中还原出原始信号 f (x)，还能得到额外的信息——位移 b real .但本文算法做到了这一点，因此必然用到了不包含在数学问题中的额外信息，否则就违反了信息论.

事实上，当选择\sum| F n |的最低点作为正确的 b trial 时，大量可能的 f (x)和 b trial 都被舍弃了.舍弃原因如 图 6 所示， 图 6(a) 、 (b) 中，经过平移、叠加，合成了一幅位移为0.3的测试数据 g (x). 图 6(c) 是以正确的尝试值从 图 6(b) 中的 g (x)恢复得到的原函数，其形状和 图 6(a) 中的 f (x)基本相同，只是多了一些噪声.而采用不正确的尝试值的 b trial 所还原出的原函数曲线 f (x)如 图 6(d) 所示，出现了幅度高达1011的高频振荡，这完全不符合通常光谱曲线的特性.正常的光谱曲线应当是比较光滑的，其频域空间低频成分比较多，高频成分很少，不会出现这样大的高频振荡.从这个角度理解，“求\sum| F n (b trial)|的最小值”这一过程实际上就是对 b trial 的不断筛选，筛选出的结果应当使得 f (x)的高频成分幅值最小.

图 6. 不同b_trial恢复出的原函数f(x)

Fig. 6. The f(x) recovered from different b_trial

下载图片查看所有图片

在这个筛选过程中，对 f (x)存在一个物理上的预期：它应当是光滑的，有一定的整体起伏波动，噪声毛刺很低，没有大的快速振荡.换句话说，频域就是 f (x)的高频成分幅值应当很小.这种预期是来源于物理学，而不在第1节所提出的数学问题之内，所以它就是“额外的”信息来源.

3.2 函数∑|F_n|曲线的更多细节

在\sum| F n |曲线的右边半部分(见 图 5)，还有一个最低点 b trial =0.65.这个点恰好和左边半部分的最低点 b trial =0.35是对称的.这一点可以用圆移的性质来解释.新设一个信号 h (x)，定义为 f (x +0.35)，容易推导出 g (x)= f (x)+ f (x +0.35)= h (x +0.65)+ h (x).这样，问题所求的位移、原函数就存在两对解：第一对{0.35, f (x)}，第二对{0.65, h (x)}.实际上，可以证明整个曲线\sum| F n |都是关于 b trial =1/2反射对称的

另外， 图 5 中有一些分母比较小的有理点用小圆圈做了标记；它们分别是局部的极大值和局部的极小值.正如第1节所计算的，当 b trial 是一个分母为小偶数的有理数，\sum| F n |的值会趋向发散，因此会产生局部极大值；但是，图中除了能看到极大值，在分母为小奇数的有理数上还存在局部极小值.这一点需要进一步的解释.发散条件式(9)要求 b trial = q / p 的分母 p 为偶数，此时在 n 取某些值的情况下会出现 nb trial = k /2，导致\sum| F n |发散.若分母 p 为奇数，无论 n 取何数值也不可能出现 nb trial = k /2.在这个条件下，| F n |取最大值时 nb trial － k /2需要取最小值.

以 b trial =1/3处为例， nb trial = n /3则 nb trial － k /2= n /3－ k /2能取到的最小值为1/6，要求是 n =(3/2) k \pm1/2.此时的| F n |为最大值| F n |=|1－(1－ ε)exp(1/6)| －1 .类似地，对于其他的 b trial = q / p (p 为奇数)，可以推导出 nb trial － k /2的最小值1/2 p ，而| F n |为最大值| F n |=|1－(1－ ε)exp(1/2 p)| －1 .这个表达式对 p 来说是单调递增函数.综上，对于分母是奇数的有理数，分母越小，| F n |的最大值越小.这可以解释为什么在分母为小奇数的有理点上会经常出现局部极小值.

3.3 高斯白噪声的影响

为了测试噪声的影响，在所有的测试数据上又加入了高斯白噪声，白噪声的范围从－40 dB到－20 dB.由于所加噪声跨越三个数量级，这里以对数单位分贝(dB)来衡量噪声的幅值.具体换算式为$N\left( {{\text{dB}}} \right){\text{ }} = 10{\text{log}}\frac{{{E_{{\text{Noise}}}}}}{{{E_{{\text{signal}}}}}}{\text{dB}} $，其中ENoise是噪声的能量，而Esignal是信号的能量.信噪比S/N和N(dB)有换算关系$S/N{\text{ }} = \sqrt {{{10}^{ - N\left( {{\text{dB}}} \right)/10}}} $.

加入白噪声后用算法处理的结果如 图 7(a) .结果表明，本文所研究的处理方法对噪声是比较敏感的：尽管当噪声 < －34 dB时，可以维持准确率在一个很高的水平，但是当噪声增长到大约－30 dB之后，准确率就急剧下降.为了提高这种处理方法在有噪声污染时的表现，本文构造了一个高斯滤波器(标准差 σ =0.02)来过滤噪声信号，然后再计算\sum| F n (b trial)|.在过滤信号之后，对大多数噪声>－30 dB的区域准确率有明显提升，如 图 7(c) .

图 7. 被高斯白噪声污染后算法的准确率、污染又滤波处理后算法的准确率，以及相应信号的频谱图

Fig. 7. The accuracies when white Gaussian noise is mixed with the original signal, the accuracies for the same signals when a Gaussian filter is applied during processing, and the corresponding frequency charts of these signals

下载图片查看所有图片

为了解释 图 7(a) 中对高斯白噪声的敏感性，这里提出一个简要的理论分析.和白噪声混合之后， F n 的表达式变为

式中， N n 是噪声的傅里叶系数.根据高斯白噪声的性质， N n 表现出近似于均匀分布的常数功率谱；然而 G n 作为一个实际信号的频谱，随着 n 的增长下降很快.因此，除了最开始的低频段， N n 都要比 G n 大很多.相应的频谱图如 图 7(b) .另一方面，固定一个 b trial = q / p (p 是偶数)，使得 X n (b trial)发散的 n 在数轴上有一个分布.这个分布是均匀的，固定间隔为 p ，并且 X n (b trial)发散处的大小也不会随 n 变化.结合这两个前提，即使噪声很小，但是由于在高频区域其功率谱比 G n 大，还乘上了发散量，还是会显著影响整个\sum| F n |曲线.考虑到噪声的频谱接近于常数，当噪声远大于信号时，函数曲线\sum| F n |的形状会和 图 2(a) 很相似.

未经污染信号和噪声污染信号的频谱如 图 7(b) .比较这两个信号的频谱，很容易得到启示：一个低通滤波器可能会有助于求解噪声污染信号的位移.滤波后两个信号的高频部分都被削去一些，如 图 7(d) ；滤波后算法的准确率如 图 7(c) .从 图 7(c) 可以看出当噪声水平高于－30 dB时，实验能够验证这个假设，但是在较低的噪声水平下，滤波之后准确率反而有所下降；尤其是位移靠近0.33处.这种准确率的下降，大部分可以归因于滤波后高频信息的损失.在0.33附近的显著变化是因为曲线在1/3处本来就有一个极小值，容易造成干扰.经检查发现，在 b real =0.33时，判断错误的测试数据，极小值基本都落在0.333 3处，证实了本节的分析判断.

3.4 泵浦光和SBS光幅值之比的影响

本研究中提出的算法要求泵浦光光谱和布里渊光光谱的幅值为严格1:1的关系，但是实验中存在偏差.实验中(装置如 图 8)，激光通过焦距为50 mm的透镜系统聚焦到布里渊池中，在聚焦透镜前放置1/4波片和格兰棱镜作为光隔离器.产生的后向SBS光和池壁的反射光等通过格兰棱镜分光后，先后经过毛玻璃散射、标准具分光、镜头会聚，最后被CCD采集.为了保证两种光进入CCD的比例为1:1，在泵浦光聚焦之前的位置增加了一片石英片(quartz plate)来反射一部分泵浦光进入CCD，并预先开发程序实时监测CCD像平面的总响应.若在SBS cell的入射窗口处用泡沫遮断，CCD采集到的只有泵浦光；反之则采集到的是泵浦光和布里渊光的总和.实验时需要调整石英片角度，使得“遮断/未遮断”两种情况下CCD总响应为1:2，则此时泵浦光光谱和布里渊光光谱的幅值比为1:1.

图 8. 测试泵浦光和SBS光幅值之比所用的实验装置

Fig. 8. Experimental configuration to measure the ratio of pump light to SBS light

下载图片查看所有图片

为了研究幅值之比的偏差造成的影响，构造了泵浦光光谱和布里渊光光谱的幅值比 k \neq1的一系列测试数据.对测试数据的算法验证结果如 图 9 .图中的数据表明，总体上，本研究提出的算法在 k < 1.05时有很高的准确率，在1.05 < k < 1.08时准确率也比较高；而当 k >1.08时基本失去作用.

图 9. 泵浦光和SBS光幅值之比不为1时算法的准确率

Fig. 9. The accuracies when ratio of pump light to SBS light is not 1

下载图片查看所有图片

在实际实验中发现，由于SBS过程本身存在不稳定性，幅值比 k 起伏变化较大.例如，在35 ℃的纯水中进行实验，泵浦光和SBS光各自采集200发脉冲，其CCD总响应的统计情况如 图 10 .泵浦光的均方差稳定性(RMS stability)为2.44%，而SBS光脉冲的均方差稳定性为24.77%.显然，在这种能量稳定性下，目前的算法还不足以满足实际需要，需要改进.

图 10. 200发激光脉冲中泵浦光和布里渊光能量稳定性

Fig. 10. The stability of pump laser and SBS singal of 200 laser pulses

下载图片查看所有图片

SBS的能量稳定性比泵浦光要差很多，可能的原因有：1)泵浦光能量在阈值附近导致不稳定.苏红等 [35] 进行的理论和实验研究表明，在 f =50 mm的透镜汇聚下，纯水SBS阈值为30 mJ.而本实验中使用的泵浦能量仅为38 mJ，超过阈值不多；此时泵浦光即便有微弱的能量起伏，反映到SBS能量上也会有较大变化.若将来实验中采用单脉冲能量更大的激光，该问题可能会得到改善. 2)激光诱导击穿(Laser Induced Breakdown, LIB)竞争导致SBS能量不稳定.水是高度极性的溶剂，击穿阈值较低，在泵浦光聚焦作用下容易击穿.水的电离本身会消耗一部分泵浦光的能量，另外LIB产生的等离子体对泵浦光的吸收也很明显.LIB产生等离子体的形状、大小具有一定的随机因素，每次都不完全相同，因此对SBS过程的能量竞争具有一定随机性，导致了SBS能量不稳定.目前已知存在一些性能稳定的介质(如SF 6 、全氟碳化物等)不容易产生LIB；另外，提高水的纯度也有助于防止LIB.但是这些措施只能用于实验室环境.SBS光谱探测的目的是用于自然水体，在防止电离这方面还需要进一步探索.

4 结论

本文提出了一种用于处理宽带布里渊光谱信号的解析算法.该算法针对的是用F-P标准具和CCD采集到的泵浦光-布里渊光混合光谱信号.首先用该信号构造出黎曼函数的一个变体，再基于黎曼函数的某些性质可以用于从混合信号 g (x)= f (x)+ f (x + b)中还原出泵浦光信号 f (x)和布里渊频移 b .为了验证该算法，从无种子注入的多纵模Q调制固体激光器上采集光谱信号，并在此基础上构造出各种频移的测试数据.无噪声的理想情况下，算法的准确率接近于100%.为了测试非理想情况下的表现，又构造了包含高斯白噪声的测试数据集及泵浦光和SBS光幅值之比不为1的测试数据集.结果表明，只有当噪声小于－30 dB，并且泵浦光和SBS光幅值之比小于1.05时，算法才能得到准确结果.另一方面，和其他关于黎曼函数应用的研究工作相比，本文中的黎曼函数完全是从连续统中自然出现的，并不涉及概率和统计.这也为黎曼函数的应用研究拓宽了思路.

本文构造的测试数据验证了该处理方法的可行性、准确率等基本要素.但是，对于未来实际外场实验而言，该算法对信噪比、能量稳定性的要求比较苛刻，目前还无法满足实际应用要求，还需要进一步的改进.

参考文献

[1] SHI Jin-wei, 石锦卫, 弓文平, GONG Wen-ping, BAI Jian-hui, 白建辉. 布里渊散射及其在激光雷达中的应用[J]. 物理, 2007, 36(10): 777-782.

[2] SHI J, OUYANG M, GONG W. A Brillouin lidar system using F-P etalon and ICCD for remote sensing of the ocean[J]. Applied Physics B, 2008, 90(3-4): 569-571.

[3] REN Xiu-yun, TIAN Zhao-shuo, ZHANG Yan-chao. Theoretical and experimental investigations on measuring underwater temperature by the coherent Brillouin scattering method[J]. Applied Optics, 2015, 54(30): 9025-9029.

[4] LIU Da-he, XU Jian-feng, LI Rong-sheng. Measurements of sound speed in the water by Brillouin scattering using pulsed Nd:YAG laser[J]. Optics Communications, 2002, 203(3): 335-340.

[5] HE Xing-dao, WEI Hong-jun, SHI Jiu-lin. Experimental measurement of bulk viscosity of water based on stimulated Brillouin scattering[J]. Optics Communications, 2012, 285(20): 4120-4124.

[6] , KUMER J, NIGHTINGALE R. Validation of CH₄ and N₂O measurements by the cryogenic limb array etalon spectrometer instrument on the upper atmosphere research satellite[J]. Journal of Geophysical Research:Atmospheres, 1996, 101(D6): 9679-9710.

[7] KUMMROW A, MENG H. Pressure dependence of stimulated Brillouin backscattering in gases[J]. Optics Communications, 1991, 83(5-6): 342-348.

[8] BRIGNON A, HUIGNARD J P. Phase conjugate laser optics:volume 9[M]. New Jersey:John Wiley & Sons, 2003.

[9] DAMZEN M, VLAD V, MOCOFANESCU A, et al. Stimulated Brillouin scattering:fundamentals and applications[M]. CRC press, 2003.

[10] BAI Zhen-xu, YUAN Hang, LIU Zhao-hong. Stimulated Brillouin scattering materials, experimental design and applications:A review[J]. Optical Materials, 2018, 75: 626-645.

[11] LIANG Kun, MA Yong, HUANG Jun. Precise measurement of Brillouin scattering spectrum in the ocean using F-P etalon and ICCD[J]. Applied Physics B, 2011, 105(2): 421.

[12] SCHORSTEIN K, POPESCU A, GÖBEL M. Remote water temperature measurements based on Brillouin scattering with a frequency doubled pulsed Yb:doped fiber amplifier[J]. Sensors, 2008, 8(9): 5820-5831.

[13] LIU Da-he, 刘大禾, KATZ J, KATZ J. 水中布里渊散射的边缘探测方法[J]. 中国激光, 1999, 26(4): 307-311.

[14] REN Xiao-bin, 任晓斌, 弓文平, GONG Wen-ping, DAI Rui, 戴瑞. 用碘分子共振吸收滤波器精确测量水中的布里渊散射频移[J]. 光学学报, 2004, 24(10): 1429-1433.

[15] BREWERR G. Stimulated Brillouin shifts by optical beats[J]. Applied Physics Letters, 1966, 9(1): 51-53.

[16] FRIMAN S. Field test of a Brillouin LIDAR f temperature profiles of the ocean[D]. Stockholm, Sweden: KTH Royal Institute of Technology. 2016.http:kth.divaptal.gsmashrecd.jsfpid=diva2%3A923913&dswid=8131

[17] LI Wen-hai, BAO Xiao-yi, LI Yun. Differential pulse-width pair BOTDA for high spatial resolution sensing[J]. Optics Express, 2008, 16(26): 21616-21625.

[18] ZHOU Deng-wang, 周登望, 王本章, WANG Ben-zhang, BA De-xin, 巴德欣. 用于动态应变测量的快速分布式布里渊光纤传感[J]. 光学学报, 2018, 38(3): 0328005.

[19] ZHANG Li-xin, 张立欣, 李永倩, LI Yong-qian, AN Qi, 安琪. 脉冲编码瑞利布里渊光时域分析温度传感技术[J]. 光学学报, 2017, 37(11): 68-75.

[20] DAI Rui, 戴瑞, 弓文平, GONG Wen-ping, SHI Jin-wei., 石锦卫. 基于布里渊散射的水下目标探测[J]. 激光与光电子学进展, 2008, 45(2): 65-69.

[21] BRIVIO S, POLLI D, CRESPI A. Observation of anomalous acoustic phonon dispersion in SrTiO₃ by broadband stimulated Brillouin scattering[J]. Applied Physics Letters, 2011, 98(21): 211907.

[22] LI Wei, HE Bin, ZHANG Chun-feng. Coherent acoustic phonons in YBa₂Cu₃O₇/La₁/3Ca₂/3MnO₃ superlattices[J]. Applied Physics Letters, 2016, 108(13): 132601.

[23] BOJAHR A, HERZOG M, MITZSCHERLING S. Brillouin scattering of visible and hard X-ray photons from optically synthesized phonon wavepackets[J]. Optics Express, 2013, 21(18): 21188-21197.

[24] TRIFONOV V, PASQUALUCCI L, DALLA-FAVERA R. Fractal-like distributions over the rational numbers in high-throughput biological and clinical data[J]. Scientific Reports, 2011, 1: 191.

[25] MEIER J, ROBERTS C, AVENT K. Fractal organization of the human T cell repertoire in health and after stem cell transplantation[J]. Biology of Blood and Marrow Transplantation, 2013, 19(3): 366-377.

[26] SAVELJEV V, KIM S K. Probability of the moiré effect in barrier and lenticular auto stereoscopic 3D displays[J]. Optics Express, 2015, 23(20): 25597-25607.

[27] LUNDHOLM D. Many-anyon trial states[J]. Physics Review A, 2017, 96(1): 012116.

[28] ROZENAS A. Detecting election fraud from irregularities in vote-share distributions[J]. Political Analysis, 2017, 25(1): 41-56.

[29] NECHAEV S. Non-euclidean geometry in nature[M]. arXiv preprint arXiv:1702.06757, 2017:61-111.

[30] YARIV A. Quantum electronics[M]. Shanghai:Shanghai Science and Technology Press, 1983.

[31] WANG Yulei. Investigation on a 100J laser system highenergy high power stimulated Brillouin scattering phase conjugation[D]. Harbin: Harbin Institute of Technology, 2007.http:cdmd.cnki.com.cnArticleCDMD102132008193233.htm王雨雷.百焦耳激光器及高能高功率受激布里渊散射相位共轭研究[D].哈尔滨: 哈尔滨工业大学, 2007.

[32] SHI Jiu-lin, CHEN Wei, MO Xiao-feng. Experimental investigation on the competition between wideband stimulated Brillouin scattering and forward stimulated Raman scattering in water[J]. Optics Letters, 2012, 37(14): 2988-2990.

[33] XIA Xu-sheng, YUAN Hong, LIU Jin-bo. Method to improve the resolution of a non-parallel Fabry-Perot etalon[J]. Applied Optics, 2018, 57(29): 8757-8765.

[34] HUANG Jun, MA Yong, ZHOU Bo. Processing method of spectral measurement using F-P etalon and ICCD[J]. Optics Express, 2012, 20(17): 18568-18578.

[35] SU Hong, 苏红, 张为俊, ZHANG Wei-jun, 高晓明, GAO Xiao-ming. 受激布里渊散射阈值性能的研究[J]. 量子电子学报, 2002, 19(3): 210-213.

夏栩笙, 蔡向龙, 李仲慧, 沈陈诚, 刘万发, 金玉奇, 桑凤亭, 郭敬为. 基于黎曼函数的宽带布里渊光谱信号解析算法[J]. 光子学报, 2020, 49(6): 0630003. Xu-sheng XIA, Xiang-long CAI, Zhong-hui LI, Chen-cheng SHEN, Wan-fa LIU, Yu-qi JIN, Feng-ting SANG, Jing-wei GUO. Algorithm for Broadband Brillouin Spectra Analysis Based on Thomae's Function[J]. ACTA PHOTONICA SINICA, 2020, 49(6): 0630003.