基于AI的淡水养殖水质溯源、优化系统方案
前言
8月24日，日本政府无视国际社会的强烈质疑和反对，单方面强行启动福岛核污染水排海。海关总署决定在当日起全面暂停了进口原产地为日本的水产品。
海鲜市场遭遇了前所未有的困境。消费者对于核污水问题的恐慌已经超越了个人口腹之欲，海鲜销售量急剧下降。同时，二级市场对核污染水排海事件的敏度同样很高。受该消息影响，核污染水排海相关多个概念股异动，出现大幅度拉升甚至涨停。
淡水水产可能会成为其替代品，从而推动淡水养殖行业的景气度提升。借助本次intel 黑客松，我们通过问题分析，提出可落地、有价值、有创新性的解决方案，同时借助于英特尔 oneAPI AI 工具分析包 实现一个AI淡水养殖水质溯源、优化系统应用程序的调试与改进。
一、核心痛点及关键需求
排海事件打击了消费者信心，推动了供给变化，导致短期市场现象级紊乱，但其背后反应的是市场对于绿色、安全、可监控、可溯源食品的需求在不断增加。对农业产业链上的企业来说，无疑是一个难得的发展机遇。
在提出方案的同时，要考虑以下痛点：
1.淡水水产品牌效应差。2.淡水水产无法完全替代海产品口感。3水产养殖模式需要规范，在节约成本的同时，经验可以复制。
接下来我们对于淡水养殖水质溯源这一个命题进行分析其价值。
1.政策引导及产业升级
新颁布的《中华人民共和国水污染防治法》严格了污染排放要求，加快推进水产养殖业绿色发展转型升级的过程中，如何推进池塘清洁生产、推进养殖尾水达标排放、发展资源节约型池塘养殖模式，推进池塘养殖节本增效富渔增收是核心问题。
基于AI技术实现的智慧水产养殖模式，具备养殖环境监测调控、投喂效果监测、病害监测预警等重要功能，可以大幅提高水产养殖的效率和收益。
2.紧跟时事及供给变化
按照洋流的流经顺序，核污染水排海，至少要在几年后才会对国内海水养殖业构成冲击。短期的恐慌过后，国内的相关产业应该好好利用窗口期，调整产业结构，减少对远洋捕捞业的依赖，多发展内陆的江河湖泊养殖，同时国内大量的淡水湖也有着相关的天然地理优势。
3.品牌打造，重拾消费者信赖
目前淡水水产的品牌影响力仍然较差。
老百姓“舌尖上的安全”有保障了，才会放心消费。
水产品可追踪溯源显得尤为重要， 养殖过程需要严格监管，通过AI分析决策系统，对水产品的生产过程、用药记录和销售记录建立健全追溯档案，是水产养殖业走向标准化生产的必然路径。
4.特色生态新模式
随着“海虾淡养”的成功，淡水养出“海味”也成为了现实，虽然技术门槛高，培育养殖成功后，一年四季都能养殖，市场需求量大，产业前景广阔。
渔业协调发展探索离不开AI的助力，通过AI水质监测优化，可以为当地开发和利用盐碱水域，构建特色的盐碱水湖塘养殖生态新模式。
二、方案设计
淡水养殖水质溯源系统，可以助力水产行业用户打造企业品牌，提高竞争力，提供有力保障，同时通过产品溯源，提供让老百姓吃的安心、吃的放心的健康有机好产品。
该方案以intel AI水质预测算法为核心，然后针对于不同场景的特征分布进行定向分析，实现水质溯源档案、重金属富集预测、污染风险监测、水华预警预测、兽药残留检测等功能，
在水产品生产的各个环节（鱼苗、投喂、水质等）如实进行采集、保存、记录追溯信息，保障水产品可溯源，实现集中化管理、规范化养殖。
1.水质溯源档案
水质溯源档案建立水产从投苗到收获的环境全生命周期监测，除了水质预测的各类特征外，基于水域的水质档案还涵盖了过程、人员、监控等辅助数据。当水域水产出现问题，可以及时通过水质回溯档案查询可能涵盖的问题，同时通过监控，精准定位可能由自然灾害、人为恶意破坏等受灾原因。为追责提供留档证据。
2.数字孪生系统
数字孪生系统是企业提高品牌竞争力的一个重要手段，它将系统数据与环境数据进行结合，在三维世界进行映射。
消费者除了可以直观的了解产品溯源，更可以增加一些知识图谱、菜系介绍、周边美食等定制化推广方案。
3.基于intel AI 水质预测算法
通过针对预测水域质量的数据集，我们从多个角度对其展开分析，包括不限于缺失值分析、异常值分析、不规则分布分析、离散/连续数据分析等。
在分析过程中，我们使用了，
英特尔® Modin 分发版提升pandas处理与探索数据；
英特尔® Extension for Scikit-learn提升基线模型（SVM等）的执行效率；
英特尔® 架构优化的 XGBoost作为我们最终的分类模型进行效果的提升与验证。
三、实践案例
1、方案简述
根据之前获取的数据集进行特征构造，用于预测目标水域质量，围绕特征与目标间的关联表现进行特征筛选，利用英特尔® Modin 分发版提升pandas处理特征集的效率，得到了50维相关特征作为最终的入模特证。最后，利用英特尔® 架构优化的XGBoost模型进行二分类。
完成模型后，使用Intel Extension for PyTorch提供的分布式训练功能，利用多台机器和多GPU进行训练。通过使用ipex.distributed.DistributedDataParallel()来包装模型，并在Tesla T4GPU上并行训练数据。
环境配置：
2、数据分析
预处理
首先，获取一份水质监测的离线数据，针对于后续可能实现的污染监测、重金属富集，加入一些随机特征进行模拟。然后将该数据集按照5：2.5：2.5的比例划分为训练集、验证集、测试集。测试集仅在模型推理阶段时使用，可以防止数据穿越等问题的发生，并能够准确对模型关于预测目标的拟合效果进行评估。
特征类型处理
1、划分类别特征, 数值特征
# 划分类别特征, 数值特征
cat_cols, float_cols = [], ['Target']
for col in data.columns:
if data[col].value_counts().count() < 50:
cat_cols.append(col)
else:
float_cols.append(col)
print('离散特征：', cat_cols)
print('连续特征：', float_cols)
离散特征： [‘Color’, ‘Source’, ‘Month’, ‘Day’, ‘Time of Day’, ‘Target’]
连续特征： [‘Target’, ‘pH’, ‘Iron’, ‘Nitrate’, ‘Chloride’, ‘Lead’, ‘Zinc’…]
2、针对离散特征，进行缺失值处理
# 查看缺失值情况
display(data[cat_cols].isna().sum())
missing=data[cat_cols].isna().sum().sum()
print("\nThere are {:,.0f} missing values in the data.".format(missing))
print('-' * 50)
其中像月份属于有序类别特征，如果单纯地映射数字会破坏掉这种有序性，对于这类数据进行单独映射。
......
# 将月份映射至数字1~12
months = ["January", "February", "March", "April", "May",
"June", "July", "August", "September", "October",
"November", "December"]
data['Month'] = data['Month'].replace(dict(zip(months, range(1, len(months) + 1))))
display(data['Month'].value_counts())
......
针对明显无序的类别特征，我们对其进行独特编码处理。对于有序类别特征，我们不对其进行额外处理。
......
# 将无序类别特征进行独特编码
color_ohe = OneHotEncoder()
X_color = color_ohe.fit_transform(data.Color.values.reshape(-1, 1)).toarray()
pickle.dump(color_ohe, open('../feat_data/color_ohe.pkl', 'wb')) # 将编码方式保存在本地
dfOneHot = pandas.DataFrame(X_color, columns=["Color_" + str(int(i)) for i in range(X_color.shape[1])])
data = pandas.concat([data, dfOneHot], axis=1)
......
由于连续变量存在数据缺失的情况，所以我们通过统计每个连续类别对应的中位数对缺失值进行填充。
display(data[float_cols].isna().sum())
missing=data[float_cols].isna().sum().sum()
print("\nThere are {:,.0f} missing values in the data.".format(missing))
# 使用中位数填充连续变量每列缺失值
fill_dict = {}
for column in list(data[float_cols].columns[data[float_cols].isnull().sum() > 0]):
tmp_list = list(data[column].dropna())
# 平均值
mean_val = sum(tmp_list) / len(tmp_list)
# 中位数
tmp_list.sort()
mid_val = tmp_list[len(tmp_list) // 2]
if len(tmp_list) % 2 == 0:
mid_val = (tmp_list[len(tmp_list) // 2 - 1] + tmp_list[len(tmp_list) // 2]) / 2
fill_val = mid_val
# 填充缺失值
data[column] = data[column].fillna(fill_val)
fill_dict[column] = fill_val
特征分布分析
# 针对每个连续变量的频率画直方图
data[float_cols].hist(bins=50,figsize=(16,12))
对于非常不规则的数据进行非线性变换。
# 针对不规则分布的变量进行非线性变换，一般进行log
log_col = ['Iron', 'Zinc', 'Turbidity', 'Copper', 'Manganese']
show_col = []
for col in log_col:
data[col + '_log'] = np.log(data[col])
show_col.append(col + '_log')
# 特殊地，由于Lead变量0值占据了一定比例，此处特殊处理Lead变量
excep_col = ['Lead']
spec_val = {}
# 先将0元素替换为非0最小值，再进行log变换
for col in excep_col:
spec_val[col] = data.loc[(data[col] != 0.0), col].min()
data.loc[(data[col] == 0.0), col] = data.loc[(data[col] != 0.0), col].min()
data[col + '_log'] = np.log(data[col])
show_col.append(col + '_log')
data[show_col].hist(bins=50,figsize=(16,12))
可以发现，上述不规则分布的变量经过log之后已经变成了模型更容易拟合的类正态分布。
3、特征构造
由于特征之间往往隐藏着很多相关的信息等待挖掘，此处进行特征关联，得到新的特征。由于时间关系，此处我们仅以类别特征进行分组，去依次统计其余特征在组内的中位数、均值、方差、最大值、最小值、计数特征，挖掘有效的交互特征。
......
# 特征交互统计
cat_interaction_dict = {}
del_feat_list = []
if 'del_feat_list' in count_fea_dict:
del_feat_list = count_fea_dict['del_feat_list']
add_feat_list = []
for cat_col1 in cat_cols:
if cat_col1 == label_col:
continue
new_col = cat_col1 + '_count'
if new_col not in data.columns and new_col not in del_feat_list:
add_feat_list.append(new_col)
temp = data.groupby(cat_col1).size()
cat_interaction_dict[new_col] = dict(temp)
temp = temp.reset_index().rename(columns={0: new_col})
data = data.merge(temp, 'left', on=cat_col1)
for cat_col2 in cat_cols:
if cat_col2 == label_col:
continue
if cat_col1 == cat_col2:
continue
new_col = cat_col1 + '_' + cat_col2 + '_count'
if new_col not in data.columns and new_col not in del_feat_list:
add_feat_list.append(new_col)
temp = data.groupby(cat_col1)[cat_col2].nunique()
cat_interaction_dict[new_col] = dict(temp)
temp = temp.reset_index().rename(columns={cat_col2: new_col})
data = data.merge(temp, 'left', on=cat_col1)
......
4、特征选择
过滤法
由于本次目标是对淡水质量是否可用进行二分类预测，所以我们针对包含构造特征的所有特征数据，对于离散特征，利用卡方检验判断特征与目标的相关性，对于连续特征，利用方差分析判断特征与目标的相关性，以此做特征集的第一轮筛选。
# 卡方检验对类别特征进行筛选
def calc_chi2(x, y):
import pandas as pd
import numpy as np
from scipy import stats
chi_value, p_value, def_free, exp_freq = -1, -1, -1, []
tab = pd.crosstab(x, y)
if tab.shape == (2, 2):
# Fisher确切概率法, 2✖️2列联表中推荐使用Fisher检验
oddsr, p_value = stats.fisher_exact(tab, alternative='two-sided')
else:
# Pearson卡方检验, 参数correction默认为True
chi_value, p_value, def_free, exp_freq = stats.chi2_contingency(tab, correction=False)
min_exp = exp_freq.min()
if min_exp >= 1 and min_exp < 5:
# Yates校正卡方检验
chi_value, p_value, def_free, exp_freq = stats.chi2_contingency(tab, correction=True)
return chi_value, p_value, def_free, exp_freq
def select_feat_by_chi2(df, cate_cols, target_col, alpha=0.05, cut=99999):
col_chi2_list = []
for col in cate_cols:
chi_value, p_value, def_free, exp_freq = calc_chi2(df[col], df[target_col])
if p_value < alpha:
col_chi2_list.append([col, chi_value])
col_chi2_list = sorted(col_chi2_list, key=lambda x: x[1], reverse=False)
rel_cols = []
for i in range(min(len(col_chi2_list), cut)):
rel_cols.append(col_chi2_list[i][0])
return rel_cols
data = {
'a': ['Ohio', 'Ohio', 'Ohio', 'Nevada', 'Nevada', 'Nevada'],
'b': [2000, 2001, 2002, 2001, 2000, 2002],
'Target': [0, 0, 1, 1, 1, 0]
}
df = pd.DataFrame(data)
cols = ['a', 'b']
select_feat_by_chi2(df, cols, 'Target')
# 方差分析法对连续特征进行筛选
def calc_fc(x, y):
import pandas as pd
import statsmodels.api as sm
from statsmodels.formula.api import ols
tab = pd.DataFrame({'group': y, 'feat': x})
tab['group'] = tab['group'].astype('str').map(lambda x: 'g_' + x)
mod = ols('feat ~ group', data=tab).fit()
ano_table = sm.stats.anova_lm(mod, typ=2)
F = ano_table.loc['group', 'F']
p_value = ano_table.loc['group', 'PR(>F)']
return F, p_value
def select_feat_by_fc(df, cate_cols, target_col, alpha=0.05, cut=99999):
from scipy import stats
col_F_list = []
for col in cate_cols:
F, p_value = calc_fc(df[col], df[target_col])
F_test = stats.f.ppf((1-alpha),
df[target_col].nunique(),
len(df[col]) - df[target_col].nunique())
# p值小于显著性水平 or MSa/MSe大于临界值，则拒绝原假设，认为各组样本不完全来自同一总体
if p_value < alpha or F > F_test:
col_F_list.append([col, F])
col_F_list = sorted(col_F_list, key=lambda x: x[1], reverse=False)
rel_cols = []
for i in range(min(len(col_F_list), cut)):
rel_cols.append(col_F_list[i][0])
return rel_cols
data = {
'a': [0.13, 0.09, 0.2, 0.11, 0.13, 0.15],
'b': [2000, 2001, 2002, 2001, 2000, 2002],
'Target': [0, 0, 1, 1, 1, 0]
}
df = pd.DataFrame(data)
cols = ['a', 'b']
select_feat_by_fc(df, cols, 'Target')
重要性排序
针对第一轮筛选剩余的150维特征，使用模型去进行拟合，并根据最后模型产出的特征重要性从大到小排序，根据重要性以及设定阈值对特征进行第二轮筛选。
# 模型训练函数
def train(data):
## Prepare Train and Test datasets ##
print("Preparing Train and Test datasets")
X_train, X_test, y_train, y_test = prepare_train_test_data(data=data,
target_col='Target',
test_size=.25)
## Initialize XGBoost model ##
ratio = float(np.sum(y_train == 0)) / np.sum(y_train == 1)
parameters = {
'scale_pos_weight': len(raw_data.loc[raw_data['Target'] == 0]) / len(raw_data.loc[raw_data['Target'] == 1]),
'objective': "binary:logistic",
'learning_rate': 0.1,
'n_estimators': 18,
'max_depth': 10,
'min_child_weight': 5,
'alpha': 4,
'seed': 1024,
}
xgb_model = XGBClassifier(**parameters)
xgb_model.fit(X_train, y_train)
print("Done!\nBest hyperparameters:", grid_search.best_params_)
print("Best cross-validation accuracy: {:.2f}%".format(grid_search.best_score_ * 100))
## Convert XGB model to daal4py ##
xgb = grid_search.best_estimator_
daal_model = d4p.get_gbt_model_from_xgboost(xgb.get_booster())
## Calculate predictions ##
daal_prob = d4p.gbt_classification_prediction(nClasses=2,
resultsToEvaluate="computeClassLabels|computeClassProbabilities",
fptype='float').compute(X_test,
daal_model).probabilities # or .predictions
xgb_pred = pd.Series(np.where(daal_prob[:, 1] > .5, 1, 0), name='Target')
xgb_acc = accuracy_score(y_test, xgb_pred)
xgb_auc = roc_auc_score(y_test, daal_prob[:, 1])
xgb_f1 = f1_score(y_test, xgb_pred)
## Plot model results ##
print("\ndaal4py Test ACC: {:.2f}%, F1 Accuracy: {:.2f}%, AUC: {:.5f}".format(xgb_acc * 100, xgb_f1 * 100, xgb_auc))
# plot_model_res(model_name='XGBoost', y_test=y_test, y_prob=daal_prob[:, 1])
importance_df = pd.DataFrame({'importance':xgb.feature_importances_,'feat_name': data.drop('Target', axis=1).columns})
importance_df = importance_df.sort_values(by='importance', ascending=False)
return importance_df
importance_df = train(data)
display(importance_df.head())
针对现有的150维特征对XGB模型进行拟合，并针对输出的特征重要性得分进行排序，将小于阈值的特征进行过滤。
5、模型训练
针对经过两轮筛选得到的50维入模特征，使用RandomizedSearchCV对XGBoost的重要参数进行搜索，同时使用StratifiedKFold对模型效果进行交叉验证。
......
import modin.pandas as pd
from modin.config import Engine
Engine.put("dask")
......
from sklearn.preprocessing import OneHotEncoder
from sklearnex import patch_sklearn
patch_sklearn()
......
from sklearn.model_selection import train_test_split, StratifiedKFold, GridSearchCV, RandomizedSearchCV
from sklearn.preprocessing import RobustScaler, StandardScaler
from sklearn.metrics import make_scorer, recall_score, precision_score, accuracy_score, roc_auc_score, f1_score
from sklearn.metrics import roc_auc_score, roc_curve, auc, accuracy_score, f1_score
from sklearn.metrics import precision_recall_curve, average_precision_score
from sklearn.svm import SVC
from xgboost import plot_importance
import pickle
......
# 模型训练函数
def train(data):
## Prepare Train and Test datasets ##
print("Preparing Train and Test datasets")
X_train, X_test, y_train, y_test = prepare_train_test_data(data=data,
target_col='Target',
test_size=.25)
## Initialize XGBoost model ##
ratio = float(np.sum(y_train == 0)) / np.sum(y_train == 1)
parameters = {
'scale_pos_weight': len(raw_data.loc[raw_data['Target'] == 0]) / len(raw_data.loc[raw_data['Target'] == 1]),
'objective': "binary:logistic",
'learning_rate': 0.1,
'n_estimators': 18,
'max_depth': 10,
'min_child_weight': 5,
'alpha': 4,
'seed': 1024,
}
xgb_model = XGBClassifier(**parameters)
## Tune hyperparameters ##
strat_kfold = StratifiedKFold(n_splits=3, shuffle=True, random_state=21)
print("\nTuning hyperparameters..")
grid = {'min_child_weight': [1, 5, 10],
'gamma': [0.5, 1, 1.5, 2, 5],
'max_depth': [15, 17, 20],
strat_kfold = StratifiedKFold(n_splits=3, shuffle=True, random_state=1024)
grid_search = RandomizedSearchCV(xgb_model, param_distributions=grid,
cv=strat_kfold, n_iter=10, scoring='accuracy',
verbose=1, n_jobs=-1, random_state=1024)
grid_search.fit(X_train, y_train)
print("Done!\nBest hyperparameters:", grid_search.best_params_)
print("Best cross-validation accuracy: {:.2f}%".format(grid_search.best_score_ * 100))
## Convert XGB model to daal4py ##
xgb = grid_search.best_estimator_
daal_model = d4p.get_gbt_model_from_xgboost(xgb.get_booster())
## Calculate predictions ##
daal_prob = d4p.gbt_classification_prediction(nClasses=2,
resultsToEvaluate="computeClassLabels|computeClassProbabilities",
fptype='float').compute(X_test,
daal_model).probabilities # or .predictions
xgb_pred = pd.Series(np.where(daal_prob[:, 1] > .5, 1, 0), name='Target')
xgb_acc = accuracy_score(y_test, xgb_pred)
xgb_auc = roc_auc_score(y_test, daal_prob[:, 1])
xgb_f1 = f1_score(y_test, xgb_pred)
## Plot model results ##
print("\ndaal4py Test ACC: {:.2f}%, F1 Accuracy: {:.2f}%, AUC: {:.5f}".format(xgb_acc * 100, xgb_f1 * 100, xgb_auc))
# plot_model_res(model_name='XGBoost', y_test=y_test, y_prob=daal_prob[:, 1])
##
importance_df = pd.DataFrame({'importance':xgb.feature_importances_,'feat_name': data.drop('Target', axis=1).columns})
importance_df = importance_df.sort_values(by='importance', ascending=False)
print('model saving...')
pickle.dump(xgb, open("../result/model_best.dat", "wb"))
return importance_df
四、Intel OneAPI AI Kit 优化
Intel ® Modin
本方案使用了Intel ® Modin 分发版提升pandas处理与探索数据；Intel 为Pandas库提供了基于Dask的并行计算功能。通过配置引擎，代码将使用Dask作为默认的并行计算引擎，在使用Pandas或Modin进行数据处理和分析时，将利用Dask的分布式计算能力来加速操作。
Intel ® Extension for Scikit-learn
本方案使用 Intel® Extension for Scikit-learn的目的是提升基线模型（SVM等）的执行效率；
其中，包含使用Intel scikit-learn的StandardScaler进行特征缩放，并使用train_test_split函数进行数据集分割。在完成分割后，函数还将scaler对象保存到文件中，并打印训练集的形状。同时使用了intel优化的NearestNeighbors、RobustScaler等类进行聚类、分类、标准化动作。
Intel daal4py
在oneAPI提供的daal4py模型的加速下，算法的推理速度能得到了更大的提升，同时，在特征选择和模型训练中，当数据集的某些特性不符合特定条件时，需要应用特定的补丁来处理这些特性。PatchingConditionsChain 类提供了一种方便的方式来创建和管理这些补丁条件链。
Intel ® XGBoost
Intel XGBoost的XGBClassifier是一个基分类器，它使用XGBoost算法进行分类任务。该分类器默认使用gbtree作为基分类器，并支持并行计算处理。XGBClassifier和XGBRegressor类提供了nthread参数，用于在训练期间可以使用的线程数。
我们通过调整模型的nthread、n_estimators、learning_rate等参数来优化模型的性能。
Intel ® Extension for PyTorch
首先，使用Intel Scikit-learn库提供的预处理函数来对数据进行归一化，以加速神经网络的训练。可以使用sklearn_ipydnn.Normalize()函数对输入数据进行归一化。
然后，使用Intel Scikit-learn库提供的优化器来替换PyTorch的优化器。例如，sklearn_ipydnn.Adam()函数来替换PyTorch的Adam优化器。这些优化器针对Intel MKL进行了优化，可以加速训练过程。
功耗更低，速度更快的Tesla T4
完成模型后，使用Intel Extension for PyTorch提供的分布式训练功能，利用多台机器和多GPU进行训练。通过使用ipex.distributed.DistributedDataParallel()来包装模型，并在Tesla T4 GPU上并行训练数据。
将数据加载到GPU内存中，以加速数据传输和处理。可以使用PyTorch的torch.cuda.is_available()函数检查是否有可用的GPU，然后使用.to(device)方法将数据和模型移动到GPU上。
总结
最终，我们经过多次训练，模型在训练集的平均准确率约91.31%，验证集上的准确率达到90.03%，针对单条数据的推理时间为0.03s。针对特殊处理出来的测试集上准确率达到90.42%，最高准确率可达92.16%。上图是其中一次训练结果的展示结果，该结果有力地验证了本模型的有效性。
上述方案我们规划并实现了一个简单的淡水养殖水质溯源系统，同时采用了Intel ® Modin 、Intel ® Extension for Scikit-learn、 Intel daal4py、Intel ® XGBoost、Intel ® Extension for PyTorch、Tesla T4加速等多种手段，大大提高了推理速度和系统鲁棒性，验证了我们提出模型的有效性、高效性、可行性。
1296
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
