泛函回歸代理及條件期望配準的機械擺動測量

鄭思凡，陳平平，蘇凱雄，吳永春

（1. 福州大學物理與信息工程學院，福建福州350116；2. 泉州黎明職業大學智能制造工程學院，福建泉州362000；3. 數字電視智能化技術國家地方聯合工程研究中心，福建福州350116）

1 引 言

機械擺動測量在精密制造與智能制造的精密閉環控制具有廣泛的應用，如五軸數控機床擺動式回轉軸角位置測量［1］，固體發動機柔性擺動噴管擺角測量［2］，起重機吊鉤擺動角測量［3］，偏擺儀指針偏擺檢測等。相應地，工程上也出現了各種基于光電技術的高精度測角儀［4］，但一般對測量相對位置與光照環境要求較高，測量對象單一，無法滿足廣泛應用需求，因而如何運用圖像運動感知技術對上述機械擺件進行逐幀提取分割，對擺角進行視覺測量則具有很大的工業實際意義。

隨著圖像采集及視頻壓縮傳輸軟硬件技術的發展，高幀率，高分辨率自帶光學防抖的CMOS 與CCD 攝像機的出現為高速的機械運動非接觸式視覺測量提供了硬件實現條件，如三星6 400 萬像素ISOCELL Bright GW1 CMOS 感光元，其單個像素面積僅0.8 μm，支持Tetracell 像素合成技術和Remosaic 還原，可直接輸出480 FPS 的1 080 pixel 圖像［5］。當前國內也有關于三維運動圖像測量分析系統的相關軟件產品，如北京現代富博科技有限公司的MIAS3.0［6］，但是產品需要手工標定跟蹤像素，在自動化方面存在提升空間。

當前 SLAM［7］（Simultaneous Localization and Mapping）地圖重建技術及其衍生品視覺里程計［8］（Visual Odometry）技術的發展為視覺測速提供了許多成熟的算法框架，二者的核心框架均在于通過相鄰幀對應目標點云配準完成對相機的位姿運動參數的估計。對于SLAM 技術的視覺測速，當前也有比較成熟的研究成果，如圖優化工具g2o（General Graph Optimization）［9］所集成的基于指數映射的李代數擾動模型的集束調整（Bundle Adjustment）［10］算法。該算法將有約束的旋轉矩陣李群目標函數，轉化為無約束的李代數矢量優化問題，從而可以用普通牛頓迭代法求解。但上述研究僅局限于相機本身運動估計，如何將點云配準結合運動分割，用于工業上機械剛體點云的運動參數估計與測量則有較大研究空間。

為實現對弱結構張量特征點光流的運動分割，T. Brox，J. Malik 等提出通過對變分稠密光流軌跡時空相似度拉氏譜圖譜聚類的運動分割算法［11-13］。此類算法是采用ε-NN 或K-NN 近鄰矩陣來構建，其超參數ε 或K的取值由數據集訓練后統一設定，無法根據實際光流稠密度不同自適應調整。為此Roberto Tron，Ren′e Vidal 等根據SFM（Structure From Motion）［14］三維重建理論為基礎，將運動分割轉化為剛體光流軌跡，聚類到2～4 維的低維超平面子空間里來完成運動分割。代表性算法有廣義主成分分析、局部子空間相似度聚類、低秩分解、低秩子空間聚類、多段學習等［15］。上述算法充分挖掘了剛體光流軌跡蘊含的子空間結構信息，克服了光流空間分散及密度不均造成的分割尺度誤差。但是各子空間的維度及個數參數仍然需要事先指定。有鑒于此，von 等于2013 年提出了稀疏子空間聚類［16］，將傳統的基于ε-NN 或K-NN 近鄰矩陣替換為L1范數圖。其中各光流軌跡為圖的節點，邊的權重則是以其他軌跡作為字典基進行稀疏學習后的自線性表達系數組成。這樣利用高維數據的稀疏性使得同一流形內部取得不為零的表達系數，而不同流形的節點其代表的軌跡互表達系數為零。從而以特征選擇的方式自適應表達了光流軌跡密度不均的流形結構，解決了ε-NN 或K 自適應調整問題。論文［17-19］等通過擴展修正不同的字典基表達系數的正則約束，改進SSC 聚類的類內一致性及稀疏性。但是，上述算法均是針對塊狀移動物體進行分割提取，在機械擺動的視覺運動分割中，因背景復雜性及運動對象的材質（如金屬件），存在大量弱結構張量特征點及被中斷為多段零散的光流軌跡，利用傳統的時空相似度譜聚類易形成碎片化的塊狀軌跡群，難以一次性完成所有軌跡的子空間聚類，而是需要分段多次進行，不僅計算量大且在光流遮擋處密度過于稀疏影響配準精度，因此SSC在機械擺動的運動分割方面存在較大的局限性。

