尹 愚 宋雪梅 堯德中（電子科技大學生命科學與技術(shù)學院，成都 60054）

2（中國科學院上海生命科學研究院，上海 200031）

引言

腦內(nèi)源性信號的光學成像（optical imaging based on intrinsic signals）是從20世紀70年代發(fā)展起來的一種腦成像技術(shù)，是研究腦的工作機制的重要手段之一。同其它研究手段相比，內(nèi)源性光學成像在時間和空間分辨率上都有較大優(yōu)勢。在哺乳動物的大腦內(nèi)部，執(zhí)行特定功能或有相同特性的細胞聚集在一起［1-2］。這些細胞發(fā)生興奮活動時，皮層神經(jīng)活動的改變會導致相應區(qū)域的血流速度、去氧血紅蛋白（Hbr）濃度等血液參數(shù)的變化，進而引起其所在皮層對特定波長光吸收率的變化。這些變化是依賴于腦活動的內(nèi)源性信號，顯示了大腦皮層的功能構(gòu)筑［3］。在活體腦組織中基于這種血氧動力學而產(chǎn)生的光學信號，其強度變化非常微弱，相對于全部反射光大約占0.1％～0.2％（605nm光照條件）［4］。

因為進行研究的腦功能區(qū)域一般在皮層以下（深度為200～800μm），光學信號在透射出皮層的過程中，大腦組織會對光產(chǎn)生散射作用。這使得記錄到的功能區(qū)域會比真實活動的皮層范圍有所擴散［4］。另外在圖像生成和傳輸過程中的各種干擾也會對圖像產(chǎn)生模糊作用，而這種模糊過程一般可以看作是系統(tǒng)點擴展函數(shù)與圖像卷積的結(jié)果。因此尋找一個對點擴展函數(shù)準確的估計，并利用反卷積方法就能對這種退化圖像進行恢復。對于內(nèi)源性光記錄信號的去噪處理，主要還是采用空間平滑濾波。本研究根據(jù)內(nèi)源性光學信號圖像二階統(tǒng)計量的分析，利用其基于生理特性的分布形式對點擴展函數(shù)進行合理的估計，并恢復出功能柱圖像。相對于恢復后的圖像，圖像噪聲得到了有效抑制，功能柱結(jié)構(gòu)清晰，而且功能區(qū)域更為集中。

1 實驗方法和數(shù)據(jù)處理

1.1 動物準備和光學成像

實驗采用體重2.0～3.0kg的成年貓3只，性別不限。肌肉注射30mg/kg鹽酸氯胺酮誘導麻醉，施行氣管插管和股靜脈插管手術(shù)。將動物俯臥固定于立體定位儀（SN-3，NARISHIGE）；持續(xù)靜脈注射3mg/（kg·h）戊巴比妥鈉和10mg/（kg·h）三碘季胺酚，以維持實驗過程中動物的麻醉和麻痹狀態(tài)；150mg/（kg·h）葡萄糖溶液為動物提供必要的營養(yǎng)物質(zhì)。人工呼吸機（CIV-100，Cloumbus Instruments）設置為20次/min，潮氣量約為12.3mL/kg。動物體溫維持儀（SS20-1，蚌埠市大江電子公司）將被試動物體溫控制在38℃～38.5℃。實驗過程中監(jiān)視動物心電（ECG）和腦電（EEG），以了解動物的麻醉狀態(tài)和生理情況。根據(jù)圖譜在Horsley-Clarke坐標系的P0-P10、L0-L6處，開顱暴露皮層17區(qū)，面積約為5mm×5mm，剝離硬腦膜。在暴露皮層的顱骨窗外圍用牙科水泥封固一個金屬小室，使用甲基硅油填充滿小室后加玻璃蓋片封閉以維持適當壓力，從而減少心跳和呼吸導致的皮層表面機械波動［5］。動物眼睛加滴1％阿托品擴瞳，5％新福林收縮瞬膜。然后佩戴接觸鏡，用以矯正視力和防止角膜干燥。用檢眼鏡將動物眼底血管投影于視覺刺激顯示屏上，以確定視盤位置。在接觸鏡外放置直徑為3mm的人工瞳孔。

