基于線陣像機(jī)的彈丸速度與姿態(tài)測量

回丙偉，趙竹新，文貢堅(jiān)

（國防科學(xué)技術(shù)大學(xué)ATR重點(diǎn)實(shí)驗(yàn)室，長沙410073）

線陣像機(jī)以及傳統(tǒng)的膠片式狹縫攝影機(jī)作為一類獨(dú)特的光學(xué)成像設(shè)備，能夠?qū)?jīng)過其視平面的目標(biāo)進(jìn)行掃描成像，記錄目標(biāo)在一個(gè)極短時(shí)間段內(nèi)的速度和姿態(tài)信息，因而在彈丸的攻角測量中具有不可替代的作用[1～5].由于早期的狹縫攝影機(jī)參數(shù)難以標(biāo)定，從而無法從攝影幾何的角度建立嚴(yán)格的數(shù)學(xué)成像模型，導(dǎo)致傳統(tǒng)測量方法通常在一定近似條件下進(jìn)行.這對提高測量的精度極為不利，同時(shí)膠片式攝影機(jī)復(fù)雜的操作過程對靶場測量人員要求較高[2].近年來隨著光電子和計(jì)算機(jī)技術(shù)的發(fā)展，一方面線陣像機(jī)成像速度越來越高，并以數(shù)字化和自動化的突出優(yōu)點(diǎn)逐步替代膠片式狹縫攝影機(jī)；另一方面研究人員對線陣像機(jī)的成像模型研究越來越多，線陣像機(jī)參數(shù)標(biāo)定技術(shù)日趨成熟[6，7]，這為在嚴(yán)格線陣像機(jī)成像方程下研究目標(biāo)的運(yùn)動與姿態(tài)參數(shù)測量方法奠定了基礎(chǔ).本文在這種背景下開展了線陣像機(jī)靶場測量方法的探索和研究.

為獲得穩(wěn)定的三維測量結(jié)果，測量中常采用交會攝影法獲取多視角圖像.根據(jù)線陣像機(jī)視場位于一個(gè)平面內(nèi)這一特性，實(shí)際中常采用鏡面反射法獲取待測彈丸的立體圖像，如圖1所示.在彈丸飛行軌跡的下方，放置一個(gè)大約與水平面成45°的平面鏡，當(dāng)彈丸從平面鏡上方飛過時(shí)，彈丸本身和彈丸在平面鏡中的像都能在該線陣像機(jī)中成像.上述方法實(shí)現(xiàn)了使用一個(gè)線陣像機(jī)從2個(gè)視角獲取目標(biāo)立體圖像的目的[8].這種測量方式涉及3個(gè)坐標(biāo)系：①測量坐標(biāo)系OXYZ，各測量參數(shù)均在該坐標(biāo)系下描述，是測量的基準(zhǔn)坐標(biāo)系；②像機(jī)坐標(biāo)系sxyz，目標(biāo)特征點(diǎn)的像點(diǎn)坐標(biāo)在該坐標(biāo)系下已知，且該坐標(biāo)系相對于測量坐標(biāo)系OXYZ的相對方位關(guān)系在經(jīng)過標(biāo)定后已知；③目標(biāo)坐標(biāo)系ox′y′z′，目標(biāo)上特征點(diǎn)的三維坐標(biāo)在該坐標(biāo)系下已知，該坐標(biāo)系相對于測量坐標(biāo)系OXYZ的方位關(guān)系在本文中是待求參數(shù).下文中，將目標(biāo)上的一個(gè)特征點(diǎn)和其對應(yīng)的圖像點(diǎn)稱為一對控制點(diǎn).

圖1 利用一個(gè)線陣像機(jī)獲取立體圖像的方法

在攝取到彈丸目標(biāo)的線陣立體圖像后，基于多對控制點(diǎn)和成像方程，聯(lián)立一組以彈丸速度、姿態(tài)和初始位置8個(gè)參數(shù)為未知數(shù)的非線性方程組，并進(jìn)行最優(yōu)化建模.在對模型求解的過程中，首先利用Powell算法對非線性的姿態(tài)參數(shù)進(jìn)行搜索，同時(shí)使用線性方程組對位置和速度參數(shù)進(jìn)行求解，并在搜索過程中，對優(yōu)化模型的目標(biāo)函數(shù)值不斷尋優(yōu)得到測量初值；然后進(jìn)一步使用Gauss－Newton法進(jìn)行精確迭代求解.

1 線陣像機(jī)測量原理與成像方程

