基于力学原理的逆向设计方法，用于控制纤维增强双层的褶皱现象

《Materials Today Communications》：Mechanics-informed inverse design for wrinkle control in fiber-reinforced bilayers

　　陈晓毅|王浩|王新宇|甘佩军|戴俊龙|舒天宇|陈宇轩|杨凡泽|刘永涛北京师范大学-香港浸会大学，广东省数据科学与应用跨学科研究重点实验室，珠海519087，中国

摘要
控制纤维增强薄膜-基底系统中的皱纹对于从抗衰老护肤到皮肤替代品工程等应用至关重要。虽然皱纹力学在正向预测方面已经发展得相当成熟，但针对微观结构控制的逆向设计策略（如纤维预应变和取向）仍然有限。在这里，我们开发了一个基于力学的逆向设计框架，以调节受皮肤启发的纤维增强双层的皱纹起始和形态。使用半解析的形态弹性板模型，我们量化了五个设计变量（薄膜厚度比h、基底刚度比、纤维角度、纤维预应变参数和纤维体积分数）如何控制临界皱纹阈值和临界模式数。这种降阶公式使得在高纤维体积分数的情况下也能进行高通量探索，而这通常是标准3D有限元模拟所面临的挑战。全局敏感性分析（Spearman相关性和XGBoost）以及局部力学扰动分析表明，纤维预应变是主要调节因素，而基底刚度和纤维取向主要控制皱纹的形成。对于逆向设计，我们实施了两种互补的方法：(i) 基于力学的正则化优化，直接利用半解析稳定性条件并实现物理指导的参数简化；(ii) 基于模型的数据训练的替代辅助遗传算法（BPNN–GA）进行可扩展的全局搜索。所得到的设计为中年和老年人的基础皮肤模型提供了与物理一致的调节策略，为个性化的、受皮肤启发的材料设计提供了定量指导。

