一種多策略改進遺傳算法的醫學圖像配準方法

秦澤青，葉志偉，王澤松，徐 川，曹 羽

1(湖北工業大學 計算機學院，武漢 430068) 2(湖北大學 計算機與信息工程學院，武漢 430062) 3(湖南大學 信息科學與工程學院，長沙 410012) E-mail：weizhiye121@163.com

1 引 言

醫學影像可以分為解剖成像和功能成像[1]，前者分辨率較高，能夠提供解剖形態信息，后者分辨率較低，但可以很好地提供功能代謝等信息.醫學圖像配準是指通過調整幾何變換的參數，從而使不同類型醫學影像的共同部分達到空間位置上一致的過程[2，3]，可以更全面地提供病灶的信息，對現代醫療有著重大的意義[4].本質上，醫學圖像配準是一種參數最優化問題[5]，可以使用優化算法進行處理.

傳統的優化方法已被應用于醫學圖像配準問題，然而，由于圖像離散化，噪聲等因素的影響，優化過程容易陷入局部最優解.近年來，基于進化算法和其他元啟發式方法的醫學圖像配準方法展現出了良好的應用潛力[6]，配準精確度高，健壯性強，在很大程度上克服了傳統方法的不足[7].

遺傳算法(Genetic Algorithm，GA)是經典的進化算法，基于進化論模擬了自然界中的染色體交叉，堿基對變異等現象，有著良好的全局搜索能力，在許多領域得到了廣泛的應用[8].但是，基于遺傳算法的醫學圖像配準方法仍然存在著局部搜索能力較差且容易早熟收斂等不足[9].

為了改善上述問題，本研究從隨機機制，種群多樣性以及遺傳算子3個方面出發，提出多種策略對遺傳算法進行改進：1)設計一種基于混沌系統與頻度記憶的隨機策略，利用Logistic混沌映射產生隨機數序列并結合頻度記憶對隨機個體的生成進行了干預，實現了隨機數和隨機個體的可再現性與可控性；2)提出一種動態的基因多樣性控制策略，通過控制進化過程中個體間基因的差異，保證了種群的多樣性；3)提出一種基于雜種優勢的進化策略，模擬三系法雜交水稻協同進化育種模式設計新的遺傳算子.最后結合多分辨率策略，利用改進的遺傳算法搜索多模態醫學圖像配準的最優參數并進行實驗仿真.結果表明本文提出的策略有效地提升了遺傳算法的搜索能力并降低了早熟收斂的概率.

2 相關理論

2.1 圖像配準

圖像配準指的是通過特定的幾何變換，從而使兩個或多個圖像的共同部分達到空間位置上一致的過程[2，3].在這個過程中，作為參照模板固定的圖像叫做參考圖像，被幾何變換調整的圖像叫做移動圖像.圖像配準可以分為基于灰度的配準方式以及基于特征的配準方式.從效率來看，前者只需要提取圖像的特征點、邊緣線等突出信息，時耗較小，而后者需要提取圖像的所有信息，效率較低；從復雜程度來看，前者需要對圖像進行特征提取等預處理，較為繁瑣，而后者無需任何預處理，十分簡潔；從應用范圍來看，基于特征的圖像配準只適用于有著明顯特征的圖片，而基于灰度的圖像配準適用于所有類型的圖像[2].本文采用基于灰度的圖像配準方式，包含5個步驟：初始化參數過程、幾何變換過程、評價過程、優化過程、輸出過程.配準流程如圖1所示.

圖1 基于灰度的圖像配準流程圖Fig.1 Intensity-based image registration flow chart

圖像配準作為一種有效的信息整合方法，已經在醫學圖像引導手術(image-guided surgery)[4]等方面得到了廣泛的應用，可以被視作參數優化問題.基于傳統優化算法的圖像配準方法已取得了一定的成效，例如鮑威爾方法(Powell′s method)，下降單純型法(downhill simplex method)，高斯牛頓法(Gauss-Newton)，梯度下降法(gradient descent)等[2]，但由于圖像的離散化，噪聲等負面因素的影響，這些方法非常容易陷入局部最優解[3].

