基于遺傳算法反求耳結(jié)構(gòu)彈性模量

程向東 姚文娟

(上海大學(xué)土木工程系，上海200072)

基于遺傳算法反求耳結(jié)構(gòu)彈性模量

程向東 姚文娟*

(上海大學(xué)土木工程系，上海200072)

本研究通過(guò)力學(xué)反問(wèn)題原理，利用已知的位移求解耳結(jié)構(gòu)彈性模量。隨機(jī)產(chǎn)生遺傳算法的初始種群，使用自編的Matlab算法程序，對(duì)初始種群進(jìn)行遺傳迭代計(jì)算，把已知的目標(biāo)位移與種群位移的均方差作為目標(biāo)函數(shù)，以目標(biāo)函數(shù)值最小控制迭代進(jìn)化的方向。通過(guò)耳結(jié)構(gòu)的鼓膜凸和鐙骨底板2個(gè)控制位移以及砧鐙關(guān)節(jié)周?chē)?個(gè)控制位移這兩個(gè)算例求解正常砧鐙關(guān)節(jié)的彈性模量，并使用耳結(jié)構(gòu)的鼓膜凸和鐙骨底板2個(gè)控制位移求解病變砧鐙關(guān)節(jié)的彈性模量。結(jié)果表明，使用基于遺傳算法的反問(wèn)題方法計(jì)算耳結(jié)構(gòu)的彈性模量是可行的，并且具有穩(wěn)定性和不受結(jié)構(gòu)力學(xué)性能影響的特點(diǎn)，相對(duì)誤差分別為0.05%和0.2%、0.03%，可為臨床病變耳提供有效的力學(xué)參數(shù)。

耳結(jié)構(gòu);反問(wèn)題;彈性模量;遺傳算法

引言

隨著科學(xué)的發(fā)展，學(xué)科之間的交叉研究已成為新的發(fā)展方向。把力學(xué)原理應(yīng)用于探索生命科學(xué)是交叉學(xué)科研究的典范。20世紀(jì)末至21世紀(jì)初，耳生物力學(xué)研究開(kāi)始起步，國(guó)內(nèi)外很多學(xué)者應(yīng)用力學(xué)原理對(duì)耳結(jié)構(gòu)進(jìn)行研究，其一是解析分析方法，如依據(jù)彈性力學(xué)原理推導(dǎo)基底膜振動(dòng)方程［1］，依據(jù)變分原理推導(dǎo)耳膜的振動(dòng)方程［2］，建立檢驗(yàn)人工聽(tīng)骨力學(xué)性質(zhì)的解析模型等［3］;其二是有限元數(shù)值模擬計(jì)算分析，如通過(guò)高分辨率計(jì)算機(jī)輔助斷層或組織切片，建立了中耳及全耳的三維有限元模型［4－5］，研究砧骨假體頭的尺寸對(duì)振動(dòng)影響［6］，研究鐙骨置換假體接頭形狀對(duì)聽(tīng)力傳導(dǎo)影響［7］。通過(guò)CT掃描建立了中耳及全耳的三維有限元模型，得到鼓膜及聽(tīng)骨鏈的最大應(yīng)力區(qū)域就是鼓膜穿孔及聽(tīng)骨鏈病變部位，同時(shí)得到人工鐙骨較為合理的形體［8－10］，得到外耳道、中耳砧鐙關(guān)節(jié)、內(nèi)耳淋巴液對(duì)傳聲特性的影響，并得到了鼓膜病變對(duì)傳聲的影響［11－12］。

應(yīng)用力學(xué)原理探索耳結(jié)構(gòu)的研究中的一個(gè)必要前提是力學(xué)參數(shù)的確定，以前研究力學(xué)參數(shù)的方法是通過(guò)對(duì)正常人耳進(jìn)行實(shí)驗(yàn)測(cè)試得到［13－14］，如彈性模量［15］及本構(gòu)關(guān)系等［16］。在耳結(jié)構(gòu)發(fā)生病變以后，其力學(xué)參數(shù)會(huì)發(fā)生相應(yīng)的改變，并且無(wú)法對(duì)臨床病變耳進(jìn)行復(fù)雜的實(shí)驗(yàn)研究以獲得力學(xué)參數(shù)，便不能應(yīng)用力學(xué)原理對(duì)其進(jìn)行探索。本研究基于遺傳算法，通過(guò)耳結(jié)構(gòu)在外力作用下的位移，反求耳結(jié)構(gòu)的力學(xué)參數(shù)。為無(wú)法通過(guò)實(shí)驗(yàn)測(cè)試的病變活體結(jié)構(gòu)力學(xué)參數(shù)的反演提供了簡(jiǎn)捷有效的方法。