靜態(tài)場景中，線陣像機(jī)成像可以看作面陣像機(jī)成像的一種特殊情況，只有位于線陣傳感器與主光軸所構(gòu)成的平面（即視平面[6]）內(nèi)的場景，能夠按中心投影原理成像.在sxyz中描述，即x≡0，y方向滿足中心投影方程[9]：

式中，（X，Y，Z）為場景中一點(diǎn)在OXYZ中的空間坐標(biāo)，y為對應(yīng)的像點(diǎn)坐標(biāo).Xs，Ys，Zs為像機(jī)攝影中心在OXYZ中的坐標(biāo)；旋轉(zhuǎn)矩陣元素al，bl，cl（l＝1，2，3）是由像機(jī)外參數(shù)姿態(tài)角φ，ω，κ決定的.y0為像主點(diǎn)位置，fy為像機(jī)主距，y0、fy統(tǒng)稱為像機(jī)內(nèi)參數(shù).

式（1）的第1式即為視平面的方程，由該式可得：

將式（2）代入式（1）中第2式，并根據(jù)旋轉(zhuǎn)矩陣的性質(zhì)進(jìn)行化簡，得到線陣像機(jī)的靜態(tài)成像方程：

將每一條成像線按時(shí)間順次排列成一幅二維圖像，即為線陣像機(jī)的動態(tài)成像結(jié)果，以下稱之為線排列圖像.通?？梢哉J(rèn)為彈體在經(jīng)過線陣像機(jī)視平面的過程中是保持固定姿態(tài)勻速運(yùn)動的，下面推導(dǎo)該過程中彈體表面點(diǎn)的動態(tài)成像方程.

一般三維剛體目標(biāo)的空間姿態(tài)需要用3個(gè)空間旋角來表達(dá).由于絕大多數(shù)的彈丸目標(biāo)都具有旋轉(zhuǎn)體的幾何結(jié)構(gòu)，因此沿彈丸中軸線的自轉(zhuǎn)角不用考慮，此時(shí)其空間姿態(tài)通常使用目標(biāo)繞2個(gè)指定坐標(biāo)軸（本文中即為目標(biāo)坐標(biāo)系的y′軸和z′軸）的2個(gè)空間旋轉(zhuǎn)角α和β來描述.并設(shè)某一時(shí)刻在OXYZ中，目標(biāo)坐標(biāo)系原點(diǎn)的坐標(biāo)為（Xo，Yo，Zo），運(yùn)動速度矢量v＝（vX vYvZ）；線陣像機(jī)的成像頻率為ν；假設(shè)目標(biāo)坐標(biāo)系中一點(diǎn)（x′，y′，z′）在線陣像機(jī)的第n次成像時(shí)被捕獲；像機(jī)相鄰兩次成像間隔t＝n/ν；并記（DX DY DZ）＝（vX/νvY/νvZ/ν）.可得目標(biāo)點(diǎn)（x′，y′，z′）在經(jīng)過線陣像機(jī)瞬間的空間坐標(biāo)（X，Y，Z）為[8]

式（4）表明，能在線陣像機(jī)中成像的彈丸表面點(diǎn)的空間坐標(biāo)由三部分確定：彈丸的姿態(tài)角α，β決定的（XR，YR，ZR）；初始位置時(shí)刻彈丸的目標(biāo)坐標(biāo)系原點(diǎn)在OXYZ中的坐標(biāo)（Xo，Yo，Zo）；彈丸的速度（vX，vY，vZ）決定的（nDX，nDY，nDZ）.因此該式包含了全部待求參數(shù).將式（4）代入式（3），并變換其中的約束條件方程得彈丸目標(biāo)的動態(tài)成像方程：

結(jié)合平面鏡成像原理，目標(biāo)點(diǎn)（x′，y′，z′）的鏡像點(diǎn)在OXYZ坐標(biāo)系中的坐標(biāo)為

同理將式（6）代入式（3）可得：

因此在上文的假設(shè)前提下，式（5）與式（7）分別描述了彈丸目標(biāo)上任意一點(diǎn)及其平面鏡中的虛像點(diǎn)與各自其圖像點(diǎn)的對應(yīng)關(guān)系.基于這樣一種嚴(yán)格幾何成像關(guān)系，本文以α，β，DX，DY，DZ，Xo，Yo，Zo為未知量實(shí)現(xiàn)彈丸各項(xiàng)參數(shù)的測量.

2 速度與姿態(tài)測量

2.1 特征點(diǎn)選取與最優(yōu)化模型

