何炎高， 徐定華， 陳瑞林

(1. 浙江省服裝工程技術(shù)中心(浙江理工大學(xué))， 浙江 杭州 310018； 2. 浙江理工大學(xué) 理學(xué)院， 浙江 杭州 310018)

?

紡織材料設(shè)計(jì)反問題的貝葉斯統(tǒng)計(jì)推斷方法

何炎高1,2， 徐定華1,2， 陳瑞林1,2

(1. 浙江省服裝工程技術(shù)中心(浙江理工大學(xué))， 浙江 杭州 310018； 2. 浙江理工大學(xué) 理學(xué)院， 浙江 杭州 310018)

針對具有不適定性紡織材料設(shè)計(jì)反問題，給出了利用貝葉斯蒙特卡洛方法求解紡織材料單參數(shù)和多參數(shù)反演問題的一種新方法。因織物穩(wěn)態(tài)熱濕傳遞模型的非線性和反問題的不適定性，基于貝葉斯統(tǒng)計(jì)推斷方法的紡織材料類型、厚度、孔隙率等參數(shù)的后驗(yàn)概率分布推斷是一種有效的方法。這種方法將參數(shù)的先驗(yàn)信息描述為先驗(yàn)概率密度，構(gòu)建了紡織材料設(shè)計(jì)反問題的數(shù)值算法。數(shù)值實(shí)驗(yàn)結(jié)果表明，與馬爾科夫鏈蒙特卡洛抽樣算法相匹配的貝葉斯推理可用來求解紡織材料設(shè)計(jì)反問題。

紡織材料； 設(shè)計(jì)； 反問題； 貝葉斯推斷； 單參數(shù)； 多參數(shù)