高靈敏度數(shù)字CCD攝像機［6］（iXon897，ANDOR）在波長550nm（半波寬為30nm）的綠光照射條件下，將鏡頭聚焦于皮層表面，拍攝出皮層表面的血管圖像。然后改為使用波長為605nm（半波寬為10nm）的紅光進行照明，并通過對CCD電動支架調(diào)節(jié)，使得鏡頭聚焦于皮層表面以下約500μm處，以減弱血管的干擾，進行腦功能光學圖像的采集。

1.2 視覺刺激和信號采集

視覺刺激通過一臺刺激服務器連接一臺顯示器給出。顯示器置于貓眼前57cm處，刺激范圍為30°×40°視角［7］。視覺刺激的平均亮度為20cd/m2，對比度為100％。運動的正弦光柵作為視覺刺激，其時間頻率為2Hz，空間頻率為0.8c/deg，刺激朝向為0°～157.5°，間隔22.5°共八種。高靈敏度CCD攝像機正對暴露皮層區(qū)域，拍攝速度為20幀/s，空間分辨率為16μm，景深約為0.1mm。每個刺激周期內(nèi)視覺刺激呈現(xiàn)前持續(xù)拍攝2s（共計拍攝40幀），刺激呈現(xiàn)后持續(xù)拍攝8s（共計拍攝160幀）。不同朝向的光柵在不同的記錄周期內(nèi)隨機分別呈現(xiàn)，實驗重復記錄10次。視覺刺激為單眼刺激，即僅暴露一只眼，遮蓋另一只眼進行刺激。

1.3 數(shù)據(jù)預處理

由于被探測皮層不可能處于理想的均勻光照下。為獲得皮層活動的圖像，首先將皮層受到刺激時獲取的圖像與基準圖像進行對比，以修正光照的不均勻性。將實驗中全部刺激呈現(xiàn)前所拍攝的圖像進行疊加平均作為皮層的基線圖像（baseline image），這被稱為空白本底（blank）。將多次刺激的全部圖像進行疊加平均的目的是減弱生物體生理活動和光照系統(tǒng)的低頻波動所引入的噪聲，多次試驗的疊加平均可以將圖像中每個像素的標準差減小為噪聲標準差的1/（K是進行疊加平均的圖像數(shù)目）。對內(nèi)源性光學信號圖像進行標準分析，將某種刺激條件下，獲得的全部拍攝圖像進行疊加平均，然后與本底圖相除。而后使用“第一幀分析”技術(shù)（first frame analysis），去除低于0.3Hz的低頻噪聲（主要是呼吸等生物體活動引入的干擾）。即將每幀圖像都減去未刺激時所采集的第一幀圖像。圖1（b）顯示了實驗動物1受朝向為0°刺激時，通過上述兩項預處理后所得到的結(jié)果。對照于圖1（a）的直接拍攝所得圖像，從中可見存在以簇狀成團方式呈現(xiàn)的功能活動區(qū)。在圖像中暗色區(qū)域是有著較強的光吸收區(qū)域（相對于本底圖），即為有著較強活動的功能柱區(qū)域。通過多次重復實驗的疊加平均，能將各種低頻周期波動所引入的噪聲進行一定程度的壓制。但是因為大腦組織對光的散射，內(nèi)源信號會發(fā)生一定程度的擴散。所以圖中所顯示的暗色區(qū)域邊緣部分并非完全是皮層的相關(guān)活動功能區(qū)，有部分是因為光學彌散而產(chǎn)生的偽跡［4］。

圖1 內(nèi)源光學信號。（a）CCD采集圖像；（b）預處理結(jié)果Fig.1 Opticalimaging based on intrinsic signals.（a）CCD image；（b）function image

1.4 二階統(tǒng)計分析

二階統(tǒng)計量既方差所表示的是偏離平均值的誤差情況，因此可以通過計算像素點（m，n）鄰域S內(nèi)的二階統(tǒng)計量來評價像素點間的差異情況。

