王 剛，肖宇峰，鄭又能，田星皓

（1.西南科技大學 信息工程學院，四川 綿陽 621010；2.中國原子能科學研究院核技術(shù)應用研究所，北京 102413）

引言

隨著我國核工業(yè)的快速發(fā)展，核設施和放射源在能源、軍事、醫(yī)療、工業(yè)等領域的應用日益廣泛。在核設施退役、核事故應急處置中，放射性區(qū)域的重建與定位能幫助作業(yè)人員直觀、準確地分析核輻射環(huán)境，提高退役去污、應急處置作業(yè)的效率。

目前對放射性區(qū)域重建主要有康普頓相機成像和γ相機成像兩種方式[1]。國外研究中，普遍采用康普頓成像儀結(jié)合深度相機或者激光雷達實現(xiàn)輻射場的三維重建[2-3]。文獻[4]提出的SDF（scene data fusion）方法使用兩個3D位置敏感高純度鍺探測器組成移動康普頓成像儀，結(jié)合同時定位與地圖創(chuàng)建（simultaneous localization and mapping，SLAM）算法能夠近乎實時重建未知輻射環(huán)境的三維場景圖，其中SLAM算法可以實時提供輻射場景的三維模型，迭代算法用于重建放射源的三維分布模型。SDF 方法實時提供了放射源的劑量分布和位置信息，極大地增強了輻射場景建圖的適用性，但是其結(jié)果未對成像效率和精度進行優(yōu)化。在此思想上，Haefner 等人使用高效多模成像儀（highefficiency multimode imager，HEMI）與Kinect相機構(gòu)成手持式移動平臺，融合γ 射線數(shù)據(jù)和視覺數(shù)據(jù)，可以實時生成精確的放射源三維分布模型，但未對放射源位置信息進行量化[5]；Lee 等人使用康普頓相機和激光雷達搭載在移動機器人上，構(gòu)建了基于體素網(wǎng)格的輻射圖，多源定位的平均誤差在0.2 m以內(nèi)[6]。

相比之下，基于γ相機三維重建技術(shù)的報道和研究不多。在國內(nèi)，中國工程物理研究院材料研究所發(fā)展了針對點源的γ相機重建技術(shù)[7-8]。γ相機是一種高效準確的在線測量儀器，利用光學成像和射線探測機理合成放射源劑量分布圖像[9]，為放射性物質(zhì)的定位、搜尋及后續(xù)處置提供依據(jù)。然而γ相機只能得到放射源的二維圖像[10]，用二維圖像重構(gòu)核污染位置三維分布，是實現(xiàn)快速測定放射性區(qū)域污染位置與分布的重要技術(shù)。對此朱寧等人采用最大似然期望最大化（maximum likelihood expectation maximization，MLEM）算法，通過模擬實驗重建出了放射性體源的三維位置和形狀分布[11]。但這一結(jié)果沒有與周圍環(huán)境信息融合，對放射性分布的描繪尚不充分。

根據(jù)上述情況，本文提出了基于Kinect與γ相機圖像信息融合的放射性區(qū)域重建與定位方法。首先，構(gòu)建Kinect與γ相機組合成像模型并完成其聯(lián)合標定；接著構(gòu)建輻射場景稠密點云地圖并得到Kinect位姿數(shù)據(jù)；然后提取γ相機圖像中放射性分布信息，根據(jù)聯(lián)合標定參數(shù)生成放射性區(qū)域的點云，并依據(jù)Kinect位姿與輻射場景點云地圖融合；最后，用最小包圍盒確定放射性區(qū)域，對中心位置進行計算。本文在真實放射性環(huán)境中開展實驗，重構(gòu)了包含放射源三維分布模型的輻射場景點云地圖，最終實現(xiàn)了放射性區(qū)域的三維可視化。測試結(jié)果證明本文方法可行，具有一定的準確度。

1 方法

1.1 放射性區(qū)域重建與定位方案

本文采用的Kinect版本是Microsoft 研發(fā)的Kinect V2[12]；γ相機是中國原子能研究院研制的全自動γ 射線掃描成像系統(tǒng)，型號為CIAE-ASC。該γ相機以等角度掃描的方式探測區(qū)域內(nèi)各個點的γ 射線強度，形成二維強度分布圖，并將放射源位置與實際環(huán)境圖像集成顯示[13]。

圖1 放射性區(qū)域重建與定位方案Fig.1 Scheme of reconstruction and localization for radioactive area

