破裂预测流程及问题（初稿）

破裂预测总体上分3步：1）数据预处理 2）模型训练 3）预测分析。
经过上一学期的工作，以上3步的程序已经写完。对原来选择的诊断信号进行了训练并预测，但效果不理想。现在正换一批信号用于训练。

破裂预测的具体工作：

一、数据预处理

这部分的工作是整个破裂预测工作中最麻烦的部分。 就是选择相应的诊断信号，从MDSplus中下载所选诊断的炮的数据，进行切片、归一化，整理成可训练的数据集。 具体点我将数据预处理分成了以下10个步骤，此10步也可以作为炮数据预处理的一般步骤。

1.从MDSplus中获取数据到客户端

将炮以HDF5格式储存到本地. 出现的问题:各种存储诊断信号数据都非常大，并且网络带宽有限，导致下载炮数据非常慢。 实际操作： 选择2019年秋季这轮2614个炮
选择诊断信号：
r'\ip', r'\Bt', r'\axuv_ca_01', r'\vs_c3_aa001', r'\vs_ha_aa001', r'\sxr_cc_049', r'\Polaris_Den_Rt01', r'\Polaris_Den_Rt02', r'\Polaris_Den_Rt03', r'\Polaris_Den_Rt04', r'\Polaris_Den_Rt05', r'\Polaris_Den_Rt06', r'\Polaris_Den_Rt07', r'\Polaris_Den_Rt08', r'\Polaris_Den_Rt09', r'\Polaris_Den_Rt10', r'\Polaris_Den_Rt11', r'\Polaris_Den_Rt12', r'\Polaris_Den_Rt13', r'\Polaris_Den_Rt14', r'\Polaris_Den_Rt15', r'\Polaris_Den_Rt16', r'\Polaris_Den_Rt17'

需要下载大概12小时，储存空间大约200G。

2.通Ip信号选出可用炮号

去掉等离子体电流没有完全爬升、没有信号的炮号。初步选择出Ip可用的正常炮和破裂炮的炮号，存储到ValidShot列表内。

出现的问题:ValidShot列表内还是存在一些不标准的正常炮和破裂炮，如下图。需要后续的手动剔除。
实际操作： 步骤1下载的2614个炮，去掉了1137个无用炮，ValidShot列表内有可用炮1477个。

3.选取可用的诊断信号

对托卡马克进行探测时，有时诊断装置、采集装置会故障，导致一些炮中的某些诊断信号缺失。在ValidShot列表中的一些炮也就缺少所选诊断（见步骤1）中的部分诊断。所有，步骤3的目的是从ValidShot列表中剔除这些诊断不完整的炮号。步骤3还需要分成两小步：

1）列出每个诊断信号的可用炮个数

通过程序，找出所选诊断中每个诊断在ValidShot列表1477炮中有效的个数。下图为输出结果。
由图可知，在ValidShot中的1477炮中没有Polaris_Den_Rt12~17这6个诊断信号，因此在所选诊断中删除。 得到了清理后的所选诊断列表：
r'\ip', r'\Bt', r'\axuv_ca_01', r'\vs_c3_aa001', r'\vs_ha_aa001', r'\sxr_cc_049', r'\Polaris_Den_Rt01', r'\Polaris_Den_Rt02', r'\Polaris_Den_Rt03', r'\Polaris_Den_Rt04', r'\Polaris_Den_Rt05', r'\Polaris_Den_Rt06', r'\Polaris_Den_Rt07', r'\Polaris_Den_Rt08', r'\Polaris_Den_Rt09', r'\Polaris_Den_Rt10', r'\Polaris_Den_Rt11'

此小步的目的：对所选择的诊断信号列表进行清洗，剔除无用的诊断信号。

2）筛选出所选信号均有效的炮号

在上一小步图中，除Polaris_Den_Rt12~17这6个诊断有效个数为0外，其它没被剔除的诊断有效个数相互之间有不相同的（Bt : 1473 , Polaris_Den_Rt01 : 1454等）。说明ValidShot中的1477炮有部分炮的部分诊断无效，需要把这些炮剔除。
实际操作： 从ValidShot中的1477炮中筛选出清理后的所选诊断均有效的炮号，然后储存到CompletionShot列表中。CompletionShot列表中有1439炮。