鄰域S范圍選擇192μm×192μm的空間大小，實驗動物1的處理結(jié)果如圖2（a）所示。其中暗色區(qū)域為f*（m，n）值較小的區(qū)域，在這些位置上像素點（m，n）與其相鄰區(qū)域S內(nèi)的像素間的差異較小。可以觀察到圖1（b）中最暗和最亮的區(qū)域都對應著圖2（a）中的最暗區(qū)域。這說明在功能圖1（b）中，因視覺刺激而激活的皮層區(qū)域內(nèi)所產(chǎn)生的光學信號的相對變化程度的一致性較高。這是由于激活的皮層區(qū)域內(nèi)大量功能相同的神經(jīng)元發(fā)生放電活動，放電活動代謝能量而引起血氧濃度變化，使得對照射光的吸收程度增強。在有著較強活動的功能柱區(qū)域內(nèi)，功能相似的神經(jīng)元的活動所產(chǎn)生的光學信號也相似，因此這些位置上圖像二階統(tǒng)計處理結(jié)果的值也較小。這樣的結(jié)果是由功能柱的結(jié)構(gòu)和生理活動所決定的。圖2（a）結(jié)果中顯示在未激活的皮層區(qū)域也同樣保持著較小的二階統(tǒng)計量。這是因為在非活動的皮層區(qū)域神經(jīng)元未被激活，其生理活動也一致的保持為較低水平，所以這個區(qū)域的皮層對照射光吸收程度的變化差異不大。因此在功能區(qū)域和非功能區(qū)域內(nèi)的二階統(tǒng)計量較為一致，其差異主要體現(xiàn)為均值的不同。白色區(qū)域表示這些區(qū)域的f*（m，n）值較大，說明像素點之間的差異較大。對應于圖1（b），這些位置是處于功能柱和未激活皮層之間。這些區(qū)域并非與刺激相關(guān)的神經(jīng)元集中的位置，因此并沒有較強的放電活動代謝能量而引起血氧濃度變化。在這些位置光學信號的變化更可能是由于附近位置血氧濃度擴散和不完全透明的腦組織造成的彌散結(jié)果。在這些部位，靠近功能柱區(qū)域的顏色偏暗，而靠近未激活皮層區(qū)域的顏色偏亮。圖2（b）為圖像二階統(tǒng)計量歸一化值的分布直方圖。

圖2 功能圖二階統(tǒng)計結(jié)果。（a）二階統(tǒng)計結(jié)果；（b）二階統(tǒng)計量分布直方圖Fig.2 Second-order statistics results.（a）the second-order statistics image；（b）the secondorder statistics histogram.

1.5 點擴展函數(shù)與估計

內(nèi)源光學信號十分微弱，而且預處理后的圖像對比度不高。為盡可能利用所采集的圖像信息，對功能柱區(qū)域的信號部分與彌散部分進行分離，需要對所采集的圖像進行恢復處理。彌散信號部分可以認為是信號在皮層中的擴散結(jié)果，這被視為是圖像的一種模糊過程，使得圖像產(chǎn)生了退化。通常使用以下卷積模型來模擬模糊過程

式中，g為模糊圖像；h為點擴展函數(shù)（PSF）；f為原圖像；n為隨機噪聲（通常被假設為具有均值為0，方差為σ2的高斯分布）；*表示卷積運算。如果可以得知點擴展函數(shù)的形式，就能通過反卷積恢復出真實的圖像。

因為無法得知點擴展函數(shù)的真實情況，為恢復出原始圖像，需要在合理的假設條件下對點擴展函數(shù)進行估計。在許多光學測量系統(tǒng)和成像系統(tǒng)模型中，一般都假定點擴展函數(shù)形式為高斯函數(shù)［8］。一種典型的高斯點擴展函數(shù)模型為

式中，σ為高斯函數(shù)的方差。要直接對二維PSF進行估計比較困難和復雜［9］，根據(jù)1.3中光記錄預處理的描述可知經(jīng)過預處理的圖像在照明上可以認為是大體上均勻的，因此該PSF可以假設為具有圓對稱特性［10］。設x方向上的線擴展函數(shù)（LSF）為h*（x），并且采用式（3）的形式，于是有

這樣將二維PSF的估計轉(zhuǎn)化為一維LSF的估計。