基于上述平臺，本文設計了如圖1所示的放射性區(qū)域重建與定位方案。Kinect 能以30幀/s的速率同時采集分辨率為1 920×1 080像素的彩色圖像和512×424像素的深度圖像，利用同一時間戳的彩色圖和深度圖構(gòu)建輻射環(huán)境地圖、估計Kinect位姿。γ相機能采集分辨率為1 920×1 080像素的γ圖像，提取γ圖像放射性分布區(qū)域，再將其根據(jù)聯(lián)合標定結(jié)果映射到與γ圖像對應時間戳的深度圖中得到其點云，并依據(jù)Kinect位姿信息與輻射場景點云融合?；谧钚“鼑泄烙嫹派湫詤^(qū)域點云中心位置，對放射源進行三維定位。

1.2 Kinect與γ相機協(xié)作基本原理

作為一類專用于核輻射環(huán)境成像的特殊設備，γ相機與常規(guī)成像設備組合工作的情形很少見。為更好解釋兩臺相機協(xié)作工作原理，基于圖2的安裝結(jié)構(gòu)，下文對相機組合成像模型、聯(lián)合標定方法進行描述。

圖2 兩臺相機的安裝結(jié)構(gòu)Fig.2 Installation structure of two cameras

1.2.1 Kinect與γ相機組合成像模型

圖3 相機組合成像模型Fig.3 Camera combined imaging model

實際工作中，Kinect 安裝在γ相機頂部，并保持兩個相機有盡可能多的共視區(qū)域?；趫D2的安裝結(jié)構(gòu)和相機針孔模型，構(gòu)建了如圖3所示的相機組合成像模型。Kinect 坐標系{D}為OdXdYdZd，世界坐標系{W}（慣性坐標系）為OwXwYwZw，γ相機坐標系{γ}為OγXγYγZγ，圖像坐標系為OXY，像素坐標系為ouv。{W}里一點Pw像素點為m，{D}下的點可通過Kinect位姿Tw，ci轉(zhuǎn)換到{W}。{γ}與{D}存在一個歐式變換關(guān)系，本文用Rγd、tγd描述γ相機相對于Kinect的旋轉(zhuǎn)和平移矩陣，可通過聯(lián)合標定獲得。

1.2.2 Kinect與γ相機聯(lián)合標定

對Kinect和γ相機聯(lián)合標定，獲取兩相機的內(nèi)參及其相對位置關(guān)系Rγd、tγd，才能實現(xiàn)γ相機圖像中放射性區(qū)域到Kinect 深度圖像的精準對齊。本文采用基于平面棋盤圖像的標定方法，實現(xiàn)了Kinect和γ相機的聯(lián)合標定。

聯(lián)合標定主要步驟如下：

1） 圖片獲取：選取內(nèi)角點為11×8的棋盤格作為靶標，如圖4所示。移動靶標，用固定好位置的Kinect和γ相機，拍攝20幅不同位置不同角度的棋盤格圖片。

圖4 棋盤格標定板Fig.4 Chessboard calibration plate

2） Kinect和γ相機的內(nèi)參獲?。悍謩e使用γ相機圖像和Kinect 強度圖檢測內(nèi)角點，將角點像素坐標與其三維坐標相對應。角點的像素坐標作為實際值，三維坐標經(jīng)相機內(nèi)外參投影得到的像素值作為觀測值，通過最小化實際值和觀測值的位置求解最大似然估計問題得到最優(yōu)參數(shù)：相機內(nèi)參K和外參Ri、ti。最大似然估計的評價函數(shù)定義為

式中：n為棋盤格圖像張數(shù)；m為棋盤格角點數(shù)；xij為第i張棋盤格圖像第j個角點的像素坐標；Xj為第j個角點的三維坐標；πs（·）為從 ?3→?2的投影函數(shù)，計算公式如（2）式：

3） 相機間相對位置估計：運用最小二乘法求解γ相機棋盤格角點在Kinect 視圖上的最小重投影誤差，得到聯(lián)合標定的最優(yōu)解，計算公式如（3）式：

式中：xd，ij為Kinect 第i張強度圖像的第j個角點坐標；xγd，ij為第i張γ相機圖像的第j個角點坐標xγ，ij在對應Kinect 強度圖上的投影坐標，計算公式如（4）式：