為此本文提出了一種以高斯隨機過程泛函回歸代理的弱結構張量機械擺動測量算法。首先獲取機械擺動件LDOF 變分光流；其次通過構建新的鄰接矩陣與相似度度量獲取過分割狀態的弧狀軌跡群；再者以時間為參變量分別對軌跡群x，y坐標進行高斯回歸，學習出該軌跡群的平均軌跡，此平均軌跡代理該軌跡群作為稀疏子空間聚類的種子樣本，一次性完成聚類；這樣即克服了碎片化導致的分段多次聚類，也通過SSC 的模型選擇得到正確的運動物體數目從而克服了過分割的局限性，最后對每幀運動分割后的機械擺件特征點云匹配關系建模為各匹配對應點為中心的GMM 分布，將該分布的數據最大似然參數估計作為運動位姿矩陣估計值，通過提取旋轉矩陣歐拉角完成擺角計算及測量。為證明所提算法的有效性，本文以理想雙搖桿模型的6 種不同照度下的車輛刮水器擺桿為對象，結合基于機器視覺的車輛安全自動化日檢工程項目于安檢現場提取了30 組數據實驗。分析比較了所提算法與當前算法對擺桿的光流軌跡聚類角位移誤差。實驗結果表明：本算法能完整學習出等長軌跡并在滿足一定的環境亮度情況下，角位移測算值與人為標定的回歸值均方差小于10%。同時解決了當前算法存在的過分割，稠密度不均導致的配準誤差及多次分段與迭代的運算量問題。

本文首先介紹測量對象四連桿機構動力學原理及其理想角位移曲線；其次介紹稠密變分光流種類及被測刮刷總成運動光流軌跡的獲?。坏谌糠纸榻B傳統塊狀譜聚類運動分割存在的過分割及碎片化兩個局限性現象及原因并通過修改傳統塊狀軌跡群的相似度度量與近鄰矩陣構建方法引出本文的弧狀超流素軌跡群概念及實現；第四部分介紹超流素泛函回歸代理參數的訓練及推斷；第五部分介紹算法的聚類性能比較；第六部分為算法角位移測量誤差效果比較；最后為總結與展望。

2 四連桿機構動力學原理

當前客運車輛較常用的刮水器總成結構由兩部分四連桿組成，為研究方便，可在solidworks建模如圖1 所示。圖中L1為左搖桿、L2為機架、L3為右搖桿、L4為耦合連桿，L1～4構成雙搖桿機構，L5為雨刮擺桿，L6為輸入連桿，L7為電機驅動曲柄，三者構成主刮曲柄擺桿系統；通過設定左右擺桿長度起到控制輸入動桿與輸出桿速度和相位的作用。

為了構建仿真算例，觀察左右搖桿角位移變化曲線，取L1=L2=75 cm，L3=L4=35 cm 并設仿真時間設置為5 s，電機轉速r=30 rmps，動畫幀率為30 FPS。算例求解完成后，得左搖桿的角速度變化曲線分別如圖2 所示。

圖1 刮水器連桿結構Fig.1 Wiper link structure

圖2 左搖桿角位移變化曲線Fig.2 Angular velocity curve of left Rocker

由圖可見四連桿刮刷運動過程呈周期性變化，且在啟動與結束處速度存在收斂的過渡過程。若按照圖像逐幀分割進行還原，則會呈現如圖3 圖樣。

圖3 刮桿逐幀運動分割圖樣Fig.3 Fame-by-frame motion segmentation of Wiper

根據矢量運算原理，同樣可以得到右搖桿的角位移曲線，因為不同擺桿其在SFM 模型對應的運動矩陣平移分量不同，所以左右擺桿處于不同仿射空間，在理想情況下，可以直接使用子空間聚類對滿幅等長軌跡一次性完成聚類分割，相應地，各幀分割形成的點云經過配準并提取位姿旋轉矩陣的歐拉角則可一次性求得角位移曲線，但在復雜背景下將呈現為多段碎片化分割區域，無法一次完成分割。