1 遺傳算法基本理論

遺傳算法(genetic algorithm，GA)是由美國(guó)Michigan大學(xué)教授 Holland及其學(xué)生于1975年創(chuàng)建，是基于對(duì)自然界中生物遺傳和進(jìn)化機(jī)理的模仿。對(duì)遺傳算法基礎(chǔ)理論的研究主要分3個(gè)方面:模式定理的拓廣和深入、遺傳算法的新模型和遺傳算法的收斂性理論。隱含并行性原理和模式定理被看作是遺傳算法的兩大基石。Holland給出了模式定理:具有低階、短的定義長(zhǎng)度、模式采樣的平均適值在種群平均適值以上的模式，在遺傳迭代過(guò)程中將按指數(shù)增長(zhǎng)率被采樣。模式定理可表達(dá)為

式中，m(H，t)表示第t代群體中包含模式H的染色體，f(H)表示在第t代群體中包含模式H的染色體的平均適應(yīng)度值，為第t代群體的所有染色體的平均適應(yīng)度值。L為染色體的長(zhǎng)度，Pc、Pm為交叉、變異概率，O(H)、δ(H)為模式位數(shù)和模式定義長(zhǎng)度。

針對(duì)不同的問(wèn)題，由不同的編碼方法和不同的遺傳算子就構(gòu)成了各種不同的遺傳算法［17］。遺傳算法有諸多優(yōu)點(diǎn):直接以目標(biāo)函數(shù)作為搜索信息，并且使用多個(gè)搜索點(diǎn)的搜索信息［18］，具有全局搜索能力，不宜陷入局部最優(yōu)等［19］。

基本遺傳算法(simple genetic algorithm)數(shù)學(xué)模型可表示為［17］

式中，C為 個(gè)體編碼方法;V為 個(gè)體適應(yīng)度評(píng)價(jià)方法;pop為初始種群;popsize為種群大小;Φ選擇算子;Γ為交叉算子;Ψ為變異算子;Τ為算法終止條件。對(duì)于右邊的參數(shù)和算子，確定好參數(shù) C、pop、popsize、Τ以后，通過(guò)算子 Φ、Γ、Ψ 來(lái)模仿生物的遺傳進(jìn)化，以V控制迭代進(jìn)化的方向。

流程圖可表示為圖1。首先確定編碼方式 C，并產(chǎn)生初始種群pop，對(duì)種群進(jìn)行適應(yīng)度值計(jì)算V，然后基于種群個(gè)體的適應(yīng)度值，對(duì)種群進(jìn)行遺傳算子Φ、Γ、Ψ操作，使種群不斷進(jìn)化，直至新產(chǎn)生種群中的個(gè)體足夠優(yōu)秀。

圖1 基本遺傳算法流程圖Fig.1 Flow chart of simple genetic-algorithm

2 求解思路

基于遺傳算法，通過(guò)已知位移 un(n=1，2，3，…，k)求解耳結(jié)構(gòu)的彈性模量。如圖1所示，基于遺傳算法的迭代計(jì)算，首先需要產(chǎn)生一個(gè)初始種群，然后通過(guò)設(shè)定的適應(yīng)度函數(shù)對(duì)種群的個(gè)體進(jìn)行適應(yīng)度的計(jì)算，基于個(gè)體的適應(yīng)度值，進(jìn)行選擇、交叉、變異算子的操作，使種群往最優(yōu)解方向進(jìn)化，每個(gè)循環(huán)結(jié)束后，評(píng)價(jià)適應(yīng)度值是否達(dá)到最優(yōu)解，如果未達(dá)到最優(yōu)解則進(jìn)行下一個(gè)循環(huán)，直到達(dá)到最優(yōu)停止。

相對(duì)于反問(wèn)題，正問(wèn)題的求解是成熟的，本文通過(guò)耳結(jié)構(gòu)有限元模型來(lái)進(jìn)行正問(wèn)題的計(jì)算，并把計(jì)算的參數(shù)和結(jié)果作為反問(wèn)題遺傳算法的初始種群，本研究用的有限元模型已經(jīng)驗(yàn)證了正確性［9，10］。

3 程序設(shè)計(jì)

整個(gè)算法的實(shí)現(xiàn)包括9個(gè)步驟。