近年來，進化計算和元啟發式算法在圖像配準問題上受到了廣泛地關注.相較于傳統優化算法，這些方法求解速度快、尋優能力強、計算成本低.Valsecchi等人[10]提出了一種改進后的散射搜索方法(Scatter Search*)，采用雙層設計(2-tier design)同時考慮了最優子集的質量和多樣性，并按照圖像配準的特征對散射搜索的基本模板進行改進，最后將改進后的散射搜索法與多分辨率進行結合，表現出較好的性能.Bermejo等人[11]首次將珊瑚礁優化算法(Coral Reef Optimization)用于醫學圖像配準問題，基于5種算子結合的進化模式為珊瑚設計了附著層，在基于灰度和特征的醫學圖像配準實驗中都表現出了較高的配準精度和較強的健壯性.Chen等人[12]將生物地理學優化算法(Biogeography-based Optimization)首次用于基于灰度的多模醫學圖像配準問題，并結合精英學習策略對其做出了改進，提升了所得解集的質量并提高了收斂的效率.

上述算法都在醫學圖像配準問題上表現出了良好的尋優能力.然而，這些方法一定程度上還存在容易陷入局部最優解的問題，收斂性也沒有得到充分證明.作為經典的最優化算法，遺傳算法具有著簡單高效的特性，其收斂性也得到了證明，本文擬通過多種策略對遺傳算法做出改進，使其更加適用于醫學圖像配準問題.

2.2 幾何變換

幾何變換分為剛體變換和非剛體變換兩類，后者又包括仿射變換以及透視變換等.利用幾何變換可以對圖像進行旋轉、平移、縮放、剪切等操作.為了減少計算代價，本文采用自由度最低的三維剛體變換.此變換包含3個旋轉參數rx、ry、rz和3個平移參數tx、ty、tz.

2.3 遺傳算法

遺傳算法是一種基于生物界自然選擇和遺傳機制構建的隨機搜索算法，由美國的Holland教授提出.這種算法來源于達爾文的進化論、魏茨曼的物種選擇學說和孟德爾的群體遺傳學說，模擬了自然界生物繁育過程中的染色體交叉和基因變異等現象，有著全局優化性、梯度信息不依賴性、簡單易施性等特點，在很多領域有著廣泛的應用[8].作為一種自適應全局優化算法，遺傳算法使用二進制編碼，其基本執行過程為：

步驟1.確定種群規模N，交叉概率Pc以及變異概率Pm，設定停止準則；

步驟2.按照某一隨機機制初始化N個個體的基因；

步驟3.計算各個體的適應度，按照選擇算子選擇親代；

步驟4.按照交叉算子遵循概率Pc執行交叉；

步驟5.按照變異算子執行概率Pm執行變異；

步驟6.計算各個體的適應度，按照選擇算子進行淘汰；

步驟7.終止檢驗.若已經滿足停止準則，則輸出目前種群中適應度最大的個體基因作為最優解；否則轉步驟3.

針對遺傳算法的研究主要包含3個部分：隨機機制(Random Mechanisms)、遺傳算子(Genetic Operators)、種群多樣性(Diversity of Population).其中遺傳算子包括了選擇、交叉和變異3種基本形式.

在醫學圖像配準問題中，遺傳算法存在著編碼方式受限、局部搜索能力較差、容易早熟收斂等不足[9].圖2表明了這些不足與本文主要研究內容之間的關系.

圖2 本文主要的研究內容Fig.2 Topics of this work

3 提出的方法

為了讓遺傳算法更加適用于基于灰度的醫學圖像配準問題，本研究從隨機機制、種群多樣性、遺傳算子3個研究方向對標準的遺傳算法提出多種策略進行改進并用于醫學圖像配準.

3.1 隨機機制

遺傳算法是一種隨機搜索算法，隨機數和隨機個體的產生機制在很大程度上決定了搜索的精度和速度.在標準遺傳算法的實現中，大多使用類如randomize()或rand()等隨機函數，存在著效率低，不可再現等問題.本文將混沌系統和頻度記憶結合，實現了隨機數和隨機個體的可再現性與可控性.

3.1.1 混沌系統

混沌系統是一種復雜的非線性動力學系統，產生的混沌序列具有著遍歷性、類隨機性、確定可再生性等特點[13].混沌系統被廣泛應用于保密通信、數據安全等領域中并取得了良好的效果.本文采用混沌系統中簡單的一維Logistic混沌映射來生成隨機數，如式(1)所示：