3 稠密變分光流定義及計算

為了解決傳統光流亮度一致性方程梯度逆陣病態帶來的孔徑問題，Horn-Schunck 等提出了數據項與正則項聯合變分的稠密光流［11］。其中數據項主要完成亮度一致性約束，具體形式如式（1）所示：

其中：T(x，y)為參考圖像，I(x，y)為當前圖像。u(x，y)，v(x，y)是圖像上每一點的偏移量，φ(x)為代價函數，可取絕對值，二次函數等一些非負對稱且正半軸單調遞增的函數。

為了對純色區域的弱結構張量像素實現幀間稠密匹配，在數據項加入一個平滑項，使得相鄰兩個特征點的偏移量相差不能太大：

為研究方便，本文采用華睿A3200CU000 面陣高速工業相機進行采集。在幀率為120 FPS，分辨率720 pixel 的設定下，在安檢臺采集客車雨刮器上升動作視頻180 幀分析。當設定空間采樣率為8、像素灰度結構張量閾值為0.4 的情況下，可以采用前述OpenCV3.4 提供的BroxOptical-Flow 類完成實時LDOF 光流計算并生成相應的flo 文件，為了保證軌跡的魯棒性，在各相鄰幀的光流連線形成軌跡前，必須對各光流進行前向與后向方向一致性檢查［14］，這里設定一致性系數0.06，一致性偏移閾值為4，可得雨刮的軌跡分布如圖4 所示。

圖4 刮水器總的光流軌跡分布Fig.4 Total optical flow trajectory distribution of wiper

4 當前運動分割局限性及原因分析

這里先考慮以軌跡時空相似度譜聚類前述光流軌跡的運動分割算法［11］，其中軌跡的相似度以共有幀的位移及速度范數距離的RBF 核來構建，此相似度進一步作為ε-NN 近鄰矩陣的邊權重A(tri，trj)即：

其中：wv(tri，trj)為軌跡tri與軌跡trj在幀t的速度相似度，定義如下：

其中：Lrw為隨機矩陣，其元素pij=wij di可以看成是節點i通過隨機行走到達節點j的概率［21］，其中wij為節點i與節點j之間相似度權值，di為連接到節點i的邊總數。這樣通過特征分解可得特征矢量如式（8）：

最后取特征值小于0.15 對應的特征矢量矩陣，構成行矢量空間進行k-mean 聚類。圖5 展示了特征矢量對應的特征值從小到大排列。

圖5 SC 特征值排列Fig.5 Eigenvalue index of Spectral Clustering

由圖可見，前5 個特征矢量特征值小于0.15，因此直接對前述光流軌跡譜聚類得到如下的聚類結果，如圖6 所示。

圖6 譜聚類運動分割結果Fig.6 Spectral clustering motion segmentation results

圖7 拉氏特征映射低維嵌入分布圖Fig.7 Low-dimensional embedding of Laplace Eigenmaps

相應地，圖7 展示了軌跡相似度在1～5 各維度的低維嵌入分布。由上述可見，直接對軌跡相似度矩陣譜聚類能夠得到較好的類內一致性。相似度矩陣也顯示出塊狀結構，各低維嵌入的同類數據也保持相鄰的拓撲關系。

但是因其構建相似度矩陣所采用統一固定的鄰域ε 參數，與尺度參數η無法反映同一擺桿不同段光流線速度存在的差異，導致其受到軌跡密度與速度分布不均的影響，將兩根擺桿過分割為5 根。因此，算法沒有利用同一剛體的軌跡處于同一子空間的線性表達關系的內在信息，得到準確的聚類超參數K，出現了過分割現象。

下面介紹分析此算法的另一個局限性即碎片化現象，為此在前述得到軌跡相似度拉氏矩陣的情況下，直接取K=24 得到超流素并以顏色區分，如圖8 所示（彩圖見期刊電子版）。由圖可以看出，傳統譜聚類運動分割在K 取較大值時，因為遮擋及軌跡的中斷將使得分割結果碎片化，無法理想將線速度相近的軌跡聚類為同一軌跡群。

