SINBAD论文笔记 | Set Features for Fine-grained Anomaly Detection

洋芋焖土豆

已于 2024-08-23 15:47:33 修改

阅读量953

点赞数 24

分类专栏：缺陷检测文章标签：论文阅读

于 2024-08-22 13:56:33 首次发布

本文链接：https://blog.csdn.net/weixin_42355902/article/details/141137917

版权

缺陷检测专栏收录该内容

6 篇文章 0 订阅

订阅专栏

文章目录

Set Features for Fine-grained Anomaly Detection

Set Features for Fine-grained Anomaly Detection

Summary

使用一个特征提取器提取特征，每张图的feature map视为一个集合，每个最小1x1的patch视为一个元素
将每个集合沿共享的随机方向投影，投影矩阵参数正太分布，然后沿着这个方向为每个集合计算一个1D直方图，
直方图描述符采用每个维度上的值的区间分布
异常分数采用test图的直方图描述符与训练集中最近的直方图描述符的带白化的马氏距离

Method(s)

3. Set Features for Anomaly Detection

3.1 一个集合不仅仅是其各部分的总和
检测逻辑图像异常或集体时间序列异常需要了解每个样本的不同部分如何相互作用，举例来说：螺丝袋类，此类中的每个正常样品均包含两个螺钉（不同长度）、两个螺母和两个垫圈，当缺少一颗螺钉或用额外的螺母代替其中一个垫圈时，可能会出现异常情况。检测此类异常需要对样本中的所有元素进行联合描述，因为每个局部元素都可能来自于正常样本

典型的聚合元素描述符特征的方法是通过平均池化（average pool）。然而，这并不总是适合于集合异常检测。在监督学习中，平均池化通常被构建到像ResNet（He等人，2016年）或DeepSets（Zaheer等人，2017年）这样的架构中，以聚合局部特征。因此，用监督损失学习到的深层特征已经被训练为池化有效。然而，对于低级特征描述符，比如我们在这里使用的这些，情况可能并非如此。正如图2所示，一组特征的平均远非集合的完整描述。这在异常检测中尤其正确，其中密度估计方法需要比监督学习所需的更具区分性的特征。
在这里插入是描述
在这里，我们希望通过模拟其元素的分布来描述一个集合，同时忽略它们之间的顺序。一种简单的做法是通过离散化的体积表示，类似于点云的3D体素。不幸的是，这种方法无法扩展到高维度，而且需要更紧凑的表示。在这项工作中，我们选择使用一组1D直方图来表示集合。每个直方图代表了集合元素在沿特定方向投影时的密度。我们在图2中提供了这个想法的说明。
沿原始轴投影集合可能不足以区分。原始轴上的直方图对应于1D边缘分布，可能会将远离的元素映射到相同的直方图箱中（见图2的说明）。另一方面，我们可以在图的底部看到，当集合元素首先沿另一个方向投影时，两个集合的直方图是不同的。这表明了一种集合描述的方法：
首先将每个集合沿共享的随机方向投影，然后沿着这个方向为每个集合计算1D直方图。通过沿多个随机方向重复此过程，我们可以获得一个更强大的描述符

PS：radon变换：Radon变换的本质是将原来的函数做了一个空间转换，即，将原来的XY平面内的点映射到AB平面上，那么原来在XY平面上的一条直线的所有的点在AB平面上都位于同一点。记录AB平面上的点的积累厚度，便可知XY平面上的线的存在性。这便是大家所公认的Radon变换的实质所在。

3.2 Preliminaries
对于一个训练集S，包含 $N_s$ 个正常的样本，在test时，对于一个测试样本 $\tilde{x}$ 。目标是学习一个模型，该模型对每个样本 $\tilde{x}$ 输出异常分数，异常分数高于预定阈值的样本被标记为异常。
方法的特点：对每个样本x视为一组元素组成，即
$x=[e_{1},e_{2}...e_{N_{E}}].$

