黃志業(yè) 陳武凡 周凌宏 徐子海 陳超敏*（南方醫(yī)科大學生物醫(yī)學工程學院，廣州 5055）

2（解放軍303醫(yī)院放射治療中心，南寧 530021）

引言

在胸腹部腫瘤的放射治療過程中，由于心跳、呼吸運動等的影響，使得腫瘤組織總是隨著它們的運動而產(chǎn)生位置的偏移，這是提高胸腹部腫瘤放射治療療效的重要障礙。為了解決呼吸運動的影響，在傳統(tǒng)的放射治療和IMRT中通常采用腫瘤靶區(qū)射野擴邊［1］、呼吸保持［2］、呼吸門控［3－4］等技術。這些技術的產(chǎn)生，對胸腹部腫瘤的放射治療有著積極的意義。但是，人們很快發(fā)現(xiàn)，這些技術并不能很好地幫助人們達到期望的目的。擴邊增加了腫瘤周圍健康組織的照射體積，對周圍的健康組織造成額外的損傷；呼吸保持方法要求病人在治療過程中保持某一呼吸狀態(tài)，對于那些有肺功能缺損的病人來說耐受性很差；門控技術要求放射線和病人呼吸周期同步，只在一個很小的呼吸期內(nèi)打開射線進行照射，因此增加了治療的時間。所以，人們正在積極地探索新的解決方法。目前，實時地跟蹤和預測腫瘤運動軌跡，是研究補償呼吸運動影響的重要方向，因為它可以克服擴邊、呼吸保持和門控等技術的限制。賽博刀（cyberknife）系統(tǒng)就是實時跟蹤技術在進行腫瘤放射治療中的成功應用。

實時跟蹤的最大挑戰(zhàn)是治療系統(tǒng)存在的延時，它必須在射線重新對準腫瘤的時間延時內(nèi)預測出腫瘤的運動軌跡，并通過控制環(huán)把腫瘤的坐標發(fā)送到位置調(diào)整系統(tǒng)，重新調(diào)整射線和腫瘤之間的位置［5］。由于治療系統(tǒng)的延時總是存在的，所以要實時地跟蹤腫瘤的位置，需要提前預測出腫瘤在將來一段時間內(nèi)的運動軌跡，以作為系統(tǒng)延時的補償。本研究采用遺傳算法（genetic algorithm，GA）和 BP（back propagation）神經(jīng)網(wǎng)絡結(jié)合，形成GA-BP神經(jīng)網(wǎng)絡算法，并把呼吸運動信號作為時間序列信號進行呼吸運動的預測。遺傳算法（GA）是模擬自然演化過程搜索全局最優(yōu)解的方法，它的群體性搜索策略使其不易陷入局部最優(yōu)；同時，遺傳算法基本上不用搜索空間的知識或其他輔助信息，而僅用適應度函數(shù)值來評估個體，在此基礎上進行遺傳操作。適應度函數(shù)不僅不受連續(xù)可微的約束，而且其定義域可以任意設定，使得遺傳算法的應用范圍大大擴展，方便與其他算法相結(jié)合［6－7］。BP（Back Propagation）神經(jīng)網(wǎng)絡算法本質(zhì)上為梯度下降法，所以容易陷入局部最優(yōu)，而且BP網(wǎng)絡上每個節(jié)點的權值和閾值都會影響網(wǎng)絡的最后輸出，存在著收斂速度慢、學習過程中常常發(fā)生振蕩等缺點。

鑒于遺傳算法的優(yōu)點和BP網(wǎng)絡存在的缺陷，有學者就提出了將兩者結(jié)合使用的設想，希望利用遺傳算法不易陷入局部最優(yōu)的特性，補償BP神經(jīng)網(wǎng)絡梯度下降易于陷入局部最優(yōu)的缺陷［8-10］。本研究基于這樣的思想，對GA-BP神經(jīng)網(wǎng)絡算法在呼吸運動預測上做了初步的探討。

1 GA-BP網(wǎng)絡預測原理

利用GA-BP網(wǎng)絡預測涉及兩個重要的過程：首先應用遺傳算法（GA），優(yōu)化BP網(wǎng)絡的權值和閾值的初始值；然后使用BP網(wǎng)絡，根據(jù)獲得的優(yōu)化初始值進行網(wǎng)絡訓練。通過利用遺傳算法能夠搜索全局最優(yōu)的優(yōu)點，結(jié)合BP網(wǎng)絡能進行任意維數(shù)輸入輸出映射的特性，達到提高網(wǎng)絡預測精度的目的。