出现的问题：步骤3的2小步均需要遍历ValidShot中1477炮的所有诊断，通过每个诊断的数据判断此诊断是否有效；并且每一炮作为一个HDF5文件储存的，所以步骤3会大量加载HDF5文件，使其执行过程较慢。 为了解决此问题，我和吴琪琪正在建一个可用诊断名的数据库。（上周组会说的那个）

4.区分正常炮和破裂炮

通过程序将CompletionShot列表中的正常炮和破裂炮区分开，对应的炮号分别储存在Normal列表和Break列表内。 实际操作：
Normal列表内967炮，Break列表内472炮。

5.通vp信号排除主动破裂炮

主动破裂炮是一般人为注气让其破裂，则注气破裂后气压迅速大幅度升高。而自然破裂炮破裂后气压不会明显变化。所以主动破裂炮的最大vp信号幅值是自然破裂炮最大vp信号幅值的3~4数量级倍数。则可选择一个阈值作为判断自然破裂炮的标准。当最大vp信号幅值小于阈值时，此炮为自然破裂炮；当最大vp信号幅值大于阈值时，此炮为主动破裂炮。
实际操作： 选择的阈值是0.015。 得到自然破裂293个，并将炮号覆盖储存到Break列表内。

6.手动剔除不良炮

由于目前还不知道那些炮是不规范炮（步骤2所说），需要将Normal列表内967炮和Break列表内293炮的Ip信号画出来。眼睛观察，将波形不规范炮的炮号剔除。如步骤2中的炮或以下的炮：

在后续降采样时，需要寻找正常炮和破裂炮的下降时刻，并且寻找下降时刻的函数程序已经完成（精确度1ms）。但是由于炮的波形噪音干扰非常厉害，找到下降点就非常困难。因此我们现在的下降时刻的函数程序的正确率不是100%（我和吴琪琪学长反复优化，正确率达到了99.6%），有些炮的下降时间找得不准确。而下降时间不准确，就会导致后续将采样和时间切片过程发生错误。因此画Ip图时也调用找下降时间函数，把下降时间点标在信号图上，观察图形，把下降时间错误的炮也剔除掉。
实际操作： 经过步骤6后，Normal列表内正常炮834个和Break列表内破裂炮238个。
出现的问题：步骤6需要人工操作，过程繁琐。

7.划分训练集和测试集

从Normal列表和Break列表内各随机取50炮，作为测试炮集（共100炮），分别储存在TestNormal列表和TestBreak列表中。
剩下的784个正常炮和188个破裂炮作为训练炮集，用于时间切片，分别储存在TrainNormal列表和TrainBreak列表中。

8.信号降采样

托卡马克的诊断信号的采样频率有10kHz、50kHz、250kHz、500kHz、2MHz、10MHz这几种。而预测的要求信号精确到1ms,并且为便于后续切片，所以需要将所选信号都降采样到1kHz。 实际操作：
我采用平均值降采样（比如10kHz降到1KHz,连续的10个样点取平均值作为样点）。并且保留 0.2s 到 炮下降时刻 这个时间段的诊断数据，降采样。
整体文件从200G降到了40M。

9.选取归一化基准值

训练前数据预处理，归一化，将数据大小缩小到[-1,1]，有助于训练。一般归一化方法有： 1）标准化

2）区间放缩法

但是，具体分析所有的诊断会发现所有的信号不能采用上述的任一方法作为通用方法来完成归一化。

如图是Ip信号图，可以观察到波形非常嘈杂，有冲击干扰。区间放缩法中用程序找到的最大值不准确，所以区间放缩法不可行。当然可以考虑先将信号滤波，信号处理好后，再用区间放缩法，这样是可行的。但是这种先滤波、后用区间放缩的方法仅对Ip、电子密度等非周期信号有效。但对于下图的axuv_ca_01信号，它的波形在0附近高频震荡，如果对其进行滤波，信号的波形会是趋于0的直线，完全失真。