步驟1:設(shè)置遺傳算法參數(shù)。編碼方式C采用實(shí)數(shù)編碼。種群的大小即種群個(gè)體的數(shù)量，用popsize表示。種群個(gè)體的大小即個(gè)體的長(zhǎng)度，用nvars表示。最大迭代次數(shù)用maxgen表示。交叉算子用pc表示;變異算子用 pm表示;繁殖因子用 pg表示。循環(huán)結(jié)束的最佳適應(yīng)度界限值條件用delt表示。

步驟2: 初始種群pop的產(chǎn)生。對(duì)耳結(jié)構(gòu)有限元模型賦予不同的彈性模量，計(jì)算出不同的位移，把計(jì)算的數(shù)據(jù)保存于 matlab工作目錄下，命名為 my_data.txt，其大小為 popsize × nvars。算法程序中導(dǎo)入文件my_data.txt作為初始種群即可。

步驟3:適應(yīng)度值計(jì)算。遺傳算子是基于種群個(gè)體適應(yīng)度值進(jìn)行的，所以先計(jì)算種群中個(gè)體的適應(yīng)度值，文中適應(yīng)度值通過(guò)下式計(jì)算［19］

種群個(gè)體位移 u'1，u'2，…，u'k與目標(biāo)位移 u1，u2，…，uk的平方差，因?yàn)楸狙芯恐惺钦业溅盏臉O小值即認(rèn)為得到最優(yōu)解，所以適應(yīng)度值需要在φ前加一個(gè)負(fù)號(hào)，即適應(yīng)度值Fit為:

計(jì)算出適應(yīng)度值Fit以后，存放在種群pop的第nvars+1列，然后依次計(jì)算每個(gè)個(gè)體得到的適應(yīng)度值Fit均存放在每個(gè)個(gè)體的 nvars+1個(gè)位置，適應(yīng)度值計(jì)算完成后，種群 pop的大小變成 popsize×nvars+1。

步驟4:保持最優(yōu)個(gè)體keepbest。保持最優(yōu)個(gè)體函數(shù) keepbest實(shí)質(zhì)上也是一個(gè)選擇算子［20］。初始種群是隨機(jī)產(chǎn)生的，初始種群中包含很多個(gè)體是具有低適應(yīng)度，不適合遺傳到下一代的個(gè)體，但是也有適應(yīng)度比較高的個(gè)體。為了保證最優(yōu)的個(gè)體不會(huì)在遺傳操作中被淘汰掉，在選擇交叉變異算子進(jìn)行操作前，先進(jìn)行保持最優(yōu)個(gè)體算子的操作。依次比較種群中每個(gè)個(gè)體的適應(yīng)度值大小，把適應(yīng)度值最大的一個(gè)個(gè)體存放于popsize+1行，保持最優(yōu)個(gè)體算子計(jì)算完成后種群pop的大小是(popsize+1)×(nvars+1)。

步驟5:選擇操作。選擇算子采用隨機(jī)競(jìng)爭(zhēng)選擇［17］。選擇操作的主要目的是為了避免有用遺傳信息的丟失，提高全局收斂性和計(jì)算效率。在選擇操作中引入繁殖因子 pg(一個(gè)操作算子)［21］，來(lái)調(diào)節(jié)競(jìng)爭(zhēng)種群規(guī)模與適應(yīng)性種群規(guī)模的比例。根據(jù)達(dá)爾文進(jìn)化論中大量繁殖的概念，將群體分成適應(yīng)性群體與競(jìng)爭(zhēng)群體兩種，競(jìng)爭(zhēng)群體的規(guī)模遠(yuǎn)大于適應(yīng)性群體的規(guī)模。這樣，除了把適應(yīng)性群體中的優(yōu)良基因型延續(xù)到下一代，還增加了一些新的基因型，實(shí)際上就是使群體包含更多的模式，從而避免陷入局部最優(yōu)解，減少未成熟收斂機(jī)會(huì)，增大搜索到全局最優(yōu)解的可能性［22］。

選擇、繁殖的操作過(guò)程:競(jìng)爭(zhēng)種群用 selpool表示，種群大小為(popsize×pg)×(nvars+1)。對(duì)適應(yīng)性種群pop的popsize個(gè)個(gè)體進(jìn)行隨機(jī)選擇，每次選擇兩個(gè)個(gè)體，然后對(duì)所選擇的兩個(gè)個(gè)體的適應(yīng)度值進(jìn)行比較，適應(yīng)度值大的一個(gè)個(gè)體就被提取到競(jìng)爭(zhēng)種群 selpool，重復(fù)操作直到填滿(mǎn)競(jìng)爭(zhēng)種群的popsize×pg個(gè)個(gè)體。