xn+1=r·xn·(1-xn)

(1)

其中：系統控制參數r∈(0，4]，x∈(0，1).

Logistic映射的動態行為與系統控制參數r密切相關，不同的r可以使Logistic映射處于周期性或混沌狀態.研究[14]表明，當r∈(3.5699456，4]時，Logistic映射進入混沌狀態，產生的混沌序列{xn|n=0，1，2，3，…}呈現出非周期性，具有良好的隨機分布特征且不收斂.在r=4時，Logistic映射的輸入輸出都在(0，1)區間內，混沌軌跡最為明顯，具有極好的遍歷特性.本文采用r=4時的Logistic混沌映射來產生隨機序列，對應的映射表達式如式(2)所示：

xn+1=4·xn·(1-xn)

(2)

3.1.2 頻度記憶

在禁忌搜索的實現中，頻度記憶(frequency-based memory)被廣泛應用于長期記憶策略[15]，這種方法可以自適應地調整隨機個體生成的概率.本文利用頻度記憶對遺傳算法中隨機個體的生成過程進行干預，使其分布更加均勻，從而提升隨機個體的可控性.基于頻度記憶的隨機個體生成策略分為3步：

步驟1.將隨機個體參數的搜索空間S從小到大劃分為4個大小相等的區間S={s1，s2，s3，s4}，并賦予每個區間一個頻率計數器freqi，初始值記為0；

步驟2.隨機生成一個屬于[0，4]的整數來選擇具體的區間.隨機整數的生成遵循Logistic混沌映射產生的混沌序列，但是其生成概率與freqi的值成反比，生成后對應的freqi值加1；

步驟3.在被選擇的區間中利用Logistic混沌映射生成隨機個體的參數值.

3.2 種群多樣性

早熟收斂是指超級個體的出現使得種群一致向其靠攏的現象.由于擁有相對良好的性狀，超級個體的基因會很快在種群中占有絕對的比例，從而導致多樣性迅速降低，整個群體喪失進化能力.為了降低這一現象發生的概率，本研究提出了一種動態的基因多樣性控制器.設個體基因的維度為n，我們要求進化過程中子代基因x與親代基因y的距離不得小于多樣性閾值Diversity_Threshold，如式(3)所示：

(3)

由于在進化的過程中群體的需求是變化的，所以多樣性控制器的閾值應該是動態的.自然界中物種的進化現象表明，群體進化初期個體之間性狀差異較大，而中后期為了適應于特定的環境，群體的性狀趨于一致[16].所以進化的早期不需要多樣性控制或者說只需要很弱的控制，過早的控制可能會丟失寶貴的臨近解；然而在進化的中后期則需要加大控制的力度來避免早熟收斂.基于上述思想，多樣性控制器的閾值Diversity_Threshold是隨著種群的進化逐漸增大的.假設種群的進化次數為n，最大差異限度為Diversity_Max.Diversity_Threshold是自變量處于[1，n]，函數值處于[0，Diversity_Max]上的線性函數.當種群進化次數為1時Diversity_Threshold=0，當進化次數為n時Diversity_Threshold=Diversity_Max.

3.3 遺傳算子

為了進一步學習自然界中優秀的繁育機制并擴大基因差異化對種群進化的正向影響，受到雜種優勢的啟發，本研究對標準的遺傳算子進行了改進.雜種優勢是指雜合體在一種或多種性狀上優于親本的一種生物學現象，在農業中得到了廣泛的應用[17].中國的三系雜交水稻[18]是雜種優勢的典型應用案例，在很大程度上解決了糧食短缺的問題.三系是指保持系、恢復系以及不育系，在繁育的過程中水稻種子被人為地劃分進這三系.位于不同系的水稻種子通過雜交、自交等方式相互學習優良的性狀并繁育后代.相較于傳統遺傳算法的雙親學習模式，中國的三系雜交水稻這種三系協同進化的模式更好地發揮群體智能的優勢并在很大程度上保證了群體基因的差異.本研究模擬了中國三系雜交水稻的育種機制式，改進了標準遺傳算法的交叉、變異以及選擇算子.具體的步驟如下：

