基于对比度增强抑制灰度的改进优化准则法

许小奎 郭宝峰 金 淼

燕山大学先进锻压成形技术与科学教育部重点实验室，秦皇岛，066004

基于对比度增强抑制灰度的改进优化准则法

许小奎 郭宝峰 金 淼

燕山大学先进锻压成形技术与科学教育部重点实验室，秦皇岛，066004

为了抑制变密度法拓扑优化结果中的灰度单元，基于对比度增强策略，提出了两种形式的对比度增强算子，引入到优化准则法迭代计算中，对优化准则法进行了改进。对比度增强算子能够拉大单元灰度差别，驱动中间密度向两极进行分化，最终使优化收敛于0/1材料分布状态。采用典型数值算例对该方法进行了验证，计算结果表明，该方法能够得到拓扑结构清晰的优化结果。

拓扑优化；对比度增强；灰度单元；优化准则法

0 引言

在机械结构设计中，结构拓扑优化为初期结构的选型提供了一种高效的技术手段，能够帮助设计者设计出创新型可靠产品。结构拓扑优化主要研究结构所用材料的最优分布问题，是当前结构优化设计领域中最具挑战性和经济效益的研究方向之一。

在连续体结构的拓扑优化研究中，变密度法[1]自被提出以来，因其求解效率高、易于实施等特点，得到了广泛的应用。基于有限元技术的变密度法在求解结构拓扑优化问题时，会不可避免地出现棋盘格问题[2]，为此，应用变密度法的拓扑优化求解中大多采取了基于图像处理的过滤技术[3-4]。过滤技术解决了棋盘格问题，但却在结构的边缘区域又形成了灰度单元，使得结构边缘扩散而不清晰，影响了结果的后续提取与加工，为直接工程化应用带来了困难。

为了解决灰度单元问题，研究人员提出了许多方法。早期一般采用简单后处理方法，如阈值控制法等，其结果一般不再满足约束条件。密度过滤结合Heaviside函数的方法包括Heaviside过滤[5]、改进Heaviside过滤[6]和体积保持的Heaviside过滤方法[7]，可以得到清晰的优化结果，但其优化需要较多的迭代步数，优化效率较低。基于形态学的方法[6]同样需要较多的迭代步数。BORRVALL等[8]采用增加约束的方法但会使优化求解过程变得复杂。GROENWOLD等[9]对优化准则进行了修改，将单元密度向单极进行压缩，实现了灰度单元的抑制。张志飞等[10]采用过滤与不过滤交替实施的优化方法，但存在棋盘格再次出现的可能。陈垂福等[11]提出了一种变过滤半径的敏度过滤方法，但失去了过滤半径控制最小尺寸的特点。还有一些其他的方法，如非线性扩散技术[12]、双边滤波方法[13]等，对过滤方法进行了修改，但都需要经验来确定合适的调节参数，从而优化得到理想的结果。

本文采用图像处理中的对比度增强技术，实现对灰度单元的抑制。对比度增强技术主要是通过扩大图像灰度差别、增强灰度对比的方法，使图像特征得到加强。本文构造了两种形式的对比度增强算子，将其引入到优化准则法(optimality criteria,OC)迭代计算中，使单元密度向两极分化，以得到边缘清晰的优化结果。该方法不需要添加附加约束条件，仅对优化准则进行了修改，实施过程简便。同时，借助OC法求解，使得优化具有较高的求解效率。

1 连续体结构拓扑优化

1.1 优化模型

基于变密度法建立优化模型，体积约束条件下求结构最小柔度的优化问题可以表示为

(1)

式中，ρ为相对密度(即设计变量)；C为结构柔度；K和U分别为结构总刚度矩阵和总位移矩阵；ke和ue分别为单元刚度矩阵和单元位移矩阵；ρe为单元相对密度；Ve为单元体积；V*为目标体积；n为结构中单元总数。