步驟6:交叉算子。交叉算子采用算術(shù)交叉［22～24］，交叉算子的目的:在生物自然進(jìn)化進(jìn)程中，交配重組是一個(gè)重要環(huán)節(jié)，兩個(gè)同源染色體通過(guò)交配而重組，形成新的染色體，從而產(chǎn)生出新的個(gè)體或物種。按照較大的概率從種群中選擇兩個(gè)個(gè)體，交換某些位產(chǎn)生子代，子代繼承了父代的基本特征。本算法采用的編碼方式是實(shí)數(shù)編碼，實(shí)數(shù)編碼對(duì)問(wèn)題的描述是自然的，一致性的，且精確度較二進(jìn)制高。算術(shù)交叉的表述:兩個(gè)基因x、y進(jìn)行算術(shù)交叉得到的子代的基因是:

式中，α為一個(gè)0～1之間的隨機(jī)數(shù)。

從競(jìng)爭(zhēng)種群的第一個(gè)個(gè)體開(kāi)始，先選定第一個(gè)個(gè)體，再在競(jìng)爭(zhēng)種群中另外隨機(jī)選擇一個(gè)個(gè)體，與選定的個(gè)體進(jìn)行算術(shù)交叉。然后選定第二個(gè)個(gè)體，進(jìn)行同樣的操作，如此依次對(duì)競(jìng)爭(zhēng)種群中的個(gè)體進(jìn)行同樣的交叉操作。

步驟7:變異算子。采用均勻變異算子［17］，均勻變異是指分別用符合某一范圍內(nèi)均勻分布的隨機(jī)數(shù)，以某一較小的概率來(lái)替換個(gè)體編碼串中各個(gè)基因座上的原有基因值。它使搜索點(diǎn)在整個(gè)搜索空間移動(dòng)，并且比非均勻變異更有效的對(duì)某一重點(diǎn)區(qū)域進(jìn)行局部搜索。變異算子的目的:(1)改善遺傳算法的局部搜索能力。(2)維持群體的多樣性，避免出現(xiàn)早熟現(xiàn)象。變異算子用新的基因值替換原有的基因值，從而可以改變個(gè)體編碼串的結(jié)構(gòu)維持群體多樣性，防止早熟，

從競(jìng)爭(zhēng)種群的第一個(gè)個(gè)體開(kāi)始，先確定變異點(diǎn)，變異點(diǎn)是一個(gè)1到nvars之間的隨機(jī)數(shù)，變異點(diǎn)mpoint=round(rand×(nvars+1)－1)，round表示對(duì)一個(gè)小數(shù)四舍五入取整，rand表示隨機(jī)產(chǎn)生的一個(gè)0到1之間的數(shù);然后確定種群的邊界，上邊界ubound和下邊界lbound;變異點(diǎn)的值就變?yōu)?lbound+rand×(ubound-lbound)，它表示種群上下邊界范圍內(nèi)的一個(gè)隨機(jī)值。如此依次對(duì)競(jìng)爭(zhēng)種群的每個(gè)個(gè)體進(jìn)行同樣的變異操作。

步驟8:適應(yīng)度值。選擇、交叉、變異算子操作的對(duì)象都是競(jìng)爭(zhēng)種群，3個(gè)算子操作結(jié)束以后，種群中產(chǎn)生了新的個(gè)體，所以對(duì)競(jìng)爭(zhēng)種群進(jìn)行適應(yīng)度值的計(jì)算。從競(jìng)爭(zhēng)種群的第一個(gè)個(gè)體開(kāi)始，進(jìn)行與步驟2一樣的操作，計(jì)算得到的適應(yīng)度值，存放于個(gè)體的nvars+1個(gè)位置，即存放于競(jìng)爭(zhēng)種群的最后一列。