步驟1.利用隨機機制初始化n個遺傳個體X的基因，并按照適應度函數值從高到低將遺傳個體等量地劃分入保持系{X1，X2，…Xm}、恢復系{Xm+1，Xm+2，…X2m}和不育系{X2m+1，X2m+2，…X3m}，每系中的個體數量m=?n/3」；

步驟2.從保持系和不育系中分別選擇親代進行雜交，對于保持系的個體Xi選擇不育系的個體Xn-i與其配對.相較于傳統遺傳算法的單點交叉模式，本研究采用BLX-α交叉算子來對親代基因進行學習.設被選中的親代為x、y，對于他們基因的的每一個維度i有xi，yi∈[ai，bi]，其中ai是維度i取值范圍的下界，bi是取值范圍的上界.子代基因zi是在區間[LEFTi，RIGHTi]上服從均勻分布的隨機數.LEFTi是區間下界，見式(4).RIGHTi是區間上界，見式(5)：

LEFTi=max{ai，min(xi，yi)-αdi}

(4)

RIGHTi=min{max(xi，yi)+αdi，bi}

(5)

其中：親代基因第i維的距離di=|xi-yi|，α為雜交率.

如果子代的適應度優于不育系親代，則后者的基因將會被直接替換；

步驟3.恢復系的個體充分學習保持系中表現最好個體的基因.在這一步，基于隨機機制選擇兩個恢復系的個體，并選擇保持系中適應度最好的個體.這些個體將進行交叉操作完成自交.設來自保持系的個體為Zbest，恢復系中的個體為x、y，且對于基因的的每一個維度i有xi，yi∈[ai，bi]，其中ai是維度i取值范圍的下界，bi是取值范圍的上界.子代基因的第i維從Fgene和Sgene兩個區間中的一個遵循隨機機制生成，Fgene和Sgene的表達式如式(6)、式(7)所示：

Fgene=[max{ai，Zbesti-d1i·β}，min{Zbesti+d1i·β，bi}](6)

Sgene=[max{ai，Zbesti-d2i·β}，min{Zbesti+d2i·β，bi}](7)

其中：d1i=|Zbesti-xi|，d2i=|Zbesti-yi|，β是自交率.

恢復系中的個體變異的概率為Pm，假設會發生基因變異的個體為x，根據高斯分布概率密度函數構造高斯變異函數，具體表達式如式(8)所示：

(8)

其中：Best是當前的最佳解，σ是高斯分布的標準偏差.由于三維剛體變換有6個參數，i是在區間[1，6]中隨機生成的整數.

受到高斯突變的干擾，更新后的個體可以逃脫最優解的引力，從而跳出局部最優值.當前解的第i個位置的值如式(9)所示：

xi=xi·GMF(x，i)

(9)

3.4 具體應用

3.4.1 遺傳個體編碼

本研究采用實數編碼的策略來代替二進制編碼.對三維剛體變換的參數進行編碼，組成遺傳個體的基因.設種群的大小為N，編碼后的遺傳個體如式(10)所示：

Xi=(xi，1，xi，2，xi，3，xi，4，xi，5，xi，6)，i=1，2，3，…，N

(10)

其中xi，j代表三維剛體變換的第j個參數.三維剛體變換共有6個參數，分別是3個旋轉參數和3個平移參數.參照文獻[10，11]，旋轉參數rx、ry、rz的取值范圍是[-90，90]，平移參數tx、ty、tz的搜索空間為[-30，30].

3.4.2 適應度函數

本研究利用歸一化互信息度量遺傳個體的適應度.互信息可以看作是一個隨機變量中包含的關于另一個隨機變量的信息量，被廣泛用于測量醫學圖像配準效果[3].設圖像X，Y的聯合分布概率為p(x，y)，p(x)和p(y)是他們的邊緣分布概率.互信息I(x，y)是聯合分布p(x，y)與乘積分布概率p(x)p(y)的相對熵，如式(11)所示：

=H(X)+H(Y)-H(X，Y)

(11)

其中：H(X)，H(Y)分別是X和Y的信息熵，H(X，Y)是他們的聯合熵.

適應度函數的計算如式(12)所示，配準效果越好時歸一化互信息值NMI越趨近于1，反之則越趨近于0.

(12)