在本文中，彈丸的速度與姿態(tài)參數(shù)的測量本質(zhì)上是利用若干對控制點(diǎn)，求解一個(gè)含有8個(gè)未知數(shù)的非線性方程組.由動態(tài)成像方程式（5）或式（7）可知，一對控制點(diǎn)可以建立2個(gè)方程.因此方程組可解的條件要求至少存在4對控制點(diǎn).由于彈丸目標(biāo)在飛行狀態(tài)時(shí)高速旋轉(zhuǎn)，其表面的特征點(diǎn)在線排列圖像中通常難以正確識別.而比較穩(wěn)定的特征點(diǎn)是彈尖和彈尾等特殊點(diǎn)，如圖2所示.

圖2 線排列圖像中彈丸上特殊點(diǎn)的選取

對于彈尖上的控制點(diǎn)，無論在圖像中還是在目標(biāo)中都是容易確定的；而對于彈尾上來說，圖像上容易識別的上下兩個(gè)端點(diǎn)（圖2中的彈尾A和彈尾B）卻在目標(biāo)坐標(biāo)系中難以確定，因此本文選擇彈尾的中心點(diǎn)作為控制點(diǎn)，在圖像中表現(xiàn)為A、B兩點(diǎn)連線的中點(diǎn).這樣從2個(gè)視角拍攝的彈丸圖像中可以確定4對控制點(diǎn)，可以滿足方程求解的基本條件.足夠數(shù)量的控制點(diǎn)是實(shí)現(xiàn)本文測量方法的必要條件.

多元非線性方程組的求解，可根據(jù)最小二乘原理轉(zhuǎn)化為非線性最優(yōu)化模型來描述.假設(shè)目標(biāo)表面有N個(gè)特征點(diǎn)，（i＝1，2，…，N，N≥4），在線陣像機(jī)中所成的像點(diǎn)為（ni，yi）.利用（x′i，y′i，z′i），根據(jù)一組參數(shù)Φ＝（αβDXDYDZXoYoZo），按式（5）或式（7）計(jì)算的投影點(diǎn) 坐 標(biāo)）.若 參 數(shù) 向 量Φ能 使（ni，yi）與滿 足 最 優(yōu) 化 目 標(biāo) 函 數(shù) 時(shí)，即 為 方 程 組 的唯一解：

對于式（8）的最優(yōu)化數(shù)學(xué)模型，Gauss－Newton法是解決這類問題的經(jīng)典方法，但對于高維的優(yōu)化問題，通常需要為其提供一個(gè)較好的初值，本節(jié)根據(jù)模型的特點(diǎn)，利用Powell算法搜索待測參數(shù)的近似初始值，然后進(jìn)一步使用Gauss－Newton迭代法精確求解.

2.2 初值搜索

對于式（5）或式（7）組成的方程，若已知參數(shù)α，β，可將其代換為

式中，lgh，mgj，ng（g＝1，2；h＝1，2，…，6；j＝1，2，3）均可根據(jù)目標(biāo)上特征點(diǎn)在ox′y′z′中的坐標(biāo)和像機(jī)參數(shù)等已知量進(jìn)行代數(shù)換算，（XR，YR，ZR）在給定α，β值的情況也為已知量，因此式（9）是關(guān)于6個(gè)待求參數(shù)DX、DY、DZ、Xo、Yo、Zo的線性方程.利用這一特性和解的唯一性，本文通過對α、β進(jìn)行搜索，并對每一給定值的α（k）、β（k），利用原方程組解算線性參數(shù)，并 尋 找 一 組 參 數(shù)在目標(biāo)函數(shù)式（8）的度量下達(dá)到最小值，此時(shí)可認(rèn)為其近似為方程組的解.因此，在式（9）的約束下，可認(rèn)為式（8）的目標(biāo)函數(shù)是關(guān)于α、β求極小值的問題.圖3顯示了某次測量中，目標(biāo)函數(shù)值與α、β的關(guān)系.為直觀地表現(xiàn)目標(biāo)函數(shù)的極小值，繪圖時(shí)對目標(biāo)函數(shù)進(jìn)行對數(shù)變換，這并不影響函數(shù)的單調(diào)性.

圖3 目標(biāo)函數(shù)值的對數(shù)與α、β的關(guān)系

下面考慮使用快速的搜索算法使目標(biāo)函數(shù)能在任意α、β初值情況下，快速收斂到極小值，該值即為待測參數(shù)的一組近似解.Powell算法是一種不涉及目標(biāo)函數(shù)導(dǎo)數(shù)的最優(yōu)化搜索算法，它是以正定二次函數(shù)為背景，以共軛方向?yàn)榛A(chǔ)的最優(yōu)化問題直接搜索方法.該方法具有良好的模塊化特性，已經(jīng)在許多工程領(lǐng)域獲得了廣泛的應(yīng)用，本文以Powell算法為基礎(chǔ)搜索方程解的近似值，具體過程如下：