隨著人民生活質(zhì)量的改善和科技水平的提高，紡織材料的功能和應(yīng)用領(lǐng)域也不斷增多，紡織服裝的熱濕舒適性要求也日益受到關(guān)注，因此，基于人體舒適性要求的紡織材料設(shè)計(jì)具有重要的現(xiàn)實(shí)意義。一般地說，在人體-服裝-環(huán)境系統(tǒng)中，基于織物的熱濕傳遞模型，根據(jù)熱濕傳遞方程和初邊值條件，給定織物的物理參數(shù)和結(jié)構(gòu)參數(shù)，從而計(jì)算人體與織物間微氣候區(qū)的溫度或濕度，該問題被稱為正問題(DP: direct problems)。根據(jù)服裝的熱濕舒適性要求來決定紡織材料的物理參數(shù)或結(jié)構(gòu)參數(shù)，該問題被稱為反問題(IP: inverse problems)。通常由于測量數(shù)據(jù)有限且?guī)в幸欢ǖ恼`差，使得反問題具有不適定性，因而在求解時(shí)存在較大的困難，針對其不適定性問題，常見的反演方法有正則化方法, 此時(shí)往往轉(zhuǎn)化為最優(yōu)化問題進(jìn)行參數(shù)估計(jì)，如 Hooke-Jevees模式搜索算法[1]、0.618法、或粒子群算法[2]來求解。近期建立在統(tǒng)計(jì)學(xué)基礎(chǔ)上的貝葉斯推理在污染源識別反問題[4]、熱傳導(dǎo)反問題[5]、熱輻射源估計(jì)[6]、熱參數(shù)的估計(jì)[7]等研究中發(fā)揮了重要的作用并可較好地對反問題進(jìn)行求解。

目前，紡織材料設(shè)計(jì)反問題研究尚處于初步階段，本文針對織物穩(wěn)態(tài)熱濕傳遞反問題建立了貝葉斯推理的反演算法，采用馬爾科夫鏈蒙特卡洛抽樣的方法對后驗(yàn)狀態(tài)空間進(jìn)行抽樣并獲得了紡織材料參數(shù)的后驗(yàn)概率分布規(guī)律及進(jìn)行相應(yīng)的估計(jì)，可為紡織材料產(chǎn)品設(shè)計(jì)或?qū)嶋H生產(chǎn)提供理論參考和實(shí)踐指導(dǎo)。

1 紡織材料熱濕傳遞模型及反問題

在如圖1所示的人體-服裝-環(huán)境系統(tǒng)示意圖中，考慮其平行圓柱孔的單層織物的穩(wěn)態(tài)熱濕傳遞模型[8-9]為

(1)

其初邊值條件為

(2)

式中：k1和k2均為與水分子質(zhì)量和氣體常數(shù)相關(guān)的常數(shù)；ε(x)為紡織品表面孔隙率，%；r(x)為纖維孔半徑，m；τ(x)為纖維孔的曲折系數(shù)；pv(x)為水蒸氣壓力，Pa；T(x)為織物溫度，K；mv(x)為水蒸氣質(zhì)量通量，kg/(m2·s)；Γ(x)為水蒸氣凝結(jié)率，kg/(m3·s)；λ為水蒸氣吸收凝結(jié)熱，J/kg；κ為織物的熱傳導(dǎo)系，W/(m·K)；T(0)為織物內(nèi)側(cè)溫度；T(L)為織物外側(cè)溫度；mv(0)為織物內(nèi)側(cè)水蒸氣質(zhì)量通量；pv(0)為織物內(nèi)側(cè)水蒸氣壓力。

圖1 人體-服裝-環(huán)境系統(tǒng)示意圖Fig.1 Schematic diagram of body-clothing- environment system

飽和水蒸氣壓力的經(jīng)驗(yàn)公式為

求解上述帶有初邊值條件的耦合常微分方程組稱為正問題(DP)，由Banach不動點(diǎn)定理可以證明其解的存在性和唯一性[9]。

服裝舒適性是紡織材料為滿足人體生理需要所必備的性能，也是紡織材料設(shè)計(jì)的核心要求。一般認(rèn)為人體皮膚與服裝內(nèi)側(cè)間的微氣候區(qū)內(nèi)溫度(32±1)℃、濕度(50±10)%、氣流(25±15)cm/s為標(biāo)準(zhǔn)服裝氣候，即為熱濕舒適性指標(biāo)[9]。因此，在一定溫度和濕度的環(huán)境下，根據(jù)服裝的熱濕舒適性要求，決定織物的物理參數(shù)(如熱傳導(dǎo)系數(shù)κ)和結(jié)構(gòu)參數(shù)(如厚度L、孔隙率ε)稱為紡織材料設(shè)計(jì)反問題(inverse problem of textile material design，IPTMD)。

單參數(shù)的決定：給定環(huán)境的溫度和濕度，根據(jù)服裝的熱濕舒適性指標(biāo)，設(shè)計(jì)織物的熱傳導(dǎo)系數(shù)κ、厚度L或孔隙率ε，分別稱為類型決定、厚度決定或孔隙率決定。

多參數(shù)的決定：給定環(huán)境的溫度和濕度，根據(jù)服裝的熱濕舒適性指標(biāo)，設(shè)計(jì)織物的熱傳導(dǎo)系數(shù)κ、厚度L、孔隙率ε中2個(gè)或全部參數(shù)。

2 正問題和反問題的數(shù)值算法

2.1 正問題DP的數(shù)值解法

將微分方程組模型(1)與(2)解耦[1]得

其中

用有限差分法離散得到以下差分方程：

當(dāng)i=2,…,N-1時(shí)，

當(dāng)i=N時(shí)，

已知T0、TN，通過插值得到T1、TN-1，這樣通過以上差分方程可以計(jì)算TN-2,…,T2，由此得到

進(jìn)而微氣候區(qū)的相對濕度(RH)的表達(dá)式可表示為

(3)

2.2 反問題IPTMD的貝葉斯推斷方法

貝葉斯推理的基礎(chǔ)是貝葉斯定理，即

(4)

從式(4)可以看出參數(shù)的所有信息都包含在后驗(yàn)分布中，一旦知道了后驗(yàn)概率密度函數(shù)的分布規(guī)律，就可以利用點(diǎn)估計(jì)的辦法，如最大后驗(yàn)估計(jì)(MAP)：

后驗(yàn)均值估計(jì)：

同樣也可以作區(qū)間估計(jì)。

通常情況下，測量數(shù)據(jù)的邊緣概率密度函數(shù)π(y)在后驗(yàn)狀態(tài)空間中可看成積分常數(shù)沒有進(jìn)行計(jì)算的必要，因而后驗(yàn)概率密度函數(shù)可以簡單地表示為

(5)

一般地，織物參數(shù)θi(i=1,2…m；m為模型參數(shù)的個(gè)數(shù))在一定的范圍θi∈[ai,bi]內(nèi)且滿足均勻分布，先驗(yàn)概率密度函數(shù)可表示為

(6)

由于各參數(shù)之間具有相互獨(dú)立性，則總的先驗(yàn)分布可表示為

(7)

測量誤差一般可以認(rèn)為是白噪聲η，其每個(gè)分量誤差均服從均值為零、標(biāo)準(zhǔn)差為σ的正態(tài)分布N(0,σ2)，似然函數(shù)可表示為

(8)

式中n為測量數(shù)據(jù)的個(gè)數(shù)。

從理論上講，利用式(5)即可求出后驗(yàn)概率密度函數(shù)，但往往由于參數(shù)維數(shù)較大或正演關(guān)系比較復(fù)雜難以得到明確的數(shù)學(xué)表達(dá)式，使得數(shù)值積分算法計(jì)算量呈指數(shù)增長，因而計(jì)算難度較大，為此，需采用特定的抽樣方法實(shí)現(xiàn)對后驗(yàn)概率密度進(jìn)行求解。

2.3 馬爾科夫鏈蒙特卡洛法

馬爾科夫鏈蒙特卡洛(Markov Chain Monte Carlo，MCMC)法是一類算法的總稱。從數(shù)學(xué)上講，其思想是產(chǎn)生一個(gè)Markov鏈，以目標(biāo)分布為平穩(wěn)分布。根據(jù)Markov鏈理論，一個(gè)Markov鏈從任意初值出發(fā)，都會收斂到其平穩(wěn)分布。MCMC就是用馬爾科夫鏈的平穩(wěn)分布，生成滿足特定分布的隨機(jī)數(shù)構(gòu)造一個(gè)目標(biāo)分布樣本。

Metropolis算法[4]是一種重要的MCMC抽樣算法，其算法可表述如下。

1)在模型參數(shù)先驗(yàn)范圍內(nèi)隨機(jī)產(chǎn)生模型參數(shù)初始點(diǎn)θ(i),i=1；

式中π(θ*)，π(θ(i))為目標(biāo)概率密度函數(shù)。

3)產(chǎn)生一個(gè)0～1之間的隨機(jī)數(shù)u；

4)如果u

5)重復(fù)步驟3)、4)直到達(dá)到迭代次數(shù)。

利用馬爾科夫蒙特卡洛方法對后驗(yàn)狀態(tài)空間進(jìn)行抽樣，在大數(shù)法則的原理下根據(jù)樣本均值依概率收斂性，可用收斂的樣本均值去估計(jì)期望值[3]。

3 數(shù)值實(shí)驗(yàn)與結(jié)果分析

式中θ=(θ1,θ2,θ3),θ1=κ,θ2=L,θ3=ε。

在2種不同低溫情形下進(jìn)行數(shù)值實(shí)驗(yàn), 微氣候區(qū)o點(diǎn)作為觀察點(diǎn)，觀察數(shù)據(jù)源于2種情形下微氣候區(qū)相對濕度的計(jì)算值(如圖2所示)，本文采用棉、羊毛、滌綸等常見的紡織材料或其混紡材料進(jìn)行數(shù)值實(shí)驗(yàn)，相關(guān)參數(shù)及先驗(yàn)信息如表1所示，針對2種不同情形運(yùn)用貝葉斯理論對反演參數(shù)進(jìn)行后驗(yàn)均值估計(jì)。

圖2 2種情形下微氣候區(qū)內(nèi)相對濕度的近似值Fig.2 Approximate RH in microclimate area at two cases.(a)Case 1; (b) Case 2

參數(shù)數(shù)值參數(shù)數(shù)值k/(W·(m·K)-1)0.1T(0)/K305.15ε0.9mv(0)/(kg·(m2·s)-1)3.083×10-5r/m1×10-5k16×10-5τ1.2k27×10-5λ/(J·kg-1)3.593×106RHe(0.4,0.9)

3.1 反演算例1：單參數(shù)κ或L的反演

針對情形1和情形2，待反演的模型參數(shù)為織物的熱傳導(dǎo)系數(shù)κ或厚度L。根據(jù)先驗(yàn)信息，待定參數(shù)的先驗(yàn)分布采用均勻分布，則其對應(yīng)的先驗(yàn)概率密度函數(shù)分別為

數(shù)值實(shí)驗(yàn)?zāi)M過程中假定觀察誤差服從白噪聲N(0,σ2),σ=0.001，則似然函數(shù)條件概率密度分別為

根據(jù)貝葉斯定理及相應(yīng)理論，織物的熱傳導(dǎo)系數(shù)和厚度后驗(yàn)概率密度函數(shù)分別為

式中，λ1、λ2為常數(shù)。

圖3 情形1下參數(shù)κ反演結(jié)果Fig.3 Inversion results of parameter κ in case 1.(a) Iteration curve graph; (b) Posterior histogram

圖4 情形2下參數(shù)L反演結(jié)果Fig.4 Inversion results of parameter L in Case 2.(a) Iteration curve graph; (b) Posterior histogram

3.2 反演算例2:參數(shù)(L,ε)同時(shí)反演

分別對2種情形下參數(shù)L和ε進(jìn)行同時(shí)反演。根據(jù)先驗(yàn)信息，待反演模型參數(shù)先驗(yàn)分布采用均勻分布，其對應(yīng)的先驗(yàn)概率密度函數(shù)為

數(shù)值實(shí)驗(yàn)?zāi)M過程中假定觀察誤差η為白噪聲N(0,σ2),σ=0.001，則似然函數(shù)條件概率密度可表示為

根據(jù)貝葉斯定理及相應(yīng)理論，參數(shù)(L,ε)后驗(yàn)概率密度函數(shù)為

式中λ為常數(shù)。

圖5 參數(shù)L和ε反演迭代曲線Fig.5 Iteration curves of parameters L and ε.(a) Case 1 (L); (b) Case 1(ε); (c) Case 2(L); (d) Case 2 (ε)

3.3 反演結(jié)果分析

從圖3～5中可看出，迭代曲線開始具有一定周期振蕩的特性，當(dāng)?shù)螖?shù)達(dá)到一定次數(shù)后趨于穩(wěn)定狀態(tài)，表明馬爾科夫鏈進(jìn)入穩(wěn)定的收斂區(qū)域。從后驗(yàn)概率直方圖可看出，單參數(shù)反演時(shí)，κ和L分別在0.1和0.008附近時(shí)的頻數(shù)最大，與預(yù)設(shè)定值完全符合，參數(shù)的后驗(yàn)統(tǒng)計(jì)結(jié)果如表2所示。2種情形下參數(shù)L、ε同時(shí)反演時(shí)分別在0.005、0.9和0.008、0.9附近時(shí)的頻數(shù)最大，與預(yù)設(shè)定值完全符合，參數(shù)的后驗(yàn)統(tǒng)計(jì)結(jié)果如表3所示。

圖6 參數(shù)L和ε后驗(yàn)概率直方圖Fig.6 Posterior histograms of parameters L and ε.(a) Case 1 (L); (b) Case 1 (ε); (c) Case 2 (L); (d) Case 2 (ε)

參數(shù)真值迭代次數(shù)后驗(yàn)均值估計(jì)均值誤差/%κ0.110000.1030413.04L0.00810000.0080480.60

表3 L和ε的后驗(yàn)統(tǒng)計(jì)結(jié)果

結(jié)果表明，采用Bayesian-MCMC抽樣算法的貝葉斯推理方法反演誤差較小、精度較高。因而，可以利用反演結(jié)果對紡織材料類型選擇和厚度設(shè)計(jì)進(jìn)行理論指導(dǎo)，后驗(yàn)估計(jì)中可剔除前面沒有達(dá)到平衡的點(diǎn)進(jìn)行分析或進(jìn)行參數(shù)反演時(shí)將Bayesian-MCMC方法與其他優(yōu)化算法結(jié)合起來。在先驗(yàn)信息很少場合下，通常取無信息的先驗(yàn)分布，常用的有Reformulation方法、不變Haar測度方法、Box-Tiao技術(shù)、Lindley方法及最大熵原理。

4 結(jié) 論

本文針對紡織材料穩(wěn)態(tài)熱濕傳遞模型中的參數(shù)決定，建立了基于貝葉斯推斷的反演數(shù)學(xué)模型及其數(shù)值算法，運(yùn)用貝葉斯推斷理論，提出了一種用Bayesian-MCMC方法決定紡織材料設(shè)計(jì)反問題的新方法。該方法可用來求解紡織材料熱傳導(dǎo)系數(shù)、厚度、空隙率等單參數(shù)或多參數(shù)設(shè)計(jì)反演問題。

Bayesian-MCMC方法是一種充分利用先驗(yàn)信息通過構(gòu)造Markov鏈進(jìn)行隨機(jī)動態(tài)模擬，得到后驗(yàn)概率密度分布，具有估計(jì)精度較高，計(jì)算速度快等特點(diǎn)，針對復(fù)雜空間上多維估值問題，MCMC方法具有很強(qiáng)的適用性，同時(shí)求解比其他方法更容易。理論上，采用Bayesian-MCMC方法只要采樣次數(shù)足夠多，馬爾科夫鏈一定收斂且可以精確得到后驗(yàn)分布。

[1] 徐定華，陳遠(yuǎn)波，程建新.低溫環(huán)境下紡織材料類型設(shè)計(jì)反問題[J].紡織學(xué)報(bào),2011,32(9):23-28. XU Dinghua, CHEN Yuanbo, CHENG Jianxin. Inverse problem of textile material design at low tempera-ture[J].Journal of Textile Research,2011,32(9):23-28.

[2] CUI Peng, XU Dinghua. Simultaneous determination of thickness and heat conductivity in fabric design: steady-state modeling and PSO algorithms[J]. Journal of Fiber Bioengineering and Informatics,2013(2):195-204.

[3] KAIPIO J,SOMERSALO E. Statistical and Computational Inverse Problems[M]. New York:Springer, 2004:91-108.

[4] 陳海洋，騰彥國，王金生，等.基于Bayesian-MCMC方法的水體污染識別反問題[J].湖南大學(xué)學(xué)報(bào):自然科學(xué)版,2012(6):74-78. CHEN Haiyang, TENG Yanguo, WANG Jinsheng, et al. Event source identification of water pollution based on Bayesian-MCMC[J]. Journal of Hunan University： Natural Science Edition, 2012(6):74-78.

[5] WANG Jingbo, ZABARAS Nicholas.A Bayesian inference approach to the inverse heat conduction problem[J]. International Journal of Heat and Mass Transfer, 2004(47):3927-3914.

[6] WANG Jingbo, ZABARAS Nicholas. Using Bayesian statistics in the estimation of heat source in radiation[J]. International Journal of Heat and Mass Transfer, 2005(48):15-29.

[7] 高思云，楊晨.利用貝葉斯模型進(jìn)行熱參數(shù)估計(jì)[J].系統(tǒng)仿真學(xué)報(bào),2006(6):1462-1465. GAO Siyun, YANG Chen. Bayesian model to parameter estimation[J].Journal of System Simulation, 2006(6): 1462-1465.

[8] XU Dinghua, CHENG Jianxin, ZHOU Xiaohong. A model of heat and moisture transfer through parallel pore Textiles[J]. Proceedings of Textile Bioengineering and Informatics Symposium, 2010(5):1139-1144.

[9] 徐定華. 紡織材料熱濕傳遞數(shù)學(xué)模型及設(shè)計(jì)反問題[M].北京:科學(xué)出版社,2014:25-29. XU Dinghua. Mathematical Model of Heat and Moisture Transfer and Inverse Problems of Textile Material Design[M]. Beijing: Science Press, 2014:25-29.

[10] XU Dinghua, GE Meibao.Thickness determination in textile meterial design:dynamic modeling and numerical algorithms[J]. Inverse Problems, 2012(28)：035011.

Bayesian statistical inference method for inverse problems of textile material design

HE Yangao1,2, XU Dinghua1,2, CHEN Ruilin1,2

(1.EngineeringResearchCenterofClothingofZhejiangProvince(ZhejiangSci-TechUniversity),Hangzhou,Zhejiang310018,China； 2.SchoolofSciences,ZhejiangSci-TechUniversity,Hangzhou,Zhejiang310018,China)

Aiming at the ill-posed of the inverse problem of textile material design(IPTMD), a new approach based on Bayesian Markov Chain Monte Carlo(Bayesian-MCMC) method was proposed for solving the problem of single and multiple parameter determination. Since the heat and moisture transfer model is non-linear and the IPTMD is ill-posed, it is proved that posterior distribution for the model parameters such as the heat conductivity, thickness and porosity of the material is an effective method. This method describes prior information of parameters as prior probability density, constructing numerical algorithms of the IPTMD. The numerical results show that Bayesian inference method agreed with Markov Chain Monte Carlo sampling algorithm can be applied to solve the IPTMD.

textile material; design; inverse problem; Bayesian inference; single parameter; multiple parameters

10.13475/j.fzxb.201501002307

2013-11-25

2014-10-08

國家自然科學(xué)基金資助項(xiàng)目(11071221，11471287)；浙江省高校重中之重紡織材料與工程一級學(xué)科和浙江省服裝工程技術(shù)研究中心開放基金項(xiàng)目(2013KF10)；浙江醫(yī)學(xué)高等?？茖W(xué)校科研項(xiàng)目(2014XZA01)

何炎高(1985—)，男，碩士生。主要研究方向?yàn)槲⒎址匠谭磫栴}的理論、計(jì)算及應(yīng)用。徐定華，通信作者，E-mail:dhxu6708@zstu.edu.cn。

TS 101.1

A