步驟9:精華模型。精華模型也是一種選擇操作，使用精華模型的目的是在對(duì)競(jìng)爭(zhēng)種群進(jìn)行遺傳算子操作以后，需要把競(jìng)爭(zhēng)種群selpool中適應(yīng)度值較高的個(gè)體提取到適應(yīng)性種群pop中，才算一次循環(huán)完成，下一次循環(huán)操作是基于這次循環(huán)完成以后新產(chǎn)生的適應(yīng)性種群pop進(jìn)行的，這樣，算法才能進(jìn)行遺傳進(jìn)化，不斷產(chǎn)生較優(yōu)的新個(gè)體，使迭代向最優(yōu)解逼近。先把競(jìng)爭(zhēng)種群個(gè)體popsize×pg中的前popsize個(gè)個(gè)體賦予給適應(yīng)性種群pop，隨之將競(jìng)爭(zhēng)種群中剩下的個(gè)體分為pg-1組，然后每組中選取適應(yīng)度值最高的個(gè)體替換pop中的最差個(gè)體，適應(yīng)度值最高的個(gè)體的選取方法是:對(duì) pg-1組個(gè)體，逐組操作，對(duì)于每一組，先對(duì)其按照適應(yīng)度值的大小進(jìn)行排序，排序后便能找出適應(yīng)度值最大的個(gè)體，用以替代適應(yīng)性種群pop中的最差值。

以上9個(gè)步驟完成以后，便結(jié)束了算法的一次迭代，第一次迭代結(jié)束以后，算法并沒(méi)有停止，算法停止的條件Τ是:迭代達(dá)到了規(guī)定的最大迭代次數(shù)或者迭代結(jié)果達(dá)到了規(guī)定誤差允許范圍以?xún)?nèi)。從第二代開(kāi)始，種群的每一代計(jì)算是從第5步到第9步，每一次循環(huán)結(jié)束后，適應(yīng)性種群pop會(huì)被更新。

4 程序框圖

圖2為程序框圖。首先對(duì)遺傳算法的初始參數(shù)進(jìn)行設(shè)置，然后產(chǎn)生初始種群，使用V對(duì)初始種群進(jìn)行適應(yīng)度值的計(jì)算;基于種群的適應(yīng)度值進(jìn)行選擇算子Φ，交叉算子Γ，變異算子Ψ的操作，產(chǎn)生新的種群;然后使用V對(duì)新種群進(jìn)行適應(yīng)度值計(jì)算，并基于個(gè)體適應(yīng)度值，選取種群中的較優(yōu)個(gè)體及上一代種群的較優(yōu)個(gè)體組成下一代種群，保證新產(chǎn)生的種群的總體優(yōu)越性。如此重復(fù)，直至達(dá)到初始設(shè)置的算法終止條件Τ。

圖2 遺傳算法反算耳結(jié)構(gòu)彈性模量程序框圖Fig.2 Program chart of genetic-algorithm for inverse derivative of elastic modulus for human ear

5 算例

算例1:變異區(qū)域?yàn)殄N砧關(guān)節(jié)，參考彈性模量為正常人耳砧鐙關(guān)節(jié)彈性模量0.6 MPa，已知其對(duì)應(yīng)的位移u1和u2。遺傳算法程序基本參數(shù)為:個(gè)體長(zhǎng)度nvars=3，位移控制點(diǎn)數(shù) k=2，一個(gè)為鼓膜凸，一個(gè)為鐙骨底板中心，在鼓膜上施加90 dB的壓力，選取在1 000 Hz下兩點(diǎn)的位移。種群大小popsize為20，即隨機(jī)選取的彈性模量為20個(gè)，那么初始種群大小為20×3的矩陣，如表1所示。最大迭代次數(shù)maxgen為500，交叉概率為 pc=0.8，變異概率為 pm=0.01，繁殖因子 pg=4，最佳適應(yīng)度界限值 delt=0.1 ×10－13。

算例1計(jì)算結(jié)果:隨著迭代次數(shù)的增加，種群彈性模量逐漸向靠近參考彈性模量的優(yōu)解方向變化。在迭代440代以后，彈性模量的值接近于真值彈性模量0.6 MPa，本算法的結(jié)束條件是迭代次數(shù)達(dá)到500代或者適應(yīng)度值的絕對(duì)值小于 delta=0.1×10－13。因?yàn)檫m應(yīng)度值是使用種群位移與真值位移的平方差來(lái)計(jì)算的，那么如果適應(yīng)度值小于一個(gè)大于零的很小的數(shù)delta，則種群位移相當(dāng)接近真值位移，便認(rèn)為此時(shí)種群位移對(duì)應(yīng)的彈性模量為最優(yōu)解。在迭代440代以后，種群趨于穩(wěn)定，適應(yīng)度值的絕對(duì)值等于4.627 1×10(15＜delta，算法停止，此時(shí)得到的種群最優(yōu)個(gè)體彈性模量為0.624 MPa。

表1 算例1遺傳算法初始種群Tab.1 Initial population of example1

算例2:變異區(qū)域和算例1相同，參考彈性模量為正常人耳砧鐙關(guān)節(jié)彈性模量0.6 MPa，已知其對(duì)應(yīng)的位移 un(n=1，2，…，8)。遺傳算法基本參數(shù)為:個(gè)體長(zhǎng)度nvars=9，位移控制點(diǎn)數(shù)k=8，取為砧鐙關(guān)節(jié)周?chē)吘壣厦娴墓?jié)點(diǎn)，在鼓膜上施加90 dB的壓力，選取在1 000 Hz下兩點(diǎn)的位移。種群大小popsize為20，即隨機(jī)選取的彈性模量為20個(gè)，那么初始種群大小為20×9，如表2所示、數(shù)量單位與表1相同。最大迭代次數(shù)maxgen為500，交叉概率為pc=0.8，變異概率為 pm=0.01，繁殖因子 pg=4，最佳適應(yīng)度界限值 delt=0.1 ×10－12。