①首先任意設(shè)定待求參數(shù)初值α（0）、β（0）；

②進(jìn)入Powell搜索，然后按照算法的搜索規(guī)則尋找目標(biāo)函數(shù)極小值，并在第k次搜索中，先將α（k）、β（k）代入式（9），利用線性最小二乘法計(jì)算

③將α（k），β（k）代入式（8）計(jì)算當(dāng)前目標(biāo)函數(shù)值ε（k）；

④重復(fù)上述步驟②和步驟③直到最終輸出α，β，DX，DY，DZ，Xo，Yo，Zo滿足精度的測量值.

雖然本節(jié)Powell搜索中使用的目標(biāo)函數(shù)是按最小二乘的形式定義的，但并不是所有的參數(shù)都統(tǒng)一參加了最小二乘平差計(jì)算，因此，其搜索結(jié)果并不是最小二乘意義下的最優(yōu)解，下一節(jié)將使用嚴(yán)格最小二乘迭代求解其精確解.

2.3 迭代求解

根據(jù)最小二乘原理，將N個(gè)特征點(diǎn)在ox′y′z′中的坐標(biāo)（x′i，y′i，z′i）視為真值，而 把 相 應(yīng) 的 圖像點(diǎn)坐標(biāo)（ni，yi）視為觀測值，對于每一個(gè)特征點(diǎn)通過將式（5）或式（7）進(jìn)行一階泰勒展開，列各點(diǎn)的誤差方程式為

式（10）用矩陣形式表示為

式中，

若有N對控制點(diǎn)，則可按式（11）列出N組誤差方程式，共同構(gòu)成總的誤差方程式：

式中，E＝（E1E2…EN）T，A＝（A1A2…AN）T，L＝（l1l2…lN）T.

根據(jù)最小二乘平差原理，可以得到待求參數(shù)近似值的改正量：

由于式（10）中的各系數(shù)取自泰勒級數(shù)的線性展開，因此迭代計(jì)算通過逐漸趨近的方法，即用原近似值與式（13）求得的改正量的和作為新的近似值，重復(fù)上述過程，直至改正量各分量的絕對值小于某一個(gè)設(shè)定的極限值為止，最后得到待求參數(shù)的解.

在得到DX，DY，DZ的精確迭代結(jié)果后，彈丸速度（vXvYvZ）＝（νDXνDYνDZ）.

3 仿真實(shí)驗(yàn)

為驗(yàn)證上述方法的測量精度，本文對不同像機(jī)參數(shù)下，不同彈丸幾何尺寸和不同的運(yùn)動姿態(tài)參數(shù)進(jìn)行了大量圖像仿真和測量實(shí)驗(yàn)，下文以一組實(shí)驗(yàn)為例進(jìn)行說明.

首先利用3dMax創(chuàng)建了一個(gè)彈長550mm，彈徑70mm的彈丸目標(biāo)三維表面模型；然后利用第1節(jié)線陣像機(jī)的成像原理對彈丸飛過像機(jī)視場時(shí)進(jìn)行8次仿真成像，每次成像中彈丸運(yùn)動與姿態(tài)參數(shù)如表1所示，像機(jī)主機(jī)像素為fy＝5 000，y0＝0，線陣像機(jī)成像速度ν＝140 000s－1.像機(jī)外參數(shù)如表2所示.

表1 仿真中彈丸目標(biāo)速度與姿態(tài)參數(shù)設(shè)定值

表2 線陣像機(jī)外參數(shù)

使用上面第1節(jié)與第2節(jié)的測量方法，分別對仿真得到的圖像進(jìn)行測量.對合速度參數(shù)v和2個(gè)姿態(tài)參數(shù)分別進(jìn)行精度比較，如表3所示，表中v為利用表1中的3個(gè)速度分量計(jì)算所得的合速度的理論值；v′，α′，β′為彈丸合速度與姿態(tài)參數(shù)的測量值.

表3 彈丸速度與姿態(tài)參數(shù)測量值

由表3的測量結(jié)果對比表1可知，合速度的相對測量誤差小于2%，姿態(tài)的測量最大誤差小于0.4°，平均偏差為0.20°.

4 結(jié)論與展望

在本文的測量方法中，彈丸目標(biāo)的初始時(shí)刻在測量坐標(biāo)系中的位置Xo、Yo、Zo，并不是待求參數(shù)，但是為了能夠從嚴(yán)格投影方程的角度來分析彈丸的運(yùn)動和姿態(tài)，該參數(shù)向量的引入是有必要的.