圖8 SC 軌跡聚類碎片化現象Fig.8 fragmentation phenomenon caused by SC

下面從相似度的距離度量定義及鄰接矩陣的構建方式兩個方面來說明上述聚類形成碎片化的原因并提出弧狀軌跡群的分割方法。

第一，在構建鄰接矩陣方面，以圖9 為例，參考文獻［11］提出以軌跡公共時間幀里兩條軌跡相應X，Y位移差最大值作為ε，則圖9 中軌跡A與B，C與D兩兩之間滿足條件構成鄰接關系，而A與C，B與D因為中間出現了軌跡中斷，沒有公共時間幀，于是鄰接關系為零，不構成鄰接關系，在后續的聚類結果里被聚類為不同類別，E，F，G，H同樣如此，從而出現了碎片化。

圖9 SC 軌跡聚類碎片化現象Fig.9 Fragmentation phenomenon caused by SC

于是本文在原鄰接關系的基礎上，將共有時間幀推廣到在軌跡終點處且時間差Δt在10 幀以內且軌跡起點在ε位移半徑以內的另一條軌跡建立軌跡鄰接關系，這樣就使得A與C，B與D分別構成鄰接關系為后面將他們聚類為同一類別創造了基礎。

其次，從公式（4）定義的時空相似度來看，圖10 中A與B的相似度大于A與C和B與D，同理E與F相似度大于E與G和F與H（塊狀軌跡群相似度＞弧狀軌跡群相似度），因此在K值較小時A與C，B與D還有可能通過某條中間橋接軌跡分別聚為同一類，但是在K取較大時A與B，C與D各自聚為同一類形成塊狀軌跡群，而AC，BD則聚為不同軌跡群，從而造成K值較大時出現碎片化。

圖10 稀疏子空間運動分割結果Fig.10 Results of Sparse subspace clustering

為了使A與C聚為同一類形成弧狀軌跡群，而A與B形成不同軌跡群，考慮到同一半徑處的像素運動軌跡不管被中斷為幾段，各段的平均曲率或線速度絕對值相似，因此這里將軌跡相似度距離度量修改為各軌跡的平均曲率，使得AC的相似度大于AB的相似度，同理，EG相似度大于EF相似度，這樣在K 取較大值時準確構建弧狀軌跡群，為后面的GP回歸及Lasso 回歸提供了正確的種子樣本。

綜合上述，這里ε取8，K取48，將各段軌跡經立方樣條spline 擬合后對函數結構求一二階導數所得平均曲率作為近鄰矩陣的邊權重A(tri，trj)后再譜聚類可得48 個弧狀軌跡群如圖11。

由圖可見，上述分割方式解決了直接譜聚類在弧向方向的碎片化問題。各弧狀軌跡群在徑向方向則形成高維的線性表達關系，這一特征可以用子空間聚類方式再進一步聚類，從而可以以無監督的機器學習方式經過兩層回歸對原始零散的軌跡群進行完整聚類。

圖11 弧狀軌跡群聚類對碎片化的改進結果Fig.11 Improved results by hierarchical clustering

5 超流素的定義及泛函回歸參數的學習與樣本推斷

5.1 超流素定義

超流素是在軌跡元（Tracklet）［20-22］基礎上經過一層聚類的軌跡群，代表一段視頻中光流軌跡時空相似度最小分割單位，將其作為視頻分割的中層表達形式由文獻［23］提出。為了與相關文獻保持名稱一致，下面以超流素的概念作為本文軌跡群的定義。

5.2 高斯回歸定義及新值統計推斷

高斯回歸代理作為一種泛函回歸［24］，是基于假定采樣點之間的距離或方向可存在變化的空間相關性基礎上，將此相關性以統計模型進行學習擬合，利用統計推斷對新測試值做出最佳預測。其本質可以看為一種貝葉斯推斷，因此如果數據X似然部分采用高斯似然作為核函數，則可以利用共軛先驗特性獲取新數據協方差∑*及平均值μ*的閉式解。由此通過貝葉斯推斷可得新測試樣本X*的泛函分布f*如下：

其中：X為樣本數據，y為樣本對應的訓練值，K為數據核函數矩陣，此矩陣可以分為樣本數據與測試數據兩部分，分別以腳標*與y標出，定義為：K=κ(X，X)；K*=κ(X，X*)；K**=κ(X*，X*)。