1. 引言
皮肤皱纹是衰老的普遍标志，从力学角度来看，可以认为是由皮肤分层结构中的刚度和预应变不匹配引起的不稳定性：相对较硬的表皮/真皮层覆盖在较软的真皮层上。与年龄相关的微观结构退化——包括胶原纤维损伤、预应变减少和真皮层变薄——在生理生长、重塑或压缩过程中进一步促进了表面皱纹的出现。这些观察结果引发了一个实际重要的逆向问题：应该针对哪些力学和几何参数来延缓皱纹的出现并减少皱纹的密度（模式数）？解决这一逆向设计问题对于指导定量抗衰老策略和设计仿生合成皮肤替代品至关重要。
在皮肤皱纹的正向建模和测量方面已经取得了实质性进展。赵等人[1]将皮肤建模为多层Neo-Hookean系统，并通过实验校准后，使用有限元模拟来量化材料和几何因素如何影响皱纹深度。Lejeune等人[2]提出了一种处理分层介质中皱纹的算法方法，通过将多层简化为等效的复合梁来获得可行的预测结果，该方法应用于表皮电子学。Kuhl及其同事还开发了针对患者的特定模拟流程，他们将临床成像与生物力学生长定律相结合，以虚拟评估手术选项，如切口布局和组织扩张器放置[3]、[4]、[5]、[6]、[7]、[8]、[9]。与此同时，体内表征也得到了改进：Machado等人[10]使用立体摄影测量法量化了皱纹深度，而Shah等人[11]和Moon等人[12]发布了大规模的面部数据集，使得群体水平分析成为可能。尽管有这些进展，大多数贡献仍然是正向的——预测给定参数下的变形和皱纹指标——而实际的逆向设计能力仍然有限，主要是因为重复评估非线性正向求解器（通常是FEM）在计算上非常昂贵。
来自双层和薄膜-基底不稳定性的力学提供了互补的基础，这些力学为皱纹阈值和模式选择提供了成熟的分析和计算工具[13]、[14]、[15]、[16]、[17]、[18]、[19]、[20]、[21]、[22]、[23]、[24]、[25]。在接近皱纹起始时，线性分岔理论可以得出临界条件（例如，临界应变；见[26]），而弱非线性分析可以捕捉振幅演变[13]。在远离皱纹起始或出现次级分岔时，通常需要完全非线性模拟。在这个范围内，Cao和Hutchinson[18]量化了Neo-Hookean双层的临界皱纹应变，并使用FEM揭示了丰富的分岔后情景。Holland等人[24]系统研究了多种加载模式和刚度比，并表明在低刚度对比下，加载模式的敏感性在使用有效度量时可以大大降低。对于纤维增强双层，Liu等人[25]结合了增量有限弹性与线性和弱非线性分析，推导出了精确的分岔条件和解析振幅方程，强调了纤维刚度、纤维取向和层间模量比如何共同控制超临界到亚临界的转变。在生物驱动的几何形状中，也使用结合分岔分析和FEM研究了生长驱动的不稳定性；例如，Ciarletta等人[21]为管状双层构建了相图。从板理论的角度来看，Audoly和Boudaoud[15]表明等轴加载可以使直线皱纹变得不稳定，形成棋盘图案，而各向异性加载可以触发向波浪条纹的次级转变。虽然这些框架在正向预测方面非常强大，但它们本身并不能提供在实际约束下选择延迟不稳定性和引导模式数的微观结构和几何参数的逆向设计方法。
数据驱动的、与物理一致的建模越来越多地用于研究软组织力学，这是由于数据集的增长以及使用纯解析理论或高保真模拟难以捕捉异质性、各向异性和多尺度非线性耦合。特别是神经网络已被嵌入到力学流程中作为替代物或本构表示；例如，Lin等人[27]开发了一个基于能量的、物理信息的神经网络，用于生长驱动的形态发生。然而，逆问题仍然具有挑战性，因为关键参数通常是潜在的，只能通过变形或表面图案间接观察到。因此，参数识别通常成为一个需要许多正向求解器的嵌套优化[28]。最近的努力包括基于PINN的逆向公式[29]、[30]，以及在设计环境中结合遗传算法和学习到的替代物的混合工作流程[31]。然而，在皮肤病学驱动的力学中，大多数干预策略仍然强调粗略的宏观描述符（例如，层厚度和体积刚度），而明确调节微观结构级别控制（如纤维结构和预应变）的逆向设计仍然相对较少探索。
在这项工作中，我们开发了一个高效的逆向设计框架，用于延缓纤维增强薄膜-基底系统中皱纹的出现并调节皱纹起始模式数，该框架模仿了皮肤层的关键方面。由于皮肤相对于其平面尺寸来说很薄，我们采用了降维的板公式来捕捉有限应变变形和不稳定，其成本远低于3D连续体FEM，从而能够进行大规模参数扫描和设计探索。我们的方法建立在Dai及其同事[32]、[33]、[34]、[35]开发的非线性、渐近一致的板理论基础上，并由Steigmann[36]进行了综述，并利用了与生长和重塑相关的形态弹性扩展。对于材料行为，我们注意到各向同性超弹性模型（Neo-Hookean、Mooney–Rivlin、Ogden）是皮肤变形的常见初始近似[37]、[38]、[39]、[40]、[41]、[42]、[43]，但真皮力学受到胶原纤维增强和各向异性的强烈影响，这促使人们开发了纤维增强超弹性框架，如HGO型模型[44]、[45]、[46]。尽管粘弹性可能很重要（例如，准线性粘弹性）[47]、[48]、[49]、[50]、[51]，但我们关注的是弹性和生长驱动不稳定性的起始。因此，我们开发了一个纤维增强形态弹性板模型，在该模型中明确引入了微观结构纤维描述符，包括纤维刚度、取向和体积分数以及纤维预应变参数（用表示）。纤维的贡献遵循标准的、成熟的纤维增强建模方法[44]、[52]、[53]，而形态弹性生长变量提供了额外的、基于生理学的不匹配来源，可以驱动不稳定性。这种组合实现了可调微观结构级别属性与用于逆向设计的起始描述符之间的力学一致联系。
我们的目标是(i) 提高皱纹阈值（延缓起始）和(ii) 降低皱纹起始模式数（减少皱纹密度）。我们在纤维增强形态弹性板框架内使用线性分岔分析推导出和。然后我们生成一个数据集，将五个几何和材料输入映射到两个输出，并使用它进行Spearman等级相关性、XGBoost回归和训练反向传播神经网络（BPNN）替代算法。通过将这些数据驱动的分析与力学模型的局部敏感性分析相结合，我们确定了最强烈影响和的参数。对于逆向设计，我们考虑了代表中年和老年组的两个基线参数集，并寻找满足预定多目标配置（例如，将增加或同时减少）。我们实施了两种互补的方法：一种是在半解析板模型上直接进行的基于力学的优化，另一种是使用替代辅助的遗传算法（BPNN–GA）。尽管对于当前的半解析设置直接优化是可行的，但我们的机器学习辅助流程旨在用于实际应用中的可扩展设计，其中分析解决方案不可用（例如，患者特定的3D几何形状或异质性损伤），并且每次评估可能依赖于昂贵的FEM或实验数据。
这项研究是我们JMPS（《固体力学与物理学杂志》）论文[54]的补充，该论文专注于纤维增强薄膜-基底皱纹的形态弹性板模型的正向公式化和验证，包括与数值模拟/实验的比较。相比之下，当前的工作针对模型的逆向使用：我们(i) 进行全局和局部敏感性分析以确定可操作的微观结构控制参数，(ii) 开发和基准测试两种互补的逆向设计方法——一种是基于力学的正则化优化，利用半解析稳定性条件；另一种是数据驱动的替代辅助GA（BPNN–GA）。为了支持可重复性和重用，我们还发布了完整的数据集生成和逆向设计脚本。
本文的其余部分组织如下。第2节介绍了纤维增强形态弹性板模型（详细信息在附录中）和数据集生成过程。第3节展示了Spearman相关性和XGBoost回归结果。第4节详细介绍了BPNN架构和训练过程。第5节讨论了参数对和的影响。第6节报告了基于力学优化和BPNN-GA获得的逆向设计结果。第7节总结了本文。