这类元素的例子包括图像的局部区域（patches）和时间序列的时间窗口（temporal windows）。我们假设存在一个强大的特征提取器 F，它将每个原始元素 e 映射到一个元素特征描述符 $f_e$ 。我们将在第4节中描述两个重要应用的特征提取的具体实现：图像和时间序列。

3.3 Set Features by Histogram of Projections
直方图描述符：我们使用每个维度上的值的直方图来描述该集合( $N_s$ , $N_e$ , $N_D$ )，我们将每个样本的每个元素中第 j 个特征的值的集合记为 $s[j] = \{f_{1}[j],f_{2}[j]..f_{N_{S}}[j]\}$ (从前面说明和代码逻辑，这应该是 $N_e$ )。我们计算所有样本中集合 s[j] 的最大值和最小值，并将它们之间的区域划分为 K 个区域。我们计算 $N_D$ 中每个维度的直方图 $H_j$ 并将它们连接成单个描述符集 h。因此，每个集合的描述符的维度为 $N_D·K$

投影：如前所述，并非所有投影方向对于描述集合的分布都具有同等的信息量。在一般情况下，不知道哪些方向对于捕获正常集和异常集之间的差异最具信息性。由于我们无法提前告知最佳投影方向，因此我们建议随机投影特征。这确保了灾难性不良投影方向的可能性较低，如图 2 中的示例所示。