將前述超流素的坐標位移容限15 看為高斯回歸的數據噪聲協方差，則可得核函數矩陣元素對應核函數如下：

5.3 高斯回歸超參數學習

設θ=(?，σ2f)則可得上述MAP 估計對超參數的梯度如下：

令此梯度為零可得超參數閉式解，而上述高斯過程對新數據的MAP 估計等價于無限維的RKHS 空間二次范數正則化泛函的線性回歸，根據里斯表示定理（Riesz Representation Theorem）可以由數據樣本張成的有限維空間表示如下：

其中：α?y為各已知樣本與測試點距離的加權系數是測試點與各樣本的核函數矢量，其元素為κ(xi，x*)，即測試數據與樣本映射到RKHS 希爾伯特空間的內積距離。

為反映各軌跡時間方面的同步與次序信息，這里對超流素軌跡的回歸以時間為中間參變量，分別對軌跡群的橫縱坐標進行插值。如前述取超參數=15，=20，?=480 為初值，利用公式（14）對每個超流素的超參數重新學習后，由公式（9）對所有時間段橫坐標值重新推斷，則可得回歸的平均橫坐標值如圖12 所示。

圖12 超流素橫坐標隨時間的高斯回歸Fig.12 Gaussian regression of Superfloxel’s x coordinate

相應地，y坐標的回歸平均軌跡如圖13所示。

圖13 超流素縱坐標隨時間的高斯回歸Fig.13 Gaussian regression of Superfloxel’s y coordinate

5.4 算法實現及時間復雜度

為了考量算法時間復雜度，這里假定樣本數為r。由公式（17）可見，為了算法穩定性，不宜直接對Ky求逆，這里采用對Ky做Cholesky 分解為Ky=LLT則算法具體實現如下：

其中步驟1 的Cholesky 分解時間復雜度為O(r3)，步驟2 求α時間復雜度為O(r2)，步驟3的矩陣乘法需要O(r)，步驟4 和步驟5 運算復雜度為O(r2)。為減少運算量，在實踐中只需要求出步驟3 所表示的坐標值為回歸軌跡坐標，步驟4～6 所求的回歸置信度可以省去。為進一步降低運算量，這里采用子集回歸的稀疏高斯回歸（SoD）［25］算法對數據核Gram 矩陣K分塊如下：

其中：Km m由m個引導點（Inducing Point）構成活動集I（Active Set）對應的分塊矩陣，則經過SVD分解后原來核Gram 矩陣K可用其前最大m個特征值及相應特征向量構成的m階降秩矩陣近似，其定義如下：

相應地，公式（17）核函數由測試點與子集樣本的核函數線性組合近似表達為：

其中：j為引導點序號，cij為線性表達系數，其構成矩陣C由式（19）確定：

其閉式解為：

Copt=Knm K-1mm，這樣新的子集回歸對應的平均值與協方差近似值如下：

這樣總的時間復雜度降為O(m2r)，這里取m值為幀數，并分別對前述48 個超流素進行坐標插值回歸，可得回歸后的軌跡形狀如圖14 所示。

圖14 超流素高斯回歸平均軌跡Fig.14 Gaussian regression of Superfloxel

6 聚類性能比較

將上述回歸的軌跡作為種子樣本，假設軌跡長度為F幀，軌跡數為N，將軌跡每個點坐標按照先橫坐標后縱坐標次序排列為矩陣Y的每一列。

則SSC 通過如下目標函數求解線性自表達系數矩陣X：

其中：E為字典重構誤差，‖E‖為誤差E 的Frobenius 范數，當假定重構誤差為高斯分布時此范數反映誤差代價，λ為權衡參數，平衡軌跡誤差與解的稀疏性及泛函的凸性，因為本文是針對回歸后的軌跡聚類，故誤差相比理想情況下的軌跡大，故對λ取較大的值0.1 以克服回歸帶來的誤差。

通過引入拉格朗日乘子Δ ∈RN×N及懲罰系數μ可將上式轉化為無約束目標函數如下：

文獻［16］通過ADMM（交替方向乘子法）算法輪流優化變量E，X，Δ 及μ將上述目標函數轉化為兩個子優化問題如下：

