唐文秀，林 虎，薛 梓，秦海濛，鐵咪咪

(1.北京信息科技大學(xué) 儀器科學(xué)與光電工程學(xué)院，北京 100192；2.中國(guó)計(jì)量科學(xué)研究院，北京 100029；3.天津大學(xué) 精密儀器與光電子工程學(xué)院，天津 300072)

1 引 言

隨著大型裝備制造業(yè)的發(fā)展，在船舶、軌道交通、航空航天等領(lǐng)域都對(duì)大尺寸測(cè)量技術(shù)及測(cè)量設(shè)備提出了越來(lái)越高的要求。激光跟蹤儀是大尺寸測(cè)量中常用的測(cè)量設(shè)備，可用于精密幾何測(cè)量和三維運(yùn)動(dòng)軌跡跟蹤，具有測(cè)量范圍大、速度快、精度高等優(yōu)點(diǎn)[1]。近些年，各種新型的激光跟蹤儀相繼被研制出來(lái)，其中德國(guó)物理技術(shù)研究院和英國(guó)國(guó)家物理實(shí)驗(yàn)室共同研制了一種新型激光跟蹤干涉儀[2]，采用球度誤差小于50 nm的標(biāo)準(zhǔn)球作為參考反射鏡，該標(biāo)準(zhǔn)球固定在熱膨脹系數(shù)極小的殷鋼支柱上，干涉測(cè)量系統(tǒng)僅圍繞該固定標(biāo)準(zhǔn)球轉(zhuǎn)動(dòng)，因此機(jī)械結(jié)構(gòu)對(duì)于精度的影響可以降低到很小，并且只利用干涉測(cè)量長(zhǎng)度進(jìn)行空間坐標(biāo)的解算，避免了角度測(cè)量引入的誤差，具有很高的精度[3,4]。其中基于多邊法的四路激光跟蹤三維坐標(biāo)測(cè)量技術(shù)是近年的研究熱點(diǎn)[5～9]。四路激光跟蹤三維坐標(biāo)測(cè)量系統(tǒng)在實(shí)際應(yīng)用前首先需要進(jìn)行系統(tǒng)多參量自標(biāo)定，確定出跟蹤干涉儀的基站坐標(biāo)。系統(tǒng)多參量的自標(biāo)定過(guò)程本質(zhì)是一個(gè)非線性最優(yōu)化問(wèn)題，求解非線性最優(yōu)化問(wèn)題常用的方法有高斯牛頓法、Levenberg-Marquart法、信賴域法等等[10]。高斯牛頓法由于避免了計(jì)算大規(guī)模的Hessian矩陣所帶來(lái)的不便，極大地簡(jiǎn)化了運(yùn)算量，實(shí)際應(yīng)用廣泛。因此，本文通過(guò)C#與MATLAB混合編程、C#語(yǔ)言編程兩種方式實(shí)現(xiàn)高斯牛頓法求解激光跟蹤干涉儀基站空間坐標(biāo)值。

2 激光跟蹤干涉儀基站空間坐標(biāo)標(biāo)定原理

如圖1所示為激光跟蹤干涉儀基站空間坐標(biāo)標(biāo)定原理圖，對(duì)于多臺(tái)激光跟蹤干涉儀組成三維坐標(biāo)測(cè)量系統(tǒng)，其每個(gè)基站坐標(biāo)標(biāo)定原理類似。在全局坐標(biāo)系統(tǒng)中，激光跟蹤干涉儀第j個(gè)基站坐標(biāo)為(x0j,y0j,z0j)(j=1,…,4)，對(duì)于任意空間目標(biāo)點(diǎn)Pi(xi,yi,zi)，激光跟蹤干涉儀可以測(cè)出其相對(duì)初始目標(biāo)點(diǎn)P0(x0,y0,z0)的相對(duì)長(zhǎng)度變化量Δlij=lij-l0j，這里需要求解基站坐標(biāo)及到初始目標(biāo)點(diǎn)的長(zhǎng)度值l0j。為了對(duì)上述參量進(jìn)行標(biāo)定和求解，需要在全局坐標(biāo)系統(tǒng)中設(shè)置一系列已知的目標(biāo)點(diǎn)Pi(xi,yi,zi)(i=1,2,…,m)，利用激光跟蹤干涉儀對(duì)上述目標(biāo)點(diǎn)進(jìn)行依次測(cè)量，根據(jù)基站坐標(biāo)、初始點(diǎn)坐標(biāo)及目標(biāo)點(diǎn)坐標(biāo)的空間幾何關(guān)系可以建立如下方程組：

