一种基于CMA‑ES优化算法设计声学超材料单元的方法与流程

本发明涉及一种声学超材料单元的设计方法，特别涉及一种基于协方差矩阵自适应进化策略优化算法设计声学超材料单元的方法。

背景技术：

新型人工声学材料是一种由不同弹性模量和质量密度的材料构成的周期性人工复合结构，具有自然界媒质所不具备的性质，可以实现负的等效弹性模量或负的等效质量密度。通过改变单元结构及其空间排列来改变声波传播路径上空间的声压，从而实现人工调控声波。声学超材料由于其非常的参数特性，因此具有许多新颖的声学传播效应，如负折射效应、平面聚焦、声隐身、超透镜效应等。

等效折射率是衡量声学超材料单元性能的重要指标，实现大折射声学单元可以将透镜做薄或者实现聚焦透镜的短焦距，大折射率在超薄声学隐身斗篷上也有很好的应用前景。在声学超材料单元领域，目前主要利用分形和卷曲空间达到增加折射率的目的，且现阶段二维声学超材料单元的等效折射率最高只能达到5左右，因此，寻求大折射率单元已成为许多学者关注的问题。

协方差矩阵自适应进化策略优化算法(简称“CMA-ES”优化算法)是启发式演化算法中的一种，它与所求问题只存在唯一的衔接点——问题的适应度值，这一特点使得算法能够在解决问题的可靠性和消耗的时间这两个量之间进行权衡，能够更加高效快速的解决问题。简单地说，CMA-ES算法就是对问题的参数空间进行高斯采样，并根据一定的选择机制对样本空间的高斯分布进行更新，继而获得新的样本空间，直至迭代次数达到预先设置的最大值或达到其他能够使算法停止运行的预设条件。CMA-ES算法在电磁优化领域有很多的应用案例，得到很多满足特殊性能的电磁超材料结构或阵列，但是在声学超材料领域却几乎没有应用，发明人将CMA-ES算法用于声学领域，形成了本

技术实现要素：
。

发明内容

发明目的：本发明的目的在于提供一种基于CMA-ES优化算法设计声学超材料单元的方法，该方法可设计出具有较高的折射率的声学超材料单元，且其阻抗失配也较小。

技术方案：本发明所述的一种基于CMA-ES优化算法设计声学超材料单元的方法，包括如下步骤：

(1)根据声学超材料单元的设计尺寸确定各个粒子的维度及取值空间，得出样本空间的初始分布，从中选择若干粒子构成父系种群；

(2)父系种群中的每个粒子均是0-1分布的阵列，根据阵列分布在仿真软件中将每个粒子均构建成一个声学单元物理模型，其中，0和1分别代表由空气或光敏树脂构成的声学单元结构的子单元；

(3)对父系种群中的每个声学单元物理模型进行提参，得到其等效折射率和阻抗值，将该值作为适应度函数的变量求得适应度函数值，根据适应度函数值判断父系种群中是否存在满足目标条件的粒子；如存在，算法停止，如不存在，更新种群，进行迭代算法，直至达到预设的停止条件；

(4)最终停止时得到满足条件的粒子，其对应的声学单元物理模型即为符合设计要求的声学超材料单元模型结构。

上述步骤(3)中，适应度函数值为满足阻抗取值范围的每个粒子的等效折射率最大值。

其中，适应度函数的表达式为：

Cost＝max(n(Xi))&(ξ(X_i)＜ξ_upper)；

对于父系种群中的粒子X_i，如果其阻抗值ξ(X_i)满足ξ(X_i)＜ξ_upper，那么其折射率n(Xi)在一定频段内的最大值max(n(Xi))即为其适应度函数值；如果其阻抗值ξ(X_i)不满足ξ(X_i)＜ξ_upper，那么这个粒子的适应度函数值为0。

较优的，步骤(3)中，先通过有限元分析方法对种群中的每个声学单元模型结构进行仿真，得到其反射系数和透射系数，然后由反射系数和透射系数求得其等效折射率和阻抗。

进一步的，等效折射率n根据下式确定：

式中，m为分支选择数，R、T分别为声学超材料单元结构的反射系数和透射系数，d为声学超材料单元在声波传播方向上的厚度；k＝2πf/c₀，其中f为频率，c₀为环境媒质中的声速。

具体的，步骤(3)中，当父系种群中不存在满足目标条件的粒子时，从父系种群中选取若干个最优的粒子构成子孙种群，得出这个种群的初始分布，获得新的满足高斯分布的样本空间并重复前述步骤、直至达到预设的停止条件。

优选的，采用八分之一对称优化算法设计声学超材料单元结构，即将预设计的声学超材料单元结构划分为八部分，通过CMA-ES优化算法优化得到其中一部分的单元结构，将该部分的单元结构翻折对称得到其余结构。采用八分之一对称优化，可只对一个单元的八分之一部分进行优化，可大大缩短优化算法的运行时间，高效的寻找到符合要求的声学单元结构。

