基于最小二乘正則相關性分析的顱骨身份識別

周明全，楊 穩，林芃樾，耿國華，劉曉寧，李 康

（1.西北大學信息科學與技術學院，陜西西安710127；2.北京師范大學信息科學與技術學院，北京100875）

1 引 言

顱骨身份識別是法醫學中一個重要的研究課題。在很多法醫學案例中，技術人員只能獲取到受害人的顱骨，并且沒有其他證據線索，這種情況下無法使用傳統的身份識別技術，顱骨識別［1］便成為受害者身份認證的重要技術手段之一。隨著CT、3D掃描等數字化技術的發展，三維顱面采集變得更加容易，計算機輔助顱骨識別成為研究熱點，內容涉及法醫學、信息學、人類學和計算機圖形學等多個領域。

目前，顱骨識別研究主要集中在顱像重合［2］和顱面復原［3］兩方面。顱像重合是將未知身源顱骨的2D投影與失蹤人2D照片相疊加，根據其相關特征找到一個合適的匹配，從而獲取到顱骨人臉疊加像，即可判斷顱骨和人臉照片是否屬于同一人。顱面復原［4］是通過分析顱骨和面皮之間的內在關系，復原出未知顱骨面貌的技術。當前主流的顱面復原方法采用統計學習技術，如統計形狀模型［5］和統計回歸模型［6］等。

顱像重合［7-11］和顱面復原［12-13］都需要利用形態學知識準確提取顱骨和面皮之間的內在關系，存在不確定性較高和識別能力較低等問題。隨著基于圖像的3D人臉建模［14］和3D數據采集技術［15］的發展，使用3D顱骨和3D面皮之間的相關性識別顱骨逐漸成為研究熱點問題。Duan［16］等人首次提出通過分析顱骨和面皮之間的形態相關度來識別未知顱骨，并使用典型關聯分析［17］提取映射后兩組相關性最大的主成分向量，用于匹配顱骨和面皮。但由于顱骨與面皮之間關系復雜，該方法的識別率有待提高。

本文提出了一種基于最小二乘正則相關性分析（Least-Squares Canonical Dependency Anal?ysis，LSCDA）的顱骨識別方法。該方法通過構建顱骨和面皮的相關性分析模型，計算未知顱骨和面皮庫中每張面皮的匹配程度，得到識別結果。與顱像重合和顱面復原不同，該方法無需準確提取顱骨和面皮之間內在的復雜關系，僅測量它們之間的相關性即可確定顱骨身份，降低了顱骨識別問題的復雜性，提高了識別率。

2 數據獲取與預處理

本研究使用了255套顱面數據，這些數據來源于西北大學可視化技術研究所，女性年齡在19～75歲，男性年齡在21～67歲。顱面CT數據由陜西中醫藥大學附屬醫院使用西門子多排探測器螺旋CT機采集，每個圖像以DICOM格式存儲，大小約為512×512×250。原始CT切片圖像在去噪后，使用Sobel算子提取顱骨和面皮輪廓，采用Marching Cube算法［17］進行重建，得到顱骨和面皮的三維模型。顱骨和面皮分別表示為大約包含220 000和150 000個頂點的三角形網格。另外，每套顱面還記錄了年齡、性別和BMI等屬性。

為了消除因位置、姿態和尺度等因素造成的數據不一致，所有樣本都被轉換在統一的法蘭克福坐標系下。法蘭克福坐標系是由左右雙側耳門上點、左眼眶下邊緣點和眉間點4個顱骨特征點確定的，分別用Lp，Rp，Mp和Vp表示。法蘭克福平面由Lp，Rp和Mp3點確定［18］。如圖1所示：以Lp Rp為法向量且過點Vp的平面與直線Lp Rp的交點記為法蘭克福平面的原點O′；Lp，Rp以及Mp3點所構成的平面為坐標系中的XO′Y平面，左耳門上點Lp到右耳門上點Rp的方向為X軸方向；過原點O′且與XO′Y平面垂直向上的方向記為Z軸；同時垂直于X軸、Z軸的直線記為坐標系的Y軸，Y軸的正方向由右手法則確定。

圖1 歸一化坐標Fig.1 Normalized coordinates