(1)

在上述方程組中，需要求解的未知參量為(x0j,y0j,z0j,l0j)，為了獲取更高的解算精度，需要更多的冗余方程，即目標(biāo)點(diǎn)的數(shù)量要大于4個(gè)，此時(shí)方程個(gè)數(shù)大于待求參量個(gè)數(shù)，構(gòu)成超定方程組。超定方程組沒有通常意義下的確定解，但是可通過(guò)數(shù)值算法解算出近似最優(yōu)解。因此，求解激光跟蹤干涉儀基站空間坐標(biāo)問(wèn)題轉(zhuǎn)化為對(duì)非線性超定方程組的最優(yōu)化求解，如何進(jìn)行最優(yōu)化求解并且通過(guò)C#編程實(shí)現(xiàn)成為本文研究的關(guān)鍵內(nèi)容。

圖1 激光跟蹤干涉儀基站空間坐標(biāo)標(biāo)定原理圖

3 高斯牛頓算法求解原理

(2)

進(jìn)一步可推導(dǎo)為：

(3)

上述非線性超定方程組具有一定的平方和結(jié)構(gòu)特殊性，屬于小殘量問(wèn)題，因此，可選用高斯牛頓法求解非線性超定方程組的最優(yōu)解。該算法對(duì)于非線性擬合具有較好的擬合精度，收斂速度較快，對(duì)初值依賴強(qiáng)，必須合理選取初值，它通過(guò)舍棄Hessian矩陣的二階偏導(dǎo)數(shù)實(shí)現(xiàn)，可用于數(shù)據(jù)擬合、參數(shù)估計(jì)和函數(shù)估計(jì)等方面，并且在實(shí)際應(yīng)用中，常用于目標(biāo)定位、單點(diǎn)定位等坐標(biāo)求解問(wèn)題[11～13]。

ri=(b4+Δlij)-

(4)

求解該最小二乘問(wèn)題即讓目標(biāo)函數(shù)最小，數(shù)學(xué)表達(dá)式為：

m>4

(5)

將f(b)按泰勒公式在bk處展開并略去高階項(xiàng)：

f(b)=f(bk)+g(bk)T(b-bk)+

(6)

函數(shù)ri(b)二次連續(xù)可微，目標(biāo)函數(shù)的梯度g(b)與Hessian矩陣G(b)分別為：

(7)

=J(b)TJ(b)+s(b)

(8)

根據(jù)高斯牛頓迭代法，對(duì)ri分別求b1、b2、b3、b4的偏導(dǎo)，其中初值為b0=[b01,b02,b03,b04]，[b01,b02,b03]為根據(jù)空間位置粗略估計(jì)激光跟蹤干涉儀初始坐標(biāo)，b04為激光跟蹤干涉儀到初始動(dòng)點(diǎn)的迭代初始值，構(gòu)成的雅克比矩陣J如下所示：

(9)

將g(b)與G(b)代入式(6)可得：

(10)

忽略Hessian矩陣的二階偏導(dǎo)s(bk)后，根據(jù)牛頓法推導(dǎo)可得高斯牛頓迭代公式為：

bk+1=bk-(JTJ)-1JTr(bk)

(11)

4 C#與MATLAB混合編程關(guān)鍵技術(shù)

為了實(shí)現(xiàn)最優(yōu)解的快速迭代求解，可直接在科學(xué)計(jì)算軟件MATLAB中調(diào)用相關(guān)函數(shù)。MATLAB具有豐富的函數(shù)庫(kù)，專門用于數(shù)值分析，調(diào)試簡(jiǎn)單，但不擅長(zhǎng)于后續(xù)測(cè)控軟件開發(fā)。為了求解出激光跟蹤干涉儀基站空間坐標(biāo)，并且顯示到用戶界面中，方便用戶進(jìn)行數(shù)據(jù)分析處理，本文考慮采用C#與MATLAB混合編程，由C#實(shí)現(xiàn)用戶界面設(shè)計(jì)，MATLAB實(shí)現(xiàn)最優(yōu)解的快速迭代求解，C#與MATLAB采用.NET組件技術(shù)實(shí)現(xiàn)混合編程，在MATLAB中調(diào)用lsqnonlin函數(shù)并編譯生成動(dòng)態(tài)鏈接庫(kù)(Dynamic Link Library，DLL)文件，在C#中添加該.DLL文件引用，即可在C#中執(zhí)行與MATLAB同樣的計(jì)算功能[14～18]。具體實(shí)現(xiàn)步驟如下：

