2.1　参数化水平集方法

2.1.1　基于水平集的结构边界隐式表达

基于水平集的结构边界隐式表达的思想是将结构边界隐式地嵌入到高一维的水平集函数中，并将结构边界视为水平集函数的零等值面，如图2-1所示。值得注意的是，这里的水平集函数应具有Lipschitz连续性[5]，[6]。假定有一个固定的欧拉空间作为参考设计域D，该设计域中包含了实体、边界和孔洞，设计域D中的各个部分可按照如下方式进行数学表达：

其中，x代表设计域D中的空间变量，即水平集网格的坐标；t为一种虚拟的时间变量；Ω表示设计域中结构的所有可能形状；Γ代表结构边界。需要注意的是Γ包含了Dirichlet边界Γ_D、Neumann边界Γ_N和无牵引力边界Γ_F：

除了将结构边界表述为零水平集外，基于水平集的结构边界隐式描述的另外一个重要问题是将结构边界几何信息映射到力学模型中。一般地，结构边界的数学表达式（2-1）可重新表达为一种包含Heaviside函数H(Φ)的格式。例如，可利用H(Φ)在结构设计域内构建泛函J的积分形式：

其中

式（2-3）中，dV代表体积积分。

在进行结构拓扑优化设计时，式（2-3）可以通过有限元方法在固定网格上进行近似计算。然而在水平集方法中，若不重新划分网格，标准的有限元方法很难准确地计算一些结构边界上的单元应变，特别是当固定的有限单元被动态的水平集边界所分割时，固定有限元网格所带来的计算误差会更大[7]，[8]。学者们通常会采用一种所谓的人工材料模型[8]，[9]来应对上述问题。该模型假设有限单元的应变或刚度等指标与该单元的材料百分比（即虚拟单元密度）成正比，且孔洞单元由一种密度极低的弱材料填充以避免刚度矩阵奇异。在该假设的基础上，泛函J的体积积分可由有限元法近似计算得到：

其中，h表示设计域D中的任一单元；D_h表示该单元的设计空间；NE代表有限单元的数量。由于精确的Heaviside函数式（2-4）不可微，在拓扑优化设计过程中需要引入一种近似的Heaviside函数（光滑且可微）[9]来替代精确的Heaviside函数：

近似Heaviside函数的偏导数，即Dirac delta函数可定义为

在式（2-6）和式（2-7）中，γ是一个足够小的正数，其大小通常等于2Δ，Δ代表固定水平集网格的边长。于是，前面提到的虚拟单元密度可被定义为

其中，τ=0.001被用于避免计算奇异问题。式（2-8）中，位于空间D_h上的体积积分可通过高斯积分法进行计算。值得注意的是，对于结构边界上的单元，需要布置大量的高斯积分点来保证计算精度[9]。类似地，下面的体积积分亦可由有限元方法近似计算得到：

通过引入时间变量t，水平集函数随时间沿特定方向动态变化从而实现结构的拓扑优化。利用链式求导法则对水平集函数Φ(x，t)=0求其关于时间变量t的偏导[10]，可得Hamilton-Jacobi偏微分方程，如下所示：

其中，v=dx/dt为结构边界上的速度场。v由法向速度场和切向速度场构成，由于仅法向速度对结构边界的形状演化起作用[11]，[12]，故式（2-10）可改写为

其中法向速度为

在水平集方法的框架下，水平集函数并没有一种显式的解析形式，因此水平集在优化过程中的运动必须通过求解Hamilton-Jacobi偏微分方程来实现。结构边界的优化过程可等价于在速度场驱动下的水平集面的演化过程，而该速度场由给定的目标函数和约束条件所决定。通过Hamilton-Jacobi偏微分方程求解合适的速度场十分困难，通常需用规则的网格将水平集函数离散化，并采用差分方法求解，例如基于Eulerian方法的显式有限差分[11]，[13]。运用有限差分方法时需谨慎处理其数值计算问题，例如速度场的扩展（velocity extension）、重新初始化（re-initialization）以及CFL条件（CFLcondition）等[14]。这些复杂的数值计算过程导致水平集方法的效率和稳定性下降，阻碍其在结构优化领域的发展和应用。为克服上述数值计算困难，并提高求解效率，下文将在传统水平集方法的基础上，讨论一种基于CSRBF和离散小波变换（discrete wavelet transform，DWT）的参数化水平集方法。

2.1.2　基于径向基函数插值的水平集函数参数化