統一坐標系后，再進行尺度歸一化。設定所有顱骨模型的Lp到Rp的距離為單位1，則對顱骨每一個頂點（x，y，z）進行尺度變換為（x/|LpRp|，y/|LpRp|，z/|LpRp|）。

由于面部特征信息主要集中于前半部分［19］，為了消除不相關數據的影響和減少算法耗時，對顱面樣本進行分割，只保留與識別有直接關系的前半部分，分割后的顱骨和面皮分別為包含約160 000和80 000個頂點的三角形網格。為了建立相同的點對應關系，選取一組分割后的顱骨和面皮作為參考模板，如圖2所示。其余所有顱骨模型和面皮模型通過非剛性配準方法［20］自動配準到參考模型中。

圖2 顱骨和面皮參考模型Fig.2 Skull and skin reference models

3 構建統計形狀模型

在進行LSCDA分析前，需要分別為顱骨和面皮建立統計形狀模型。統計形狀模型［21］是一個參數化的可變形模板，每個顱面都有一個參數化的表示，即每個顱面對應著參數空間中的一個點。主成分分析（Principal Components Analysis，PCA）是建立顱面統計形狀模型的有力工具。PCA能夠反映訓練樣本數據集真實分布的一組正交軸，這組正交軸原點設在數據集的中心，表示訓練數據中具有最大方差的一組方向，數據能夠以最小誤差重構，從而達到降維的效果。

通過連接所有頂點的坐標，顱骨和面皮可以分別表示為一個高維度的向量。因此，得到一個顱骨的訓練數據集和一個面皮的訓練數據集其中N表示訓練樣本個數，m和n分別為顱骨和面皮的頂點數，每個坐標索引標記訓練集上的對應點。假定P1，P2，...，Pp是均值歸一化顱骨協方差矩陣的單位正交特征向量是樣本顱骨向量的均值。PCA變換可以將每一個顱骨向量Si通過變換投影到以均值為原點，以這些主方向為坐標軸的坐標系中。那么，顱骨統計形狀模型可以表示為：

其中：a=(a1，a2，...，ap)Τ是組合參數，ai滿足高斯分布，在模型中只要給定組合參數a即可產生顱骨參數模型S(a)。

所有顱骨數據和面皮數據均進行了非剛性配準操作，具有相同的點對應關系，可以通過PCA變換來確定形狀模型參數。對于每一個顱骨向量Si，變換為：

其中：P=[P1，P2，...，Pp]，ai∈Rp稱為顱骨向量Si的主成分向量。類似地，可通過相同的處理過程來構建面皮統計形狀模型：

因此，通過構建統計形狀模型，將顱骨和面皮數據分別投影到其低維形狀參數空間ap和bq中。

4 顱骨和面皮的主相關信息提取

通過構建統計形狀模型，每個顱骨可以表示為p維特征向量，每個面皮表示為q維特征向量。LSCDA旨在通過平方損失互信息［22］（Squaredloss Mutual Information，SMI）找到ap和bq之間的統計相關關系，它可以在線性投影下捕獲兩個變量的主相關信息。

因此，用LSCDA找到ap和bq的低維子空間，投影的相關性通過線性降維達到最大值，即有：

其中U和V是轉換矩陣。

其中p和q分別是u和v的投影維數，Ip是p維單位矩陣，Iq是q維單位矩陣。

由此可知，尋找轉換矩陣U和V使得u和v之間的相關關系最大化。這里通過SMI來測量其相關性：

其中：pu，v(u，v)是聯合概率分布密度，pu(u)和pv(v)分別是u和v的邊緣概率分布密度。

然而，pu，v(u，v)，pu(u)和pv(v)通常是未知的，不能直接最大化SMI。因此，本文采用最小二乘互信息（Least-Squares Mutual Information，LSMI）方法直接估計密度比，而不用估計每個密度，即有：

密度比函數g?(u，v)由線性組合建模表示為g(u，v)，具體如下：

其中α=(α1，α2，...，αb)Τ是從樣本中學習的參數，φ(u，v)=(φ1(u，v)，φ2(u，v)，...，φb(u，v))Τ是一個基函數，并且對于所有的u和v，φ(u，v)≥0b，0b是一個b維零向量。

參數α是通過最小化平方誤差δ來估計的，即：