算例2計(jì)算結(jié)果:當(dāng)算法的終止條件設(shè)置為適應(yīng)度值絕對(duì)值達(dá)到界限值delt=0.1×10－12時(shí)，算法迭代了4代便停止了，因?yàn)榻?jīng)過(guò)第4次迭代計(jì)算，種群中已經(jīng)有個(gè)體的適應(yīng)度值達(dá)到了delt以?xún)?nèi)，此時(shí)最佳適應(yīng)個(gè)體的適應(yīng)度值是 －5.49×10－14，其絕對(duì)值在界限delt之內(nèi)，最佳個(gè)體中對(duì)應(yīng)的彈性模量值是0.64 MPa。當(dāng)算法終止的界限值設(shè)置為delt=0.1×10－13時(shí)，迭代次數(shù)在達(dá)到260次以后，種群不再隨著迭代次數(shù)的增加而變化，直到到達(dá)對(duì)大迭代次數(shù)500代，種群最優(yōu)中個(gè)體的彈性模量為0.639 MPa。

算例3:變異區(qū)域?yàn)殄N砧關(guān)節(jié)，參考彈性模量為病變?nèi)硕桤嬯P(guān)節(jié)彈性模量0.8 MPa，已知其對(duì)應(yīng)的位移u1和u2。遺傳算法程序基本參數(shù)為:個(gè)體長(zhǎng)度nvars=3，位移控制點(diǎn)數(shù) k=2，一個(gè)為鼓膜凸位移，一個(gè)為鐙骨底板中心位移。種群大小popsize為20，如表3所示。最大迭代次數(shù)maxgen為500，交叉概率為pc=0.8，變異概率為pm=0.01，繁殖因子 pg=4，最佳適應(yīng)度界限值 delt=0.1 ×10－13。

表2 算例2遺傳算法初始種群Tab.2 Initial population of example2

表3 算例3遺傳算法初始種群Tab.3 Initial population of example 3

算例3計(jì)算結(jié)果:算法迭代450代以后，種群中最優(yōu)個(gè)體的適應(yīng)度值小于初始設(shè)定的 delt，種群個(gè)體區(qū)域穩(wěn)定，得到彈性模量0.784 MPa，接近于參考彈性模量 0.8 MPa。

6 分析討論

算例1中，從初始離散的種群，經(jīng)過(guò)遺傳迭代440代以后便找到一個(gè)解0.624 MPa，真值彈性模量為0.6 MPa。所采用的初始種群的上限為48 MPa，下限為0.01 MPa，那么在區(qū)域 U=0.01～48的范圍內(nèi)求解，得到最優(yōu)解0.6 MPa的相對(duì)誤差為:δ1=(0.624－0.6)÷(48－0.01)=0.05%

由于中耳有復(fù)雜的力學(xué)機(jī)理，而文中遺傳算法計(jì)算反問(wèn)題的程序是可以排除力學(xué)機(jī)理等的影響的，為了驗(yàn)證本程序不受力學(xué)機(jī)理的影響，在算例2中選取的控制位移點(diǎn)是待求砧鐙關(guān)節(jié)周?chē)墓?jié)點(diǎn)，以排除聽(tīng)骨鏈力學(xué)的影響。迭代計(jì)算求得彈性模量0.639 MPa與參考彈性模量0.6 MPa也非常接近，相對(duì)誤差為:

為了更進(jìn)一步驗(yàn)證程序的有效性及穩(wěn)定性，算例3中取參考彈性模量為病變的砧鐙關(guān)節(jié)彈性模量0.8 MPa，采用與算例1中一樣的控制位移及初始種群，設(shè)定特定的適應(yīng)度函數(shù)，控制種群進(jìn)化的方向，對(duì)種群進(jìn)行迭代計(jì)算，迭代450代以后得到彈性模量值為0.784 MPa，與彈性模量0.8 MPa的相對(duì)誤差為

3個(gè)算例的計(jì)算結(jié)果與真實(shí)結(jié)果的相對(duì)誤差都非常小，說(shuō)明本算法是可行的，算例1和算例2的控制位移的不同說(shuō)明了算法不受結(jié)構(gòu)力學(xué)機(jī)理的影響，算例1和算例3參考彈性模量的不同說(shuō)明了算法的穩(wěn)定性及有效性。根據(jù)本研究的結(jié)果，學(xué)者在對(duì)臨床病變耳進(jìn)行研究時(shí)，可先通過(guò)所提出算法求解出結(jié)構(gòu)的力學(xué)參數(shù)，然后通過(guò)力學(xué)原理對(duì)結(jié)構(gòu)進(jìn)行探索，為臨床醫(yī)學(xué)提供參考。本算法具有重現(xiàn)性，但是其穩(wěn)定性不是非常好，以后的研究中需對(duì)算法進(jìn)行改進(jìn)，提高算法的穩(wěn)定性。

7 結(jié)論

1)使用遺傳算法求解耳結(jié)構(gòu)的彈性模量是可行的，所需要的只是不同的彈性模量對(duì)應(yīng)的結(jié)構(gòu)的位移響應(yīng)。對(duì)于其它的力學(xué)參數(shù)，只要能用正問(wèn)題求解出相應(yīng)的位移，便可使用遺傳算法通過(guò)位移求解出相應(yīng)的參數(shù)。

2)用遺傳算法反求結(jié)構(gòu)的材料屬性，不受結(jié)構(gòu)力學(xué)性能的影響。算例1中所取的位移是鼓膜凸位移和鐙骨足板的位移，兩者之間是一個(gè)完整的聽(tīng)骨鏈，具有復(fù)雜的力學(xué)性能;算例2中所取的位移是錘砧關(guān)節(jié)周?chē)狞c(diǎn)，沒(méi)有復(fù)雜的力學(xué)性能，但是兩個(gè)算例的結(jié)果非常接近，即不受聽(tīng)骨鏈力學(xué)性能影響。

3)本研究為無(wú)法通過(guò)實(shí)驗(yàn)測(cè)試的病變活體結(jié)構(gòu)力學(xué)參數(shù)的反演提供了簡(jiǎn)捷有效的方法。

［1］ Pozrikidis C. Boundary-integral modeling of cochlear hydrodynamics［J］.Journal of Fluids and Structures，2008，24(3):336－365.

［2］ 姚文娟，李武，黃新生，等.耳膜振動(dòng)方程的建立與求解［J］.振動(dòng)與沖擊，2008，27(3):465 －472.

［3］ 姚文娟，李武，李曉青，檢驗(yàn)人工聽(tīng)骨力學(xué)性質(zhì)的解析方法［J］.力學(xué)學(xué)報(bào)，2009，41(2):216 －220.

［4］ Lee GF，Chen PR，Lee WJ，et al.Computer aided threedimensional reconstruction and modeling of middle ear biomechanics by high-resolution computed tomography and finite element analysis［J］.Biomedical Engineering，2006，18(5):214－221.

［5］ Gan RZ，F(xiàn)eng B，Sun Q，Three-Dimensional Finite Element Modeling of Human Ear for Sound Transmission［J］.Annals of Biomedical Engineering，2004，32(6):847 －859.

［6］ Gan RZ， ReevesBP， WangXuelin. Modelingofsound transmission from ear canal to cochlea［J］.Annals of Biomedical Engineering，2007，35(12):2180 －2195.

［7］ Tange RA，Grolman W，An analysis of the air-bone gap closure obtained by a crimping and a non-crimping titanium stapes prosthesis in otosclerosis［J］.Auris Nasus Larynx，2008，35(2):181－184.

［8］ Yao Wenjuan， Huang Xinsheng， Fu Lijie， Transmitting Vibration of Artificial Ossicle［J］.International Journal of Nonlinear Sciences and Numerical Simulation，2008，9(2):131－139.