2. 力学框架和数据生成
2.1. 形态弹性板模型的发展
我们将在本小节简要介绍形态弹性板模型和主要结果（更详细的推导在附录A中给出）。
我们考虑在生长加载条件下薄膜基底系统的皱纹（见图1）。薄膜是具有两种纤维族的Neo-Hookean材料，基底被视为Winkler基础。变形梯度的分解如下[55]（1），弹性变形的不可压缩约束为（2）。名义应力由（3）给出。
下载：下载高分辨率图像（200KB）
下载：下载全尺寸图像
图1. 薄膜-基底系统的示意图。薄膜被建模为放置在Winkler基础上的纤维增强板。
三维平衡方程和边界条件为（4），其中该薄膜-基底系统在两端受到约束，薄膜在其表面由Winkler基础支撑，并且表面无牵引力。
我们考虑在平面应变条件下的系统：。为了推导形态弹性板模型，我们假设位置和拉格朗日乘数足够平滑，可以将它们展开为泰勒级数（5）。将这些展开式代入运动学关系、本构方程和3D控制系统中，并进行适当的代数操作，我们得到了展开系数的递归关系以及简化后的板方程。递归关系为（6），其中是名义应力张量的展开系数，表示的元素，表示的展开系数。这些递归关系用于消除和的更高阶展开系数（即）。
板方程为（7）。将递归关系代入板方程后，我们得到了主要未知数的常微分方程（ODEs）。边界条件为（8）。
（7）、（8）的线性分岔分析得出了生长加载的临界条件，由（9）给出，从中我们得到（10）。
给定模式数，我们可以确定相应的临界生长因子。在所有可能的模式中，具有最小临界生长因子的模式首先被激发，从而定义了实际的临界生长因子（表示为）和相应的临界模式数（表示为）。对进行最小化得到临界模式数为（11）。
将（11）代入（9），我们得到满足（12）的相应临界生长因子。
从（11）、（12）中，我们预测了和关于参数的理论曲线。我们在之前的论文[54]中已经表明，我们的理论预测通过数值模拟和实验得到了很好的验证，见图2。参数的解释见表1。
下载：下载高分辨率图像（425KB）
下载：下载全尺寸图像
图2.关键生长因子和关键模态数作为纤维角度的函数。结果在理论预测、数值模拟和实验实现之间进行了比较。摘自我们的配套论文[54]。表1. 本研究五个输入和两个输出的解释。参数 物理解释 薄膜厚度与薄膜长度的比率 Winkler基底刚度与薄膜刚度的比率 纤维对齐角度 纤维预张力参数 纤维体积分数 皱纹阈值 皱纹模态数通过我们的模型，之前研究[54]的敏感性分析表明，是控制皱纹阈值和分叉后振幅的主要参数，其敏感性至少比、和高出一个数量级。我们还发现，最大化阈值的参数集几乎与最小化皱纹振幅的参数集相同。因此，在本研究的逆问题中，我们只关注优化皱纹阈值和皱纹模态数。与我们之前的工作[54]相比，其中纤维体积分数通常在到之间，当前研究考虑了一个更宽的范围，这是由于皮肤组织的特性，纤维体积分数可能取更高的值（高达）。在密集纤维区域，观察到的敏感性趋势是否仍然存在仍然是一个未解决的问题。在这种高纤维区域，有限元方法（FEM）面临重大挑战：严重的网格变形和计算奇异性使得标准的3D数值模拟在计算上变得不可行或无法收敛。因此，一个稳健的半解析框架对于探索这种纤维密集区域中的薄膜-基底系统的皱纹行为至关重要，它填补了数值方法失效的空白。为了展示实际应用价值，我们进行了逆参数设计，以探索延迟皱纹出现和减少皱纹模态数的策略。该模型故意保持在概念验证阶段，以便快速探索设计空间；在必要时，可以在未来的工作中加入更复杂的本构特性。我们关注两个年龄组，即中年组和老年组。这些组的基线参数集来自文献，并在表2中总结。以基线值为参考点，我们扩展了表3中列出的范围，以便进行系统性的研究。所选范围受到文献[1]、[56]、[57]、[58]、[59]、[60]中报告的生理相关值的指导。表2. 中年组和老年组的基线参数。组 中年组 0.01 0.22 50.93 0.8 老年组 0.008 0.17 50.96 0.6 表3. 数据集构建的输入参数范围。对于每个参数集，我们进行分岔分析以确定皱纹阈值和皱纹出现时的模态数。我们的逆向设计目标是增加皱纹阈值，同时减少模态数。为了解决这个逆问题，我们采用了多种方法，包括（i）基于机械模型的正则化优化，以及（ii）基于深度学习的替代建模。对于数据驱动的方法，我们通过采样参数空间并使用板模型的半解析预测来构建训练数据集。备注：根据我们之前的基于力学的分析，在感兴趣的参数范围内，无量纲纤维刚度参数对关键载荷和选定的模态数的影响可以忽略不计。此外，对不稳定性响应的影响与厚度参数的影响大体一致；因此，将和视为独立的设计变量会引入冗余。为了减少设计空间的有效维度，我们在整个研究中固定和，并专注于随后的参数采样和替代建模。2.2. 数据集构建基于第2.1节中得出的半解析结果，我们高效地生成了一个数据集，将关键模型参数映射到两个皱纹描述符：皱纹阈值和相应的模态数。每个数据记录存储为。数据集构建在两个方面利用了机械模型：1. 输入参数的选择。根据我们的敏感性分析，我们筛选出对皱纹影响可以忽略的参数，只保留最有影响力的变量作为输入。具体来说，输入集是，其中参数的解释列在表1中。2. 输出参数的选择。原始的设计目标包括提高皱纹阈值（），减少模态数（），以及减小皱纹振幅（）。机械模型表明，在考虑的参数范围内，最大化几乎同时发生与最小化。因此，我们选择和作为主要输出，这减少了学习和优化的复杂性，同时不损害本研究的设计目标。我们使用表3中列出的参数范围来确定数据集的大小。在五维输入空间中采用了两种采样策略：均匀网格采样和随机采样。均匀网格采样产生了有效的数据记录，包括：，，，和。均匀网格总共包含组合；然而，由于自动化批量求解器偶尔无法产生有效的（，），因此丢弃了组合。随机采样是通过Mathematica生成的，提供了表3中范围内的记录。此外，为了提高两个年龄组的替代精度和逆向设计性能，我们在两组的基础上进行了局部随机采样，增加了记录。总共构建了形式的数据记录。这个数据集用于两个目的。首先，我们将统计分析和机器学习模型与基于力学的敏感性分析结合起来，以确定控制和的最有影响力参数。其次，我们利用这些结果进行逆向设计，使用（i）基于力学的优化和（ii）基于神经网络模型的替代辅助优化，结合遗传算法（BPNN–GA）。为了进一步阐明数据集的高维分布特征，采用了t-分布随机邻域嵌入（t-SNE）将五维输入空间投影到二维流形上。如图3所示，对于“网格”和“随机”数据集，t-SNE算法有效地保留了局部邻域结构，揭示了不同的拓扑模式。在网格数据集中，空间排列显示出高度有序的几何规律性，反映了参数采样的系统性，如图3(a)和(b)所示。相反，随机数据集显示出更随机的分布，同时保持了全局聚类，如图3(c)和(d)所示。下载：下载高分辨率图像（849KB）下载：下载全尺寸图像图3. 网格和随机数据集在两个性能指标上的t-SNE。（a）网格数据集在。（b）网格数据集在。（c）随机数据集在。（d）随机数据集在。最后，关于计算成本，数据集的生成是使用标准的个人计算设施完成的。作为一个代表性的基准，在配备Apple M2 Max处理器的机器上，评估皱纹阈值和模态数的半解析解每个参数集大约需要到2秒。整个数据集的生成，包括超过10,000个参数组合，分布在几台标准个人计算机上完成。在初始计算之后，进行了严格的验证步骤，以过滤或重新计算任何非物理结果，从而确保最终数据集的高质量和可靠性。3. 统计分析和机器学习方法3.1. 斯皮尔曼相关性我们使用斯皮尔曼相关性分析来量化五个输入变量和两个输出之间的单调关联。斯皮尔曼相关系数范围从到：表示正相关，表示负相关，衡量关联的强度，对应于完美的单调关系。图4显示了斯皮尔曼相关热图。输入变量显示出弱的成对相关性，表明采样参数空间中的多重共线性有限。此外，输入-输出相关性是可以解释的，并且与我们的模型预测的机制趋势大体一致，支持使用该数据集进行后续的替代建模。下载：下载高分辨率图像（249KB）下载：下载全尺寸图像图4. 输入参数和皱纹输出的斯皮尔曼相关系数。对于关键载荷，主导因素是纤维预张力参数，它与显示出强烈的负相关。因此，在考虑的范围内，增加预张力（减小）是提高的最有效方法。对于关键模态数，两个因素最有影响力：纤维角度与强烈正相关，而厚度比率与强烈负相关。因此，为了减小，可以增加薄膜厚度（更大）和/或减小纤维取向角度（更小）。3.2. XGBoost回归我们使用XGBoost回归（极端梯度提升）来建模输入参数和目标输出之间的关系。XGBoost是一种基于梯度提升的树基集成方法：它通过依次添加浅层决策树来构建强大的预测器，每个新树都训练用于减少当前集成的预测误差。与简单的线性模型相比，XGBoost可以捕捉到力学和材料相关数据集中常见的非线性参数-属性关系和特征交互。它还结合了正则化和高效的数值实现，有助于减轻过拟合并提高计算效率。在这项研究中，我们训练了两个独立的XGBoost回归器，一个用于，另一个用于。我们使用完整数据集的训练-测试分割来评估预测性能。性能使用均方根误差（RMSE）和决定系数（）来报告，分别由（13）和（14）定义，其中是样本数量，和表示第个样本的真实值和预测值（其中代表或），是样本均值。RMSE和的结果总结在表4中。总体而言，XGBoost替代模型在训练集和测试集上都取得了高精度。这是预期的，因为数据集是从半解析机械模型生成的，并且基本上没有噪声，产生了可以通过灵活的回归器很好地近似的平滑输入-输出映射。因此，测试集指标提供了在采样参数范围内的泛化能力的有意义度量。第5节提供了进一步的分析和讨论。表4. XGBoost在训练集和测试集上的预测性能，分别针对和。变量 训练RMSE 测试RMSE 0.9999 740.048 0.048 0.580.9912 290.827 3390.9999 500.0028 110.9934 880.033 438. 使用BPNN进行的前向预测本工作中使用的反向传播神经网络（BPNN）是一个使用反向传播算法[61]训练的前馈多层感知器。通过将信息从前向传播到互连的神经元，并在训练期间将预测误差向后传播，BPNN学习输入向量与输出向量之间的非线性映射：（15）其中输入向量是，输出向量是。网络架构如图5(a)所示。输入层包含五个神经元，对应于五个输入参数，输出层包含两个神经元，对应于和。在输入层和输出层之间放置了三个全连接的隐藏层，每个隐藏层有神经元。连接仅存在于相邻层之间（全连接），同一层内没有连接。下载：下载高分辨率图像（391KB）下载：下载全尺寸图像图5. BPNN的结构。（a）本工作中使用的BPNN模型，包括输入层、输出层和隐藏层。（b）神经元是来自前一层的权重、偏置和神经元的线性组合的非线性激活函数。如图5(b)所示，每个神经元首先对前一层的输入进行仿射变换（加权和加上偏置），然后对结果应用非线性激活函数（16）其中是第层中神经元的激活函数，是连接层中神经元到层中神经元的权重，是偏置项。在这项研究中，采用Sigmoid函数作为激活函数。训练旨在通过迭代更新权重和偏置来最小化预测输出和真实输出之间的差异。数据集包含11,640个样本，并随机分为训练、验证和测试子集，比例分别为70%、15%和15%。对于这个回归任务，使用均方误差（MSE）损失：（17）其中是样本数量，表示欧几里得范数。决定系数（）也用于评估预测性能。由于当前问题有两个输出，是在二维输出向量上联合计算的：（18）其中表示数据集上真实输出向量的均值。MSE和的定义在（13）、（14）中，输出是向量形式，而对于和，输出是标量。MSE和的结果将在第6.5节给出。参数的特征重要性基于第3节观察到的相关性趋势，本节进一步量化了中年组和老年组参数的相对重要性。我们采用了两种互补的方法：一种是基于力学模型的局部逐一（OAT）敏感性分析，另一种是基于XGBoost替代模型的全局特征重要性分析。重要的是，OAT研究通过评估特定基线周围的扰动来关注局部敏感性。因此，得出的参数排名可能与第3节中的Spearman分析和XGBoost模型（第5.2小节）提供的全球洞察不同。