其中：θ=‖Y‖?X(Xk，Ek，Δk，μk)為目標函數L對X的偏導函數，k為迭代次數，則公式（25）變為標準LASSO 問題，可以采用通用迭代收縮閾值算法（Generalized Iterated Shrinkage Algorithm）［26］求解。另外通過引入輔助變量，即假設B=YX(k)-Y-Δ(k)/μ(k)，則公式（27）可改寫為：

一樣轉化為針對每個數據εj的標準LASSO問題從而同上求解。由式（26），式（27）可見ADMM 算法每進行一次迭代都需要求N 階逆陣，其運算復雜度為O(N3)，在N取1 000，m取100，超流素K取50 時，則r為2 000，可見ADMM 的單步復雜度O(N3)大于稀疏高斯回歸O(m2r)，但是高斯回歸不需要迭代，因此本文算法可減少一階運算量。

最后將所求的矩陣X元素作為L1 范數圖節點間的權值并預先取K=2 進行譜聚類可得如圖15 所示結果。

圖15 代理軌跡的SSC 聚類結果Fig.15 SSC clustering result of the surrogate trajectories

由圖可見，回歸軌跡代理的SSC 聚類（RSSSC）準確標識了左右擺桿，但上述的聚類數目為事先取定，依據前述譜圖理論，聚類個數的值可由特征值為零的個數決定，相應圖5，這里將特征矢量對應的特征值從小到大排列如圖16所示。

相應的48 條種子軌跡相似度拉氏圖矩陣及1～2 維的拉氏映射嵌入分別如圖17～圖18 所示。由圖可見由線性自表達系數構成L1 范數圖構成的軌跡相似度圖相比原來的相似度矩陣，在特征值分布和低維嵌入數據的類間距離，更接近超參數K真實值2。

圖16 回歸軌跡代理的SSC 特征值排列Fig.16 Eigenvalue index of the regression surrogate SSC

從而說明了經過超流素回歸軌跡代理的SSC 聚類在K-mean 聚類前可以通過矩陣干擾理論或模型選擇等方法獲取最佳聚類數目，克服了前述譜聚類過分割缺陷，同時這一優勢使得本算法可以從光流軌跡中發現真正擺桿數目，在機械故障診斷中（如檢測因疲勞而停止擺動的賈卡針數目）有一定實際意義。

圖17 回歸軌跡代理的SSC 的相似度矩陣Fig.17 Affinity matrix of the regression surrogate SSC

圖18 回歸軌跡代理的SSC 低維嵌入分布圖Fig.18 Low-dimensional embedding of RS-SSC

7 角位移曲線的測量及比較

除了類內一致性性能，下面從聚類稠密度提高的優勢說明經過代理回歸后SSC 聚類相比直接稀疏子空間聚類可提高后續擺件分割點云配準得到的角位移精度。

7.1 非種子軌跡的分割及點云提取

代理回歸SSC 得到各幀種子軌跡聚類標簽后，將種子軌跡標簽在超流素內部擴散至回歸前的所有軌跡樣本，使得運動分割后各種分割點云稠密度達到光流采樣的密度. 在第125 幀將擴散后非種子軌跡的像素點云及所屬聚類標簽用顏色標出如圖19 所示（彩圖見期刊電子版）。

圖19 回歸軌跡代理的SSC 聚類運動分割結果Fig.19 Motion segmentation results of the RS-SSC

7.2 最大條件期望點云配準及角位移測量

由前述可以看出，因為光流阻斷，各幀分割的運動擺桿點云個數不同，不能直接采用ICP 等配準算法求位姿，故這里采用論文［27］基于條件期望的最大似然參數估計來求點云間的位姿變換矩陣，為減少匹配離群點個數，這里取個數最大的點云為模型數據，表示為X={Xi}1≤i≤n，而待匹配各幀點云為Y={Yj}1≤j≤m，模型經過位姿參數配準后的點云為：μ(Xi；Θ)=RXi+t，Θ：={R，t}。假定數據點云每個點以其匹配的相應模型點云中的點為中心的高斯分布，則配準問題轉化為以匹配關系矩陣Z為隱含變量并以觀察到的數據為條件的GMM 最大似然參數估計問題，被估計的參數為Ψ={Θ，Σ1，…，Σn}，其中Σn為各點分布協方差，Θ：={R，t}為位姿參數。