這里，令：

于是有：

近似估計H和h的值如下：

將問題轉換為：

通過計算上述目標函數的導數并將它設為零，可以得到：

其中Ib是b階單位矩陣。

那么SMI的近似值表示如下：

綜上所述，顱骨和面皮之間的主相關關系提取算法過程如下：

首先，最合適的變換矩陣U和V是以迭代的方式確定的，具體步驟如下：

Step1：根據式（4）初始化U和V；

Step2：根據式（16）對參數α進行優化；

Step3：采用序列二次規劃技術來更新U和V，重復此過程直至收斂。

然后，通過線性投影將顱骨形狀參數轉換成LSCDA空間參數，即：

同理，通過線性投影將面皮形狀參數轉換為LSCDA空間參數，即：

最后，計算庫中的未知顱骨與每張面皮的匹配分數，匹配得分最高的面皮就是識別結果。定義顱骨和面皮之間的匹配程度如下：

其中表示內積。

5 基于LSCDA的顱骨識別方法

未知顱骨與面皮之間的識別將分別從整體和區域兩個方面進行。首先，構建整體相關性分析模型，從全局來分析顱骨和面皮之間的關系；其次，由于顱面形態非常復雜，不同顱骨和面皮不同區域的相關度不同。例如，前額區域的相關性較強，因為該區域的軟組織厚度較薄；而鼻子、眼睛和嘴巴等區域的相關性則較弱，因為這些區域的結構復雜。強關聯區域具有較高的識別能力。然而，整體相關性分析模型只反映了全局意義上顱骨和面皮之間的相關性，并沒有充分利用它們的特性信息。為了更加可靠地測量顱骨與面皮之間的相關性，提高識別能力，提出了一種基于區域的局部化識別方法。識別流程如圖3所示。

圖3 基于區域的顱骨識別流程Fig.3 Flow chart of skull recognition based on regional fusion method

5.1 基于LSCDA的整體顱骨識別

基于LSCDA整體顱骨識別的目的是從三維面皮庫中找到一個與未知顱骨最可能匹配的面皮。因此，需要計算未知顱骨與每個3D面皮之間的匹配分數。整體顱骨識別的框架如圖4所示，識別過程包括兩個階段：訓練階段和識別階段。其中，S表示顱骨訓練集，F表示面皮訓練集，u和v分別表示LSCDA提取后的子空間向量。

訓練階段，建立相關性分析模型的步驟如下：

Step1：使用所有訓練數據構建顱骨和面皮的統計形狀模型，如式（1）和式（3）所述；

Step2：根據式（2）將所有訓練數據投影到形狀參數空間中；

Step3：針對訓練顱骨和面皮的形狀參數特征執行LSCDA，并獲得兩個映射向量u和v。

識別階段，使用相關性分析模型識別未知顱骨，步驟如下：

Step1：通過上述的統計形狀模型將未知顱骨和面皮庫中的面皮分別投影到形狀參數空間中；

Step2：根據式（4），將它們的形狀參數特征投影到LSCDA子空間中。

圖4 整體顱骨識別框架Fig.4 Flow chart of skull recognition based on holistic approach

Step3：根據式（20）計算未知顱骨與面皮庫中每個面皮之間的匹配分數，具有最高匹配分數的面皮即為識別結果。

5.2 基于LSCDA的區域顱骨識別

如圖5所示，將顱骨和面皮分成幾個相應的生理特征區域，即前額、眼睛、鼻子、嘴巴和輪廓。訓練集之間的對應關系已經建立，只需要手動分割一個樣本，其他樣本可以通過點對應關系實現自動分割。基于LSCDA的區域顱骨識別框架與上述整體顱骨識別框架大致相同，采用該框架分別為5個區域構建相關性分析模型。在識別階段，可以單區域相關性分析進行未知顱骨的識別；還可以將5個相關性分析模型進行融合，進一步提高未知顱骨的識別率。

對于一個顱骨和面皮匹配對，可以使用5個相關性分析模型獲得5個區域匹配分數，最終匹配分數通過加權和計算，如式（21）所示：

圖5 顱骨和面皮的各個區域Fig.5 Various areas of skull and face