式中：Kγ、Kd分別是γ相機和Kinect的內(nèi)參矩陣；I為單位矩陣。

按照上述聯(lián)合標定步驟，得到Kinect和γ相機的內(nèi)參結(jié)果如表1所示，聯(lián)合標定的結(jié)果如表2所示，單張圖像的重投影誤差結(jié)果如圖5所示，平均誤差為0.154 63個像素，滿足后續(xù)圖像融合要求。

表1 單個相機標定相機內(nèi)參及重投影誤差Table1 Intrinsic parameter and re-projection error in single camera calibration

表2 聯(lián)合標定的兩個相機位置關(guān)系及重投影誤差Table2 Position relation and re-projection error in joint calibration between two cameras

1.3 基于Kinect 的放射性環(huán)境點云地圖構(gòu)建

現(xiàn)有的ORB-SLAM2[14]方法可用于構(gòu)建環(huán)境的三維稀疏特征地圖，其關(guān)鍵幀保存了對應的彩色圖像、深度圖像和相機位姿，但因為特征稀疏不利于放射性區(qū)域的完整重建。在該方法的基礎上，本文對關(guān)鍵幀的所有像素進行點云重建，構(gòu)建放射性場景的稠密點云地圖。

圖5 聯(lián)合標定重投影誤差Fig.5 Re-projection error of joint calibration

在圖3相機組合成像模型中，設{W}空間一點Pw（x，y，z）的二維像素坐標為m（u，v）。根據(jù)二維彩色圖像和深度圖像計算三維點云的公式如（5）式：

式中：fx、fy分別為u軸和v軸上的歸一化焦距；（uo，vo）為圖像主點坐標；d為Kinect 測得目標點的距離（單位為mm）；λ為實際距離與測得距離的比例系數(shù)，取值為1000。

從{D}到{W}的點云變換公式為

式中：Xci，j為第i個關(guān)鍵幀坐標系上的點云，Tw，ci為第i個關(guān)鍵幀的位姿，Xw，j是Xci，j變換后在{W}上的點云。

當Kinect 運動產(chǎn)生新的關(guān)鍵幀時，構(gòu)建新的局部地圖。待系統(tǒng)運行結(jié)束，將所有的局部地圖整合為完整的地圖。至此可以得到全局一致性的放射性環(huán)境稠密點云地圖和精確Kinect 相機位姿。

1.4 多相機圖像信息融合的放射性區(qū)域點云生成

由γ相機成像原理可知，放射源圖像是疊加在其光學圖像上的，所以可將γ圖像中的放射性區(qū)域根據(jù)聯(lián)合標定參數(shù)投影到對應Kinect 深度圖中得到其點云，最后依據(jù)對應的關(guān)鍵幀位姿將其融合到輻射場景點云地圖中。

1.4.1 γ圖像中放射性分布信息提取

γ相機提供包含放射源的圖像和背景圖像，如圖6（a）、6（b）所示。于是采用背景減除法用包含放射源信息的圖像減去背景圖像，得到放射性區(qū)域，即：

式中：Iγ表示包含放射源信息的伽瑪圖像；Ib表示背景圖像；I表示放射性區(qū)域。

圖6 放射性區(qū)域輪廓提取Fig.6 Contour extraction of radioactive area

對放射性區(qū)域進行灰度化，然后進行閾值處理，得到放射性區(qū)域的二值圖像，如圖6（c）所示。

其中：I（x，y）表示像素點（x，y）的灰度值，T為自適應灰度閾值。得到放射性區(qū)域的二值圖像后，對其進行輪廓提取，結(jié)果如圖6（d）所示。

1.4.2 融合放射性分布信息的點云生成方法

設放射性分布區(qū)域輪廓內(nèi)一點像素坐標為Pγ=（uγ，vγ，1）T，對應深度圖中的坐標為Pd=（ud，vd，1）T，對應三維空間點坐標為Pw=（Xw，Yw，Zw）T，由相機針孔模型得：

其中外參Rγ和tγ是{γ}相對于{W}的旋轉(zhuǎn)和平移矩陣，同理對于Kinect的對應點有：

其中外參Rd和td是{D}相對于{W}的旋轉(zhuǎn)和平移矩陣，聯(lián)立式子（9）和（10）得：

記

則（11）式可以簡化為

將聯(lián)合標定結(jié)果Rγd、tγd代入（13）式可以計算出γ相機視野下的深度值Zγ，進而得到該點點云。