变密度法中对中间密度的材料属性进行了描述，本文采用改进的SIMP(solid isotropic microstructures with penalization)对材料属性进行描述，单元的弹性模量表示为

(2)

式中，E0为实体单元的弹性模量；Emin为空洞单元的弹性模量，为了防止刚度矩阵奇异，通常设置为大于0的值；p为惩罚系数。

柔度和体积的敏度分别为

(3)

(4)

式中，ke0为单位弹性模量的单元刚度矩阵。

1.2 基于OC法的问题求解

通过OC法进行优化求解时，迭代公式可以表达为

(5)

(6)

式中，m和η分别为移动极限和阻尼系数，是保证迭代稳定和收敛的控制参数，本文中取m=0.2，η=0.5；λ为拉格朗日乘子，在每次迭代中根据体积约束条件采用二分法对其进行求解。

直接采用式(5)进行迭代求解，会出现棋盘格问题，需要进行过滤处理。

1.3 过滤处理

过滤处理方法如图1所示，对于任意单元e，以其中心为圆心，以半径rmin作圆，圆形区域Ωe是单元e的过滤影响范围。敏度的过滤处理表达式为

(7)

Hi=rmin-d(i,e)

(8)

i=1,2,…，Ne

其中，Ne为过滤范围内单元的数量；Hi为距离权重；d(i,e)为过滤范围内单元i和中心单元e的距离。为防止ρe为0时用作分母，ρe取为max(ρe,10-3)。

图1 网格过滤示意图Fig.1 The filtering scheme

将过滤之后的柔度敏度(式(7))替代迭代公式(式(5))中的柔度敏度，可以有效消除棋盘格问题，但由于过滤方法均化作用的影响，又使得结构的边缘区域出现大量灰度单元，边缘扩散而不清晰。图2所示是应用过滤处理的优化结果，结构边缘模糊且呈现为斜坡型边缘，其剖切平面如图3a所示，而理想的结构边缘应是阶跃型边缘(图3b)。结构边缘的剖切平面是沿结构边缘的法向进行剖切所得到的平面。过滤半径rmin越大，优化结果中的灰度单元就会越多，斜坡越缓。

图2 过滤优化结果Fig.2 The result with filtering treatment

(a)斜坡型边缘 (b)阶跃型边缘图3 结构边缘剖切平面图Fig.3 The sectional drawing of structure edge

2 基于对比度增强的改进优化准则法

对比度增强是一种通过改变灰度值来改变图元的对比度，从而改善图像质量的图像处理方法。通常根据需求将原有的灰度值通过函数变换进行重新映射，达到增强图像特征的目的。

下面借助对比度增强方法，实现抑制拓扑优化结果中的灰度单元、增强结构的边缘特征的目的，进而得到理想阶跃型的结构边缘。在实施过程中，构造了对比度增强算子，在过滤处理之后，将其引入到准则迭代计算中，对准则迭代中的相对密度进行修改，表达式为

(9)

由于优化准则法采用二分法求解拉格朗日乘子以使结果满足体积约束，故应用修改之后的迭代公式进行优化求解能够保证体积约束条件得到满足。

对比度增强算子应起到抑制中间密度，使单元密度两极分化，增强边缘的作用。为此，本文提出了两种分段线性变换的对比度增强算子，如图4所示。分段线性变换的优点是可以拉伸感兴趣的灰度区间，抑制不感兴趣的灰度级，使得细节特征得到突出。

(a)算子1 (b)算子2图4 对比度增强算子Fig.4 The contrast enhancement operator

算子1的表达式为

(10)

式中，c1为斜率控制参数，c1≥1。

如图4a所示，算子1将低密度边缘[0,ρa]调整为0、将高密度边缘[ρb, 1]调整为1、将中间密度范围(ρa,ρb)线性拉伸到(0,1)，即算子1通过截取一定范围内的中间密度并拉伸至整个密度范围，以提高单元密度对比度，达到增强边缘的目的。随着c1的增大，对比度增大越明显，当c1→∞时，算子1退化为阈值函数。

