课程报告分享|机器学习与声子晶体分析报告
1 问题描述
由于卷积神经网络可以有效的学习到某些潜在的对称性特征，因此推测其可以用来解决力学中周期性结构的特征值求解问题。其中，弹性超材料就是一种典型的周期性结构，它可以用做周期隔震基础和周期屏障，乃至于地震波的衰减。而评价弹性超材料特性的指标之一就是其带隙特性。但带隙特性需通过频散曲线的计算来表征，而这其实就是特征值求解的过程。下面通过对二维Bravais格子的描述来进一步说明如上问题。
图1-1 二维Bravais格子示意图（格矢、倒格矢和直角坐标系单位矢）
考虑直角坐标系下的二维波动问题。其中，体积力分量为
表示，
为材料弹性系数张量，则二维线弹性介质的动力问题可用以下三类控制方程来描述。
几何方程（梯度方程）:
(1.1)
运动方程（散度方程）:
(1.2)
本构关系（物理方程）:
(1.3)
将式（2.2）带入式（2.3），得到用位移表示的应力，然后再带入式（2.1），控制方程简化为用三个位移分量表示的运动方程（纳维方程）：
(1.4)
当忽略体积力时，有：
(1.5)
引入线性微分算子
后，式（2.5）可改写为如下矢量形式：
； 
(1.6)
其中，
且此处位移矢量已经过分离变量处理，因此之后公式推导中位移仅为位置的函数而与时间无关。回看图1-1所示二维Bravias格子示意图，对于二维平面问题，任一格点可由格矢
和直角坐标系单位矢量
唯一确定：
。而
为倒易空间中的格矢（倒格矢），与正格矢满足关系：
。周期结构的材料属性同样有如下周期性：
； 
(1.7)
其中
为整数。由于结构的周期特性，位移场解形式可设为：
，其中
为Bloch波矢，是倒格矢的线性组合：
，将位移解形式带入式（1.6），忽略体力，则该二维周期问题的动力学方程如下所示：
(1.8)
其中：
(1.9)
且
为倒易空间中的矢量（倒格矢的线性组合）。
对于给定波矢，求解式（1.8）得到特征值
，就可获得对应的频率，从而的到该周期结构的色散曲线。因此在第一段中说本实验是要解决力学中周期性结构特征值求解问题。在之前的设计工作中，一般采用有限元方法，利用大型商用软件COMSOL进行求解，本实验用于探索利用神经网络取代方程求解过程的可行性，最终要得到的结果是对特定结构声子晶体特征频率的预测值，是一个数组。
2 形式化
详细的问题转化和目的在第一节中已经进行了阐述，以一维问题为例，如图2-1(a)所示即为一维层状周期结构示意图。晶格常数为a，且有
，则倒格矢为：
。下考察其带隙结构。当平面纵波沿
方向传播，则周期结构仅沿
方向有位移，不为零的应力分量为
，则运动方程和本构方程为：
； 
(2.1)
对于一维二组分周期结构，通过传递矩阵法可以获得其频散方程如下：
(2.2)
其中
是第i层的厚度，
是该层密度，
则是该层纵波速度，让波矢扫描第一布里渊区，就得到了如图2-1(b)所示带隙图。因此，当给定波数值，其对应的频率是与声子晶体的特征值有关的。而图2-1(b)是通过求方程（2.2）的数值解获得的，其中，所用的材料参数值如表2-1所示，其中A是混凝土，B是橡胶。而本实验的目的是通过对神经网络的训练（对于图2-1(a)所示一维模型，样本就是不同杨氏模量、密度的材料A、B构成的一维层状周期结构）来替代求解方程（2.2）的过程。
(a) (b)
图2-1 (a)一维层状周期结构示意图 (b)图(a)所示声子晶体带隙图
表2-1 材料A（混凝土）和材料B（橡胶）的材料参数
材料 | 杨氏模量 E(Gpa) | 泊松比 v | 密度
|
混凝土 | 30 | 0.2 | 2500 |
橡胶 | 1.37×10-4 | 0.463 | 1300 |
故而上述一维声子晶体特征值问题可以替代为输入—输出关系，其形式为如下格式：
(2.3)
其中，
是通过非线性函数组
（即不同样本）获得的由特征值组成的向量。其中是随空间位置变化的材料参数（
）和所选波数的函数。本实验目的在于通过单胞材料属性来获得其特征值的近似值。故而首先对单胞长度进行归一化处理，然后将归一化后的长度
分成N份。给N小份分别定义材料参数以作为输入变量，从而输入—输出关系变为如下形式：
(2.4)
其中
是对连续性输入—输出关系的近似。在之后的分析中，取
，并且从训练数据集中通过神经网络学习和预测。
此外，因为本实验的最终目的是对声子晶体特征频率进行一个预测，因此是一个回归任务。
3 数据集
首先，我们通过一组训练数据集学习和预测。训练数据由输入、输出数据集构成且通过第2节中所述的一维示例解析解（即求解式（2.1）（2.2）（2.3））获得。其中，输入数据由随机生成的材料参数向量
、
构成，输出数据则为归一化处理后的所选波数
相应的前两阶特征频率，其中值在[0,0.5]范围内选取十个点进行计算。因此输出向量长度实际上为20，输入向量长度实际上为200，其中100个为弹性模量值，剩余100个为密度值。
对于卷积神经网络，输入数据同样遵循上述物理规则生成。其中，材料属性向量和通过杨氏模量和密度的均等概率分布生成。其中，概率分布通过杨氏模量上下界（100MPa~300GPa）和密度上下界（800kg/m3~8000 kg/m3）生成。随后，随机生成100个单胞
并用杨氏模量
和密度
来表征。对于上述生成的单胞，计算前20个特征值构成输出向量
。这个数据集被称作数据集A用作之后参考，总样本数为30万个。在各自的数据集中，一小部分用于模型的训练、验证和初始阶段的网络测试。且由于数据集A中的杨氏模量和密度是从两个独立的概率分布分别采样的，任一给定单元的材料属性可能不对应于任何一种真实材料。
4 算法
4.1 多层感知机
MLP是全连接的多层神经网络，它的目标是寻求输入
和输出
之间的近似函数关系（也是本实验的任务）：
(4.1)
对于如图4-1所示的单隐藏层MLP，因为隐藏层和输出层都是全连接的，故而有隐藏层权重
和输出层权重
，那么（4.1）所述输入—输出关系可表示为：
(4.2)
权重初始值通过随机生成，之后再经优化。网络训练方法就是将已知输入—输出相联系，并且通过优化算法调整权重值使输出误差最小化。
图4-1 单隐藏层MLP示意图
4.2 卷积神经网络
卷积神经网络与传统神经网络有三个主要区别：使用卷积运算；有效激活函数（ReLU）,以及系统化空间不变性。首先一点，卷积神经网络区别于全连接MLP的最明显的特征是使用卷积运算代替标准的矩阵乘法运算。对于图4-1所示的情况，中间隐藏层的值由权重和输入作简单的矩阵乘法运算获得；而对于卷积神经网络，输入将与卷积核做卷积运算得到特征图
，假设输入是长为N，深为d的2D矢量，将其表示为大小为
的三维矩阵。卷积操作在输入空间中由卷积窗口定义好维度的局部区域内进行。核向量可类似的定义成大小为
的矩阵。因此，特征图中的元素可以通过爱因斯坦求和公式：
来计算，其中l指卷积窗口所在当前位置。然后，卷积窗口从当前位置经预定好的步幅前进至下一位置，进行下一特征图的计算。在本实验中，步幅明确为1。上述过程将一直重复直到输出在长度维度上与输入相同并且完成特征图的计算。采用符号表示：
(4.3)
其中
表示卷积运算。例如，在章节2所述1D声子晶体特征值问题中，输入是由深度为2的输入向量
表示的材料属性值的空间有序序列。沿深度方向，第一个元素为所给单元的杨氏模量，第二个元素则对应于密度。卷积层对每个输入矩阵用k个卷积核
进行卷积并得到总共k个特征核。每个特征图
由下式进行计算：
(4.4)
其中
是偏差。CNNs通过调整卷积核和偏差参数来获得所需的输入—输出关系。这里有几点值得注意：首先，由于杨氏模量和密度所对应的特征图是由相同的卷积核运算得到的，因此二者并不完全独立。其次，卷积核使特征图可以表达局部间的相互作用关系，这种稀疏性使得CNNs在学习依赖于空间和事件结构的输入—输出关系使较MLPs更有效。
在计算了预激活特征图之后，将非线性激活函数ReLU作用其上。ReLU作为分段线性激活函数，当输入为负时输出为零，当输入为正时输出等于输入，已被证明有很高的的计算效率。用公式表达为如下形式：
(4.5)
其中
和分别代表ReLU的输出和输入。
当卷积运算和ReLU函数一起被使用后，输出通常经过一个参数减少层，即最大池化层。这一层的输出是p个相邻单元中最大的单元值。在上述一维问题中，沿着转化后的杨氏模量和密度的数据轴进行池化。根据计算机视觉和语言识别中的最优化CNN架构，本实验中池化操作仅在两个连续的卷积层之后操作一次。
4.3 优化算法
此外，不论是对于MLPs还是对于CNNs，都使用了带动量的随机梯度下降法来进行网络训练。训练是小批量进行的，算法通过沿负梯度方向上的小步长来更新模型参数：
(4.6)
其中
表示迭代次数为l时的可优化模型参数，是学习率，是当前小批量的损失函数，是上一梯度步骤当前迭代的贡献。所采用的具体算法为Adam。
4.4 损失函数
前人已证明基于平方误差函数的神经网络可以较准确估计有限样本容量下的后验概率。其中，均方误差损失函数定义如下：
(4.7)
其中
是预测值，是真实值（在本实验中就是特征值），对所有输出和小批量中的数据点进行求和，m是求和项的总数。
4.5 数据归一化
所有输入的材料属性值都通过平均值进行归一化处理：
(4.8)
其中
是输入向量的平均值。此外，特征值也将通过参考最大值进行归一化处理。
4.6 模型架构
在经过超参数迭代后，确定对于上述声子特征值问题的模型有如下架构：对于一维情况，如图4-2所示，100个单胞有100个输入节点，每个输入节点有两条通道，分别对应于单胞的密度和杨氏模量值。在输入节点之后是两个以ReLU函数做激活函数的卷积层，卷积窗口在输入向量大小的维度上大小为3。最大池化层卷积窗口维度为
。随后是两个以ReLU函数为激活函数的全连接层。最后输出（特征值）通过一个线性/高斯连接的全连接层。对于一维问题，基于28000个样本的训练所耗时间小于13分钟（单GPU）。当网络完成训练后，预测时间少于一秒。此外，多层感知机的架构最终确定精度较高的架构含6个隐藏层，每层有1024个计算单元。
图4-2 所用卷积神经网络架构
5 代码
具体代码见源代码压缩包。
6 评价
下面定义预测的平均绝对准确度如下：
(6.1)
其中求和是对所有特征值及数据集进行，p为求和总数。结果表明，卷积神经网络在特征值预测精度方面优于MLP。如图6-1所示可以发现，随着训练样本数的增加，不论是CNN还是MLP预测精度都得到了提高。然而，CNN仅在20000个样本的基础上进行训练，预测准确度已高于95%，而同等情况下，MLP预测准确度仅高于70%。且在使用10000个样本进行训练时，MLP仅能达到70%的准确度而CNN准确度已高于98%。值得一提的是，在30000~40000样本之间，两个网络的平均绝对误差百分比都下降了。但总的来说，对于力学中 特征值预测问题，CNN比MLP更为适合。
图6-1 CNN和MLP预测精度
7 参考文献
[1]Finol D, Lu Y, Mahadevan V, Srivastava A. Deep convolutional neural networks for eigenvalue problems inmechanics. Int J NumerMethods Eng. 2019;118:258–275.
[2]Tan P, Steinbach M, Kumar V. Introduction to Data Mining[M].人民邮电出版社.