径向基函数（radial basis function，RBF）的特点是能够利用一元函数来描述多元函数，其最为常见的应用场景是用于近似或者拟合大规模散乱数据的几何形状。径向基函数运算简单、存储方便及计算效率高，在学术界和工程界得到了非常广泛的关注和应用。理论上，当非线性等式在插值节点处构成的插值矩阵可逆时，RBF可确保近似曲面或超曲面的光滑性。因此与水平集方法结合，RBF能够快捷地逼近任何可能的形状与拓扑，同时保留水平集方法隐式边界描述的优点[13]，[15]。RBF方法的一些独特优点，例如插值系统解的唯一性、插值计算的效率、RBF的光滑性和收敛性等，使该方法在拓扑优化领域的应用很有吸引力[16]。

利用具有N个插值点的一组RBF，任意函数f(x)：Rd→R(d≥1)的近似形式可定义为

其中，x_i为第i个插值控制点的位置坐标；α为插值扩展系数；φ为径向基函数；‖·‖代表d维空间上的欧氏距离。

RBF插值的精度与效率很大程度上由所选取的插值核函数决定，按照插值过程中形成矩阵的稀疏性来区分，可将RBF分为全局支撑的径向基函数（GSRBF）和紧支撑径向基函数（CSRBF）两类[17]。GSRBF具有精度高、光滑性好的特点，但其在近似水平集面的过程中会产生元素全部非零的插值矩阵，导致计算复杂度达到O（N2）甚至O（N3）[18]，难以适应大规模结构优化问题。此外，GSRBF插值精度很大程度上取决于所谓的自由形状参数（free shape parameter），目前并没有可靠的方法来确定不同拓扑优化问题中自由形状参数的取值。相比之下，CSRBF具有严格正定、插值系数矩阵稀疏、参数对全局计算精度影响小，以及插值操作能够继承CSRBF的连续性等优点，使其更加适合结构拓扑优化问题。虽然CSRBF插值精度稍低于GSRBF，但针对常见的拓扑优化设计问题，两类基函数均能够满足水平集方法的精度要求。

本书将采用CSRBF对水平集函数进行插值，以此将水平集方法转化为一种广义的参数化形式。Wendland[19]所提出的系列CSRBF可表示为

对于二维问题，r为定义在欧式空间内的支撑半径（radiusof support）：

其中，d_I表示在CSRBF的支撑半径范围内，当前样本点（x，y）到插值控制节点（x_i，y_i）间的距离。d_mI反映CSRBF在节点（x_i，y_i）处的影响范围大小。要保证插值的稳定性和效率，需要选择大小合适的支撑半径，过小的支撑半径会导致插值奇异，过大的支撑半径显著增加计算成本。在本书中采用文献[20，21]给出的支撑半径：令参数d_mI=d_max×C_I，其中尺度因子d_max取值为2.0~4.0，C_I决定当前插值控制点的邻域范围，以保证在邻域内搜索到足够多的控制点。

图2-2对比了具有C2阶、C4阶和C6阶连续性CSRBF的形状。可见，不同阶的CSRBF均具有良好的光滑性，表明若空间内插值点布置合理，3种CSRBF都能够很好地对水平集函数进行逼近，并保证插值所得函数的光滑性和完备性。

利用CSRBF插值，水平集函数可被近似表达为一组CSRBF与其扩展系数的乘积形式：

其中，α（t）=[α₁（t），α₂（t），…，α_N（t）]T为CSRBF插值的扩展系数，φ（x）=[φ₁（x），φ₂(x)，…，φ_N(x)]包含了全部的CSRBF。可见，CSRBF仅与空间变量相关，扩展系数仅与时间变量相关，水平集函数则实现了时空变量解耦。Gaussian RBF通过考虑设计域内的全部插值点信息来获取当前样本点的水平集函数值，因此具有极高的插值精度，但这也势必将带来高昂的计算成本。为方便论述，现将式（2-18）改写为如下矩阵运算形式：

其中，A为可逆矩阵。不难发现，A的稀疏性对插值效率有重要影响。

将式（2-18）代入到经典的Hamilton-Jacobi偏微分方程中[7]，可将该偏微分方程转化为一系列易解的常微分方程：

其中，α₀代表在时间t=0时的扩展系数向量，可以在初始迭代阶段通过求解式（2-19）计算获得。于是，水平集的速度场V_n可进一步表示为

至此，标准的水平集方法被转化为一种参数化形式，能够在保持传统水平集方法有点的基础上，克服其应用缺陷。通过式（2-22）可知，速度场V_n在设计域内的全部节点上进行计算，说明水平集的速度场被自然地扩展到了整个设计域上，从而无需额外引入所谓的速度扩展策略[22]。由于Hamilton-Jacobi偏微分方程被转化为常微分方程，水平集函数的更新也不再受限于CFL条件[23]，优化效率将被显著提升。此外，此时的设计变量不再是水平集函数，而是CSRBF插值的扩展系数，表明结构拓扑优化设计问题被简化为一种广义的“尺寸”优化问题，而该问题能够被高效的梯度型算法[24]，[25]求解，这将进一步提升水平集方法的适用性。