算子2的表达式为

(11)

式中，c2为斜率控制参数，0≤c2≤1。

如图4b所示，算子2将[0,0.5)内的密度向下进行压缩，将[0.5,1]内的密度向上进行压缩，增大了两个区间密度之间的差异，因而能够达到增强边缘的目的。随着c2的减小，对比度增大越明显，当c2→0时，算子2退化为阈值函数。

在拓扑优化计算过程中，对比度增强算子的参数大小应当选择合适，以保证迭代过程稳定快速收敛。迭代计算的初期，是主体框架的成形期，大的对比度增强算子会引起对单元密度改变过多从而导致删除单元过多，进而影响随后的最优结果的搜索，因此，应该采用较小的对比度增强效果或不采用增强效果。在迭代计算的后期，拓扑结果演变缓慢且变化幅度较小，此时应采用较大的对比度增强效果，使迭代快速收敛。

根据上述思路，设计对比度增强算子的参数取值。对于大多数结构的优化计算，当设计变量的最大变化量小于0.07时，形成了拓扑结构的主体框架，此时的迭代步数记为k(c<0.07)。在拓扑结构的主体框架形成过程中，不对单元密度进行对比度增强处理。当主体框架形成之后，采用逐步增强的对比度增强算子进行优化求解。

在实施过程中，设定c1的变化表达式为

(12)

c2的变化表达式为

(13)

式中，k为迭代次数；Δ1和Δ2为变化幅度。

c1和c2的变化曲线如图5所示。在迭代步数达到k(c<0.07)以前，c1和c2都保持不变；在迭代步数k(c<0.07)以后，c1逐渐增大，c2逐渐减小，Δ1和Δ2分别控制c1和c2变化的幅度。在本文优化算例中取Δ1=0.03，Δ2=0.02。

(a)c1

(b)c2图5 c1和c2的变化曲线Fig.5 The varying curves of c1 and c2

为了保证迭代过程的稳定，还应该保持每次迭代中单元密度的变化在移动极限内。

本文采用的两种对比度增强算子驱动中间密度向两极进行分化，与文献[5]中的向单极进行压缩相比，更有利于优化得到最优拓扑结构。

3 算例与分析

采用下面典型的连续体结构拓扑优化数值算例对本文方法进行验证。算例中的结构均采用四边形单元进行离散，设材料的弹性模量E0=1，Emin=10-9，泊松比为0.3，惩罚因子p为3。优化迭代终止的判定条件为max|ρ(k+1)-ρ(k)|≤0.01。

为了衡量优化结果中灰度单元的比例，采用灰度单元占比δ进行描述，表达式为

(14)

δ的值越大，表示存在的灰度单元越多。当结果中的所有单元都是灰度单元时，δ为1；当结果中的单元密度非0即1时，δ为0。

算例1 图6所示为一个悬臂梁，设计域大小为100 mm×50 mm，厚度为1 mm。左侧固定，右侧中间作用竖直向下的载荷F，大小为1 N。设计域离散为100×50个单元，设计域体积的50%作为目标体积，过滤半径取为3.5 mm。算例1的优化结果构型如图7所示，图7a为应用OC方法优化的结果，图7b、图7c分别为引入对比度增强算子1、2的优化结果，图7d为采用文献[5]中方法的优化结果。图8所示为优化结构边缘的剖切平面密度变化图，剖切线的位置如图7a和图7b中黑线所示。各优化结果的柔度、灰度占比和迭代次数如表1所示。

图6 算例1设计域及边界条件Fig.6 Design domain and boundary conditions of example 1

(a)OC方法 (b)算子1

(c)算子2 (d)文献[5]图7 算例1的拓扑优化结果Fig.7 The topology results for example 1

图8 优化结构边缘剖切平面图Fig.8 The sectional drawing of structure edge for topology results

柔度C(N·mm)灰度占比δ(%)迭代次数OC方法68.929.873算子162.2067算子262.4073文献[5]62.34.0480