5.1 来自力学模型的敏感性分析
为了量化每个设计参数如何影响皱纹的出现和模式数，我们使用基于力学的模型对中年组和老年组进行了局部（逐一）敏感性分析。从特定年龄的基线参数集开始（见表2），我们在保持其他参数不变的情况下，将每个输入变量扰动一个单位，并计算皱纹阈值和临界模式数的相对变化（以百分比表示）。当扰动值超过表3中规定的上限或下限时，将其剪切到相应的边界值。结果总结在表5和表6中。

表5. 中年组的敏感性分析：当每个参数增加或减少时，皱纹阈值和皱纹模式数的相对变化。红色表示最具影响力的参数，蓝色表示次要影响的因素。

表6. 老年组的敏感性分析：当每个参数增加或减少时，皱纹阈值和皱纹模式数的相对变化。红色表示最具影响力的参数，蓝色表示次要影响的因素。

中年组
表5表明，纤维预张力参数是调节皱纹阈值的主要因素。特别是，减少纤维预张力会导致皱纹阈值的最大增加。在目前的定义下，减少纤维预张力是提高皱纹阈值和延迟皱纹出现的最有效策略。下一个最有影响力的因素是纤维排列角度，其扰动产生的响应明显大于其他参数，但与纤维预张力相比仍然处于次要地位。

