基于自适应阈值的拓扑图像密度单元优化

王 冲， 茅 健， 郑 武

（1 上海工程技术大学 机械与汽车工程学院， 上海 201620； 2 华融普瑞北京科技有限公司， 北京 100089）

0 引 言

拓扑优化设计是在满足一定约束条件下寻求给定材料的指定性能最优的材料分布，与其他尺寸、形状等常规优化方法相比，设计应用范围更大且可以最大化结构性能［1］。 随着数学规划、计算机科学等学科发展，其已成为结构设计领域重要方法。 然而拓扑优化所得结果为非参数化模型，在拓扑边缘存在大量灰度单元、导致边缘不清晰、可制造性较差等问题。

针对无实际制造意义的灰度单元问题，基于图像处理技术，Wang［2］等人提出一种非线性双边过滤器，可以将中间密度单元有效逼近0-1 边界；Sigmund［3］提出一种基于形态学的图像处理算子，通过算子“膨胀”和“收缩”有效提升了优化结构的清晰度；Guest［4］利用线性投影Heaviside 阶跃函数实现中间密度单元的0／1 逼近；Fuchs［5］通过引入倒变量及添加SRV（Sum of Reciprocal Variables）约束条件减少了结构中灰度单元；Zhao［6］提出一种平滑类Heaviside 函数，帮助SIMP 模型获得锐化的0／1 分布；Lee［7］通过在SIMP 中加入可移动和正则化的Heaviside 函数，可以得到收敛解且提高了材料之间边界清晰度；Liang［8］通过序列整数规划和正则松弛算法实现了对具有离散变量灵敏度的拓扑优化问题的求解，黑白度得到提升。

在数字图像领域，图像的边缘锐化及降噪处理可以有效退化部分图像特征，实现可利用部分图像显著增强，提升整体结构清晰度。 赵艳妮［9］等基于模糊理论，通过改进变分因子计算得到边缘熵和背景熵，实现图像的局部自适应增强；娄联堂［10］等通过图像灰度变换和最大类间方差的联系对比，提出一种OTSU 阈值优化结合分段灰度变换方法，实现图像分割；胡涛［11］等基于L2范数设计正则项模型，结合非局部均值算法实现图像灰度部分去噪；牛为华［12］等基于RiemannLiouville 分数阶微分，构造8方向微分卷积模板，实现多项图像质量指标的提升；刘巧红［13］等通过引入各向异性扩散张量，结合二阶总广义变分（Total Generalized Variation， TGV）正则项平衡，实现了图像边缘结构的降噪。

本文在研究基于变密度法SIMP 模型拓扑结构基础上，针对其存在的灰度单元问题，基于数字图像处理的相关研究思路，提出一种自适应阈值空域灰度变换算法，通过提取拓扑结构灰度单元，依据histogram 函数确定均匀宽度值，根据整体结构灰度值约束条件建立自适应阈值计算公式，为方便计算，将连续可积黎曼积分转换为离散单元面积计算，求得自适应整体阈值算子，在满足体积约束条件下将图像增强算法应用于拓扑图像。 与其他数字图像处理OTSU 最大类间方差法、动态范围压缩、阶梯量化等算法对比，自适应阈值方法整体效果更为优秀，同时可以满足拓扑结构体积约束情况下实现灰度单元有效过滤，偏离值仅为0.01%，收敛快速且优化结构清晰。

1 连续体拓扑图像生成方法

变密度法舍弃了均匀化法中微结构单胞的概念，通过引入密度可变材料，单元密度为ρ，ρ ＝0 代表单元为空材料，ρ ＝1 代表单元存在实体材料。 并设定该单元的材料属性（如弹性模量） 与单元密度存在对应关系。

基于变密度法单元的设计变量由均匀化法的3个设计变量即洞孔的长度、宽度以及角度变更为1个设计变量即单元密度，优化效率得到提高。 然而优化所得的中间密度无法应用于生产实际的工程制造，为此BendsØe［14］和Zhou［15］提出了固体各向同性材料惩罚模型，即SIMP，将单元密度变量的取值范围由集合｛0，1｝放松到连续区间［0，1］，引入惩罚因子p，以幂指数惩罚的方式使得优化结果向两端逼近，最终使结果收敛于0／1 材料分布，有效减少了中间密度。 此时材料杨氏模量与单元密度之间的关系模型如式（1）所示：

其中，Ei和ρi分别是第i个单元的杨氏模量和密度；E0是材料的杨氏模量；N为优化单元的数量；ρi的取值范围为｛0，1｝ 。

基于上述给定限制体积约束下的SIMP 材料插值模型，整体结构柔顺度最小的优化问题可表示为式（2）：