設第i張γ圖像放射性分布區(qū)域輪廓內(nèi)第j個像素點的點云為Xγ，ij，在{W}下坐標是Xwγ，j，則有如下計算公式：

其中Tw，ci表示第i張γ圖像對應的關(guān)鍵幀位姿，可由1.3節(jié)得到。遍歷γ圖像中放射性分布區(qū)域內(nèi)所有像素點，得到該張γ圖像在世界坐標系下的點云。

由于有測量誤差以及深度配準誤差，需對點云進行濾波處理：以點云質(zhì)心為中心點，選取放射性分布區(qū)域的質(zhì)心和輪廓點集中距離質(zhì)心最遠的點，以該兩點距離為半徑r，剔除大于r的點云。

對每張γ圖像采取上述操作，至此可以建立包含放射源三維分布重建模型的輻射場景點云地圖。

1.5 基于最小包圍盒的放射性區(qū)域中心位置估計

主元分析法（principal component analysis，PCA）是求點集最小包圍盒的常用方法，由于包圍球的緊密性不夠好，因此本文采用長方體包圍盒確定放射性區(qū)域點云的最小區(qū)域。

設世界坐標下放射性區(qū)域的點云由Xwγ={Xwγ，1，Xwγ，2，…，Xwγ，j，…，Xwγ，J}組成，令放射性區(qū)域的中心點為

其中J表示放射性區(qū)域點云數(shù)目。放射性區(qū)域的協(xié)方差矩陣可以表示為

由于協(xié)方差矩陣是實對稱矩陣，故其有3個非負特征值λ1、λ2、λ3，對應特征向量為u、v、w。將每個特征向量進行正交標準化結(jié)合質(zhì)心坐標組成坐標變換矩陣C。設Xwγ經(jīng)過坐標變換矩陣轉(zhuǎn)化到主方向坐標系下的點云為Xcγ，其計算公式如（17）式：

為計算Xwγ的最小包圍盒頂點坐標，需先計算主方向坐標系下Xcγ的軸向包圍盒。取Xcγ中各坐標分量的最小值和最大值，構(gòu)造點和以Pcγmin和Pcγmax為對角頂點建立軸向包圍盒。根據(jù)軸向包圍盒的各個頂點信息，利用（18）式將各個頂點轉(zhuǎn)化為{W}下的點，最終獲得點云Xwγ的最小包圍盒。

其中Pw，k，Pc，k分別表示{W}下和主方向坐標系下的包圍盒頂點，k=1，2，…，8。

2 實驗

利用C++語言、OpenCV3.4.0函數(shù)庫和PCL點云函數(shù)庫實現(xiàn)了上述算法，并開展了實驗測試。計算機采用Intel的3.20 GHz（4核）處理器，內(nèi)存容量為4 GB，操作系統(tǒng)為Ubuntu 16.04。

由于γ相機是基于實時掃描的系統(tǒng)，在采集γ相機圖像期間暫停點云地圖構(gòu)建，等γ相機采集完該位置之后就復位到暫停前的位置。這樣的實驗策略可確保同一時間戳的γ相機和Kinect 實現(xiàn)數(shù)據(jù)同步，確保生成的放射性區(qū)域的點云與環(huán)境地圖能夠?qū)崿F(xiàn)空間對齊。

2.1 完整的輻射場景點云效果

在8×12 m2大小的實驗室環(huán)境中，基于Kinect構(gòu)建的完整輻射場景稠密點云地圖如圖7所示。從圖7可以看出，點云地圖和實際三維環(huán)境相符，沒有出現(xiàn)地圖重疊和變形等情況，并且地圖中的暖氣管、椅子、紙箱等雜物清晰可見。其中紅色點云代表Kinect 運動期間光心坐標在全局地圖中的位置，即運動軌跡。在整個實驗中，系統(tǒng)共處理了2 510張視頻幀，關(guān)鍵幀數(shù)目為164幀，F(xiàn)PS（frames per second）為19幀/s。

圖7 輻射場景重建效果Fig.7 Radioactive scene reconstruction effect

目前，針對重建點云地圖的精度評定方式?jīng)]有固定原則。本文以中誤差作為評價點云地圖精度的標準[15]。選取門、顯示屏、紙箱3個規(guī)則形狀特征物，利用鋼尺測量所選特征物的實際尺寸，用Cloud Compare 軟件在點云中量取相應特征物的尺寸，以兩者差值作為Δ 計算標準。從表3可以看出特征物的差值都在厘米級，經(jīng)過實驗對比，最終得到的特征物尺寸中誤差為0.016 m。由此可見，重建的輻射場景點云地圖精度較高。