(1)在MATLAB中新建一個(gè).m文件，運(yùn)行編譯工具(deploytool)，生成動(dòng)態(tài)鏈接庫(kù)。

(2)打開C#，添加該引用以及程序集文件“MWArray.dll”并加入命名空間，MWArray是用于C#與MATLAB之間的數(shù)據(jù)交換類。

(3)在C#中進(jìn)行數(shù)據(jù)類型轉(zhuǎn)換，調(diào)用動(dòng)態(tài)鏈接庫(kù)開始計(jì)算，并將結(jié)果返回到界面中顯示。

C#與MATLAB混合編程的方式可充分發(fā)揮兩種語(yǔ)言與平臺(tái)的優(yōu)勢(shì)，特別是調(diào)用MATLAB的數(shù)學(xué)函數(shù)能保證迭代算法的正確性，但是混合編程的方式也存在其缺點(diǎn)：一是運(yùn)行程序之前必須提前安裝整個(gè)MATLAB軟件，準(zhǔn)備工作比較繁瑣，而最終只用到其中的幾個(gè)函數(shù)，當(dāng)程序需要商業(yè)化時(shí)，其經(jīng)濟(jì)性就凸顯不足；二是每次運(yùn)算時(shí)，數(shù)據(jù)需要頻繁地在兩個(gè)平臺(tái)之間傳輸，當(dāng)測(cè)量點(diǎn)數(shù)量較多時(shí)，程序的運(yùn)行速度大大減慢，也不利于軟件界面的開發(fā)。為解決上述不足，本文進(jìn)一步研究了僅依靠C#語(yǔ)言編程來(lái)實(shí)現(xiàn)高斯牛頓法求解激光跟蹤干涉儀基站空間坐標(biāo)，并通過(guò)與MATLAB混合編程的結(jié)果比較，驗(yàn)證C#編程算法的準(zhǔn)確性。

5 C#實(shí)現(xiàn)高斯牛頓算法

5.1 數(shù)據(jù)定義

定義已知的m個(gè)目標(biāo)點(diǎn)的空間坐標(biāo)，l為激光跟蹤干涉儀的測(cè)距讀數(shù)，b為待求參量(x0j,y0j,z0j,l0j)T，設(shè)置迭代次數(shù)為20次，精度因子為10-8。

5.2 設(shè)計(jì)類進(jìn)行矩陣運(yùn)算

設(shè)計(jì)一個(gè)類為Matrix用于矩陣的運(yùn)算，包含矩陣的加法、減法、乘法、轉(zhuǎn)置以及求逆，并且在計(jì)算之前進(jìn)行是否符合矩陣運(yùn)算的規(guī)則判斷。

5.3 重復(fù)迭代計(jì)算

根據(jù)激光跟蹤干涉儀在三維空間的位置估計(jì)，給出一個(gè)粗略的迭代初值，開始迭代循環(huán)計(jì)算，并最終得到滿足精度因子的解。

第1步：按式(4)定義殘差公式數(shù)學(xué)模型；

第2步：求出雅克比矩陣J及其轉(zhuǎn)置矩陣JT，其中ri對(duì)b1、b2、b3、b4的偏導(dǎo)分別為：

(12)

第3步：求出的矩陣JTJ，以及殘差ri；

第4步：計(jì)算(JTJ)-1JTri，代入迭代公式bk+1=bk-(JTJ)-1JTr(bk)；

6 實(shí)驗(yàn)驗(yàn)證