對于PSF的估計需要選取合理的考察區(qū)域，因為功能柱圖中沒有明確的邊緣特征，傳統(tǒng)的階躍邊緣法［9，11］和邊緣檢測法［12］等利用邊界信息進行PSF估計的方法并不適用。根據(jù)上文對內(nèi)源性光學圖像的二階統(tǒng)計分析，在圖像中共有三種活動模式：功能柱區(qū)、彌散區(qū)和非功能區(qū)。假設活動模式是具有高斯分布的形式而進行的，并且認為在預處理中其它噪聲已被得到有效的處理，可以暫不考慮其影響。所以將每個區(qū)域活動形式與LSF相卷積可以得到

式中，fa（x）為功能柱區(qū)的活動函數(shù)，fu（x）為非功能區(qū)的活動函數(shù)，fd（x）為彌散區(qū)的活動函數(shù)，其形式均為高斯函數(shù)。根據(jù)高斯函數(shù)的卷積規(guī)則，由式（5）直接可得

式中，σa、σu分別是功能柱區(qū)、和非功能區(qū)的真實活動方差。由于彌散區(qū)中既有功能柱區(qū)的活動也有非功能區(qū)的活動，反映的是兩個活動之和，所以可以用兩者的高斯函數(shù)相加來描述。根據(jù)高斯函數(shù)相加的原則，可知相加后的函數(shù)方差應為（σa+σu）2。σh為需要估計的LSF的方差。σga、σgd、σgu為模糊圖像中這三個區(qū)域的方差。從功能圖像中這三個區(qū)域取若干樣本，即可以由式（6）解出σh。實際上此選取不同區(qū)域的樣本有所不同，得到的σh也不會完全相同?？梢愿鶕?jù)如下準則：在線擴展函數(shù)h*（x）定義域的任何子集Ω中，整合后h*（x）的能量應該與所有樣本得到的（x）的能量的平均值相同［13］，即

基于上述準則，可以從整個圖像的全部樣本點中估計出與圖像唯一對應的h*（x）。

2 結(jié)果

2.1 點擴展函數(shù)估計結(jié)果

利用對內(nèi)源信號功能圖二階統(tǒng)計分析的結(jié)果，結(jié)合預處理功能圖的信號與非信號區(qū)域，針對實驗動物1，共選取200組特征差異明顯的像素點作為樣本。通過式（6）和式（7）得到對該組數(shù)據(jù)的PSF的估計，結(jié)果如圖3所示。

圖3 點擴展函數(shù)估計結(jié)果Fig.3 Estimation of point spread function

2.2 反卷積恢復

根據(jù)對點擴展函數(shù)的估計結(jié)果，對實驗動物1的預處理圖像結(jié)果圖1（b）進行反卷積圖像恢復，圖4（a）為功能圖像恢復結(jié)果。相比預處理圖像，圖中的噪聲得到了有效的抑制，同時保留的圖像細節(jié)。更為重要的是二階統(tǒng)計量分布直方圖發(fā)生了明顯變化。預處理圖像的二階統(tǒng)計量的分布（見圖2（b））沒有明顯的峰值，很難從中選取一個閾值以劃分二階統(tǒng)計量較小區(qū)域和較大區(qū)域。在反卷積恢復后圖像的二階統(tǒng)計量分布圖中（見圖4（c）），其分布形式發(fā)生明顯變化，左端的分布變得更為集中?？梢娫诙A統(tǒng)計值的較小端形成明顯的峰值（圖中用“*”標出），并且在峰的左右分布形狀對稱。而在其附近同時存在另外一個較為微弱的分布峰。相對于恢復前圖像的直方圖，圖4（c）的主峰明顯，之前未見的次峰也能可見。根據(jù)圖像閾值劃分的原理，可以在這兩個分布峰之間的谷值處選取一個值以作為二階統(tǒng)計量的大小的劃分閾值，如圖中“↓”標出。因為從預處理功能圖和功能柱的生理機制可知，功能柱區(qū)域是一個功能相似的神經(jīng)元集中區(qū)域，根據(jù)1.4節(jié)可知，共同相似的神經(jīng)功能活動所引起的光學信號在二階統(tǒng)計量上應表現(xiàn)為較小和分布集中。因此有效的功能區(qū)域在二階統(tǒng)計量的分布上對應于其二階統(tǒng)計值較小端的峰值。