1.1 遺傳算法優(yōu)化BP神經(jīng)網(wǎng)絡的權值和閾值

在遺傳算法中，首先要進行基因表述，即對遺傳算法中的基因進行編碼，因為基因包含了待求問題的所有參數(shù)。在遺傳算法優(yōu)化BP神經(jīng)網(wǎng)絡的權值和閾值中，遺傳算法的基因包括了 BP神經(jīng)網(wǎng)絡各層間的權值和閾值。以3層BP網(wǎng)絡的權值和閾值為例，可以用式子來表示遺傳算法的一個基因，即式中，w和W分別表示神經(jīng)網(wǎng)絡輸入層到隱層和隱層到輸出層的權重，b和B分別表示神經(jīng)網(wǎng)絡隱層和輸出層的閾值。

其次，設計合理的適應度函數(shù)，這是進行遺傳算法選擇運算操作的重要依據(jù)。進化論中的適應度，是表示某一個體對環(huán)境的適應能力，也表示該個體繁殖后代的能力。遺傳算法的適應度函數(shù)也叫評價函數(shù)，是用來判斷群體中個體優(yōu)劣程度的指標，可根據(jù)所求問題的目標函數(shù)來進行評估。適應度函數(shù)的選取應具有規(guī)范性、合理性和通用性，可采用公式來加以定義，有

式中，ei表示第i個個體樣本輸入得到的實際輸出和期望輸出差的平方和。Tkj、okj分別表示第k個樣本在第j個輸出神經(jīng)元的期望輸出和實際輸出，l和p分別表示訓練樣本數(shù)目和輸出神經(jīng)元數(shù)目，fi表示第i個個體的適應度值。

最后，實現(xiàn)遺傳算法的選擇、交叉和變異3個基本遺傳操作。在簡單遺傳算法（SGA）中，選擇的策略采用輪盤賭的方法，雖然每一個個體都可以獲得復制的機會傳到下一代，但并不能很好地體現(xiàn)優(yōu)秀個體的競爭力和遺傳算法優(yōu)勝劣汰的基本準則。采用有條件的最佳保留策略，可以更好地體現(xiàn)優(yōu)秀個體的競爭力，同時也能有效地防止算法的“早熟收斂”。

交叉操作是遺傳算法中產(chǎn)生新個體的主要方式。在遺傳算法中，交叉操作決定著算法的全局搜索性能。在算法搜索的前期，由于大部分個體分散在問題域中，應該強化算法的全局搜索求最優(yōu)解，即在前期給予較大的概率Pc，從群體中選擇個體進行交叉，形成新的個體，有利于算法的全局搜索，但是在算法的后期，假設個體的平均適應度達到了某個值，即認為群體中的大部分個體都聚集在全局最優(yōu)解的附近，這時候如果仍然采用大概率進行交叉操作，則很可能破壞優(yōu)秀的基因結(jié)構，不利于算法的穩(wěn)定性。所以，根據(jù)具體的情況，自適應地改變交叉概率Pc的大小，即

式中，Pc表示交叉概率；fmax和favg分別表示群體個體的最大適應度值和平均適應度值。

變異操作的作用是保持群體中基因的多樣性，偶然（Pm的概率很?。┑淖儺愒谶z傳算法中起著輔助和微調(diào)的作用，遺傳算法中的變異操作決定著算法的局部搜索性能?；谝陨系募僭O，在算法搜索的前期，由于大部分個體分散在問題域中，應該強化算法的全局搜索來求最優(yōu)解，所以只需要很小的Pm對個體進行局部的微調(diào)；在算法的后期，假設個體的平均適應度達到了某個值，群體中的大部分個體都聚集在全局最優(yōu)解的附近，這時候可以加大局部搜索的概率，提高局部搜索的效率。在局部進行個體的調(diào)整，使得在最優(yōu)解附近的個體進一步靠近。所以，根據(jù)具體的情況，自適應地改變變異概率Pm的大小，有式中，Pm表示變異概率，fmax和favg分別表示群體個體的最大適應度值和平均適應度值。