对于模式数，敏感性排名与中年组不同。在中年组中，模式数主要受微观结构方向和厚度比的控制。特别是，减少厚度比是降低模式数的最有效方法，而增加纤维预张力则提供了稳健的次要手段。相比之下，在研究的范围内，纤维角度和纤维排列角度对模式数的贡献相对较小。这种角色分离为中年组提供了实际的设计逻辑：使用纤维预张力来改变皱纹出现的时机（增加预张力），并通过调整纤维排列角度来塑造皱纹的形态。

老年人组
表6中的老年人组敏感性分析显示，纤维预张力在控制皱纹出现方面具有更强的主导作用。具体来说，通过扰动纤维预张力引起的变化远大于其他参数引起的变化，这证实了与预张力相关的效应是随着年龄增长延迟皱纹出现的主要机制。与中年组一样，纤维排列角度是影响模式数的第二大因素，而其他参数对出现阈值的影响较小。因此，控制皱纹出现的策略在两个年龄组中是一致的：首先也是最有效的行动是减少纤维预张力（增加预张力），然后调整纤维排列角度以进一步优化皱纹形态。

老年人组
对于模式数，老年人组的主导因素与中年组有所不同。表6显示，增加厚度比可以最大程度地减少模式数，这使得厚度比成为减少皱纹数量的最有效参数。第二有效的调整是降低刚性比：减少刚性比比同等幅度的纤维角度扰动更能有效降低模式数。相比之下，纤维角度和纤维排列角度在调整模式数方面的贡献相对较小。这种角色分离为中年组提供了实际的设计逻辑：使用纤维预张力来改变皱纹出现的时机（增加预张力），并通过调整纤维排列角度来塑造皱纹的形态。