圖4 反卷積處理。（a）功能圖恢復；（b）二階統(tǒng)計量劃分模板；（c）恢復圖像的二階統(tǒng)計量分布直方圖Fig.4 Deconvolution results.（a）function image restoration；（b）edge template based on second-order statistics；（c）second-order statistics histogram of deconvolution result.

圖4（b）中顯示了利用該選取的閾值從反卷積恢復圖像中得到分割模板結(jié)果。因為功能柱區(qū)域和非功能柱區(qū)域均屬于二階統(tǒng)計量值小的區(qū)域，所以該模板圖中既包含功能柱區(qū)域，又包含非功能柱區(qū)域。結(jié)合反卷積恢復圖可以準確劃分激活功能柱區(qū)域。

2.3 功能柱區(qū)域分割

激活區(qū)域應滿足兩點：位于由二階統(tǒng)計分析得到的分割模板的白色區(qū)域（二階統(tǒng)計量小的區(qū)域）；位于反卷積恢復圖中灰度低于全圖平均值區(qū)域（激活區(qū)域）。據(jù)此標準對激活功能柱區(qū)域進行了劃分，3只實驗動物通過該方法所得結(jié)果分別表示于圖5中。

圖5 功能柱區(qū)域劃分。（a）～（c）分別為3只實驗動物結(jié)果Fig.5 Edge segmentation.（a）～（c）the results of the three cats，respectively

3 討論

3.1 反卷積圖像增強

光學信號在傳遞、采集過程中都會受到各種因素的干擾，使得圖像產(chǎn)生模糊。這種模糊可以看作是系統(tǒng)點擴展函數(shù)與圖像卷積的結(jié)果。因而如果可以預先知道系統(tǒng)的點擴展函數(shù)，就可以通過傳統(tǒng)的圖像復原算法對模糊圖像進行逆運算，得到較清晰的復原圖像。腦內(nèi)源性光記錄信號不僅非常微弱而且噪聲成分復雜。傳統(tǒng)處理手段主要以空間濾波為主，雖然處理結(jié)果使得圖像更加平滑，圖像的顆粒噪聲也能得到抑制。內(nèi)源信號光學成像去噪使用的平滑濾波是利用濾波核對圖像進行卷積，實質(zhì)是個模糊圖像的過程。腦內(nèi)源性光記錄信號又具有彌散性的特點，生理結(jié)構(gòu)的組織上沒有明顯邊界。用空間濾波去增強圖像會使得原本彌散的結(jié)果更加模糊，對功能柱區(qū)域劃分產(chǎn)生干擾。利用腦內(nèi)源性光記錄信號的生理機制特點，以其二階統(tǒng)計量的分布關(guān)系為基礎所估計的點擴展函數(shù)，反卷積圖像增強的結(jié)果良好。功能柱區(qū)域不再是不可分辨的一團模糊的暗色區(qū)域，而是形成具有類似等高線結(jié)果的階梯區(qū)域。這樣的結(jié)果更準確的反應了功能柱區(qū)域發(fā)生活動最強烈的位置。由于反卷積圖像增強是去除圖像的模糊過程，在恢復結(jié)果中使得功能柱區(qū)域更為集中，使彌散區(qū)域能產(chǎn)生一定程度的復原。

3.2 內(nèi)源光學信號功能柱定位

