基于自適應閾值的拓撲圖像密度單元優化

王 沖， 茅 健， 鄭 武

（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%情況下實現拓撲灰度圖像有效過濾，整體優化結構清晰。