5.2 来自XGBoost的特征重要性分析
我们还利用XGBoost从全局角度识别出控制输出结果的最重要参数。该分析提供了关于哪些参数在调节皱纹出现和模式数方面最有效的补充数据集级视图。

5.2.1 来自XGBoost的特征重要性分析
从图6可以看出，对于皱纹阈值，两个最重要的因素是纤维预张力参数和纤维角度；对于模式数，主要因素是薄膜厚度比和纤维角度。在研究的范围内，其他参数的贡献较小。

5.2.2 两参数预测等高线图的构建
根据图6中建立的参数排名，图7展示了通过改变最具影响力的参数对并在保持其他参数在其数据集平均值的情况下生成的二维等高线图。这些图是通过在选定的平面上对训练有素的替代模型进行评估来构建的。具体来说，图7(a)展示了皱纹阈值对纤维预张力的响应，而图7(b)展示了模式数对纤维角度的响应。数学上，这些可视化表示了完整五维替代模型的两参数截面（或切片）：

图7. 预测的两参数等高线图。(a) 皱纹阈值作为纤维预张力的函数，(b) 模式数作为纤维角度的函数。其他参数固定在其数据集平均值。点表示中年组，星号表示老年人组。箭头指示增加皱纹阈值和减少模式数的最有效方向。

5.2.3 第一次设计
两参数截面可以用来指导“第一次”设计选择：(i) 为了提高皱纹阈值（延迟皱纹出现），图7(a)建议在等高线水平增加的区域选择纤维预张力和纤维角度，即朝较小的纤维预张力和最大化等高线的角度带移动。(ii) 为了减少临界模式数，图7(b)建议增加厚度比并选择在最小化模式数的等高线谷中调整纤维角度。特别是，如果纤维角度通过设计固定，增加纤维预张力提供了减少模式数的稳健途径；如果纤维角度受到限制，沿着低等高线带调整纤维角度可以提供最有效的减少效果。这些启发式方法为快速、定性地调整皱纹出现和形态提供了指导。我们在等高线图上标记了中年组和老年组的基线参数值。具体来说，点表示中年组，星号表示老年人组。此外，图7中的箭头指示增加皱纹阈值和减少模式数的最有效方向。

在同时调整控制参数以增加皱纹阈值和减少模式数时，存在权衡。例如，减少纤维预张力（即增加纤维预张力）会提高皱纹阈值，但也可能增加皱纹数量。同样，降低刚性比可以减少皱纹数量，但也可能降低皱纹阈值。因此，需要平衡的多参数调整来实现皱纹阈值和模式数的目标值。尽管XGBoost提供了强大且易于解释的基线，但对于当前的正向映射，它的准确性不如BPNN。因此，在逆向设计框架中使用BPNN作为替代模型，而XGBoost主要用于补充性的解释性分析。逆向设计在第6.6节中给出。

5.3 逆向设计
我们使用两种互补的方法进行逆向设计：(i) 直接基于半解析稳定性条件的正则化力学优化，以及(ii) 可扩展的替代辅助遗传算法（BPNN–GA）。第一种方法提供了可解释性，并允许系统地进行参数减少研究，而第二种方法即使在半解析条件不可用且正向模型计算成本较高时也适用。

6.1 来自力学模型的逆向设计
通过解决第2节中的两个方程(12)和(11)，可以获得临界皱纹阈值和临界模式数：

(20)
其中系数取决于以及设计参数向量（见附录B）。给定这些条件，方程(20)通过数值方法求解以获得皱纹阈值和模式数。

6.2 基于正则化的力学优化
逆向问题寻求这样的参数组合，使得皱纹阈值增加且模式数减少，以达到预设的目标。因为在可行的设计空间内可能不存在精确解，我们解决了一个正则化最小二乘问题：

(21)
其中和下标“base”指的是基线参数值，是用户定义的正则化权重。这种分阶段的正则化允许我们逐步“冻结”某些参数（例如较大的参数值），并检查减少的控制参数集是否仍能满足目标。

我们在表7和表8中报告了中年组和老年组的逆向设计结果。总体而言，正则化有效地限制了解决方案接近基线值，同时满足目标变化。请注意，表中的参数值是“设计努力”，表示所需的参数调整百分比形式。

表7. 使用五个和三个参数的正则化力学优化。中年组的逆向设计结果，目标分别为(a) 和 (b) 。蓝色表示五参数逆向设计和三参数逆向设计。

表8. 使用五个和三个参数的正则化力学优化。老年组的逆向设计结果，目标分别为(a) 和 (b) 。蓝色表示典型的五参数逆向设计和三参数逆向设计。逆向设计结果的验证（目标为20%）是通过对完整3D模型进行线性分岔分析来进行的，详见附录C。

6.3 从五参数控制到三参数控制
对于中年组，在目标下，将控制参数从五个减少到三个（例如，将纤维预张力减少到接近基线值）只会使皱纹阈值的相对变化增加约%，而和的额外变化保持在较低水平。因此，在低目标水平下，使用三个参数的逆向设计是一种实用的简化方法。对于目标，虽然减少仍然是可行的，但所需的额外变化较大。对于目标，的额外变化约为，而的变化约为，的变化较小。这些结果表明，对于温和的目标，减少参数控制是合适的；而对于激进的目标，则更倾向于使用完整的五参数控制以避免某些参数的过度偏差。