腦內(nèi)源性光記錄研究工作需要對腦活動的功能柱區(qū)域進行定位和提取。特別是針對某些需要利用功能柱邊界或者明確的功能柱區(qū)域信息所進行的研究。這類研究希望能通過內(nèi)源性光學記錄結(jié)果的引導，迅速準確的獲取相關(guān)區(qū)域的信息。傳統(tǒng)上多采取圖像灰度閾值劃分，但由于圖像對比度低，不同情況下采集數(shù)據(jù)結(jié)果均有不同，而閾值選擇主觀性比較強。另外更是因為功能柱區(qū)域的彌散，生理結(jié)構(gòu)無明顯邊界，這些都使得直接采用圖像灰度閾值劃分的結(jié)果無法準確定量，很難準確分割出功能柱區(qū)域。本研究中，利用功能柱活動的生理結(jié)構(gòu)基礎，根據(jù)相關(guān)活動區(qū)域的二階方差應當較小的這一結(jié)果，通過反卷積恢復后圖像彌散區(qū)域被有效去除的條件下，給出閾值選取的合理標準。最后結(jié)果顯示，劃分結(jié)果清楚，對血管等干擾的屏蔽作用強，基于這個方法下的功能柱劃分是具有實用意義。

4 結(jié)論

本研究基于估計點擴展函數(shù)對內(nèi)源性光學成像信號的提取，在對本實驗室采集的內(nèi)源性光記錄信號的處理應用中得到了較好結(jié)果，證明了此方法適用于內(nèi)源性光記錄信號的圖像增強和功能柱區(qū)域提取。并在一定程度上比直接采用圖像灰度閾值劃分更為合理。此方法對于內(nèi)源性光記錄信號在腦功能成像的研究，不僅提升了圖像質(zhì)量，更為明確的確定了功能區(qū)域。對于利用內(nèi)源性光記錄信號在神經(jīng)活動的探索中具有進行挖掘信號的實用價值。在此后進一步研究中，計劃結(jié)合電生理記錄對本方法進行的功能柱區(qū)域劃分結(jié)果進行單細胞記錄的驗證。

（致謝：感謝李朝義院士的大力支持與指導。感謝徐杏珍老師在動物實驗中的認真指導與熱情幫助。）

［1］Mountcastle VB.Modality and topographic properties of single neurons of cat＇s somatic sensory cortex［J］.Journalof Neurophysiology，1957，20（4）：408-434.

［2］HubelDH，Wiesel TN.Receptive fields and functional architecture in two nonstriate visual areas（18 and 19）of the cat［J］.J Neurophysiol，1965，28：229-289.

［3］Grinvald A.Functional architecture of cortex revealed by optical image of intrinsic signals［J］.Nature，1986，324：361-364.

［4］Grinvald A.In-vivo optical imaging of cortical architecture and dynamics［M］//Windhorst U，Johansson H，eds.Modern Techniques in Neuroscience Research.Berlin：Springer-Verlag，1999：893-969.

［5］俞洪波，邢大軍，壽天德.腦內(nèi)源信號光學成像術(shù)：貓視皮質(zhì)方位功能柱的活體顯示［J］.中國神經(jīng)科學雜志，2000，16（4）：355-359.

［6］White BR，Bauer AQ，Snyder AZ，et al.Imaging of functional connectivity in the mouse brain［J］.PloS I，2011，6（1）：e16322.

［7］俞洪波，壽天德.用腦光學成像術(shù)研究不同空間拓撲位置貓初級視皮層的空間頻率反應特性［J］.生理學報，2000，52：411-4151.

［8］呂成淮，何小海，陶青川，等.圖像復原中高斯點擴展函數(shù)參數(shù)估計算法研究［J］.計算機工程與應用，2007，43（10）：31-34.

［9］李東興，趙剡，許東.基于參數(shù)估計的降晰函數(shù)辨識及圖像復原算法［J］.紅外與激光工程，2010，1：166-172.

［10］Claxton CD，Staunton RC.Measurement of the point spread function of a noisy imaging system［J］.Journal of the Optical Society of America A：Optics，Image Science，and Vision，2008，25（1）：1592170.

［11］Vishwakumara K，Martens J B.Estimation of perceived image blur using edge features［J］.International Journal of Imaging Systems and Technology，1996，7：102-109.

［12］Elder JH，Zucker SW.Local scale control for edge detection and blur estimation［J］.IEEE Trans on Pattern Analysis and Machine Intelligence，1998，20（7）：699-716.

［13］吳俊芳，劉桂雄，韋寧.散焦含噪圖像的點擴散函數(shù)估計與邊緣檢測［J］.華南理工大學學報（自然科學版），2010，38（6）：118-121.