2.1.3　基于离散小波分解的插值矩阵压缩技术

面对复杂的大规模优化问题，基于CSRBF的参数化水平集法的计算效率虽然由于插值矩阵A的稀疏性得到了一定提高，但其具备进一步提升的空间[16]。本章引入离散小波变换（discrete wavelet transform，DWT）技术对由CSRBF构建的插值矩阵A进行再压缩，从而获得一种极其稀疏的插值系统，减少优化过程的处理时间及其所需的计算机存储空间。

DWT是一种输入信号或信息的多分辨率分解技术，近几十年来被广泛用于信号处理（signal processing）、图像压缩（image gompression）、去噪（denoising）等方面[26]。因其具备利用极少元素捕获大规模信息集合中关键信息的能力，使DWT的应用范围得以扩展到大规模线性系统或方程组中矩阵和向量的压缩[27]，[28]。若将式（2-19）给出的水平集函数插值操作视为规模与插值控制节点数量相关的线性系统，那么基于DWT的矩阵压缩技术则可作为一种“黑箱”操作嵌入到参数化水平集模型中，这意味着能在不改变参数化水平集模型的前提下，通过增加几次代价较小的小波变换，即可起到减少插值计算成本的效果。简言之，DWT能够在几乎不损失计算精度的前提下将稠密插值矩阵压缩为一种极其稀疏的矩阵。

针对式（2-19）中的水平集函数插值矩阵A，DWT首先将其转化为相同规模的小波基矩阵 。依据该小波基矩阵 的数据分布，能够分辨出其中的关键元素和噪声元素。因此，接下来可以通过阈值过滤技术，清除掉 中相当部分的无用元素，从而构建一个极其稀疏的新矩阵 _s。最终，水平集函数能够通过该稀疏矩阵 _s计算得到。

其中，α(0)表示原始向量α中的任一元素。 和 均为小波变换后的新向量 中的元素，其中k=1，2，…，N/2。需要注意的是，为保证一级小波变换具备完整的小波阶，通常需要向量的长度N为偶数。若N为奇数，可采用一种简易的处理方式，即在小波变化过程中在原始向量α中额外增加一项“0”元素，在重构原始向量α时再去除掉额外增加的一项“0”元素。[h₁，h₂]=[ /2， /2]为高通滤波器，[g₁，g₂]=[- /2， /2]为低通滤波器[29]。二维插值矩阵A的分解过程本质上为多次一维DWT过程的叠加。典型的二维分解过程如图2-3所示[26-30]。为方便论述，引入N×N卷积矩阵W[31]描述该塔式算法：

由式（2-23）可以看出，W为正交矩阵。利用该W矩阵，经小波变换后的水平集函数 、扩展系数向量 和插值矩阵 可分别表达为

在式（2-19）两端分别左乘W可得

其中，WT与W为单位矩阵。

通过比较式（2-24）~式（2-26）可以得到新的插值系统：

接下来，利用下述阈值过滤算法来去除小波基矩阵 中的无用元素，从而构建极其稀疏的插值矩阵 _s：

其中，q为矩阵 中的任意元素； 为矩阵 中全部元素的平均绝对值；κ为阈值调整参数。经过阈值过滤后，可以将式（2-27）的插值过程重建为一种极其稀疏的插值系统：

下面简述引入DWT后的优化过程。特别地，在第一次迭代中，将t=0时的水平集函数记为Φ₀，其小波变换后的形式记为 ，式（2-21）中定义的初始扩展系数向量α₀可由逆过程 计算得出。而α₀的小波变换形式 ₀则可由 计算得到。在后续迭代中，首先可利用基于梯度的优化算法更新设计变量α。接着可通过式（2-24）计算获得 ，并通过式（2-29）更新水平集函数的小波基形式 。最后，通过重构过程 得到水平集函数的原始基形式。

如上所述，本章将DWT与CSRBF结合起来构建了一种新的参数化水平集方法以提高求解大规模优化问题时的计算效率。与基于CSRBF的参数化水平集方法相比，上述方法在插值矩阵中含有更多零元素，求解更快速。值得指出的是，在上述方法中稀疏的插值矩阵 _s仅需在整个优化过程中计算一次。可见，上述参数化水平集方法仅仅是在插值过程中增加了一次小波变换操作 =W α和一次重构操作 ，而这两个操作均利用极其稀疏的矩阵W，因此该过程所带来的额外计算成本可以忽略不计。