利用 ML 构建水质评价模型|颜色特征

🎲 本文章所有代码基于 Python 完成，所需依赖包如下：

问题背景

有经验的从事渔业生产的从业者可通过观察水色变化调控水质，以维持养殖水体生态系统中浮游植物、微生物类、浮游动物等合理的动态平衡。
由于这些多是通过经验和肉眼观察进行判断的，存在主观性引起的观察性偏倚，使观察结果的可比性、可重复性降低，不易推广应用。

当前，数字图像处理技术为计算机监控技术在水产养殖业的应用提供更大的空间。在水质在线监测方面，数字图像处理技术是基于计算机视觉的，以专家经验为基础，对池塘水色进行优劣分级,实现对池塘水色的准确快速判别。

我们给出了某地区的多个罗非鱼池塘水样的数据（见文末），包含水产专家按水色判断水质分类的数据以及用数码相机按照标准进行水色采集的数据，每个水质图片命名规则为 类别_编号.jpg ，如 1_1.jpg 说明当前图片属于第 1 类的样本。水质类别由下表给出：

水色	浅绿色（清水或浊水）	灰蓝色	黄褐色	茶褐色（姜黄、茶褐、红褐、褐中带绿等）	绿色（黄绿、油绿、蓝绿、墨绿、绿中带褐等）
水质类别	1	2	3	4	5

请根据这些数据，利用图像处理技术，通过水色图像实现水质的自动评价。

问题分析

根据题目背景，我们需要利用机器学习相关算法实现对给定水样照片正确地对其水质情况进行划分（分类）。期间因为运用到了专家经验和大量历史图片及其正确的分类标签，所以是一种有监督的分类问题。

由于输入是一张 $2D$ 图像，因此我们应当对图像进行特征提取，再根据所提取的特征选取合适的分类模型进行训练。

一般来说，对图像特征提取和识别问题，大名鼎鼎的 CNN（卷积神经网络）算法可以很好的完成任务。但是对于本问题来说，较为直观也是最有用的特征无疑是水样的颜色特征。而对于颜色特征，有着颜色直方图、颜色矩这些特征，它们可以充分利用起来作为主要特征，因此无需搭建卷积神经网络模型，利用常见的分类模型或许也能完成任务。

综上所述，我们需要完成以下任务：

1.        对水样图片进行切割，提取水样图片中的特征（这里选择各阶颜色矩）；
2.        基于提取的特征数据，选择合适的分类模型以构建水质评价模型；
3.        对构建的模型进行评价，评价模型对于水色的识别效率；
4.        比较不同分类模型对该问题的性能表现，并分析优缺点。

数据预处理

图像切割

对于图像切割，在图像处理和机器视觉领域已经有足够成熟的算法可供使用，比如常见的单阈值分割、 $k\text{ - }means$ 分割、 $Otus$ 算法等。
我们只需将图像中，有水质颜色的部分提取出来，然后再剪切出颜色分布较为均匀的 $100\times 100$ 像素大小的部分，对其进行颜色特征的计算即可。

不过对于本问题来说，大可不必如此复杂，不难发现所有给定的图像数据中，我们感兴趣的颜色部位基本都在整个图像的正中央，这也符合对水质检验拍照时的一般习惯，（我们拍照时通常把主要目标位于镜头中央）。

因此前期图像切割部分我们只需截取每一张图像中心 $100\times 100$ 像素的部分即可进行下一步操作。

此处我们引入 opencv-python 库编程实现对图像进行切割

颜色矩

我们将彩色图像的 RGB 三个通道分离，分别计算一阶、二阶和三阶颜色矩，将这些特征作为计算的主要特征。

一阶颜色矩
一阶颜色矩采用的是一阶原点矩，反映了图像整体的明暗程度。从算式上看就是像素点灰度值的平均值（mean）：

$E=\frac 1 N\sum_{i=1}^Np_{i}$

式中， $p_i$ 表示某个颜色通道中第 $i$ 个像素的灰度值。

二阶颜色矩
二阶颜色矩采用的是二阶中心矩的平方根，反映了图像颜色的分布范围。从算式上看就是像素点灰度值的有偏估计标准差（Standard Deviation）：

$\sigma=\sqrt{\frac 1 N\sum_{i=1}^N(p_i-E)^2}$

三阶颜色矩
三阶颜色矩采用的是三阶中心矩的立方根，反映了图像颜色分布的对称性。从算式上看就是像素点灰度值的有偏估计偏度（skewness）：

$s=\sqrt[3]{\frac 1 N\sum_{i=1}^N(p_i-E)^3}$

以上各阶颜色矩均可以直接利用 numpy 库进行计算，并采用 pandas库保存为 .csv文件以便后续处理

import os
import re
import cv2
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt

# 根据指定图像计算其各阶颜色矩
def calculate(img):
    # 将图像数据转为一维数组以 numpy.array 存放
    data = np.array(img).flatten()/255
    tmp = np.mean((data - data.mean()) ** 3)
    return [data.mean(), data.std(), np.sign(tmp)*np.abs(tmp)**(1/3)]