实际操作： 我选用的归一化方法是：归一化值 = 信号幅值 / 基准值
基准值：对一具体信号归一化时，我先随机画出一定数量的信号波形图，然后观察信号绝对值的最大值，在最大值的基础上留下一定裕度估选一个数作为基准值。
某种信号基准值找到后，可将其保存在数据库里，以后一直使用。 出现的问题：当需要换数据库里没有的新信号训练时，就要将新信号的波形图画出来，观察并估选基准值。对换信号不太友善。

10.时间窗切片并形成训练矩阵

预测时是以时间窗进行训练的。我选择的时间窗长度32ms(看论文选取的，程序中用全局变量赋值设置，易更改) 。经降采样后，信号采样率为1kHz，即每个诊断信号32个采样点为一个时间窗。将每炮的每时刻对应的各所选诊断的时间窗拼接在一起作为一个实例，特征有 32*len(所选诊断列表) 个。一个炮可以被切成多个时间窗切片实例。 打标签：将步骤8中降采样数据保存的时间范围 0.2s 到 炮下降时刻。
将破裂炮下降时刻前的30ms~5ms这25个时间窗的实例打标签为-1。 将破裂炮下降时刻前的30ms之前和正常炮的所有时间窗的实例打标签为1。 实际操作：
将步骤7中已经分成训练炮集进行切片，作为训练集，即TrainNormal列表和TrainBreak列表中的炮。将标签为1的实例储存为TrainPositive.npy，标签为-1的实例储存为TrainNegative.npy。
将步骤7中已经分成测试炮集进行切片，作为测试集，即TestNormal列表和TestBreak列表中的炮。将标签为1的实例储存为TestPositive.npy，标签为-1的实例储存为TestNegative.npy。

####个人想法： 在换炮号训练、换信号训练时，均需要对数据进行预处理，将上述步骤跑一遍。每都要等待下载、画和浏览波形图、手动剔除炮。个人觉得数据预处理最麻烦。

#二、模型训练
采用 SVM 模型
rbf核中可调超参数：区分间隔gamma 、 惩罚系数 C linear核中可调超参数：惩罚系数 C sigmoid核中可调超参数：区分间隔gamma 、 惩罚系数 C 用时间窗切片训练集来训练模型，寻找最佳超参数，模型训练完成后，储存模型。
用切片测试集来测试模型的预测切片精确性，在未过拟合下精确度最高的模型对应的超参数为最佳超参数。 选取最佳超参数预测模型用于对炮的预测。

####个人想法: 破裂预测中，我感觉最花时间、最麻烦是数据预处理部分。虽然寻找最佳超参数模型训练花费时间长（现在12天），但这个过程现在自动化实现，等结果即可。训练的这几天。可以干些其他的事情。
如果换其它模型也好办，在程序上仅改几行代码即可(当然了解模型原理也仅需12天)。 #三、预测分析
将预测的那一炮时间窗切片，在通过最佳切片预测模型预测这炮每个切片，返回预测标签1或-1，将切片预测标签按对应时间拍列。 理想情况下，预测破裂炮时，开始的一段时间对应的时间窗切片的预测标签是1；在接近炮下降时刻的一段时间对应的时间窗切片预测标签都会是-1。
理想情况下破裂炮预测标签时间图：
图中，下降沿对应的时刻即使提前预测时刻。

理想情况下，预测正常炮时，此炮所有时刻对应的时间窗切片预测标签都会是1。
理想情况下正常炮预测标签时间图：

出现的问题：对炮的预测分析都是基于对时间窗切片预测精确度。实际中对切片预测的精确度不高时，预测标签图的波形会变得非常嘈杂。如下图： 当要预测的炮非常多时，不可能人工筛选炮，需要用程序筛选。当对切片预测的精确度不高，预测标签图的波形嘈杂时，正常炮和破裂炮的特点将难以用程序语言描述，不能准确的判断是正常炮还是破裂炮。所有必须尽量提高对时间窗切片的预测模型的精度。

####个人想法： 对炮预测分析这个环节的程序已经写好，每次换炮号或换信号时都不用大改，输出结果自动保存，不会花费太多时间。

---

本文章使用limfx的vsocde插件快速发布