我们进一步检查是否可以使用两个参数进行逆向设计。我们考虑了三种代表性情况，其中中的变量相对于中的变量被强烈正则化（几乎固定）：

(22)
我们设置或，以便中的参数偏离基线的程度小于，从而隔离了两个参数的影响。结果总结在表9和表10中。基于正则化的力学优化方法，采用两个参数。中年组的逆向设计结果分别为（a）和（b）目标。最后一行是从表7（a）和（b）中提取的三参数逆向设计结果。蓝色表示三参数逆向设计；红色表示不理想的结果，表明与基线参数或目标有较大偏差。表10. 基于正则化的力学优化方法，采用两个参数。老年组的逆向设计结果分别为（a）和（b）目标。最后一行是从表8（a）和（b）中提取的三参数逆向设计结果。蓝色表示三参数逆向设计；红色表示不理想的结果，表明与基线参数或目标有较大偏差。对于中年组，固定（情况3）无法达到目标，表明是不可或缺的。固定（情况1）或固定（情况2）仍然可以达到目标，但需要在剩余的自由参数上进行较大的补偿性调整，这对于实际应用来说是不理想的。对于老年组，固定（情况2）是可行的，对和只有轻微的影响，而固定或（情况1或情况3）则无法达到目标。因此，通常不建议将中年组的控制参数从三个参数减少到两个参数，而对于老年组，使用两个参数的设计可以取得与三个参数设计相当的性能。

关于两组抗皱策略的参数减少洞察总结如下：
1. 中年组：对于低或中等目标，优先选择三参数设计；对于高目标，使用五参数设计以避免参数过度变化。
2. 老年组：对于低或中等目标，优先选择两参数设计；对于高目标，推荐使用五参数设计。

6.2. 使用BPNN的逆向设计
为了识别能够（i）提高抗皱阈值和（ii）减少关键皱纹模式数的设计参数，我们构建了一个逆向设计问题，并使用基于替代模型的遗传算法（BPNN-GA）来解决它。具体来说，训练有素的反向传播神经网络（BPNN）被用作快速的前向映射替代品，而遗传算法（GA）在五维设计空间中进行全局搜索（见图8）。

GA程序实现如下：
- 表示每个个体代表一个候选设计。染色体包含五个基因，对应于设计参数。
- 初始化：初始种群由11,640个个体组成，这些个体直接来自我们数据集中的样本参数集，从而使搜索从设计空间的物理上合理的区域开始。
- 进化：在每一代中，应用交叉、变异和精英选择：
- 交叉：通过配对父母节点并执行单点交叉来生成11,640个后代，交叉点在五个基因位置中均匀选择。
- 变异：以一定的概率应用变异，以生成额外的1,164个个体，以保持多样性并防止过早收敛。
- 选择：将父代种群、交叉后代和变异个体合并成24,444个候选者。适应度最高的11,640个个体被保留到下一代（精英选择）。
- 终止：当最佳适应度连续100代没有改进时，搜索停止。相应的参数向量被报告为逆向设计解。

适应度函数
每个候选者通过BPNN预测与规定目标之间的绝对百分比误差（APE）来进行评估。GA最小化以下公式：
$$APE = \frac{|BPNN_{\text{预测}} - \text{目标}|}{\text{目标}}$$
其中，$BPNN_{\text{预测}}$是第i个替代模型的预测值，$\text{目标}$是相应的目标值。我们采用APE是因为$BPNN_{\text{预测}}$和$\text{目标}$具有不同的数值尺度；这种归一化可以防止一个目标主导搜索，并促进平衡的多目标匹配。

然后我们使用训练有素的替代模型（见第4节）进行特定年龄组的逆向设计。目标是增加抗皱能力同时减少皱纹形成。

使用完整数据集的逆向设计（初次运行）
我们首先使用在整个数据集上训练的替代模型运行BPNN–GA（如第2节所述）。尽管获得了可行的设计，但使用机械模型进行验证后发现了系统性的偏差，特别是在老年组中，无论是还是目标（见表11(b)）。

表11. 使用完整数据集的BPNN-GA逆向设计结果。结果是（a）中年组和（b）老年组。红色表示与目标的较大偏差。

诊断和数据集整理
我们从机械模型中观察到，当非常小时，和之间的关系会出现急剧变化。此外，该区域明显超出了两个年龄组的基线值。因此，为了提高预测准确性，我们移除了这些样本并重新训练了BPNN。这种方法提高了模型性能，如图9和图10所示。

使用整理后的数据集的逆向设计（最终结果）
使用重新训练的替代模型，BPNN–GA一致地满足了多目标。如表12所示，验证结果与期望的和变化非常吻合。

表12. 使用整理后的数据集的BPNN-GA逆向设计结果。结果是（a）中年组和（b）老年组。注意，和的目标都得到了持续满足。

（a）中年组
基线目标 目标
（b）老年组
基线目标 目标

这些逆向设计趋势与敏感性分析一致：增加主要是通过减少（即增加预张力）来实现的，而减少主要是通过增加薄膜厚度和/或调整纤维角度来驱动的。

