自适应最优核时频分析与CNN在气液两相流流型识别中的应用

简介：一份基于自适应最优核（AOK）与卷积神经网络（CNN）相结合的气液两相流流型识别研究论文，面向流体力学、测控技术与人工智能交叉领域的研究人员和研究生，重点解决传统流型识别中特征提取代表性不足、识别率不高的问题。压缩包内共1个PDF文件，大小约3.21MB，即《化工学报》2018年发表的期刊论文全文，包含中文摘要、英文摘要、方法原理、实验装置、结果分析与参考文献。目前已有179人学习下载。读者可从中获得完整的研究思路：利用V锥流量计和动态差压传感器采集波动信号，通过自适应最优核算法生成时频谱图，再由CNN自动提取特征并用Softmax分类器实现流型识别；论文对常见流型进行了试验分析，并讨论了方法在油气井测井、流体机械设计等领域的应用前景，对从事两相流检测和深度学习应用的研究者具有直接参考价值。

1. 从一维波动到二维图像，流型识别为什么绕不开时频变换

气液两相流系统里最常见的烦恼不是测不准流量，而是你根本不知道管道里此刻是什么状态。泡状流、弹状流、环状流的压差信号和电导信号在统计特征上往往只差几个百分点的方差，现场用阈值逻辑判断流型，工况一切换就频繁误报。把信号直接喂给卷积神经网络，看似绕开了特征工程，但一维时序波形里其实塞满了噪声和测量伪迹，卷积核在原始采样点上什么也学不出来。

这类问题的标准解法是换一个视角看数据：先用时频分析把一段非平稳信号摊成二维时频平面，再把这张图交给CNN去分类。这里的关键在于时频图的质量直接决定识别上限。短时傅里叶变换窗口一短频率糊，窗口一长时间轴又失去局部性；Wigner-Ville分布分辨率虽高却被交叉项污染得一塌糊涂。自适应最优核方法恰好卡在两者之间——它在模糊域里用一个随信号旋转的高斯核把交叉项滤掉，同时保住自项的时频聚集性。整个逻辑链就是：信号 → 模糊域 → 自适应最优核 → 时频分布 → CNN分类。这篇就按这个链路走一遍，从核函数的数学原理讲到能落地的PyTorch实现，再把参数调整和现场验证的坑逐一说明。

2. 自适应最优核为什么比固定核更适合气液两相流信号

2.1 固定核函数的局限：STFT与WVD的两种尴尬

时频分析的核心矛盾是海森堡不确定性原理：你不可能同时获得任意高的时间分辨率和频率分辨率。短时傅里叶变换用固定窗长做局部截断，窗口选定后整个时频平面上的分辨率就锁死了；一旦气液两相流的表观流速发生变化，泡径和段塞频率随之漂移，固定窗长的适应性马上出问题。

Wigner-Ville分布（WVD）用双线性变换摆脱了窗长的束缚，对单分量线性调频信号可以获得理论上最尖锐的能量聚集。但问题在于，气液两相流传感器信号的能量分布在时频平面上通常有多个分量同时存在，而这些分量之间的交叉项会产生虚假的能量峰，幅值有时比真实成分还高。工程上最常见的错误是在WVD上直接取峰值当作流型特征，结果弹状流的段塞周期被交叉项干扰得完全不可信。

2.2 模糊域视角与核函数的旋转思想

解决交叉项问题的自然思路是把信号变换到模糊域（Ambiguity Domain），因为WVD的交叉项在模糊域里总是远离原点分布的，而自项倾向于集中在原点附近。于是问题变成了设计一个二维滤波器，让核函数在模糊域里只保留自项的分布区域，滤掉远处的交叉项。固定核方法比如平滑伪WVD，用固定的矩形或椭圆窗做低通滤波，但流型变化时交叉项出现的位置和方向也变了，固定核不可能覆盖所有工况。

自适应最优核的核心思想在于：核函数不是人为预先设计的，而是针对当前这一段信号的模糊域分布动态求解出来的。这个优化目标的数学形式是：在给定体积约束的前提下，让核函数与模糊域的乘积的二维傅里叶变换结果能量最大。直观理解就是，把核函数当成一个可以在模糊域里自由伸缩的覆盖层，让它在约束条件下尽量罩住信号自身的能量分布，同时自然地避让开零散分布的交叉项。