上述遺傳算法根據(jù)特定的網(wǎng)絡結(jié)構編碼，選擇相應的適應度函數(shù)，自適應地改變交叉概率（Pc）和變異概率（Pm）大小，經(jīng)過反復迭代地操作，最后輸出優(yōu)化的BP網(wǎng)絡權值和閾值，并將此輸出作為下一步BP網(wǎng)絡訓練的初始條件。

1.2 BP神經(jīng)網(wǎng)絡結(jié)構和預測

用遺傳算法得到的輸出作為BP神經(jīng)網(wǎng)絡的初始權重和閾值，因為遺傳算法全局搜索的優(yōu)點使得到的輸出可以保證是全局最優(yōu)解或接近最優(yōu)解，而不同于直接使用BP神經(jīng)網(wǎng)絡，初始的權值和閾值是隨機獲得的。優(yōu)化的初始權值和閾值可以幫助克服BP神經(jīng)網(wǎng)絡梯度下降容易陷入局部最優(yōu)的缺點，使得BP網(wǎng)絡依據(jù)優(yōu)化初始權值和閾值在最優(yōu)解附近進行搜索，以此提高網(wǎng)絡的全局性能和網(wǎng)絡的預測精度。

本研究采用3層BP神經(jīng)網(wǎng)絡，其拓撲結(jié)構如圖1所示。該結(jié)構有3個輸入神經(jīng)元、1個輸出神經(jīng)元，由于BP神經(jīng)網(wǎng)絡對隱層神經(jīng)元數(shù)目的確定尚無理論上的指導，所以按經(jīng)驗取6。圖1所示的BP神經(jīng)網(wǎng)絡拓撲結(jié)構，首先把呼吸運動位移數(shù)據(jù)D分解成每三個連續(xù)數(shù)據(jù)作為一組神經(jīng)網(wǎng)絡輸入序列S，然后利用S的每一列呼吸運動數(shù)據(jù)通過如圖1所示的網(wǎng)絡，在對應預期呼吸運動位移輸出T的指導下進行網(wǎng)絡權值和閾值的學習，最后形成BP神經(jīng)網(wǎng)絡輸入輸出的映射關系為

式中，λ表示預測超前的步數(shù)，采用30Hz的數(shù)據(jù)采樣頻率，在試驗中分別用到100、300和500 ms的延時序列，所以在這些不同的預測延時序列中，λ分別取3、9和15。

圖1 BP神經(jīng)網(wǎng)絡拓撲結(jié)構Fig.1 Topology of BP neural network

在圖1中，S為式（6）表示的呼吸運動位移輸入序列，圓圈和里面的符號表示各層網(wǎng)絡節(jié)點和序號，w表示輸入層神經(jīng)元到隱層神經(jīng)元的權重，W表示隱層神經(jīng)元到輸出神經(jīng)元的權重，b和B分別表示隱層和輸出層神經(jīng)元的閾值，［DtDt－1Dt－2］表示一組 BP 神經(jīng)網(wǎng)絡輸入，Dt+λ表示輸入［DtDt－1Dt－2］序列得到的預期運動位移T中的輸出值，λ表示預測超前的步數(shù)。

BP神經(jīng)網(wǎng)絡預測的基本思想是利用各個神經(jīng)元的權值和閾值記憶序列的發(fā)展模式。在合適的樣本下通過訓練得到的權值和閾值在下次有相似的序列輸入，可以得到對應預期的輸出。在神經(jīng)網(wǎng)絡訓練和預測中，樣本的選取對最后的結(jié)果有非常重要的影響，一般要求訓練的樣本必須大于一個周期，并且能夠體現(xiàn)總體數(shù)據(jù)的基本特征。訓練樣本數(shù)據(jù)選擇越合理，訓練后的神經(jīng)網(wǎng)絡權值和閾值包含越準確的信息，預測的精度也越高。

2 實驗結(jié)果與分析

實驗是在Windows XP操作平臺、VC6.0和Matlab 7.6軟件系統(tǒng)下進行的，用VC6.0對胸外標記物隨呼吸產(chǎn)生運動位移的視頻進行檢測、跟蹤，并記錄標記物的運動坐標數(shù)據(jù)；利用Matlab 7.6對獲得的運動坐標數(shù)據(jù)進行預處理，并使用了Matlab 7.6提供的神經(jīng)網(wǎng)絡工具箱（NNET Toolbox）進行GA-BP、BP神經(jīng)網(wǎng)絡的訓練和仿真預測。圖2（a）是所獲得的其中一個樣本數(shù)據(jù)歸一化后的波形，由于在視頻采集的過程中存在著一定的噪聲，所以在實驗的過程中還應該對實驗數(shù)據(jù)進行平滑濾波，以獲得更佳的實驗效果，如圖2（b）所示。