需要满足的约束条件为式（3）：

其中，xi为材料的单元相对密度；ui为材料单元列向量；k0为材料密度为“1”时的单元刚度矩阵；xmin为最小材料单元密度，目的是避免有限元计算时出现刚度矩阵奇异；vi为材料单元的相对体积。

2 数字图像处理方法

2.1 最大类间方差法（OTSU）

最大类间方差法，即使用最大类间方差自动确定全局阈值，以方便后续进行图像的分割处理［16］。关键思想在于利用数字图像的灰度阈值将图像的目标和背景分离开。 在目标与背景像素值差异较大时，类间方差也越大，分割阈值的选取更为精确，图像分割效率越高。

OTSU 最大类间方差法原理：令［0，255］为数字图像灰度级取值范围，设ni为图像中灰度为i的像素的数量，则图像中像素为i的概率为：Pi ＝ni／（n0＋n1＋…＋n255）。 设阈值t将图像的像素分为两簇，A类簇取值范围为［0，t］，B类簇取值范围为［t ＋1，255］。

设阈值为t时，像素点在A类簇的概率为PA（k）、分配到A簇时像素平均灰度值为mA（k）。同理，像素点在B簇的概率为PB（k）、 分配到B簇时像素的平均灰度值为mB（k），对于灰度值k，其累加均值记为m，数字图像的灰度为mG，则有式（4）和式（5）：

通过求解方差，可得式（6）：

满足式（7）：

通过等效代入可得式（8）：

所得的最佳阈值为使得等式（8）取得最大值的灰度级阈值t。

可选目标区域在整体图像的占比达30%以上时，基于最大类间方差法的分割效果比较理想，而在整体区域占比过小时，优化效果迅速降低。 对于拓扑图像的优化，OTSU 最大类间方差法可以实现灰度单元的过滤处理，提升整结构0，1 数据集占比，但是对于拓扑约束无法满足，即产生一定偏差值，需要后续修正。

2.2 动态范围压缩算法

动态范围压缩（DRC）主要用来实现高动态范围图像HDR（High Dynamic Range）到低动态范围图像LDR（Low Dynamic Range）的映射，由于自然界中真实场景的亮度颜色区间非常广泛，而人工硬件设备无法显示全部范围，因此基于DRC 算法进行压缩实现图像显示。

常见的DRC 算法在图像单元的处理方式上可以分为对全局单元映射以及局部灰度映射两大类。灰度全局映射不考虑局部细节，统一降低图像分辨率，实现对比度的降低［17］。 而局部映射复杂度较高，优势是可以保留图像的局部细节。 可实现动态范围压缩的算法包括：线性移位算法、对数映射以及分段函数映射算法等等。

线性移位算法原理是将HDR 高比特图像直接右移固定位数变为可显示LDR 低比特图像，容易降低数值比较集中的数据段，导致失真情况。 常用对数映射即以2 为底，通过将高动态范围对数化为低动态范围，然后线性放缩至可选区间，实现图像的压缩显示。 优化后图像单元取值g（x），式（9）：

其中，f（x） 为优化前图像单元值，系数λ取值需保证图像g（x） 值符合区间范围。

2.3 自适应阈值算法

传统阈值分割算法主要在于提取待处理图像中目标区域，将其与背景区域区分，其阈值的确定有多种方式，如基于最大熵值，最大类间方差值或者迭代式阈值等。 此类阈值确定方法可以实现一般图像的处理，由于拓扑结构图像的产生需要满足给定的体积分数及其他约束条件，即图像本身的灰度值需要满足限制条件，传统阈值分割算法在处理拓扑图像时会产生偏离值，导致最终选择的目标结构不符合设计条件。 本文针对拓扑结构的限制条件提出一种自适应阈值灰度变换算法，在满足整体图像灰度值不变下建立阈值计算公式，整体流程如下：

（1）提取拓扑结构单元并数值化处理，得到数据集Pm*n，其中m为横向单元数，n为纵向单元数，数据集P（i，j） 为单元值，取值区间为［0，1］；

（2）基于OTSU 最大类间方差法确定初始阈值，减少后续计算次数；

（3）区间离散化灰度分布函数，根据histogram函数确定均匀宽度值，实现整体数据集的最优划分；

（4）依据整体灰度值约束建立自适应阈值计算公式，求解得出阈值算子t；

（5）通过全局阈值算子t对灰度图像数据作处理，即在灰度级数X高于算子t时置为目标区间上限b，在灰度级数X低于算子t时置为目标区间下限a。根据拓扑结构实际灰度范围选择a ＝0，b ＝1。