这个求解过程可以用径向高斯核来参数化，核函数沿模糊域角度方向的截面是高斯型曲线，高斯宽度由信号决定。弹状流和气液两相流中段塞边界的突变特征会造成瞬时频率跳变，自适应核能自动选择窄的核宽度来保持时频分辨率，而在平稳段自动变宽来压制噪声。这种性质是STFT的固定窗完全不具备的。

2.2.1 约束优化问题的数学表达与求解思路

工程实现上不需要读者从头推导凸优化理论，但需要理解求解过程的关键约束。常见做法的优化目标是最大化如下泛函：

[ \max_{\phi} \int_{0}^{2\pi}\int_{0}^{\infty} |P(r,\psi) \cdot \Phi(r,\psi)|^2 , r, dr, d\psi ]

其中(P(r,\psi))是极坐标下的模糊函数，(\Phi(r,\psi))就是要求解的径向高斯核。约束条件是(\Phi(r,\psi))的总体积不超过某个阈值(\alpha)，且(\Phi(r,\psi))的径向截面单调递减。阈值的物理意义是控制核函数覆盖范围与噪声抑制程度之间的平衡，通常取2到5之间，数值越小抑制交叉项越激进，但过大时真实信号的高频跳变也会被抹平，实际调参时(\alpha=3)是多数气液两相流信号的安全起点。

求解方法一般用两步迭代：先用当前信号模糊函数估算能量集中方向，再用带拉格朗日乘子的牛顿法更新核函数的径向高斯宽度。这个过程并不需要复杂的深度学习框架，用NumPy加SciPy的优化器就能在单段信号上完成计算，时间成本在毫秒到几十毫秒的范围内，完全能够支撑在线识别的需求。

2.3 气液两相流信号的时频特征与其核自适应行为

把电导传感器采集到的波动信号拿来分析，泡状流的特征是能量沿时间轴均匀散布、频率集中在低频段；弹状流则表现为周期性的能量带——每个液段通过时产生一次低频高幅值脉冲；环状流的能量集中在气芯携带液滴形成的高频扰动上。这些差异在模糊域中的表现非常不同：弹状流的自项分布在沿着径向的长条形区域，环状流的自项则分布更广且分散。

自适应最优核在这几种形态之间表现出的行为差异，恰好就是下一步CNN分类的判别基础。核函数在弹状流信号上会被拉成细长形，在环状流信号上则趋于椭圆形，这个核的形态本身就携带着流型信息。

3. 把自适应最优核的数学求解写成可运行的Python代码

3.1 信号预处理与模糊函数的数值计算

实际项目中用的传感器采样率常见是1kHz到10kHz这个范围，分析窗长取512或1024个采样点比较合适。窗长过长会让段塞的瞬态特征被平均掉，过短则频率分辨率不够。首先把原始信号做零均值化和幅值归一化，然后用一个汉明窗截取分析片段：

import numpy as np from scipy.signal import windows, hilbert def preprocess_segment(raw, win_len=1024): """ 截取分析片段并生成解析信号 raw: 原始采样序列 win_len: 分析窗长度，必须为2的幂 """ seg = raw[:win_len].copy() # 去趋势：减去线性拟合的均值趋势，避免后续时频平面出现虚假的零频成分 t = np.arange(win_len) coef = np.polyfit(t, seg, 1) seg = seg - np.polyval(coef, t) # 幅值归一化，防止不同工况下信号能量差异主导分类结果 seg = seg / (np.std(seg) + 1e-8) # 用解析信号计算模糊函数，抑制实信号负频率成分带来的干扰 analytic = hilbert(seg) return analytic

逻辑说明：解析信号通过希尔伯特变换构造，把实信号的频谱限制在正半轴，这样后续模糊函数计算时自项更集中，交叉项的特征也更稳定。标准差归一化这一步容易被忽略，但对CNN训练的收敛速度影响非常大，因为不同流型下传感器信号的绝对幅值差异可能达到十倍以上。

