系列文章目錄

“光晰本質(zhì)，譜見(jiàn)不同”，光譜作為物質(zhì)的指紋，被廣泛應(yīng)用于成分分析中。伴隨微型光譜儀/光譜成像儀的發(fā)展與普及，基于光譜的分析技術(shù)將不只停留于工業(yè)和實(shí)驗(yàn)室，即將走入生活，實(shí)現(xiàn)萬(wàn)物感知，見(jiàn)微知著。本系列文章致力于光譜分析技術(shù)的科普和應(yīng)用

前言

典型的光譜分析模型(以近紅外光譜作為示意，可見(jiàn)光、中遠(yuǎn)紅外、熒光、拉曼、高光譜等分析流程亦相似)建立流程如下所示，在建立過(guò)程中，需要使用算法對(duì)訓(xùn)練樣本進(jìn)行選擇，然后使用預(yù)處理算法對(duì)光譜進(jìn)行預(yù)處理，或?qū)庾V的特征進(jìn)行提取，再構(gòu)建校正模型實(shí)現(xiàn)定量分析，最后針對(duì)不同測(cè)量?jī)x器或環(huán)境，進(jìn)行模型轉(zhuǎn)移或傳遞。因此訓(xùn)練樣本的選擇、光譜的預(yù)處理、波長(zhǎng)篩選、校正模型、模型傳遞以及上述算法的參數(shù)都影響著模型的應(yīng)用效果。

建模流程.jpg

針對(duì)光譜分析流程所涉及的常見(jiàn)的訓(xùn)練樣本的劃分、光譜的預(yù)處理、波長(zhǎng)篩選、校正模型算法建立了完整的算法庫(kù)，名為OpenSA(OpenSpectrumAnalysis)。整套算法庫(kù)的架構(gòu)如下所示。

OpenSA.jpg

樣本劃分模塊提供隨機(jī)劃分、SPXY劃分、KS劃分三種數(shù)據(jù)集劃分方法，光譜預(yù)處理模塊提供常見(jiàn)光譜預(yù)處理，波長(zhǎng)篩選模塊提供Spa、Cars、Lars、Uve、Pca等特征降維方法，分析模塊由光譜相似度計(jì)算、聚類(lèi)、分類(lèi)(定性分析)、回歸(定量分析)構(gòu)建，光譜相似度子模塊計(jì)算提供SAM、SID、MSSIM、MPSNR等相似計(jì)算方法，聚類(lèi)子模塊提供KMeans、FCM等聚類(lèi)方法，分類(lèi)子模塊提供ANN、SVM、PLS_DA、RF等經(jīng)典化學(xué)計(jì)量學(xué)方法，亦提供CNN、AE、Transformer等前沿深度學(xué)習(xí)方法，回歸子模塊提供ANN、SVR、PLS等經(jīng)典化學(xué)計(jì)量學(xué)定量分析方法，亦提供CNN、AE、Transformer等前沿深度學(xué)習(xí)定量分析方法。模型評(píng)估模塊提供常見(jiàn)的評(píng)價(jià)指標(biāo)，用于模型評(píng)估。自動(dòng)參數(shù)優(yōu)化模塊用于自動(dòng)進(jìn)行最佳的模型設(shè)置參數(shù)尋找，提供網(wǎng)格搜索、遺傳算法、貝葉斯概率三種最優(yōu)參數(shù)尋找方法?？梢暬K提供全程的分析可視化，可為科研繪圖，模型選擇提供視覺(jué)信息。可通過(guò)幾行代碼快速實(shí)現(xiàn)完整的光譜分析及應(yīng)用（注: 自動(dòng)參數(shù)優(yōu)化模塊和可視化模塊暫不開(kāi)源，等畢業(yè)后再說(shuō))
本篇針對(duì)OpenSA的光譜預(yù)處理模塊進(jìn)行代碼開(kāi)源和使用示意。

一、光譜數(shù)據(jù)讀入

提供兩個(gè)開(kāi)源數(shù)據(jù)作為實(shí)列，一個(gè)為公開(kāi)定量分析數(shù)據(jù)集，一個(gè)為公開(kāi)定性分析數(shù)據(jù)集，本章僅以公開(kāi)定量分析數(shù)據(jù)集作為演示。

1.1 光譜數(shù)據(jù)讀入