3.4.3 算法流程

本文算法的執行流程如圖3所示.

圖3 改進后的遺傳算法的流程圖Fig.3 Flow chart of improved genetic algorithm

4 實 驗

本節通過一系列實驗來探究改進后的遺傳算法(以下標記為GA*)在醫學圖像配準問題上的搜索能力以及早熟收斂的狀況.

4.1 數據介紹和多分辨率策略

本研究選用來自Retrospective Image Registration Evaluation Project (RIRE)(1)https://www.insight-journal.org/rire/index.php提供的公開數據集中6號患者腦部的CT成像與MR-T1成像、MR-T2成像以及MR-PD成像進行配準.其中CT成像的大小為512×512×30voxels，MR成像的大小為256×256×26voxels.圖4(a)、(d)、(g)分別是MR-PD、MR-T1、以及MR-T2成像，作為配準中的參考圖像(Reference image)；圖4 (b)、(e)、(h)是CT成像，作為配準中的浮動圖像(Moving image).圖4(c)、(f)、(i)顯示了這3組圖像初步對齊的效果.

圖4 數據集可視化Fig.4 Visualization of datasets

本研究采用多分辨策略降低優化過程陷入局部最優解的概率.在進行配準之前首先對原圖像進行下采樣和高斯平滑處理，從而構建雙層圖像金字塔.金字塔的第1層是低分辨率的圖像，第2層是未經處理的原始圖像.本研究首先對第1層的低信息圖像進行配準，記錄最優解，并將最優解用于第2層配準參數的初始化.為了防止算法陷入局部最優解，優化算法需要在低分辨率圖像上運行多次.

4.2 實驗設計與參數設置

研究[6，7]表明，相較于傳統的優化算法，進化計算和元啟發式優化算法在圖像配準問題上有著更佳的效果.本文選取5種元啟發式優化算法與GA*進行對比，分別為改進后的散射搜索算法(SS*)[10]、珊瑚礁優化算法 (CRO)和帶有附著層的珊瑚礁優化算法(CRO-SL)[11]、生物地理學優化算法(BBO)和基于精英學習策略改進的生物地理學優化算法(BBO-EL)[12].

本研究的配準流程如圖5所示.

圖5 配準執行流程Fig.5 Image registration process

所有的優化算法的參數設定嚴格與相關研究[10-12]保持一致.為了保證公平性，參照文獻[11]，在圖像金子塔第1層上的配準時間限制為180秒，在第2層的時間限制為280秒.第1層最優解集的大小由算法具體執行的效率決定.本研究將多樣性控制閾值Diversity_Max設置為0.015.在基于灰度的圖像配準實驗中，文獻[11]將所有優化算法的交叉概率和混合因子分別設置為0.5和0.3，所以GA*的雜交率α=0.5，自交率β=0.3.參照文獻[9]，GA*的變異率Pm=0.08.

4.3 實驗結果與分析

4.3.1 固定時間下的配準效果研究

本研究在給定的時間限制下對所給圖像進行了配準.將MR-T1成像作為參考圖像，CT成像作為浮動圖像，采用GA*、SS*、CRO、CRO-SL、BBO以及BBO-EL這6種優化算法各自獨立地進行了20次配準，并將歸一化互信息值(NMI)作為衡量指標，實驗結果如表1所示.

表1 多模態醫學圖像配準20次結果Table 1 Results of multi-model image registration

實驗結果的描述性統計如表2所示.

表2 配準結果的描述性統計Table 2 Statistical description of registration results

從表1可以看出，相較于其他算法，GA*在所給多模醫學圖像配準的場景上總體表現最佳，在20次實驗中有11次精度最高，占比55%.對表2的數據進行分析，從20次配準的NMI平均值來看，GA*達到了0.1839，而其余算法都在0.183以下.這說明GA*在所給數據集上的搜索能力較強，在給定的時間內尋優精度最高.除此之外，本研究參照所有算法20次配準的最差值對是否錯誤配準進行評判.去除CRO的極端值，其余算法最差值的平均值為0.1767，本研究將配準結果低于0.1767的視為誤配.GA*的誤配率為0%，和BBO-EL同時達到了最低，遠小于CRO的50%.這說明GA*早熟收斂的概率相對較低.