我们通过从高斯分布 N(0, 1) 的每个维度采样值来生成随机投影矩阵 P ∈ $R^{(N_D,N_P)}$ 。我们投影 x 的每个元素的特征 f，产生投影特征 $f^{'} = P f$
我们在投影特征上运行上面描述的直方图描述符过程。因此，最终的集合描述符 $h_x$ 变成了 $N_P$ 直方图的串联，得到了 $N_P$ ·K维，最后这个集合描述为( $N_s$ ,K* $N_D$ )

3.4 Anomaly scoring
我们现在希望使用直方图特征 $h_x$ 来检测异常样本。我们使用高斯密度估计器估计正态数据特征{ $h_x$ }∈S的概率密度函数(PDF)。我们为f中的分布假设提供了理论证明。我们计算了所有集合的描述符的平均值µ和协方差Σ。估计PDF由下表给出:
$p(h)=\mathcal{N}(h|\mu,\Sigma)$
由于异常评分与样本的负似然相关，我们将异常评分a(h)定义为负对数似然。忽略常数项，这就是马氏距离:
$a(h)=(h-\mu)^T\Sigma^{-1}(h-\mu)$
在实践中，我们发现到kNN(带白化变换的kNN)的马氏距离比到µ的马氏距离要好一些，因此我们用它来计算我们的结果。

白化：白化处理后相比于原始数据有两个特点：数据的不同维度去相关；数据每个维度的方差为1。白化的目的就是降低输入的冗余性；更正式的说，我们希望通过白化过程使得学习算法的输入具有如下性质：(i)特征之间相关性较低；(ii)所有特征具有相同的方差。

在这里插入图片描述

文中表中的总异常分数ours是针对ResNet layer 3 、layer 4和原始像素集获得的异常分数的平均值分别由以下因子 (1, 1, 0.1) 加权。No pixels：使用两个 ResNet 块的结果，但没有原始像素级别平均异常分数，代码中在处理测试集在各block的anomaly map归一化是，用对应验证集在block上的平均值作为缩放因子

Application to Image and Time Series Anomaly Detection

4.1 图像作为集合
图像可以被视为由不同粒度级别的元素组成的集合。这从像素到小局部区域，再到低级元素如线条或角点，一直到高级元素如对象。对于异常检测，我们通常事先不知道正确的粒度级别，以便在正常和异常样本之间进行区分。这取决于异常，而在训练期间是未知的。相反，我们为不同的粒度级别执行异常检测并将分数组合起来。这些粒度级别对应于不同大小的局部区域。

在实践中，我们使用预训练的ResNet（He等人，2016）中间块的表示。由于ResNet网络同时嵌入了每个图像的许多局部区域，我们通过网络编码器传递图像，并从不同残差块末端的中间激活中提取我们的表示（见图3）。我们将激活图中的每个空间位置定义为一个元素。请注意，由于不同的块具有不同的分辨率，它们产生不同数量的元素。我们以每个残差块末端的元素作为我们的集合。
在这里插入图片描述

4.2 时间序列作为集合
时间序列数据可以被视为时间窗口的集合。与图像类似，通常事先不知道哪个时间尺度与检测异常相关；即包含语义现象的窗口持续时间。受到Rocket Dempster等人（2020）的启发，我们将时间序列的基本元素定义为时间窗口金字塔的集合。每个金字塔包含L个窗口。金字塔中的所有窗口都以相同的时间步为中心，每个窗口包含τ个样本（见图3）。第一级窗口包括τ个元素，步长为1；第二级窗口包括τ个元素，步长为2，依此类推。这样的窗口金字塔为时间序列中的每个时间步计算，整个序列被表示为其金字塔元素的集合。更多的实现细节描述在第C.2节中。

代码分析

# 使用wideresnet获取featur map
    def get_super_patch_embeddings(im_from_load):
        with torch.no_grad():
            # x.unfold(窗口的方向, 窗口大小, 窗口的步长)，crop操作
            super_patches = im_from_load.unfold(2, crop_size, stride_size).unfold(3, crop_size, stride_size) # # 1, 3, 1024, 1024 ->torch.Size([1, 3, 1, 1024, 1022])-> torch.Size([1, 3, 1, 1, 1022, 1022])
            super_patches = torch.reshape(super_patches,(3,super_patches.shape[2]*super_patches.shape[3],super_patches.shape[4],super_patches.shape[5])) # ->torch.Size([3, 1, 1022, 1022])
            super_patches = torch.transpose(super_patches, 0, 1)# ->torch.Size([1, 3, 1022, 1022])
            super_patches = F.interpolate(super_patches, (224, 224)) # ->torch.Size([1, 3, 224, 224])

            if args.pyramid_level > 56:
                trans_clr = transforms.Compose([transforms.Resize(int(args.pyramid_level))])
                super_patches = trans_clr(super_patches)
                X = torch.reshape(super_patches, (super_patches.shape[0], super_patches.shape[1],-1))
                token = 0
            else:
                x, x77, x1414, x2828, x5656 = net(super_patches.cuda()) # (1,2048), (1,2048, 7,7 ), (1,1024, 14,14 ), (1,512,28,28),(1,256,56,56)

                if args.pyramid_level == 7:
                    X = torch.reshape(x77, (x77.shape[0], x77.shape[1], x77.shape[2]*x77.shape[3])) # (1,2048, 49 )
                if args.pyramid_level == 14:
                    X = torch.reshape(x1414, (x1414.shape[0], x1414.shape[1], x1414.shape[2]*x1414.shape[3]))
                if args.pyramid_level == 28:
                    X = torch.reshape(x2828, (x2828.shape[0], x2828.shape[1], x2828.shape[2]*x2828.shape[3]))
                if args.pyramid_level == 56:
                    X = torch.reshape(x5656, (x5656.shape[0], x5656.shape[1], x5656.shape[2]*x5656.shape[3]))

                token = x

            return X, token
#######################################
# 从 feature map中提取直方图描述符
   def extract_radon_feaures(radon_extractor, local_mini_patches_emb, is_projection = True):
        ind_c = 0
        if is_projection:
            radon_feaures = np.zeros((len(local_mini_patches_emb), args.n_projections*args.n_quantiles))
        else:
            radon_feaures = np.zeros((len(local_mini_patches_emb), n_channels*args.n_quantiles)) # (69,5000)
        # np.ceil(a) np.floor(a) : 计算各元素的ceiling 值， floor值（ceiling向上取整，floor向下取整）
        for batch_idx in range(int(np.ceil(len(local_mini_patches_emb)/batch_size))):
            batch_local_mini_patches_emb = local_mini_patches_emb[ind_c:int(ind_c+batch_size)] #torch.Size([32, 2048, 49])

            if ind_c == 0:
                radon_extractor.fit(batch_local_mini_patches_emb)
            batch_train_radon, _ = radon_extractor.forward(batch_local_mini_patches_emb) # (32,5000)5000:N_P*k ,(32,5000,49)
            radon_feaures[ind_c:int(ind_c+batch_size)] = batch_train_radon # 填充

            ind_c = int(ind_c + batch_size)

        return radon_feaures
	#radon_extractor的定义
class CumulativeSetFeatures(torch.nn.Module):
    def __init__(self, n_channels, n_projections=100, n_quantiles=20, is_projection=True):
        self.n_channels = n_channels 
        self.n_projections = n_projections # 投影的数量，默认为100。这通常用于数据降维或特征提取
        self.n_quantiles = n_quantiles # 分位数的数量，分位数用于将数据划分为多个区间。
        self.projections = torch.randn(self.n_projections,  self.n_channels, 1) #从标准正态分布中随机初始化的1D卷积核，用于当 is_projection 为 True 时对输入数据进行投影。
        self.is_projection = is_projection # 指示是否使用投影步骤来转换输入数据。

    def fit(self, X): #用于计算训练数据的最小值和最大值，这些值用于后续的分位数计算。
        if self.is_projection: # conv1d(input, weight, bias=None, stride=1, padding=0, dilation=1, groups=1），weight是一维卷积核，其形状为(Out_Channel,In_Channel/Group,Kernel_Size),当Group=1时，每一层输出由所有输入分别与卷积核卷积的累加得到
            a = F.conv1d(X, self.projections).permute((0, 2, 1)) # 用1D卷积对输入数据 X 进行投影，torch.Size([32, 2048, 49])->torch.Size([32, 1000, 49])->torch.Size([32, 49, 1000])
            a = a.reshape((-1, self.n_projections)) # 展平，torch.Size([1568, 1000])
        else:
            a = X.permute((0, 2, 1))
            a = a.reshape((a.shape[0]*a.shape[1], -1))
        self.min_vals = torch.quantile(a, 0.01, dim=0) #计算展平后数据的1%和99%分位数 （1000）
        self.max_vals = torch.quantile(a, 0.99, dim=0)


    def forward(self, X):
        if self.is_projection:
            a = F.conv1d(X, self.projections)
        else:
            a = X
        cdf = torch.zeros((a.shape[0], a.shape[1], self.n_quantiles)) #torch.Size([32, 1000, 5])
        set = torch.zeros((a.shape[0], a.shape[1], X.shape[-1], self.n_quantiles,)) # torch.Size([32, 1000, 49, 5])
        for q in range(self.n_quantiles):
            threshold = self.min_vals + (self.max_vals - self.min_vals) * (q + 1) / (self.n_quantiles + 1) # 循环遍历每个分位数，计算每个分位数对应的阈值
            set[:, :, :, q] = (a < threshold.unsqueeze(0).unsqueeze(2)).float()
            cdf[:, :, q] = set[:, :, :, q].mean(2)

        set = torch.transpose(set, 2, 1)
        set = set.reshape((X.shape[0], X.shape[-1], -1))
        set = torch.transpose(set, 2, 1).numpy()
        cdf = cdf.reshape((X.shape[0], -1)).numpy()

        return cdf, set # 描述符(32,5000) , 集合(32,5000,49)


#######################################
# 使用knn获取最近的点
        kNN_index = kNN_shrunk(train_radon, args.n_score, args.is_cpu,  is_whitening, is_vector = True, shrinkage_factor = args.shrinkage_factor)
        pred_score = kNN_index.score(test_radon)
        pred_score_valid = kNN_index.score(valid_radon)
        pred_score = pred_score[:,0]
        pred_score_valid = pred_score_valid[:,0]

        anom_map[i] = pred_score
        anom_map_valid[i] = pred_score_valid