接下来计算模糊函数。数值上有一种高效近似方式：对解析信号做短时傅里叶变换，再对得到的三维复频谱矩阵沿频率方向做自相关。这样把直接计算二维模糊函数的复杂度从(O(N^4))降到(O(N^2\log N))：

def compute_ambiguity(analytic, nfft=1024): """ 通过STFT的自相关近似计算模糊函数 analytic: 解析信号（复数） nfft: FFT点数，和窗长保持一致 """ # 用较短的窗做STFT，窗长取nfft//4，保证时频平面上有足够的帧数 win = windows.hann(nfft // 4) stft_matrix = [] hop = nfft // 8 for start in range(0, len(analytic) - nfft // 4, hop): frame = analytic[start:start + nfft // 4] * win stft_matrix.append(np.fft.fft(frame, nfft)) S = np.array(stft_matrix) # 形状: (帧数, 频点数) # 对频率轴做自相关 amb = np.zeros((nfft, nfft), dtype=complex) for tau in range(-nfft // 2 + 1, nfft // 2): shifted = np.roll(S, tau, axis=1) corr = S * np.conj(shifted) amb[tau + nfft // 2, :] = corr.sum(axis=0) return amb

逻辑说明：模糊域的横轴对应时延，纵轴对应频偏，核心操作是计算每个频点上的信号与其时延版本的相关性。这个实现在精度上比直接公式略低，但在工程上足够用，计算速度比完整模糊函数快一到两个数量级，适合批处理。

3.2 核函数的迭代优化：从初始椭圆到自适应形状

有了模糊函数矩阵，接下来用迭代优化求解径向高斯核。不需要引入复杂的优化库，只要实现一个朴素的梯度上升即可：

def solve_optimal_kernel(amb, alpha=3.0, max_iter=50, lr=1e-2): """ 求解自适应最优核的径向高斯宽度 amb: 模糊函数矩阵 alpha: 体积约束参数，控制核的覆盖范围 """ nfft = amb.shape[0] # 将直角坐标转换成极坐标索引，预先计算好半径和角度网格 y, x = np.mgrid[-nfft//2:nfft//2, -nfft//2:nfft//2] r = np.sqrt(x**2 + y**2) theta = np.arctan2(y, x) r = r.flatten() theta = theta.flatten() amb_flat = amb.flatten() # 沿角度方向将模糊函数分成36个扇区，每个扇区独立优化高斯宽度 sigma = np.full(36, 15.0, dtype=float) # 有效半径只取到nfft//4，远处的模糊函数值基本是噪声 mask = r < nfft // 4 for _ in range(max_iter): grad = np.zeros_like(sigma) for i in range(36): sector = (np.mod(theta - i * 10 * np.pi / 180, 2 * np.pi) < np.pi / 18) idx = mask & sector if idx.sum() == 0: continue # 径向高斯核表达式：exp(-r^2/(2*sigma^2)) kernel_val = np.exp(-r[idx]**2 / (2 * sigma[i]**2)) # 目标函数对sigma的梯度 k = (r[idx]**2 / sigma[i]**3) * kernel_val weighted_energy = np.abs(amb_flat[idx])**2 grad[i] = np.sum(k * weighted_energy) sigma += lr * grad # 投影约束：sigma的上下界限制核的伸缩范围 sigma = np.clip(sigma, 5.0, 60.0) return sigma

参数说明：lr是学习率，取太大核形参数会振荡，太小收敛慢；50次迭代通常足以让sigma稳定。sigma的初始值15对应一个中等宽度的径向高斯，这个值对最终结果影响不大，因为迭代会把它推向数据对应的方向。alpha在这里不是直接以显式约束形式体现的，而是通过sigma的上下界来控制核函数的体积——这实际是带约束问题的一个简化处理，对工业现场的批处理场景已经足够稳健。

3.3 从核函数到时频分布，生成CNN的输入张量

核函数求出之后，把它和模糊函数逐点相乘，然后做二维傅里叶逆变换，得到自适应最优核时频分布：