圖2 樣本數(shù)據(jù)濾波前后的波形。（a）濾波前的波形；（b）濾波后的波形Fig.2 Sample data waveform before and after filtering.（a）waveform before filtering；（b）waveform filtered

2.1 實驗的基本步驟

步驟1：在受試者胸左下部貼一標記物，平躺在正常呼吸的情況下；設置攝像機采集頻率為30Hz，拍攝一段標記物隨呼吸產(chǎn)生運動位移的視頻。然后，根據(jù)所得到的視頻，在 VC6.0下利用幀間差分法的運動物體檢測跟蹤方法［11－12］，跟蹤并記錄呼吸過程中胸外標記物的運動坐標數(shù)據(jù)D。

步驟2：在Matlab7.6下，對獲取的標記物運動坐標數(shù)據(jù)D進行預處理，包括歸一化、平滑濾波。

步驟3：取已經(jīng)過預處理的某一段數(shù)據(jù)（大于一個周期），并將其分解成如式（6）左邊S所示的序列，作為GA-BP和BP網(wǎng)絡（神經(jīng)網(wǎng)絡的拓撲結(jié)構如圖1所示）的訓練輸入樣本。在實驗中，用100、300、500ms延時的預期輸出如式（6）右邊T所示，在采樣頻率為30Hz下 分別為3、9和15。對訓練輸入樣本S進行訓練，并獲得相應的訓練網(wǎng)絡。

步驟4：利用所獲得的數(shù)據(jù)和訓練網(wǎng)絡，分別進行GA-BP網(wǎng)絡和純BP網(wǎng)絡的仿真預測。

步驟5：對比GA-BP網(wǎng)絡和純BP網(wǎng)絡的預測結(jié)果。

2.2 實驗結(jié)果及分析

實驗對9個樣本數(shù)據(jù)進行了測試，為了使比較有統(tǒng)一的標準，采用預測結(jié)果的均方誤差（MSE）作為比較預測結(jié)果的準則。實驗分別對GA-BP網(wǎng)絡和BP網(wǎng)絡在不同延時條件下的9個樣本數(shù)據(jù)進行了預測，結(jié)果如表1所示。從表中可以看出，在相同的延時條件下，用GA-BP網(wǎng)絡預測的結(jié)果比BP網(wǎng)絡預測的結(jié)果更佳。而在同一種方法中，較小的延時具有更小的預測誤差，這是因為預測的點與最相鄰的點相關性最高，隨著點之間的距離加大，它們的相關性也會逐漸減?。?3］，致使長時間預測精度下降。

表1 GA-BP網(wǎng)絡和BP網(wǎng)絡的預測結(jié)果比較Tab.1 Comparison of GA-BP network's and BP network's predict error

圖3給出了兩個樣本（樣本1和樣本5）在不同延時下兩種預測算法的預測曲線。其中，（a）和（b）是在100ms下兩個樣本分別用兩種預測算法得到的曲線，（c）和（d）是在300 ms下兩個樣本分別用兩種預測算法得到的曲線，（e）和（f）是在500ms下兩個樣本分別用兩種預測算法得到的曲線。

圖4給出了兩種算法在不同延時下樣本平均預測誤差的曲線。在圖中可以更清楚地看出，GA-BP網(wǎng)絡在不同的延時條件下對預測精度都有一定的提高，而且隨著預測延時的變大，預測誤差也逐漸變大。

圖3 不同延時下的預測曲線。（a）樣本1在100 ms時；（b）樣本5在100 ms時；（c）樣本1在300 ms時；（d）樣本5在300 ms時；（e）樣本 1在 500 ms時；（f）樣本 5在 500 ms時Fig.3 Prediction curve under different delay condition.（a）sample 1 with 100 ms；（b）sample 5 with 100 ms；（c）sample 1 with 300 ms；（d）sample 5 with 300 ms；（e）sample 1 with 500 ms；（f）sample 5 with 500 ms

圖4 不同延時條件下的樣本平均預測誤差。（a）100 ms；（b）300 100 ms；（c）500 100 msFig.4 Average prediction error underdifferent delay condition.（a）100 ms；（b）300 ms；（c）500 ms