表3 特征物體尺寸對比Table3 Size comparison of feature objects

2.2 測試放射性區(qū)域點云生成效果

為了生成放射性區(qū)域點云，實驗中，對活度為3 mCi的137Cs 放射源，γ相機在不同位置以不同角度采集了8 張γ圖像，如圖8所示。

圖8 γ圖像數(shù)據(jù)Fig.8 Data of γ images

本節(jié)選取如圖9（a）所示的一張γ圖像，分析其點云生成過程。圖9（b）為圖9（a）中的放射性區(qū)域所對應的二值圖，其中白色區(qū)域是放射性區(qū)域；圖9（c）為將γ圖像中的放射性區(qū)域投影到已融合深度信息的彩色圖像，并用線條描繪出了放射性區(qū)域輪廓。對輪廓里的點根據(jù)（13）和（14）式生成放射性區(qū)域的點云，圖9（d）是圖9（a）中的放射性區(qū)域所對應的點云。

圖9 放射性區(qū)域點云生成過程Fig.9 Point cloud generation process of radioactive area

將8 張γ相機圖像對應的放射性區(qū)域點云根據(jù)各自關(guān)鍵幀位姿信息融合到環(huán)境地圖中，得到如圖10所示的融合放射性區(qū)域點云的場景圖，其中黃色區(qū)域代表放射性區(qū)域的點云模型，8個綠色的點代表γ相機采集γ圖像時在輻射場景點云模型中的位置。

圖10 融和放射性區(qū)域點云的場景圖Fig.10 Scene diagram of point cloud in fusion radioactive area

2.3 測試放射性區(qū)域定位準確性和穩(wěn)定性

為確定放射性區(qū)域點云的三維空間分布區(qū)域，采用基于主元分析法的最小包圍盒算法測試放射性區(qū)域。如圖10所示，箭頭表示放射源的位置，右上角矩形區(qū)域表示放大的包圍盒。從實驗效果圖可以看出，放射性區(qū)域點云分布在放射源周圍，最小包圍盒可以包圍所有的放射性區(qū)域的點云，沒有漏掉，證明了最小包圍盒算法的可行性。

在同一室內(nèi)環(huán)境中，采取相同的實驗策略進行3組測試，采用實際值和估計值的距離作為誤差度量最小包圍盒的定位精度。3次實驗的真實位置、估計位置以及誤差見表4，誤差控制在0.15 m以內(nèi)。采用真實位置和估計位置的RMSE（root mean squared error）度量放射源的定位精度，均方根誤差計算公式為

式中：ti、分別為估計位置和真實位置；n=3；均方根誤差為0.11 m。由此得出，最小包圍盒的中心坐標可以準確估計放射源的位置，證明了最小包圍盒算法對于放射性區(qū)域定位的準確性和穩(wěn)定性。

傳統(tǒng)基于最小二乘法[16]、最大似然估計法[17]、融合-迭代[18]的輻射源定位算法已取得較好的二維定位效果，本文提出的方法可對放射源三維定位。相比文獻[6]采用康普頓相機結(jié)合激光雷達實現(xiàn)的放射源三維定位方法，本文方法更加經(jīng)濟實惠。取文獻[6]中137Cs 相關(guān)實驗數(shù)據(jù)，誤差對比如表5所示。從表5得出，文獻[6]在小場景中取得較好定位效果，在相對更大場景中卻不適用；而本文方法在較大場景中定位精度較高。

表5 兩種不同方法定位誤差對比Table5 Comparison of error by 2 different methods

3 結(jié)論

本文應用Kinect 創(chuàng)建室內(nèi)環(huán)境地圖并獲取位姿信息，同時通過γ相機采集放射源的二維分布圖像，將這兩種視覺傳感器集成到數(shù)據(jù)采集和算法處理框架中，完成對放射性區(qū)域的重建與定位。從實驗數(shù)據(jù)的重建結(jié)果來看，可以在少量γ相機圖像的情況下重構(gòu)放射源三維分布模型，并且能與輻射場景地圖融合，利用最小包圍盒算法確定放射源分布區(qū)域。在獲取放射性區(qū)域重建模型后，融合放射源強度信息將會是今后研究的重點。