def aok_tfr(analytic, sigma, nfft=1024): """ 用求得的核参数生成时频分布图像 """ amb = compute_ambiguity(analytic, nfft) y, x = np.mgrid[-nfft//2:nfft//2, -nfft//2:nfft//2] r = np.sqrt(x**2 + y**2) theta = np.arctan2(y, x) # 根据每个采样点的角度索引到对应的sigma值 kernel = np.zeros_like(r, dtype=float) for i in range(36): sector = (np.mod(theta - i * 10 * np.pi / 180, 2 * np.pi) < np.pi / 18) kernel[sector] = np.exp(-r[sector]**2 / (2 * sigma[i]**2)) filtered = amb * kernel tfr = np.fft.fftshift(np.fft.ifft2(np.fft.fftshift(filtered))) return np.abs(tfr).T

这段代码输出的二维矩阵直接可以归一化后存成灰度图或者伪彩色图。但需要注意的是，时频分布的动态范围非常大——强能量区域和噪声区域的幅度差可能超过三个数量级。如果直接把原始幅值作为CNN的输入，梯度会被高能量区域完全主导。常见做法是取对数缩放到0到255的范围：

# 对数压缩 + 百分位截断，避免极端值影响图像对比度 def to_image(tfr, clip_percentile=99): log_tfr = np.log1p(tfr) vmax = np.percentile(log_tfr, clip_percentile) img = np.clip(log_tfr / vmax, 0, 1) return (img * 255).astype(np.uint8)

clip_percentile=99的含义是允许1%的最亮像素被截断，这样做是为了让CNN看到的图像对比度不至于被一两个极亮点完全压暗。这个参数值得在使用时微调，实验中发现对弹状流这类有脉冲冲击的信号，截断到98或99效果最好，而泡状流用95就足够了。

4. 卷积神经网络结构设计：怎么让网络学会读时频图

4.1 输入尺寸与网络的匹配策略

时频图生成之后，输入尺寸的选择是第一个要确定的事。直接把完整的1024×1024矩阵送入网络，计算量太大，也没必要。224×224这个尺寸是Imagenet预训练网络的标准输入，但不建议在时频图上直接用预训练权重，因为时频图的低频纹理特征与自然图像的高频语义特征差异太大。

完整时频图的频率轴比时间轴信息密度更高，因为流型特征主要体现在频率分布形态上。裁剪时优先保留频率轴的完整范围，时间轴可以适当截段。一个具体经验是：把输出图像固定为256×256，其中频率轴覆盖原始的全频段，时间轴覆盖约四分之一窗长的信号内容。这样既能保留流型的频率特征差异，又让不同信号的图像在时间维度上保持可比性。

4.2 轻量级CNN结构的选择与参数设定

气液两相流时频图的类别数通常只有四到六种，数据总量可能只有几千到几万张，在这个规模下去堆ResNet或者EfficientNet这类深层网络容易过拟合，而且推理延迟在工业在线场景里也扛不住。更合理的做法是采用一个四层卷积的浅层网络结构，配合全局平均池化替代全连接层来降低参数量：

import torch import torch.nn as nn class FlowPatternCNN(nn.Module): """ 用于流型识别的轻量卷积神经网络 输入: (B, 1, 256, 256) 的时频图灰度图 输出: (B, C) 的类别对数概率 """ def __init__(self, num_classes=4): super().__init__() self.features = nn.Sequential( # 第一层：提取基本时频形态 nn.Conv2d(1, 16, kernel_size=5, stride=2, padding=2), nn.BatchNorm2d(16), nn.ReLU(inplace=True), nn.MaxPool2d(2), # 第二层：捕捉局部纹理变化 nn.Conv2d(16, 32, kernel_size=3, stride=1, padding=1), nn.BatchNorm2d(32), nn.ReLU(inplace=True), nn.MaxPool2d(2), # 第三层：扩大感受野，覆盖段塞周期特征 nn.Conv2d(32, 64, kernel_size=3, stride=1, padding=1), nn.BatchNorm2d(64), nn.ReLU(inplace=True), nn.MaxPool2d(2), # 第四层：全局结构编码 nn.Conv2d(64, 128, kernel_size=3, stride=1, padding=1), nn.BatchNorm2d(128), nn.ReLU(inplace=True), ) self.global_pool = nn.AdaptiveAvgPool2d((1, 1)) self.classifier = nn.Linear(128, num_classes) def forward(self, x): x = self.features(x) x = self.global_pool(x) x = torch.flatten(x, 1) return self.classifier(x)

结构选型思路：第一个卷积层用5×5的大核且步长为2，因为时频图上最关键的能量聚集区域通常面积较大，小卷积核的初始感受野不足以覆盖。中间层全部用3×3配padding保持尺寸，让特征图尺寸逐层减半。BatchNorm加在每个卷积之后，对时频图这种分布随工况变化明显的输入特别重要，它能把每层的输入分布拉回标准范围。全局平均池化替代Flatten加全连接的做法能大幅减少过拟合风险，因为擦除了空间位置信息，只保留通道级的语义。

4.2.1 通道数放大比例与FLOPs的权衡

从16到128的通道数放大比例是经验性的：每层翻倍能给高层特征足够的表达容量，同时参数量控制在约15万左右，单张256×256图像在CPU上的推理时间能控制在一百毫秒量级。如果现场有GPU，可以把首个卷积层的输出通道调到32，识别精度能提升一到两个百分点，但部署在嵌入式设备上时保持16通道起步更稳妥。

4.3 数据增强操作与时频图的物理约束

时频图不同于自然图像，增强操作不能随意做。水平翻转相当于时间反转，这在流型识别上是物理可实现的——因为时间反转后的流型序列仍然是合法的流型信号，对CNN来说是一种有效的等效样本扩充。垂直翻转则完全不可行，它会把低频能量翻到高频区，生成的样本在物理上不存在于任何真实流型中。

更实用的增强手段是频率轴的局部扰动和整体缩放。频率缩放模拟的是流速变化——流速提高时整个时频图的频率成分按比例上移。可以用PyTorch的torch.nn.functional.interpolate插值实现：

def augment_tfr(img, scale_range=(0.85, 1.15)): """ 时频图的频域缩放增强 img: (H, W) 的张量，H轴为频率轴 """ scale = np.random.uniform(*scale_range) new_h = int(img.shape[0] * scale) resized = torch.nn.functional.interpolate( img.unsqueeze(0).unsqueeze(0), size=(new_h, img.shape[1]), mode='bilinear', align_corners=False ) # 缩放后裁剪或填充回原始尺寸 if new_h > img.shape[0]: start = (new_h - img.shape[0]) // 2 return resized[:, :, start:start + img.shape[0], :] else: pad = (img.shape[0] - new_h) // 2 return torch.nn.functional.pad(resized, (0, 0, pad, img.shape[0] - new_h - pad))

频率缩放的系数范围取0.85到1.15，对应流速波动大约正负百分之十五的工况范围。插值方式必须用bilinear，nearest会在缩放倍数不为整数时产生严重的方块状伪影，这些伪影会被CNN误认为时频结构特征。

5. 参数调优与训练策略：让模型收敛且不踩过拟合的坑

5.1 训练参数速查表与初始值推荐

整理一份训练参数的推荐初始值表格，实际使用中大部分参数不需要大动，真正影响精度的是学习率调度方式：

训练时的学习率调度使用torch.optim.lr_scheduler.CosineAnnealingLR，把训练轮数作为周期参数。线性衰减在时频图分类任务上效果不稳定，因为训练到中后期损失面变得平坦，线性学习率可能来回震荡。

5.2 损失函数与类别不均衡的处理策略

气液两相流实验数据天然存在类别不均衡——泡状流容易出现在宽工况范围，收集到的样本数量远多于段塞频率高的弹状流。如果直接用交叉熵损失，模型会偏向多数类而忽略少数类。简单加权交叉熵可以缓解，但对噪声标注的鲁棒性较差。

一个工程上更有效的技巧是同时使用标签平滑和类别频率加权的组合。权重计算方式采用中位数频率平衡：取所有类别样本数目的中位数除以每个类别的样本数，并截断在0.5到10的范围内。这样不会让少数类的权重膨胀到主导整个训练的程度。

损失函数的PyTorch实现可以直接用内置的CrossEntropyLoss(weight=class_weight, label_smoothing=0.05)。需要注意的是，标签平滑会轻微降低训练集上的准确率，但对测试集的泛化能力几乎总是正向的提升，尤其是当测试数据来自新的工况条件时。

5.3 训练曲线判读与过早停止策略

训练时经常遇到的陷阱是：验证集准确率在训练早期就已经接近顶点，随后进入平台期甚至缓慢下降，但训练集准确率还在持续上升。这是过拟合时频图的明确信号。解决手段不是调低模型容量，而是加强频率缩放增强并提前触发早停。

ReduceLROnPlateau在这种任务上比固定轮数训练更有效。把验证损失作为监控指标，连续5个epoch没有下降时把学习率降低到原来的十分之一，连续10个epoch不下降就终止训练。这个策略在所有流型识别项目中都比固定训练轮数平均高出两到三个百分点的测试精度，原因在于时频图特征的分布差异使得模型泛化能力在训练过程中存在一个显著的峰值区段，靠轮数猜测不如靠验证集信号直接捕捉。

训练中最好顶住只看总体准确率的诱惑，需要把混淆矩阵打印出来看特定流型对的混淆方式。气液两相流里弹状流和段塞流之间的误分经常是一个方向性的问题——弹状流被误判为段塞流的比例远高于反向。如果出现这种情况，需要检查时频图的频率轴是否被压缩了或者数据增强的缩放范围是否符合物理工况，单纯调整网络结构通常解决不了这个层面的问题。

6. 验证协议的进阶技巧：把流型识别结果拆回时频域段做审计

模型训练完成后，常规做法是把测试集整体准确率报一个数字。这个数字在流型识别任务上参考价值有限，因为不同流型的持续时间差别极大，弹状流一段可能只占窗口的四分之一时长，样本层面的准确率会被时间占比稀释。更可靠的验证方式是把测试集按流型类别分开统计F1分数，并且贴着时间窗口的图谱逐条回溯预测错误的样本。

具体的审计技巧是做类别激活映射。四层卷积网络最后一层输出的128个特征图，每张都对应一个特定的时频模式响应。把这些特征图按分类权重加权求和，叠加到输入时频图上，可以直接看出模型在做决策时到底在看哪个区域。弹状流样本激活区域通常集中在低频高能量的条带附近，如果激活区域漂移到高频区，说明模型学到的是噪声伪影而非物理特征。

# 简化的类别激活映射计算流程 # 1. 取最后一个卷积层的输出特征图 A，形状为 (128, H', W') # 2. 取分类层权重 w，形状为 (128, C) # 3. 对类别 c，计算加权和 saliency = sum_i(w[i, c] * A[i]) # 4. 插值到输入尺寸并归一化显示

对激活图检查之后，如果发现了明显的问题，常见措施是把出现频率异常高的激发区域所对应的频率带进行滤波预处理。比如传感器噪声固定在50Hz附近，而激活图恰好在该频带形成热点，可以在时频图预处理阶段直接做一个陷波滤波。这个手段比单纯增加数据增强更直接，因为根本不需要让网络去学一个"这频带是噪音"的映射，直接从输入层面排除掉。

针对部署现场替换和型号变更的情况，还有一个值得加上的验证维度——用不同直径管道采集的少量数据做跨域测试。时频图上的绝对频率值会随管径和流速基准变化而整体漂移，如果模型是在8mm管径的数据上训练的，直接部署到15mm管道上性能跌幅通常在10个百分点以上。简单迁移方案是在新数据上对模型做低学习率微调，只用新工况数据训练最后两层卷积和分类头；若目标是在完全无标注的目标域上工作，则需要引入对抗域适应策略让特征提取器输出在源域和目标域之间不可区分。无论如何，迁移场景下的流型识别性能验证必须放在时频图分布对齐的前提下进行，这不是一个只被模型结构强弱主导的问题。

最后的落地验证技巧是设置运行时的时频图质量门控——计算每张时频图的能量集中度指标，即前百分之十像素占整体能量的比例，该值偏离训练集均值超过两个标准差时，自动标记该样本为低置信度并触发重采样或人工复核。这个指标不需要任何额外标注，实现成本极低，但能让模型从"给出一个结果"提升到"知道自己什么时候结果不可靠"。对这种安全敏感性高的两相流测量场合，一个能主动承认不确定性的模型，远比一个经常自信出错的模型有用。

本文还有配套的精品资源，点击获取