基于堆栈自编码网络的的近红外光谱分类

系列文章目录

近红外光谱分析技术属于交叉领域，需要化学、计算机科学、生物科学等多领域的合作。为此，在（北京邮电大学杨辉华老师团队）指导下，近期准备开源传统的PLS，SVM，ANN，RF等经典算和SG，MSC，一阶导，二阶导等预处理以及GA等波长选择算法以及CNN、AE等最新深度学习算法，以帮助其他专业的更容易建立具有良好预测能力和鲁棒性的近红外光谱模型。

---

---

前言

NIRS是介于可见光和中红外光之间的电磁波，其波长范围为（1100∼2526 nm。 由于近红外光谱区与有机分子中含氢基团（OH、NH、CH、SH）振动的合频和 各级倍频的吸收区一致，通过扫描样品的近红外光谱，可以得到样品中有机分子含氢 基团的特征信息，常被作为获取样本信息的一种有效的载体。 基于NIRS的检测方法具有方便、高效、准确、成本低、可现场检测、不 破坏样品等优势，被广泛应用于各类检测领域。但 近红外光谱存在谱带宽、重叠较严重、吸收信号弱、信息解析复杂等问题，与常用的 化学分析方法不同，仅能作为一种间接测量方法，无法直接分析出被测样本的含量或 类别，它依赖于化学计量学方法，在样品待测属性值与近红外光谱数据之间建立一个 关联模型(或称校正模型，Calibration Model) ，再通过模型对未知样品的近红外光谱 进行预测来得到各性质成分的预测值。现有近红外建模方法主要为经典建模 （预处理+波长筛选进行特征降维和突出，再通过pls、svm算法进行建模）以及深度学习方法（端到端的建模，对预处理、波长选择等依赖性很低）

---

本篇主要讲述堆栈自编码(SAE)的的近红外光谱定性分析建模方法

一、数据来源

使用药品数据，共310个样本，每条样本404个变量，根据活性成分，分成4类
图片如下：