估計的目標函數即條件數據最大似然通過隱含變量Z可分解為：

其中，對應各數據的隱含變量矩陣為Z={Zj}，1≤j≤m，矩陣元素取值為該點所匹配對應模型點的序號，可定義該隨機變量取值分布如下：

αji=P(Zj=i∣Yj)，引入克羅內克離散取值定義后可將上述隨機變量取值視為條件期望：

從而可以采用E-M 架構對目標函數求解，其中E 步為求解該匹配關系隱含變量的后驗估計如式（31）：

其中?3D為失配離群點集，定義如下：

其次在M 步分別對另外兩個運動參數估計如下：

為保證全局最優，求解前，先對擺桿點云進行RANSAC 擬合得到的轉角結合平移分量為零構成的旋轉矩陣作為參數R初值。

7.3 角位移測量誤差比較

為方便與其他兩種算法比較，將視頻分割為4 段長度為35 的視頻共140 幀數據，分別進行SSC 運動分割以保證每段有一定數量的種子軌跡. 為簡單起見，點云所在平面為XOY平面，Z坐標為0，采用右手系確定Z軸方向后，通過反正切公式提取旋轉矩陣R的Z軸歐拉角為擺桿擺角測量值，基準值采用人工標定，標定方法為在matlab2017 里figure 對象的鼠標點擊回調函數里實現了取點及對雨刮根部主擺杠RANSAC 擬合與反正切求角，最后可得三種算法各幀擺角測量結果如圖20 所示（彩圖見期刊電子版）。圖中不同顏色點為SSC，SC，RS-SSC 各算法的測量值，曲線為測量值相應的高斯回歸平均值，圖中SSC，SC，RS-SSC 對標定回歸曲線的均方差分別為16.96，11.71，8.79。因此，RS-SSC 最接近標定值，這是因為SSC 只能選長度大于35 的光流，而RS-SSC 可以將回歸前超流素所包含的所有的被阻斷的光流軌跡經過標簽擴散后一同參與點云配準，其像素的稠密度大于前者，而SC 在擺桿接近處軌跡會出現混疊導致部分異常點降低了一部分精度。

圖20 SC，SSC，RS-SSC 各算法的擺角曲線檢測結果Fig.20 Angle curve detection of SC，SSC，RS-SSC

為了在工程現場考量算法的測量精度，結合客運車輛的機器視覺自動安檢項目，在現場采集了6 種不同環境亮度同一車型客車雨刮擺桿的5組角位移曲線作為測量對象，將人為標定回歸值作為基準則測量誤差平均值如下：

由表1 可知，隨著環境照度的提高，角位移檢測精確率均有所提高。這是因為環境亮度的提高增加了擺桿前景對比度，這樣減少了軌跡中斷現象使得種子軌跡數目增加，同時因前景特征點結構張量增加也增加了總體軌跡的數目與稠密度?？傮w上看，當環境亮度大于700 lux 時，代理回歸的SSC 的擺角測量誤差不超過10%，工程上可以作為此安檢項目中刮桿卡扣及其驅動電機的速度是否正常無故障的判斷依據。

表1 各算法在不同照度下的擺角檢測誤差Tab.1 Angle detection error under different algorithm and illuminance

8 結 論

本文在利用圖像視覺對單色的機械擺件的運動參數測量中，針對當前兩類運動分割的缺陷，提出了一種通過超流素回歸平均軌跡代理的SSC 聚類的運動分割算法。該算法因采用軌跡自表達系數構成的L1 范數圖取代傳統的變分光流的時空相似度拉式譜圖進行聚類，從而克服了選擇超參數ε、K的困難并以無監督機器學習方式通過子空間結構的信息挖掘出運動物體數目K，為運動分割用于機械擺動故障檢測提供了算法依據。

另一方面，本算法利用高斯泛函回歸學習出因背景復雜而被阻斷的軌跡并代理超流素軌跡群一次性完成SSC 運動分割，相比不回歸直接分時間片段的SSC 算法降低了運算量，增加了子空間聚類種子樣本的稠密度，結合SLAM 點云配準技術提高了逐幀角位移測量精度. 所提出的算法提供了一個將圖像學與圖形學有機結合的機械運動測量框架，具有廣闊的研究空間及工業機器視覺測量應用價值。