3 結(jié)論

本研究對GA-BP網(wǎng)絡結(jié)合的原理和優(yōu)勢進行了探討。利用遺傳算法（GA）良好的全局搜索能力，優(yōu)化BP網(wǎng)絡權值和閾值的初始值，從而克服BP網(wǎng)絡容易陷入局部最優(yōu)、收斂速度慢和學習過程中容易發(fā)生振蕩的不足，使得BP網(wǎng)絡在權值和閾值的訓練學習中能更好地達到全局最優(yōu)解。通過以上GA-BP網(wǎng)絡預測和BP網(wǎng)絡預測的實驗結(jié)果分析，在相同的延時條件下，GA-BP神經(jīng)網(wǎng)絡預測呼吸運動相比于純BP神經(jīng)網(wǎng)絡精度有所提高，特別是在較小延時條件下，GA-BP網(wǎng)絡的預測結(jié)果優(yōu)勢更加明顯。

上述預測是在獲取外標記物運動坐標的基礎上進行的，外標記物運動是隨胸腹部在呼吸時運動而發(fā)生，與內(nèi)部腫瘤的實際運動必然存在著差異。內(nèi)部腫瘤的運動非常復雜而且差異性很大，其運動規(guī)律還需要進一步研究和探索。本研究對外標記物隨呼吸過程的運動預測取得了相對滿意的結(jié)果，如果在以后的工作中可以建立起外標記物運動和腫瘤運動之間的精確關系，那么上述方法在胸腹部腫瘤放射治療中將會有更加實際的用途。

［1］Allen AM，Siracuse KM，Hayman JA，et al.Evaluation of the influence of breathing on the movement and modeling of lung tumors［ J］.International Journal of Radiation Oncology，Biology，Physics，2004，58（4）：1251 －1257.

［2］Treea A，Brocka J，McNaira H，et al.Improvement in tumor control probability with active breathing controland dose escalation：A modeling study［J］.Radiotherapy and Oncology，2009，91（3）：325－329.

［3］Minohara S，Kanai T，Endo M，et al.Respiratory gated irradiation system for heavy-ion radiotherapy［J］.International Journal of Radiation Oncology，2000，47（4）：1097－1103.

［4］Emery R，Rodriguez L，Barsa J，et al.Clinical experience using respiratory gated radiation therapy：comparison of freebreathing and breath-hold techniques［J］.International Journal of Radiation Oncology，2004，60（2）：419－426.

［5］Tchoupo G，Docef A.Nonlinear set membership time series prediction ofbreathing ［A］.In：Luciano Sbaiz， Olivier Cuisenaire，eds.Proceedings of 16th European Signal Processing Conference（EUSIPCO 2008）［C］.Lausanne：EURASIP，2008.25－29.

［6］雷英杰，張善文，李續(xù)武，等.MATLAB遺傳算法工具箱及應用［M］，西安：西安科技大學出版社，2009.30－65.

［7］Man KF，Tang KS，Wong SK.Genetic algorithms：concepts and designs［M］.London：Springer，2001.

［8］Yang Yang，Li Kaiyang.Neural network based on GA-BP algorithm and its application in the protein secondary structure prediction ［J］.Chinese Journal of Biomedical Engineering，2006，15（1）：1－9.

［9］Belew RK，McInerney J，Schraudolph NN.Evolving networks：Using the genetic algorithm with connectionist learning［R］.CS90 －174，1990.

［10］Kitano H.Empirical studies on the speed of convergence of neural network training using genetic algorithms［R］.AAAI－90－118，1990.

［11］吳大鵬，程衛(wèi)平，于盛林.基于幀間差分法和運動估計的Camshift目標跟蹤算法［J］.光電工程，2010，37（1）：55－60.

［12］Kim Changick，Hwang Jenq-Neng.Fast and Automatic Video Object Segmentation and Tracking for Content-Based Applications［J］.IEEE Transactions on Circuits and Systems for Video Technology，2002，12（2）：122－129.

［13］Wu Huanmei， Betty Salzberg， Gregory C Sharp， etal.Subsequence matching on structured time series data［A］.In：Widom J，Chirkova R，eds.Proceedings of the ACM SIGMOD International Conference on Management of Data［C］.New York：ACM，2005，682－693.