其中：rp，rf，re，rm和rn分別代表輪廓、前額、眼睛、嘴巴和鼻子的區域匹配分數；wp，wf，we，wm和wn分別表示各個區域相對應的權值。

如何選擇式（21）中每個區域的權值對顱骨識別起到決定性的作用。由于顱骨識別依賴的是顱骨與面皮之間的相關性，因此應該將高權重分配給高相關的區域，如輪廓和前額。相比之下，顱骨和面皮之間的關系在鼻子，眼睛和嘴巴等相關性較弱的區域是不確定的，但是這些區域可以為人臉識別提供更多的判別信息。

6 實驗結果與分析

本文使用第2節介紹的三維顱面數據進行實驗，并采用五折交叉驗證來評估本文方法的有效性。將數據集分為5組，每組包括51個三維顱骨和面皮對。對于每一次折疊，選擇其中一組51個三維顱骨面皮對作為測試集，剩余的204個顱骨和面皮構成訓練集，也就是說，在5次折疊中共有255個三維顱骨作為測試集。對于每一個測試顱骨，所有的255個面皮構成查詢集。整個實驗的正確識別率是5次驗證的平均值。

6.1 整體相關性分析

在PCA中大特征值對應的主成分向量表示訓練數據集的變化模式，而小特征值對應的主成分向量則表示噪聲方向。在這里，本文首先討論了統計形狀模型的方差貢獻率對顱骨識別率的影響。由式（1）所述，模式數量由從累積特征值計算的方差貢獻率確定，當方差貢獻率從95%變化到99.5%時，整體模型的正確識別率如表1所示。從表1中可以看出，隨著方差貢獻率的增加，識別率逐漸提高；在方差貢獻率達到99%時，識別率達到最高；當方差貢獻率再繼續增加時，識別率反而會下降。

表1 不同方差貢獻率的整體模型的正確識別率Tab.1 Correct recognition rates of overall model with different variance contribution rates

對于每次折疊，構建一個整體統計形狀模型，當方差貢獻率為99%時，顱骨和面皮數據投影到低維形狀參數空間中，分別表示為548和356維特征向量。基于LSCDA的全局顱骨識別通過提取基向量，正確識別率達85.2%，由此表明基于顱骨和面皮相關性分析的顱骨識別方法是有效的。進一步分析識別錯誤發現，誤匹配面皮與真實面皮之間的全局形態相似性較高，從而導致了誤匹配。圖6為兩個面皮之間的幾何距離比較，從誤差分析圖可以看出，相似人臉與真實人臉大部分區域的幾何距離都很小。由于兩個面皮的全局形態相似性程度高，因此具有相似的形狀模型參數，使得整體模型產生錯誤匹配。

圖6 相似人臉之間距離的比較Fig.6 Comparison of distances between similar faces

6.2 區域相關性分析

為了驗證區域相關性分析模型，分別為每個區域構建區域統計形狀模型，然后進行區域LSCDA映射，完成匹配。當方差貢獻率為99%時，與基于整體相關性分析的LSCDA顱骨識別相對比，可以得到區域相關性分析的顱骨識別率，如表2所示。

表2 顱骨單一區域識別率Tab.2 Recognition rate of skull single area（%）

通過對比發現，輪廓、前額和整體模型的識別率相當，這是因為輪廓和前額在整體模型相關性測量中起著決定性作用。鼻子主要由軟骨組成，嘴唇與口腔區域的牙齒分離，因此與嘴巴、鼻子區域相比，眼睛的形態對眼眶骨形狀的依賴性更強，從而識別率也較高。

表3 為區域相關性分析的正確識別率，圖7顯示了不同方差貢獻率的局部模型識別率的變化趨勢。

表3 不同方差貢獻率的局部模型的正確識別率Tab.3 Correct recognition rates of local models with dif?ferent variance contribution rates （%）

圖7 不同方差貢獻率的局部模型識別率變化趨勢Fig.7 Variation of local model recognition rate with vari?ance contribution rate

從表3和圖7可以看出，當方差貢獻率增加時，前額和輪廓的識別率變化很快，尤其是前額區域，這是因為較多地舍棄了變化模式中包含的一些特性信息。由圖7可知，當輪廓、前額、眼睛、鼻子和嘴巴區域的方差貢獻率分別為99%，99.5%，98%，97%和98%時，這些區域的識別率最高。在使用基于區域相關性分析進行顱骨識別時，對不同區域選擇不同的方差貢獻率來構建局部統計形狀模型，可以提高顱骨識別的準確率。