由图7可以看出，应用对比度增强算子的两种方法抑制了灰度单元的出现，得到了拓扑结构清晰的优化结果。由图8可以看出，与OC方法相比，应用对比度增强算子之后，得到了阶跃型的结构边缘。由表1可知，本文方法与OC方法相比，优化结果的柔度值更小，优化效果得到了提升；与文献[5]中的灰度单元抑制方法相比，所得结果的灰度单元占比更小，优化结果更接近于0/1分布状态，并且在优化求解效率方面占有明显优势。

算例2 图9所示为一个简支梁，尺寸为240 mm×60 mm，厚度为1 mm。上端面中间作用竖直向下的载荷F，大小为1 N。设计域体积的50%作为目标体积。鉴于是对称结构，选用右侧1/2模型进行优化计算，模型离散的数目为120×60。优化计算中采用了两种过滤半径进行过滤处理。图10为优化结果的拓扑构型，其中图10a、图10c、图10e采用的过滤半径是2.5 mm，分别为OC方法、引入算子1和算子2的优化结果；图10b、图10d、图10f采用的过滤半径是4.5 mm，分别为OC方法、引入算子1和算子2的优化结果。各优化结果的柔度、灰度占比和迭代次数如表2所示。

图9 算例2设计域及边界条件Fig.9 Design domain and boundary conditions of example 2

(a)OC方法(rmin=2.5 mm)(b)OC方法(rmin=4.5 mm)

(c)算子1(rmin=2.5 mm)(d)算子1(rmin=4.5 mm)

(e)算子2(rmin=2.5 mm)(f)算子2(rmin=4.5 mm) 图10 算例2的拓扑优化结果Fig.10 The topology results for example 2

柔度C(N·mm)灰度占比δ(%)迭代次数图10a81.922.5264图10b86.429.861图10c77.20.087图10d78.00.068图10e77.40.082图10f78.20.078

由图10可以看出，基于不同的过滤半径，应用本文方法，都得到了结构边缘清晰的优化结果，所得结果的拓扑形式与采用OC方法的一致，过滤半径越大，所得结果中最小结构尺寸越大。该方法保持了过滤半径控制优化结果中结构最小尺寸的优点。从表2中可以看出，在灰度方面，增大过滤半径，会使OC法结果中的灰度增大，对本文方法没有影响；在柔度方面，增大过滤半径，会使优化结果的柔度值增大。

上述算例的结果显示，所提出的两种对比度增强算子得到的拓扑形态结果基本一致。这主要是因为对比度增强算子在拓扑结构主体框架形成之后才开始发挥作用，对比度增强算子对优化结果的拓扑形式几乎不产生影响，主要作用是消除灰度单元、增强结构边缘。算子1是逐步将靠近0-1两极的密度进行增强，算子2是将中间密度向0-1两极进行压缩。在优化进程中，算子1稳步减少中间密度，在拓扑结构复杂的情况下往往能表现出更好的稳定性。

4 结论

本文提出了两种形式的对比度增强算子，对优化准则法进行了修改，驱使中间密度向两极分化，以抑制变密度法拓扑优化结果中灰度单元的出现。典型数值算例的计算结果显示，该方法消除了过滤处理导致的边缘扩散，得到了拓扑结构清晰的优化结果。

[1] BENDS∅E M P. Optimal Shape Design as a Material Distribution Problem[J]. Structural Optimization, 1989, 1(4): 193-202.

[2] SIGMUND O, PETERSSON J. Numerical Instabilities in Topology Optimization: a Survey on Procedures Dealing with Checkerboards, Mesh-dependencies and Local Minima[J]. Structural Optimization, 1998, 16(1): 68-75.

[3] SIGMUND O. Design of Material Structures Using Topology Optimization[D]. Lyngby: Technical University of Denmark, 1994.