6.3. 讨论
事实证明，基于力学的优化和BPNN-GA都可以为五个参数提供合理的逆向设计，同时满足和的目标。然而，逆向设计的努力程度不同。基于力学的优化在减少设计空间的维度方面具有明显优势。根据第5节中的敏感性分析，纤维硬度和体积分数对皱纹形成的影响可以忽略不计，这种方法策略性地冻结了这些不太关键的变量。通过仅关注高影响力参数（例如薄膜厚度和纤维角度），它在简化的、物理上重要的子空间内通过适度调整实现了目标模式的减少。相比之下，BPNN-GA在没有这种先验特征选择的情况下探索全局设计空间，以找到更激进的配置。例如，在中年组中，针对目标（见表12(a)），遗传算法通过大幅减少纤维角度并结合预张力调整实现了目标。特定参数策略的差异突显了这些方法的互补性。基于力学的方法通过物理信息驱动的特征选择来最小化活跃设计变量的数量，对于临床或美容干预特别有价值。在这些情况下，保持参数集尽可能接近患者的自然生理基线至关重要，该方法的效率确保只改变必要的机械属性。相反，BPNN-GA方法在探索合成材料设计的更广泛领域时表现出色，在那里约束更宽松，非直观的参数组合——例如老年组观察到的硬度和角度的大幅度变化——可能提供超出传统分析范围的可行工程解决方案。

我们观察到，对于目标为的老年组，设计工作量很大（见表12(b)）。因此，我们采用了修改后的神经网络架构来重新训练模型。具体来说，通过使用ReLU激活函数，我们获得了显著改进的结果（见表13）。

表13. 使用新的激活函数（用ReLU替换Sigmoid）的BPNN-GA逆向设计。对于目标为的老年组，设计工作量得到了改进。这种逆向设计的验证是通过完整3D模型的线性分岔分析来进行的，详见附录C。

老年组目标 改进
基线目标 改进

备注：我们使用具有单一标量适应度函数的GA来搜索设计空间，并识别满足两个目标的候选者。由于两个目标可以在可行范围内同时实现，因此没有观察到明显的权衡，这不需要进行帕累托前沿分析。在获得可行候选者后，我们根据设计工作量作为次要标准选择最终设计。

7. 结论
在这项研究中，我们建立了一个基于力学的逆向设计框架，用于调节纤维增强双层系统的表面皱纹，具体应用于皮肤老化和仿生材料工程。通过将半解析的纤维增强形态弹性板模型与机器学习技术相结合——具体来说是用于敏感性分析的XGBoost和用于优化的反向传播神经网络结合遗传算法（BPNN-GA）——我们成功地映射了五个生理参数（薄膜厚度、基础硬度、纤维方向、预张力和体积分数）与皱纹行为（抗皱阈值和模式数）之间的复杂关系。我们的结果证明了逆向工程材料微结构的可行性，以实现目标延迟皱纹形成（增加）和减少皱纹模式数（减少）对于中年和老年皮肤模型。

这项工作的主要贡献和亮点有三点：
1. 高纤维区域的高效探索：与在高纤维体积分数下计算成本高昂的标准有限元方法不同，我们的半解析板模型能够在整个生理相关纤维密度范围内高效、物理一致地生成数据集。这促进了在纯数值研究中经常缺失的严格全局敏感性分析。
2. 控制机制的解耦：我们明确了微观参数之间的不同作用。全局和局部敏感性分析一致表明，纤维预张力是皱纹阈值的主要调节器，而皱纹形态主要由几何和结构参数（薄膜厚度和纤维方向）控制。这为设计提供了清晰的、可操作的逻辑：首先优先调整预张力以延迟皱纹形成，然后进行几何调整以定制图案。
3. 互补的逆向设计策略：我们展示了两种不同优化途径的实用性。基于力学的正则化优化提供了可解释性和维度减少，表明它可能非常适合于需要保持基线生理状态的保守临床干预。相反，BPNN-GA框架为合成材料设计提供了强大的全局搜索能力，能够识别复杂、非凸景观中的非直观参数组合。尽管在这个框架中取得了进展，但仍有一些局限性指出了未来研究的关键方向：
- 建模限制：当前模型有几个局限性需要注意。首先，将真皮层建模为Winkler基础忽略了真皮-真皮界面处的剪切耦合。其次，采用纯超弹性本构律忽略了皮肤的固有粘弹性。第三，没有明确建模胶原纤维与周围基质的微观力学相互作用，这可能会影响预测的准确性。可以参考[62]来了解如何建模这种相互作用。正如Horgan和Murphy [62]指出的，这种相互作用对于捕捉生物材料的独特Poynting效应（即在扭转下的缩短）至关重要。最后，忽略了模型参数之间的潜在耦合效应，这可能无法捕捉生物材料的实际行为（例如复杂的泊松效应[63]）。这些建模简化可能会影响最终的皱纹阈值和皱纹模式数，这需要在未来的工作中进一步探索。
- 实验挑战：尽管我们的逆向设计算法产生了特定的最优配置，但由于制造限制，实验验证仍然具有挑战性。特别是，制造具有如此高纤维密度的物理模型已被证明是困难的，正如我们使用水凝胶样本进行的初步尝试所证明的。
- 临床转化：理论软组织力学与体内实施之间仍存在差距。尽管该模型为参数的个性化逆向设计提供了定量指导，但未来的严格临床试验是将这些机械见解转化为实际抗衰老治疗的必要条件。