此外，通過實驗發現，各個區域權重經過融合之后將達到最優的識別準確率。區域權重如下：

其中i=p，f，e，n，m，對應輪廓、前額、眼睛、鼻子和嘴，ri表示對應區域模型的正確識別率。因此，融合權重可以根據式（22）在實際情況中具體設定。

表4 區域融合權重Tab.4 Regional fusion weights

表4 顯示了本實驗中選擇不同方差貢獻率進行顱骨識別時各區域的融合權重。基于得到的融合權重，通過區域融合識別顱骨。在測試階段，第一次驗證時，51個測試顱骨中出現錯誤識別個數為2，識別準確率為96.0%；第二次驗證時，51個測試顱骨中出現錯誤識別個數為1，識別準確率為98.0%；第三次驗證時，51個測試顱骨中出現錯誤識別個數為3，識別準確率為94.1%；第四次驗證時，51個測試顱骨中出現錯誤識別個數為2，識別準確率為96.0%；第五次驗證時，51個測試顱骨中出現錯誤識別個數為4，識別準確率為92.1%。因此，通過融合5個區域建立的相關性分析模型得到，識別準確率為（96.0%+98.0%+94.1%+96.0%+92.1%）/5=95.2%。

為了驗證本文方法的優越性，我們分別與文獻［16］和文獻［23］中的方法進行對比。3種方法的對比結果如表5所示。

表5 不同方法識別準確率的對比Tab.5 Comparison of recognition accuracy of different methods

從表5可以看出，本文方法的識別準確率最高，文獻［16］方法次之，文獻［23］方法的識別準確率最低。文獻［23］中使用嵌入式隱馬爾可夫模型，雖然充分利用了顱骨的形狀和整體結構信息，但是使用顱骨二維圖像，丟失了顱骨的空間信息；而文獻［16］和本文均使用顱骨三維模型，并且兩者都是利用顱骨和面皮之間的形態關系進行識別。文獻［16］首次提出通過分析顱骨和面皮之間的形態相關度來識別未知顱骨，并使用典型關聯分析提取映射后兩組相關性最大的主成分向量，用于匹配顱骨和面皮。與文獻［16］相比，本文采用LSCDA方法分析顱骨和面皮的相關性，通過兩個對變量投影的LSMI測量相關性，并通過線性變換使相關性達到最大。該方法不僅能通過線性投影提取兩個變量的最相關信息，還可以消除無關信息，因此識別準確率更高。

7 結 論

本文提出了一種基于三維顱骨與三維人臉形態相關性度量的未知顱骨識別技術。使用三維顱骨數據作為探針，三維人臉數據作為匹配庫，通過探針和匹配庫之間的相關性度量，將未知顱骨與已配準的三維人臉進行匹配，從而實現顱骨身份識別。首先，構建顱骨和面皮的統計形狀模型，將高維顱骨和面皮映射到低維形狀參數空間；然后，基于LSCDA提取顱骨和面皮的主相關信息，構建整體相關性分析模型；最后，計算顱骨與面皮庫中每個面皮之間的匹配分數，具有最高匹配分數的面皮即為正確的識別結果。為了提高匹配的準確性，將顱骨和面皮劃分為5個生理特征區域，分別為前額、眼睛、鼻子、嘴巴和輪廓，建立5個區域的相關性分析模型，并通過區域融合進行顱骨識別。與顱像重合和顱面復原相比，本文方法不需要顱骨和面皮之間的精確關系，就可以實現未知顱骨身份的確定。本研究證明了顱骨和面皮之間確實存在著較為緊密的聯系，為顱像重合和顱面復原的研究提供了理論支持。此外，在提取顱骨和面皮之間的關系方面，基于區域的策略要優于整體策略，這也為基于區域的顱面復原研究提供了參考。

隨著基于圖像的三維人臉建模技術和三維數據采集技術的發展，三維顱面模型獲取變得越來越容易，擴展三維顱面數據集從而進一步提升本文方法的識別率，使它能夠更為廣泛地用于實際案例中，是下一步的研究重點。