与现有技术相比，本发明的显著优点为：(1)本发明的声学超材料单元设计方法可设计出具有较高折射率的声学超材料单元，其折射率远高于现有的二维声学超材料单元；(2)本发明的设计方法设计出的声学超材料单元结构，在工作频点上不仅折射率较高，其阻抗匹配也较理想；(3)本发明提供的设计方法基于协方差矩阵自适应进化策略优化算法，只需预设尺寸维度，即能得到预设计的声学超材料单元，制备方法智能高效；(4)本发明提供的设计方法同时具有可拓展性，能够进一步的实现具有其他参数特性的声学单元。

附图说明

图1为本发明的基于CMA-ES优化算法设计声学超材料单元的方法流程图；

图2为基于八分之一对称算法将声学超材料单元分为呈翻转对称的八部分的结构示意图，图中，标示“I”的部分代表空气，标示“II”的部分代表光敏树脂材料；

图3(a)为本发明中制得的声学超材料单元的结构示意图；

图3(b)为图3(a)的声学超材料单元的简化结构示意图；

图3(c)为图3(a)和图3(b)的声学超材料单元结构的等效折射率曲线；

图3(d)为图3(a)和图3(b)的声学超材料单元结构的等效阻抗曲线；

图4为入射波垂直打入本发明中设计的45°三棱台的声场分布图；

图5(a)为本发明中设计的离散偏折透镜在不同位置处的折射率分布曲线及1530Hz时对应的声学超材料单元结构的伸缩率；

图5(b)为本发明中设计的离散偏折透镜的原型示意图；

图6(a)为声波垂直入射在本发明设计的理想偏折透镜上形成的声场分布图；

图6(b)为声波垂直入射在本发明设计的离散偏折透镜上形成的声场分布图。

具体实施方式

下面结合附图对本发明的技术方案做进一步说明。

本发明的声学超材料单元的设计方法是基于协方差矩阵自适应进化策略(CMA-ES)优化算法设计实现的。声学超材料单元结构一般由若干子单元构成，这些子单元构成了一个阵列整体单元结构；每个子单元由光敏树脂或者空气构成，可以用数字“1”和“0”表示子单元填充或者不填充材料，在设计中，利用CMA-ES优化算法可以对0-1排布的阵列进行优化，每一个尺寸维度的0-1阵列，都可以对应一种声学超材料单元结构，因此可利用CMA-ES优化算法进行编码优化来设计得到声学超材料单元的结构。

在优化过程中，每一个声学超材料单元结构的参数特性可以通过有限元分析方法计算提取出来，作为优化算法中适应度函数的变量；通过对适应度函数的值进行优化，最终可以得到最优的、满足设计要求的声学超材料单元结构。

具体的，如图1，本发明的基于CMA-ES优化算法设计声学超材料单元的方法，包括如下步骤：

(A)根据声学超材料单元的具体尺寸确定各个粒子的维度和各个维度的取值空间，从而计算出样本空间的初始分布(或者直接给出样本的初始分布)；

以设计尺寸为15mm*15mm的声学超材料单元为例，每个子单元为大小为0.5mm*0.5mm的正方形，使用的材料光敏树脂的密度为1300kg/m³，声速为716m/s；同时，空气的密度和声速分别为1.29kg/m³和343m/s。

由于优化设计的时间与单元结构的维数，即子单元的个数正相关，因此减少一个单元结构的维数可以有效的降低仿真运算时间，因此，优选采用八分之一对称优化算法，在确定每个粒子的维数时，八分之一算法可以将一个粒子的维数缩小为原来的八分之一左右，减少程序运行时间。如图2，左上角加粗三角形内的部分即为需要优化的部分，剩下部分的结构均可通过该三角形内的结构翻折对称得到。本设计采用八分之一对称优化算法，此时一个粒子的维数为120，每个粒子的每一维的取值为0或1。

给定粒子的维数和每一维度取值范围后，可以计算出这个种群的初始分布，即期望值X和方差σ。

(B)在满足高斯分布的初始样本空间中随机选择一定数量的粒子(λ个粒子，用X_i表示，i＝1,2,…,λ)，其中反映高斯分布幅度的方差由各个维度的协方差矩阵构成，这些粒子构成父系种群，父系种群的粒子期望值和方差满足初始分布给定的X，σ；

(C)父系种群中的每个粒子均是0-1分布的阵列，每个0-1阵列均对应着一个声学超材料单元的结构，在仿真软件中根据每个粒子的阵列分布构建对应的声学单元物理模型；其中，0代表由空气构成的声学超材料的子单元，1代表由光敏树脂构成的子单元；

(D)对种群中的每一个物理模型求其反射系数和透射系数，由反射系数和透射系数得到其等效折射率和阻抗；

具体的，每个粒子对应的声学超材料单元的等效折射率n由下式确定：

其中，m为分支选择数，R为声学超材料单元结构的反射系数，T为声学超材料单元结构的透射系数，k＝2πf/c₀，其中f为频率，c₀为环境媒质中的声速，在本设计中，环境为空气，空气中的声速c₀为343m/s,d为声学超材料单元在声波传播方向上的厚度，在本设计中为15mm。