如圖2所示，利用三坐標(biāo)測(cè)量機(jī)及激光跟蹤干涉儀搭建實(shí)驗(yàn)系統(tǒng)，激光目標(biāo)靶鏡裝在三坐標(biāo)測(cè)量機(jī)的測(cè)頭座上。在該實(shí)驗(yàn)系統(tǒng)中，三坐標(biāo)測(cè)量機(jī)的機(jī)器坐標(biāo)系為全局空間坐標(biāo)系，激光跟蹤干涉儀的基站空間坐標(biāo)是未知的，三坐標(biāo)測(cè)量機(jī)的坐標(biāo)讀數(shù)為激光目標(biāo)靶鏡的坐標(biāo)值，當(dāng)激光目標(biāo)靶鏡在全局空間坐標(biāo)系內(nèi)運(yùn)動(dòng)時(shí)，激光跟蹤干涉儀能夠跟蹤目標(biāo)靶鏡并進(jìn)行長(zhǎng)度測(cè)量。

圖2 激光跟蹤干涉儀基站空間坐標(biāo)標(biāo)定實(shí)驗(yàn)

利用上述實(shí)驗(yàn)系統(tǒng)對(duì)激光跟蹤干涉儀開展基站空間坐標(biāo)標(biāo)定實(shí)驗(yàn)。首先將激光跟蹤干涉儀置于第1個(gè)站位LT1上，令激光目標(biāo)靶鏡在全局坐標(biāo)系中運(yùn)動(dòng)至設(shè)定好的6個(gè)目標(biāo)點(diǎn)，測(cè)量每個(gè)目標(biāo)點(diǎn)相對(duì)初始目標(biāo)點(diǎn)的相對(duì)長(zhǎng)度變化值；然后將激光跟蹤干涉儀陸續(xù)移動(dòng)至第2個(gè)站位LT2和第3個(gè)站位LT3，在每個(gè)站位下令激光目標(biāo)靶鏡運(yùn)動(dòng)至不同的6個(gè)目標(biāo)點(diǎn)，分別采集激光跟蹤干涉儀的數(shù)據(jù)如表1所示。

表1 空間目標(biāo)點(diǎn)坐標(biāo)及激光跟蹤干涉儀測(cè)距讀數(shù)

將表1中數(shù)據(jù)分別導(dǎo)入到C#的兩個(gè)程序代碼中，第一個(gè)程序代碼為C#與MATLAB混合編程，第二個(gè)為C#編程，編程軟件平臺(tái)為微軟Visual Studio 2017。

在兩個(gè)程序中分別對(duì)激光跟蹤干涉儀的基站空間坐標(biāo)進(jìn)行優(yōu)化迭代計(jì)算和求解，迭代初值為每個(gè)基站站位的估計(jì)值，基站站位的精確計(jì)算結(jié)果以及兩種方式的運(yùn)行時(shí)間如表2所示。

表2 C#編程與C#調(diào)用MATLAB計(jì)算結(jié)果

從表2列出的激光跟蹤干涉儀3個(gè)基站空間坐標(biāo)計(jì)算結(jié)果可以看到：采用C#編程的高斯牛頓法在迭代5或6次后即與C#調(diào)用MATLAB的計(jì)算結(jié)果一致，差值在10-7數(shù)量級(jí)；并且從運(yùn)行時(shí)間來(lái)看，用時(shí)更少，效率更高。實(shí)驗(yàn)結(jié)果表明C#編程實(shí)現(xiàn)了高斯牛頓算法對(duì)非線性超定方程組最優(yōu)解的快速求解，驗(yàn)證了利用C#編程求解激光跟蹤干涉儀基站空間坐標(biāo)的準(zhǔn)確性以及可行性。

7 結(jié) 論

本文將激光跟蹤干涉儀空間坐標(biāo)標(biāo)定問(wèn)題轉(zhuǎn)化為非線性最小二乘問(wèn)題，采用C#編程實(shí)現(xiàn)了高斯牛頓算法。通過(guò)實(shí)驗(yàn)將C#編程計(jì)算結(jié)果與C#調(diào)用MATLAB的計(jì)算結(jié)果相比，差值在10-7數(shù)量級(jí)，并且效率更高，驗(yàn)證了C#編程算法的正確性以及可行性，說(shuō)明了高斯牛頓法進(jìn)行非線性擬合具有較好的擬合精度，收斂速度較快，為后續(xù)進(jìn)行多個(gè)站位下的聯(lián)合測(cè)量以及利用多邊法原理反求動(dòng)點(diǎn)運(yùn)動(dòng)軌跡研究打下了基礎(chǔ)。