［9］ 姚文娟，李武，付黎杰，等.中耳結(jié)構(gòu)數(shù)值模擬及傳導(dǎo)振動(dòng)分析［J］.系統(tǒng)仿真學(xué)報(bào)，2009，21(3):651 －654.

［10］ 姚文娟，李曉青，李武，等.中耳病變及人工鐙骨形體研究［J］.醫(yī)用生物力學(xué)，2009，(24)2:118 －122.

［11］ 劉迎曦，李生，孫秀珍.人耳傳聲數(shù)值模型［J］.力學(xué)學(xué)報(bào)，2008，40(1):107 －113.

［12］ 劉迎曦，李生，孫秀珍.人耳鼓膜病變數(shù)值分析［J］.醫(yī)用生物力學(xué)，2008，23(4):275－278.

［13］ Herrmann G，Liebowitz H.Mechanics of bone fractures［M］//Fracture:An Advanced Treatise.New York:Academic Press，1972:772－840.

［14］ Békésy GV.Experiments in hearing［M］.New York:McGraw-Hill，1960.

［15］ Kurokawa H，Richard L.Sound pressure gain produced by the human middle ear［J］.Otolaryngol Head Neck Surg，1995，113(4):349－355.

［16］ Tao C，Rong Z.Mechanical properties of stapedial tendon in human middle ear［J］.Journal of Biomechanical Engineering，2007，129:913－918.

［17］ 雷英杰，張善文，李續(xù)武，等.MATLAB遺傳算法工具箱及其應(yīng)用［M］.陜西:西安電子科技大學(xué)出版社，2005.

［18］ 徐波，陳曉江.拋物型方程反問(wèn)題的遺傳算法［J］.空軍雷達(dá)學(xué)院學(xué)報(bào)，2008，22(04):293 －295.

［19］ 蔡傳寶，湯文成.遺傳算法的生物組織彈性模量反求研究［J］.應(yīng)用科學(xué)學(xué)報(bào)，2008，26(03):295 －300.

［20］ 曳永芳，杜永清，行小帥.一種抑制早熟收斂的改進(jìn)遺傳算法［J］.山西師范大學(xué)學(xué)報(bào)，2010，24(02):24 －28.

［21］ 田延碩，劉曉云.一種提高局部搜索能力的混合遺傳算法［J］.電子科技大學(xué)學(xué)報(bào)，2006，35(02):232 －235.

［22］ 黃衛(wèi)華，許小勇，范建坤.實(shí)數(shù)編碼遺傳算法中常用交叉變異算子的Matlab實(shí)現(xiàn)及應(yīng)用［J］.廣西輕工業(yè)，2007，(1):77－79.

［23］ 梁昔明，秦浩宇.一種求解約束優(yōu)化問(wèn)題的遺傳算法［J］.計(jì)算機(jī)工程，2010，36(14):147 －149.

Inverse Derivative of Elastic Modulus For Human Ear Using Genetic-Algorithm

CHENG Xiang-Dong YAO Wen-Juan*
(Department of Civil Engineering，Shanghai University，Shanghai 200072，China)

The objective of this paper is to calculate the elastic modulus of the ear structure by the known displacement based on inverse problem theory of mechanics.The initial population of the genetic algorithm was generated in random.Genetic iterative calculation was carried out on the initial population by a self－compiling program using MATLAB language.The target function was the mean-square deviation between the known displacement and the displacement of the population.We completed the following work by the calculation example:The elastic modulus of the normal incudostapedial joint was obtained by two controlled displacements of umbomembranetympani and stapesfootplateofhuman ear.Theelasticmodulusofthenormal incudostapedial joint was obtained by eight controlled displacements around incudostapedial joint.The elastic modulus of the diseased incudostapedial joint was obtained by two controlled displacements of umbo membrane tympani and stapes footplate of human ear.Results showed that it was feasible to calculate the elastic modulus of the ear structure by the inverse problem method based on genetic algorithm，which has the advantage of stability and independence of structure mechanical property.The relative errors were 0.05% ，0.2% ，and 0.03%respectively，which supplied valid mechanical parameters for diseased human ear in clinic.

ear structure;inverse problem;elastic modulus;genetic algorithm

R318

A

0258-8021(2012)04-0487-07

10.3969/j.issn.0258-8021.2012.04.000

2012-04-02，錄用日期:2012-06-08

國(guó)家自然科學(xué)基金(11072143);上海市科委基礎(chǔ)研究重點(diǎn)項(xiàng)目(08jc1404700)

* 通信作者。E-mail:wenjuan@mail.shu.edu.cn