每个粒子对应的声学超材料单元的阻抗

(E)将每个粒子的等效折射率和阻抗值作为适应度函数的变量，根据适应度函数值判断种群中是否存在粒子的适应度值满足目标条件，由判断结果决定是否进行下一步的迭代，若不需要迭代算法停止；

在本发明中，适应度函数为Cost＝max(n(Xi))&(ξ(X_i)＜ξ_upper)，即父系种群中满足阻抗取值范围的每个粒子的折射率的最大值。对于父子种群中的粒子X_i，如果其阻抗值ξ(X_i)满足ξ(X_i)＜ξ_upper，那么其折射率n(Xi)在一定频段内的最大值max(n(Xi))则作为是适应度函数值；如果其阻抗ξ(X_i)不满足ξ(X_i)＜ξ_upper，那么这个粒子的适应度函数值为0。

设置本设计中该算法的停止条件为适应度函数值达到7，即如果存在阻抗小于ξ_upper，折射率的最大值max(n)>7，那么算法停止，这个满足条件的粒子则作为最优解，对应的声学超材料单元即为本设计所需的单元。

(F)若需要迭代，从父系种群中选取前面若干个相对最优的粒子构成子孙种群，得出新的期望值与方差，获得新的满足高斯分布的样本空间并重复上述过程至最终达到预设的停止条件；

(G)算法停止运行。

经过优化设计，最终得到的声学超材料单元的结构如图3(a)，其中，空白部分代表空气，阴影部分代表光敏树脂材料部分，经测试，该单元结构的等效折射率在1530Hz时达到7.46。

为了美观和阻抗匹配，将这个结构简化，在简化单元结构的过程中，秉持折射率基本不变，阻抗值尽量降低的原则。首先将图3(a)中四个拐角处离散的子单元删除，再提取它的等效折射率和阻抗，等效折射率的最大值从7.42降低到7.37，和图3(a)中的单元结构相比降低了0.7％；再者，考虑到加工上的工艺和美观的效果，将图3(a)中间离散杂乱的子单元用一个正方形环状结构代替，此时折射率的最大值达到7.56，比原始结构的折射率增加了1.9％；最后，将光敏树脂材料和空气的交界面上的毛刺去除，即采用平滑的弧形代替原来的梯形边界，这时折射率的最大值为7.46，和原始结构的最大折射率相比，差值不到0.6％，而其等效阻抗值却在工作频率1530Hz时降低了很多。

经过简化，最终得到图3(b)中的声学超材料单元结构，图中，“I”代表空气部分，“II”代表光敏树脂材料部分；简化后的声学超材料单元结构的尺寸为：a＝15mm，a0＝14.5mm，L1＝3.5mm，L2＝6.4mm，L3＝3mm，W＝1.5mm，R＝5.2mm；其中，a为一个声学超材料单元的厚度，a0为一个声学超材料单元中光敏树脂材料的最大长度，L1为声学超材料单元中间空气方环的内正方形长度，L2为声学超材料单元中间空气方环的外正方形长度，L3为声学超材料单元中光敏树脂材料在四边的中点处挖除的长方形的长，W为声学超材料单元中光敏树脂材料在四边的中点处挖除的长方形的宽，R为光敏树脂和空气的分界面处弧形边界的半径。

对图3(a)和图3(b)的单元结构进行等效折射率和阻抗测试，如图3(c)和3(d)，可以看出，简化后的声学超材料单元结构的等效折射率与简化前无限接近，在1530Hz时折射率达到7.46；可以证明，简化后的单元结构完全具备原单元结构的折射性能，满足设计要求，而且在工作频点1530Hz时，简化后的单元结构的阻抗失配性较原单元结构小很多。

为进一步验证设计的正确性，本发明设计了两个实验。

实验一：本发明设计了一个45°的劈尖三角形，其以图3(b)中简化后的声学超材料单元结构组成，如图4，入射声波从左边垂直入射到三角形的直角边上，通过结构折射到介质环境中，介质环境设置为等效折射率为7.46的理想媒质，可以看到，声波在三角形的斜面上没有发生偏折，间接地证明了图3(b)的声学超材料单元的等效折射率在工作频点上是7.46。

实验二：本发明设计了一个离散的偏析透镜，如图5(a)，根据设计的偏析透镜在不同位置处的折射率分布，确定不同位置处对应的声学超材料单元结构的伸缩率，然后在透镜的不同位置处按照对应的折射率大小放置对应伸缩率的声学超材料单元结构，那么即可构成所需要的偏折30度的偏折透镜，如图5(b)；分别对得到的离散偏析透镜和理想的偏析透镜进行测试，入射声波从左面入射，垂直入射到偏折透镜上，在透镜与空气的交界面上偏折30°，如图6(a)和6(b)，比较理想透镜和离散透镜的仿真图，两者在误差允许的范围内一致，这也说明本发明的声学大折射率超材料单元的设计的正确性。