本文利用成像方程和從目標(biāo)表面選取的若干個(gè)特征點(diǎn)建立了以目標(biāo)的初始位置、速度和姿態(tài)參數(shù)共8個(gè)未知量的非線性方程組.通過將方程組的求解問題進(jìn)行最優(yōu)化建模和求解，實(shí)現(xiàn)了各參數(shù)的測量.在求解過程中，首先利用了方程組中大量參數(shù)是線性的這一特性，以Powell搜索為基礎(chǔ)，同時(shí)利用原方程對線性參數(shù)進(jìn)行求解，得到了方程組的一組近似解；然后進(jìn)一步使用Gauss－Newton迭代法對待求參數(shù)進(jìn)行精確求解.本文方法的優(yōu)勢體現(xiàn)在：①將線陣像機(jī)的嚴(yán)格成像方程引入彈丸姿態(tài)參數(shù)的測量中，使測量方法有嚴(yán)格的理論依據(jù)；②在嚴(yán)格成像模型下，給出了基于特征點(diǎn)坐標(biāo)的彈丸目標(biāo)初始位置、速度和姿態(tài)參數(shù)的解算方法.

考慮到像機(jī)鏡頭的光學(xué)畸變，在圖像預(yù)處理的過程中根據(jù)像機(jī)標(biāo)定所得的徑向與切向畸變系數(shù)進(jìn)行預(yù)先矯正，因此本文在對像機(jī)成像模型的討論中并未涉及畸變模型.

在具體的靶場實(shí)踐中采用本文測量方法，可能為測量引入較大的誤差的因素主要有：①目標(biāo)在經(jīng)過線陣像機(jī)視平面過程中的勻速特性；②像機(jī)成像中的積分時(shí)間帶來的圖像點(diǎn)模糊；③人工拾取特征點(diǎn)的不確定性.量化研究并盡量消除上述因素對測量精度造成的影響是下一步需要研究的重點(diǎn)內(nèi)容.此外，本文測量方法只針對具有旋轉(zhuǎn)體結(jié)構(gòu)的彈丸目標(biāo)而設(shè)計(jì)，而對于一般結(jié)構(gòu)的三維目標(biāo)，由于其空間姿態(tài)需要3個(gè)空間旋轉(zhuǎn)進(jìn)行描述，在使用本文測量基本原理的情況下，需要為之設(shè)計(jì)新的解算方法.

[1]高昕，王建軍，湯陽.利用狹縫相機(jī)實(shí)現(xiàn)彈丸攻角高精度測量研究[J].飛行器測控學(xué)報(bào)，2004，23（1）：69－72.GAO Xin，WANG Jian－jun，TANG Yang.Precision measurement of nutation angle using slit cameras[J].Journal of Spacecraft TT ＆ C Technology，2004，23（1）：69－72.（in Chinese）

[2]張三喜，薛以輝，盧宇.狹縫攝影膠片圖像運(yùn)動參數(shù)測量和處理[J].光子學(xué)報(bào)，1999，28（12）：1 117－1 121.ZHANG San－xi，XUE Yi－h(huán)ui，LU Yu.The measuring and processing of streak camera’s image motion parameters[M].ACTA Photonica Sinica，1999，28（12）：1 117－1 121.（in Chinese）

[3]LI J K，CHEN L Y.A slit photography system based on linear CCD[J].Laser ＆Infrared，2009，39（3）：300－303.

[4]SONG W D.Software design for measurement of bullet attitude based on linear CCD[C].9th International Conference on Electronic Measurement and Instruments.Xi’an：CIE，2007：521－524.

[5]HUGHETT P.Projectile velocity and spin rate by image processing of synchro－ballistic photography[C].Ultrahigh－and High－Speed Photography，Videography，Photonics and Velocimetry.San Diego：SPIE，1990：237－248.

[6]RADU H，ROGER M，BOGUSLAW L.On Single－scanline camera calibration[J].IEEE Transactions on Robobtics and Automatic，1993，9（1）：71－75.

[7]LUNA C A，MAZO M，LAZARO J L.Calibration of line－scan camera[J].IEEE Transactions on Instruementation and Measurement，2010，59（8）：2 185－2 190.

[8]馮文灝.近景攝影測量——物體外形與運(yùn)動狀態(tài)的攝影法測定[M].武漢：武漢大學(xué)出版社，2007：72－73.FENG Wen－h(huán)ao.Close－range photogrammetry（measurement of object’s facade and motion）[M].Wuhan：Wuhan University Press，2007：72－73.（in Chinese）

[9]WANG Zhi－zhuo.Principles of photogrammetry[M].Press of Wuhan Technical University of Surveying and Mapping ＆Publishing House of Surveying and Mapping of Beijing，1990：18－19.