[4] BRUNS T E, TORTORELLI D A. Topology Optimization of Non-linear Elastic Structures and Compliant Mechanisms[J]. Computer Methods in Applied Mechanics and Engineering, 2001, 190(26/27): 3443-3459.

[5] GUEST J K, PRÉVOST J H, BELYTSCHKO T. Achieving Minimum Length Scale in Topology Optimization Using Nodal Design Variables and Projection Functions[J]. International Journal for Numerical Methods in Engineering, 2004, 61(2):238-254.

[6] SIGMUND O. Morphology-based Black and White Filters for Topology Optimization[J]. Structural and Multidisciplinary Optimization, 2007, 33(4): 401-424.

[7] XU S L, CAI Y W, CHENG G D. Volume Preserving Nonlinear Density Filter Based on Heaviside Functions[J]. Structural and Multidisciplinary Optimization, 2010, 41(4): 495-505.

[8] BORRVALL T, PERTERSSON J. Topology Optimization Using Regularized Intermediate Density Control[J]. Computer Methods in Applied Mechanics and Engineering, 2001, 190(37/38): 4911-4928.

[9] GROENWOLD A A, ETMAN L F P. A Simple Heuristic for Gray-scale Suppression in Optimality Criterion-based Topology Optimization[J]. Structural and Multidisciplinary Optimization, 2009, 39(2): 217-225.

[10] 张志飞, 徐伟, 徐中明, 等. 抑制拓扑优化中灰度单元的双重SIMP方法[J]. 农业机械学报, 2015, 46(11): 405-410.

ZHANG Zhifei, XU Wei, XU Zhongming, et al. Double-SIMP Method for Gray-scale Elements Suppression in Topology Optimization[J]. Transactions of The Chinese Society of Agricultural Machinery, 2015, 46(11): 405-410.

[11] 陈垂福, 杨晓翔. 一种考虑密度补偿的变过滤半径敏度过滤方法[J]. 中国机械工程, 2017, 28(6): 669-675.

CHEN Chuifu, YANG Xiaoxiang. Variable Filter Radius Sensitivity Filtering Method Considering Density Compensation[J]. China Mechanical Engineering, 2017, 28(6): 669-675.

[12] WANG M Y, ZHOU S, DING H. Nonlinear Diffusions in Topology Optimization[J]. Structural and Multidisciplinary Optimization, 2004, 28(4): 262-276.

[13] WANG M Y, WANG S. Bilateral Filtering for Structural Topology Optimization[J]. International Journal for Numerical Methods in Engineering, 2005, 63(13): 1911-1938.

AModifiedOptimalityCriterionMethodBasedonContrastEnhancementforGrayScaleSuppression

XU Xiaokui GUO Baofeng JIN Miao

Key Laboratory of Advanced Forging & Stamping Technology and Science,Yanshan University, Qinhuangdao,Hebei,066004

In order to suppress gray scale elements in density-based topology optimization, two kinds of contrast enhancement operators were proposed based on contrast enhancement strategy and introduced into the iterative calculation to modify the optimality criteria method. The contrast enhancement operator might enlarge the gray scale differences of the elements, drive the intermediate densities to the two ends of the range, and finally make the results convergence to the 0/1 material distribution. The effects of the proposed method were investigated with classical numerical examples. The results show that the optimization results with crisp boundaries are obtained by applying the proposed method.

topology optimization; contrast enhancement; gray element; optimality criteria method

2017-08-21

国家自然科学基金资助项目(51575474)；河北省自然科学基金资助项目(E2015203220)

TH122

10.3969/j.issn.1004-132X.2017.24.008

(编辑王艳丽)

许小奎，男，1990年生。燕山大学机械工程学院博士研究生。主要研究方向为成型设备结构分析与优化设计。郭宝峰(通信作者)，男，1958年生。燕山大学机械工程学院教授、博士研究生导师。E-mail: guobaofengysu@126.com。金淼，男，1968年生。燕山大学机械工程学院教授、博士研究生导师。