表2的數據顯示，GA*標準差為0.0017，明顯低于大多數算法，僅高于BBO_EL的標準差0.0001，表現出了很強的魯棒性.在20次迭代中，BBO_EL雖最高值低于GA*，但其最差值略高于GA*，解質量的穩定性略優.

4.3.2 固定代數下的配準效果研究

為了探究不同算法進化過程之間的差異，本研究進一步在固定的迭代次數下利用不同優化算法對所給數據集進行了配準實驗.

4.3.1的實驗已經表明，CRO-SL的配準效果要優于CRO且BBO-EL的配準效果要優于BBO，本實驗主要針對SS*、CRO-SL、BBO-EL以及GA*展開進一步分析.實驗嚴格按照圖5所示的配準流程執行，在圖像金子塔第1層上的配準時間限制為180秒，在第2層上取消了時間限制，迭代次數固定為100代，其余參數保持不變.本實驗在所給數據集的3個多模配準場景上展開，4種算法隨著迭代次數增加在原圖像上的配準效果如圖6所示.

圖6 固定代數下的配準Fig.6 Image registration in fixed iteration

由于多分辨率策略的存在，不同的算法在圖像金子塔第1層上表現的差異導致了在原圖像上配準初始值的不同.從圖6所示的3個場景可以看出，GA*的初始值均為最高，甚至高于CRO-SL迭代100次后的最終配準值.這說明多分辨率策略在很大程度上提高了本文算法的搜索能力，GA*在低分辨率圖像上的表現最好，局部搜索能力較強，為高分辨率圖像上的配準打下了很好的基礎.

除此之外，GA*的全局收斂能力佳，可以在較少的迭代次數內接近或達到全局最優的配準效果，圖6(a)和圖6(c)中GA*在第10代的配準效果就已經超過了其他算法在迭代100次后所能達到的效果.最后，GA*的搜索不易陷入局部最優解，圖6(a)中GA*在第10代已經達到了很好的配準效果，在第30代卻仍然能夠出現上升趨勢，突破局部最優解.

表3給出了4種算法在3個場景下最終的配準參數、配準結果的NMI以及總耗時.從最終的配準結果來看，GA*在3個場景上的總體表現最好，時間性能略優.從總耗時來看，CRO-SL在最短的時間內完成了3個場景上的一百次迭代，其次是BBO-EL，GA*每一代的搜索的時間相對較長，有進一步提升的空間.

表3 4種算法的配準結果Table 3 Image registration results of four algorithms

以上實驗表明，本研究提出的基于混沌映射和頻度記憶的隨機策略以及基于雜種優勢的進化策略增強了遺傳算法的搜索能力；基因多樣性控制策略保證了種群的多樣性，從而降低了遺傳算法陷入局部最優解和早熟收斂的概率.相較于SS*、CRO、CRO-SL、BBO以及BBO-EL這5種優化算法，GA*在所給數據集上總體表現最佳.在給定時間限制的前提下GA*的性能有著較好的優勢，而在固定進化代數的情況下GA*在早期低分辨率圖像上的搜索效果較好，但在后期進化幅度不大.這說明相較于所對比的優化算法，GA*的主要優勢在于可以在較短的時間內達到或接近全局最優解，適用于時間優先的應用場景.

總的來說，本文提出的基于多策略改進的遺傳算法誤配率低，搜索能力強，不易早熟收斂且有著很好的魯棒性.

5 結 論

本文從隨機機制、種群多樣性以及遺傳算子3個方面出發，提出多種策略來解決遺傳算法在醫學圖像配準問題上局部搜索能力較差，容易早熟收斂等不足.并將多分辨率策略與改進后的遺傳算法相結合，以NMI作為配準效果的測度函數，實現了醫學圖像的精確配準.實驗結果表明，與SS*、CRO、CRO-SL、BBO以及BBO-EL相比，本文的方法可以較好地搜索到質量最高的解且不易早熟收斂，在很大程度上克服了標準遺傳算法的不足，應用到醫學圖像配準問題上配準精度高、誤配率低，有著良好的應用前景.未來我們將嘗試把改進后的遺傳算法與其他優化算法相結合，在不同的分辨率上發揮各自的優勢，進一步提高醫學圖像配準的精度和效率.