# 根据指定路径和命名规则读取图像
def LoadandProcessImage(path):

    # 存储 5 类水质图片的颜色矩数据
    ImageDatas = { 'id':[], # 图像编号
                   'category':[], # 所属类别
                   'Rmean':[], 'Rstd':[], 'Rskewness':[], # R通道的三颜色矩
                   'Gmean':[], 'Gstd':[], 'Gskewness':[], # G通道
                   'Bmean':[], 'Bstd':[], 'Bskewness':[], # B通道
                   }

    namelist = os.listdir(path = path)
    namelist.sort() # 按照文件编号排序
    for name in namelist:
        ans = re.findall('(\d)_\d+.jpg',name)
        if ans: 
            # 读取图像并截取中心 101*101 的矩形区域作为兴趣点
            img = cv2.imread(path + name) 
            M, N = img.shape[0:2]
            img = img[M//2-50:M//2+50,N//2-50:N//2+50,:]

            # 观察图像
            # cv2.imshow('sample',img)
            # cv2.waitKey(0)

            # 分离三通道并计算各阶矩
            b, g, r = cv2.split(img)
            BB = calculate(b)
            GG = calculate(g)
            RR = calculate(r)

            # 存储数据
            ImageDatas['id'].append(name)
            ImageDatas['category'].append(ans[0])
            for key in ['B', 'G', 'R']:
                exec(f"ImageDatas['{key}mean'].append({key+key}[0])")
                exec(f"ImageDatas['{key}std'].append({key+key}[1])")
                exec(f"ImageDatas['{key}skewness'].append({key+key}[2])")
            
    return ImageDatas

if __name__ == "__main__":
    path = "images/"
    ImageDatas = LoadandProcessImage(path)
    pd.DataFrame(ImageDatas).to_csv('image_datas.csv', index = False) # 导出 csv文件

模型构建

此处我们采用以下四种分类模型对该问题进行模型的构建：

决策树模型（Decision Tree, DT）
多层感知机模型（Multi-Layer Perceptron, MLP）
AdaBoost 模型
支持向量机模型（Support Vector Machine, SVM）

Decision Tree

首先读取数据，并将其按照 80% 训练集和 20% 测试集的比例进行划分：

datas = pd.read_csv('image_datas.csv')
# print(datas.head())
features = datas.iloc[:,2:11].values
target = datas['category'].values
X_train, X_test, y_train, y_test = train_test_split(features, target, 
                                                    test_size=0.2, # 二八分割数据
                                                    random_state=1) #程序每次运行结果是否不变

然后通过 sklearn.tree.DecisionTreeClassifier 构建模型，并且利用 sklearn.model_selection.GridSearchCV 进行参数寻优：

def tree_train():
    print("-"*20+" Descision Tree "+"-"*20)

    # 网格搜索选取最优参数
    param = [{'criterion':['gini'],
              'max_depth':[30,50,60,100],
              'min_samples_leaf':[2,3,5,10],
              'min_impurity_decrease':[0.1,0.2,0.5]},
             {'criterion':['gini','entropy']},
             {'max_depth': [30,60,100], 
              'min_impurity_decrease':[0.1,0.2,0.5]}]
    
    model = GridSearchCV(DecisionTreeClassifier(), param_grid = param, cv = 3)
    model.fit(X_train*30, y_train)
    print('最优参数:', model.best_params_, '最优分数:', model.best_score_)
    # 构建模型并预测
    descision_tree = model.best_estimator_
    pred = descision_tree.predict(X_test*30)
    # 混淆矩阵
    cm = confusion_matrix(y_test, pred)
    print('混淆矩阵为\n', cm)
    # 准确率
    acc = accuracy_score(y_test, pred)
    print('准确率为\n', acc)
    # 保存模型
    pickle.dump(descision_tree, open('tree.model','wb'))

MLP

def mlp_train():
    print("-"*20+" Multi-Layer Perceptron "+"-"*20)

    param = {"hidden_layer_sizes": [(81,24)],#, (50, 15)],
             "solver": ['lbfgs'],#,'sgd', 'adam'],
             "max_iter": [4000],
             "verbose": [False]
            }
    model = GridSearchCV(MLPClassifier(), param_grid = param, cv = 5)
    model.fit(X_train*30, y_train)
    print('最优参数:', model.best_params_, '最优分数:', model.best_score_)
    
    mlp = model.best_estimator_
    pred = mlp.predict(X_test*30)
    pickle.dump(mlp, open('mlp.model','wb'))

AdaBoost

def boost_train():
    print("-"*20+" AdaBoost Classifier "+"-"*20)

    param = [{'n_estimators':range(40,60,5),'learning_rate':np.arange(0.5,1.0,0.1)}]
    model = GridSearchCV(AdaBoostClassifier(), param_grid = param, cv = 5)
    model.fit(X_train*30, y_train)
    print('最优参数:', model.best_params_, '最优分数:', model.best_score_)

    boost = model.best_estimator_
    pred = boost.predict(X_test*30)
    pickle.dump(boost, open('boost.model','wb'))

SVM

def svm_train():
    print("-"*20+" SVM "+"-"*20)
    model = svm.SVC(kernel = 'linear',C = 1).fit(X_train*30, y_train) # linear,c=1 最优
    pred = model.predict(X_test*30)
    pickle.dump(model, open('svm.model','wb'))