对于连续型灰度直方图分布函数，满足灰度约束的计算公式（10）：

其中：x为灰度图像单位像素所在的取值范围；t值为自适应阈值取值；f为拓扑优化体积约束条件。

由于实际图像存在最小化像素点，因此根据数据图像实际单元灰度像素阶跃性，对灰度直方图分布函数做离散化处理，经过离散化处理后，满足灰度约束的计算公式（11）：

其中：r为离散灰度函数满足约束最近点；bin值为选取的离散集合均匀宽度，依据histogram 函数取得近似值；t值为满足实际约束的自适应阈值取值。

3 仿真实验

3.1 二维悬臂梁优化

采用单侧受不变载荷悬臂梁作为数值计算模型，其结构如图1 所示。 不失一般性，其左端固定，整体梁身为优化设计区域，设计域离散单元数为240*80，右侧中心部分承受竖直向下力F。

图1 悬臂梁结构示意图Fig.1 Schematic diagram of cantilever beam structure

取体积约束f ＝0.5，惩罚系数p ＝3，过滤半径rmin＝2.5，约束条件为体积原结构50%，优化目标为满足约束条件下最小化柔度。 悬臂梁优化结果对比（V ＝0.5）如图2 所示。

图2 悬臂梁优化结果对比（V＝0.5）Fig.2 Comparison of optimization results of cantilever beam （V＝0.5）

经过迭代得到的拓扑结构如图2（a）所示；对拓扑结构采用OTSU 最大类间方差法优化，所得结果如图2（b）所示；采用动态范围压缩方法优化所得结果如图2（c）所示；采用阶梯量化方法所得结果如图2（d）所示；采用普通阈值分切方法所得结果如图2（e）所示；采用本文所述自适应阈值算法所得优化对比实验的结果如图2（f）所示。

3.2 实验结果比较与分析

对基于SIMP 悬臂梁模型拓扑图采用不同方法优化数据见表1，统计数据包括不同数值优化方法的总灰度集，优化后数据的偏离值以及整体数值偏离程度。 本文设定拓扑结构离散单元数为240*80，优化选用的体积约束f为0.5，初始图像灰度单元总数为9 600。

表1 悬臂梁不同方法优化对比Tab.1 Optimization comparison of cantilever beam with different methods

分析表1 数据可知，OTSU 最大类间方差法由于其阈值的确定充分考量整体数据的方差值，准确度较高，所得数据集合与约束值为下偏差，整体数值较为接近，偏离度为-1.04%，仍需寻求偏离度较小结构。

动态范围压缩DRC 法由于数映射将多数中间图像密度单元向上边缘逼近，导致整体数据偏上限输出。 因此所得数据呈上偏差，整体偏离程度为＋2.11%，仍存在中间密度单元，优化效果不理想。

阶梯量化算法由于阶梯值为均匀选取，因此对中间灰度的放缩效应较差，导致图像较为模糊。 所得数据呈下偏差，且偏离度达到-3.73%，整体优化效果不明显，且存在较多灰度单元。

普通全局阈值分切算法较好的区分图像结构，有效抑制中间密度单元，图像较为清晰，但是整体数据偏差＋1.73%，不能满足优化后结构偏离度低的要求。

自适应阈值算法与普通阈值相比，在灰度级整体偏移程度上有较大提升，几乎完全覆盖约束条件，且优化图像黑白分明，可以有效完成中间密度单元向0和1 边界元素的放缩，整体数据偏离度仅为0.01%。

为验证自适应阈值算法针对不同拓扑结构图像具有普适性，通过不同网格规模的对比实验，验证算法的可靠性。 对不同结构采用自适应阈值优化后，对比数据见表2。

由表2 可知，对不同网格规模的拓扑结构进行优化计算，所得灰度集在标准值附近上下浮动，整体差异不大，偏离度均保持在0.02%以内，可以视为完全复原。 经过自适应阈值方法过滤所得结构在不改变整体体积约束条件下，中间密度单元明显得到过滤，完全实现图像黑白单元化。

4 结束语

本文针对拓扑优化变密度法中运用SIMP 模型得到优化结构中存在大量灰度单元，且中间密度单元对整体结构清晰度影响较大的问题，通过数值化提取拓扑结构灰度单元，依据最大类间方差法确定初始阈值，接着通过histogram 函数确定均匀宽度值，依据灰度值约束条件建立自适应阈值计算公式，得到自适应阈值算子，并通过典型算例验证方法有效性，与多种数字图像增强算法对比，本文提出算法在拓扑约束偏离值仅为0.01%情况下实现拓扑灰度图像有效过滤，整体优化结构清晰。