# 分別使用一個(gè)回歸、一個(gè)分類(lèi)的公開(kāi)數(shù)據(jù)集做為example
def LoadNirtest(type):

    if type == "Rgs":
        CDataPath1 = './/Data//Rgs//Cdata1.csv'
        VDataPath1 = './/Data//Rgs//Vdata1.csv'
        TDataPath1 = './/Data//Rgs//Tdata1.csv'

        Cdata1 = np.loadtxt(open(CDataPath1, 'rb'), dtype=np.float64, delimiter=',', skiprows=0)
        Vdata1 = np.loadtxt(open(VDataPath1, 'rb'), dtype=np.float64, delimiter=',', skiprows=0)
        Tdata1 = np.loadtxt(open(TDataPath1, 'rb'), dtype=np.float64, delimiter=',', skiprows=0)

        Nirdata1 = np.concatenate((Cdata1, Vdata1))
        Nirdata = np.concatenate((Nirdata1, Tdata1))
        data = Nirdata[:, :-4]
        label = Nirdata[:, -1]

    elif type == "Cls":
        path = './/Data//Cls//table.csv'
        Nirdata = np.loadtxt(open(path, 'rb'), dtype=np.float64, delimiter=',', skiprows=0)
        data = Nirdata[:, :-1]
        label = Nirdata[:, -1]

    return data, label

1.2 光譜可視化

    #載入原始數(shù)據(jù)并可視化
    data, label = LoadNirtest('Rgs')
    plotspc(data, "raw specturm")

采用的開(kāi)源光譜如圖所示:

image.png

二、光譜預(yù)處理

2.1 光譜預(yù)處理模塊

將常見(jiàn)的光譜進(jìn)行了封裝，使用者僅需要改變名字，即可選擇對(duì)應(yīng)的光譜分析，下面是光譜預(yù)處理模塊的核心代碼

"""
    -*- coding: utf-8 -*-
    @Time   :2022/04/12 17:10
    @Author : Pengyou FU
    @blogs  : https://blog.csdn.net/Echo_Code?spm=1000.2115.3001.5343
    @github :
    @WeChat : Fu_siry
    @License：

"""
import numpy as np
from scipy import signal
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import MinMaxScaler, StandardScaler
from copy import deepcopy
import pandas as pd
import pywt


# 最大最小值歸一化
def MMS(data):
    """
       :param data: raw spectrum data, shape (n_samples, n_features)
       :return: data after MinMaxScaler :(n_samples, n_features)
       """
    return MinMaxScaler().fit_transform(data)


# 標(biāo)準(zhǔn)化
def SS(data):
    """
        :param data: raw spectrum data, shape (n_samples, n_features)
       :return: data after StandScaler :(n_samples, n_features)
       """
    return StandardScaler().fit_transform(data)


# 均值中心化
def CT(data):
    """
       :param data: raw spectrum data, shape (n_samples, n_features)
       :return: data after MeanScaler :(n_samples, n_features)
       """
    for i in range(data.shape[0]):
        MEAN = np.mean(data[i])
        data[i] = data[i] - MEAN
    return data


# 標(biāo)準(zhǔn)正態(tài)變換
def SNV(data):
    """
        :param data: raw spectrum data, shape (n_samples, n_features)
       :return: data after SNV :(n_samples, n_features)
    """
    m = data.shape[0]
    n = data.shape[1]
    print(m, n)  #
    # 求標(biāo)準(zhǔn)差
    data_std = np.std(data, axis=1)  # 每條光譜的標(biāo)準(zhǔn)差
    # 求平均值
    data_average = np.mean(data, axis=1)  # 每條光譜的平均值
    # SNV計(jì)算
    data_snv = [[((data[i][j] - data_average[i]) / data_std[i]) for j in range(n)] for i in range(m)]
    return  data_snv



# 移動(dòng)平均平滑
def MA(data, WSZ=11):
    """
       :param data: raw spectrum data, shape (n_samples, n_features)
       :param WSZ: int
       :return: data after MA :(n_samples, n_features)
    """

    for i in range(data.shape[0]):
        out0 = np.convolve(data[i], np.ones(WSZ, dtype=int), 'valid') / WSZ # WSZ是窗口寬度，是奇數(shù)
        r = np.arange(1, WSZ - 1, 2)
        start = np.cumsum(data[i, :WSZ - 1])[::2] / r
        stop = (np.cumsum(data[i, :-WSZ:-1])[::2] / r)[::-1]
        data[i] = np.concatenate((start, out0, stop))
    return data


# Savitzky-Golay平滑濾波
def SG(data, w=11, p=2):
    """
       :param data: raw spectrum data, shape (n_samples, n_features)
       :param w: int
       :param p: int
       :return: data after SG :(n_samples, n_features)
    """
    return signal.savgol_filter(data, w, p)


# 一階導(dǎo)數(shù)
def D1(data):
    """
       :param data: raw spectrum data, shape (n_samples, n_features)
       :return: data after First derivative :(n_samples, n_features)
    """
    n, p = data.shape
    Di = np.ones((n, p - 1))
    for i in range(n):
        Di[i] = np.diff(data[i])
    return Di


# 二階導(dǎo)數(shù)
def D2(data):
    """
       :param data: raw spectrum data, shape (n_samples, n_features)
       :return: data after second derivative :(n_samples, n_features)
    """
    data = deepcopy(data)
    if isinstance(data, pd.DataFrame):
        data = data.values
    temp2 = (pd.DataFrame(data)).diff(axis=1)
    temp3 = np.delete(temp2.values, 0, axis=1)
    temp4 = (pd.DataFrame(temp3)).diff(axis=1)
    spec_D2 = np.delete(temp4.values, 0, axis=1)
    return spec_D2