![药品数据光谱](https://img-blog.csdnimg.cn/7d76d4b1ed4d4f71a150298bdae3da74.png?x-oss-process=image/watermark,type_ZHJvaWRzYW5zZmFsbGJhY2s,shadow_50,text_Q1NETiBARWNob19Db2Rl,size_20,color_FFFFFF,t_70,g_se,x_16

二、代码解读

2.1 读取数据

#载入数据
def TableDataLoad(tp, test_ratio, start, end, seed):

    # global data_x
    data_path = './/Data//table.csv'
    Rawdata = np.loadtxt(open(data_path, 'rb'), dtype=np.float64, delimiter=',', skiprows=0)
    table_random_state = seed

    if tp =='raw':
        data_x = Rawdata[0:, start:end]
    if tp =='SG':
        SGdata_path = './/Data//TableSG.csv'
        data = np.loadtxt(open(SGdata_path, 'rb'), dtype=np.float64, delimiter=',', skiprows=0)
        data_x = data[0:, start:end]
    if tp =='SNV':
        SNVata_path = './/Data//TableSNV.csv'
        data = np.loadtxt(open(SNVata_path, 'rb'), dtype=np.float64, delimiter=',', skiprows=0)
        data_x = data[0:, start:end]
    if tp == 'MSC':
        MSCdata_path = './/Data//TableMSC.csv'
        data = np.loadtxt(open(MSCdata_path, 'rb'), dtype=np.float64, delimiter=',', skiprows=0)
        data_x = data[0:, start:end]
    data_y = Rawdata[0:, -1]
    x_col = np.linspace(7400, 10507, 404)
    plotspc(x_col, data_x[:, :], tp=0)

    x_data = np.array(data_x)
    y_data = np.array(data_y)
    
   

    print('训练集规模：{}'.format(len(X_train[:,0])))
    print('测试集规模：{}'.format(len(X_test[:,0])))

2 .2 划分数据集

 X_train, X_test, y_train, y_test = train_test_split(x_data, y_data, test_size=test_ratio,random_state=table_random_state)

2 .3 返回训练集和测试集并绘制光谱图片

    data_train, data_test = ZspPocess(X_train, X_test,y_train,y_test,need=True) #for cnn :false only used in proseesing comparsion
    return data_train, data_test
    

2.4 建立SAE模型

Stacked Autoencoder（SAE）模型是一个由多层稀疏自编码器组成的深度神经网络模型，如下图所示。其前一层自编码器的输出作为其后一层自编码器的输入，最后一层是个分类器（logistic分类器或者softmax分类器）。

因此，首先建立单层自编码器

lass AutoEncoder(nn.Module):

    def __init__(self, inputDim, hiddenDim):
        super().__init__()
        self.inputDim = inputDim
        self.hiddenDim = hiddenDim
        self.encoder = nn.Linear(inputDim, hiddenDim, bias=True)
        self.decoder = nn.Linear(hiddenDim, inputDim, bias=True)
        self.act = F.relu

    def forward(self, x, rep=False):

        hidden = self.encoder(x)
        hidden = self.act(hidden)
        if rep == False:
            out = self.decoder(hidden)
            #out = self.act(out)
            return out
        else:
            return hidden

然后构建SAE单层

class SAE(nn.Module):

    def __init__(self, encoderList):

        super().__init__()

        self.encoderList = encoderList
        self.en1 = encoderList[0]
        self.en2 = encoderList[1]
        # self.en3 = encoderList[2]

        self.fc = nn.Linear(128, 4, bias=True)

    def forward(self, x):

        out = x
        out = self.en1(out, rep=True)
        out = self.en2(out, rep=True)
        #out = self.en3(out, rep=True)
        out = self.fc(out)
        # out = F.log_softmax(out)

        return out

最后构建完整的SAE网络预训练和分类器

class SAE_net(object):
    def __init__(self, AE_epoch = 200, SAE_epoch = 200,
                 input_dim = 404, hidden1_dim = 512,
                 hidden2_dim = 128, output_dim = 10,
                 batch_size = 128):
        self.AE_epoch = AE_epoch
        self.SAE_epoch = SAE_epoch
        self.input_dim = input_dim
        self.hidden1_dim = hidden1_dim
        self.hidden2_dim = hidden2_dim
        self.output_dim = output_dim
        self.batch_size = batch_size
        self.train_loader = None

        encoder1 = AutoEncoder(self.input_dim, self.hidden1_dim)
        encoder2 = AutoEncoder(self.hidden1_dim, self.hidden2_dim)
        self.encoder_list = [encoder1, encoder2]


    def trainAE(self, x_train, y_train, encoderList, trainLayer, batchSize, epoch, useCuda=False):
        if useCuda:
            for i in range(len(encoderList)):
                encoderList[i].to(device)

        optimizer = optim.Adam(encoderList[trainLayer].parameters())
        ceriation = nn.MSELoss()

        data_train = MyDataset(x_train, y_train)
        self.train_loader = torch.utils.data.DataLoader(data_train, batch_size=batchSize, shuffle=True)

        for i in range(epoch):
            sum_loss = 0
            if trainLayer != 0:  # 单独处理第0层，因为第一个编码器之前没有前驱的编码器了
                for i in range(trainLayer):  # 冻结要训练前面的所有参数
                    for param in encoderList[i].parameters():
                        param.requires_grad = False

            for batch_idx, (x, target) in enumerate(self.train_loader):
                optimizer.zero_grad()
                if useCuda:
                    x, target = x.to(device), target.to(device)
                x, target = Variable(x).type(torch.FloatTensor), Variable(target).type(torch.LongTensor)
                # x = x.view(-1, 404)
                x = x.view(x.size(0), -1)
                # 产生需要训练层的输入数据
                # inputs = Variable(inputs).type(torch.FloatTensor).to(device)  # batch x
                # labels = Variable(labels).type(torch.LongTensor).to(device)  # batch y
                out = x
                if trainLayer != 0:
                    for i in range(trainLayer):
                        out = encoderList[i](out, rep=True)

                # 训练指定的自编码器
                pred = encoderList[trainLayer](out, rep=False).cpu()

                loss = ceriation(pred, out)
                sum_loss += loss.item()
                loss.backward()
                optimizer.step()

    def trainClassifier(self, model, epoch, useCuda=False):
        if useCuda:
            model = model.to(device)

        # 解锁参数
        for param in model.parameters():
            param.requires_grad = True

        optimizer = optim.Adam(model.parameters())
        ceriation = nn.CrossEntropyLoss()

        for i in range(epoch):
            # trainning
            sum_loss = 0
            for batch_idx, (x, target) in enumerate(self.train_loader):
                optimizer.zero_grad()
                if useCuda:
                    x, target = x.to(device), target.to(device)
                x, target = Variable(x).type(torch.FloatTensor), Variable(target).type(torch.LongTensor)
                x = x.view(-1, 404)

                out = model(x)

                loss = ceriation(out, target)
                sum_loss += loss.item()
                loss.backward()
                optimizer.step()
        self.model = model

    def fit(self, x_train = None, y_train = None):
        x_train = x_train[:, np.newaxis, :]
        x_train = torch.from_numpy(x_train)
        x_train = x_train.float()

        # pre-train
        for i in range(2):
            self.trainAE(x_train=x_train, y_train=y_train,
                        encoderList = self.encoder_list, trainLayer=i, batchSize=self.batch_size,
                         epoch = self.AE_epoch)
        model = SAE(encoderList=self.encoder_list)
        self.trainClassifier(model=model, epoch=self.SAE_epoch)

    def predict_proba(self, x_test):
        x_test = torch.from_numpy(x_test)
        x_test = x_test.float()
        x_test = x_test[:, np.newaxis, :]
        x_test = Variable(x_test)
        x_test = x_test.view(-1, 404)

        out = self.model(x_test)
        outdata = out.data
        self.y_proba = outdata
        y_proba = outdata.numpy()
        return y_proba

    def predict(self, x_test):
        _, y_out = torch.max(self.y_proba, 1)
        y_pred = []
        for i in y_out:
            y_pred.append(i)
        return y_pred

在进行函数的分装，建立模型与数据集间的关系，输入数据集，输出预测精确度和，训练时间以及预测时间

def SAE1(X_train, y_train, X_test, y_test):

    clf = SAE_net()

    start_time = time.time()
    clf.fit(X_train, y_train)
    train_time = time.time()
    y_proba = clf.predict_proba(X_test)
    y_pred = clf.predict(X_test)
    # ACC
    acc = accuracy_score(y_test, y_pred)
    test_time = time.time()

    return acc, (train_time - start_time), (test_time - train_time)

总结

建立了一个SAE模型，实现对药品光谱的预测，预测准确率96.3%，低于CNN，目前来看，CNN网络在大多数情况下，还是优于自编码网络的，也有相关文献佐证，完整代码后续整理上传GitHub仓库 如果对您有用，请点赞！
代码仅供学术使用，如有问题