# 趨勢(shì)校正(DT)
def DT(data):
    """
       :param data: raw spectrum data, shape (n_samples, n_features)
       :return: data after DT :(n_samples, n_features)
    """
    lenth = data.shape[1]
    x = np.asarray(range(lenth), dtype=np.float32)
    out = np.array(data)
    l = LinearRegression()
    for i in range(out.shape[0]):
        l.fit(x.reshape(-1, 1), out[i].reshape(-1, 1))
        k = l.coef_
        b = l.intercept_
        for j in range(out.shape[1]):
            out[i][j] = out[i][j] - (j * k + b)

    return out


# 多元散射校正
def MSC(data):
    """
       :param data: raw spectrum data, shape (n_samples, n_features)
       :return: data after MSC :(n_samples, n_features)
    """
    n, p = data.shape
    msc = np.ones((n, p))

    for j in range(n):
        mean = np.mean(data, axis=0)

    # 線(xiàn)性擬合
    for i in range(n):
        y = data[i, :]
        l = LinearRegression()
        l.fit(mean.reshape(-1, 1), y.reshape(-1, 1))
        k = l.coef_
        b = l.intercept_
        msc[i, :] = (y - b) / k
    return msc

# 小波變換
def wave(data):
    """
       :param data: raw spectrum data, shape (n_samples, n_features)
       :return: data after wave :(n_samples, n_features)
    """
    data = deepcopy(data)
    if isinstance(data, pd.DataFrame):
        data = data.values
    def wave_(data):
        w = pywt.Wavelet('db8')  # 選用Daubechies8小波
        maxlev = pywt.dwt_max_level(len(data), w.dec_len)
        coeffs = pywt.wavedec(data, 'db8', level=maxlev)
        threshold = 0.04
        for i in range(1, len(coeffs)):
            coeffs[i] = pywt.threshold(coeffs[i], threshold * max(coeffs[i]))
        datarec = pywt.waverec(coeffs, 'db8')
        return datarec

    tmp = None
    for i in range(data.shape[0]):
        if (i == 0):
            tmp = wave_(data[i])
        else:
            tmp = np.vstack((tmp, wave_(data[i])))

    return tmp

def Preprocessing(method, data):

    if method == "None":
        data = data
    elif method == 'MMS':
        data = MMS(data)
    elif method == 'SS':
        data = SS(data)
    elif method == 'CT':
        data = CT(data)
    elif method == 'SNV':
        data = SNV(data)
    elif method == 'MA':
        data = MA(data)
    elif method == 'SG':
        data = SG(data)
    elif method == 'MSC':
        data = MSC(data)
    elif method == 'D1':
        data = D1(data)
    elif method == 'D2':
        data = D2(data)
    elif method == 'DT':
        data = DT(data)
    elif method == 'WVAE':
        data = wave(data)
    else:
        print("no this method of preprocessing!")

    return data

2 .2 光譜預(yù)處理的使用

在example.py文件中，提供了光譜預(yù)處理模塊的使用方法，具體如下，僅需要兩行代碼即可實(shí)現(xiàn)所有常見(jiàn)的光譜預(yù)處理。
示意1：利用OpenSA實(shí)現(xiàn)MSC多元散射校正

 #載入原始數(shù)據(jù)并可視化
    data, label = LoadNirtest('Rgs')
    plotspc(data, "raw specturm")
    #光譜預(yù)處理并可視化
    method = "MSC"
    Preprocessingdata = Preprocessing(method, data)
    plotspc(Preprocessingdata, method)

預(yù)處理后的光譜數(shù)據(jù)如圖所示:

image.png

示意2：利用OpenSA實(shí)現(xiàn)SNV預(yù)處理

    #載入原始數(shù)據(jù)并可視化
    data, label = LoadNirtest('Rgs')
    plotspc(data, "raw specturm")
    #光譜預(yù)處理并可視化
    method = "SNV"
    Preprocessingdata = Preprocessing(method, data)
    plotspc(Preprocessingdata, method)

預(yù)處理后的光譜數(shù)據(jù)如圖所示:

image.png

總結(jié)

利用OpenSA可以非常簡(jiǎn)單的實(shí)現(xiàn)對(duì)光譜的預(yù)處理，完整代碼可從獲得GitHub倉(cāng)庫(kù) 如果對(duì)您有用，請(qǐng)點(diǎn)贊！
代碼現(xiàn)僅供學(xué)術(shù)使用，若對(duì)您的學(xué)術(shù)研究有幫助，請(qǐng)引用本人的論文，同時(shí)，未經(jīng)許可不得用于商業(yè)化應(yīng)用，歡迎大家繼續(xù)補(bǔ)充OpenSA中所涉及到的算法，如有問(wèn)題，微信：Fu_siry