回丙偉，趙竹新，文貢堅(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)行.這對(duì)提高測量的精度極為不利，同時(shí)膠片式攝影機(jī)復(fù)雜的操作過程對(duì)靶場測量人員要求較高[2].近年來隨著光電子和計(jì)算機(jī)技術(shù)的發(fā)展，一方面線陣像機(jī)成像速度越來越高，并以數(shù)字化和自動(dòng)化的突出優(yōu)點(diǎn)逐步替代膠片式狹縫攝影機(jī)；另一方面研究人員對(duì)線陣像機(jī)的成像模型研究越來越多，線陣像機(jī)參數(shù)標(biāo)定技術(shù)日趨成熟[6，7]，這為在嚴(yán)格線陣像機(jī)成像方程下研究目標(biāo)的運(yùn)動(dòng)與姿態(tài)參數(shù)測量方法奠定了基礎(chǔ).本文在這種背景下開展了線陣像機(jī)靶場測量方法的探索和研究.

為獲得穩(wěn)定的三維測量結(jié)果，測量中常采用交會(huì)攝影法獲取多視角圖像.根據(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)系相對(duì)于測量坐標(biāo)系OXYZ的相對(duì)方位關(guān)系在經(jīng)過標(biāo)定后已知；③目標(biāo)坐標(biāo)系ox′y′z′，目標(biāo)上特征點(diǎn)的三維坐標(biāo)在該坐標(biāo)系下已知，該坐標(biāo)系相對(duì)于測量坐標(biāo)系OXYZ的方位關(guān)系在本文中是待求參數(shù).下文中，將目標(biāo)上的一個(gè)特征點(diǎn)和其對(duì)應(yīng)的圖像點(diǎn)稱為一對(duì)控制點(diǎn).

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

在攝取到彈丸目標(biāo)的線陣立體圖像后，基于多對(duì)控制點(diǎn)和成像方程，聯(lián)立一組以彈丸速度、姿態(tài)和初始位置8個(gè)參數(shù)為未知數(shù)的非線性方程組，并進(jìn)行最優(yōu)化建模.在對(duì)模型求解的過程中，首先利用Powell算法對(duì)非線性的姿態(tài)參數(shù)進(jìn)行搜索，同時(shí)使用線性方程組對(duì)位置和速度參數(shù)進(jìn)行求解，并在搜索過程中，對(duì)優(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為對(duì)應(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ī)的動(dòng)態(tài)成像結(jié)果，以下稱之為線排列圖像.通?？梢哉J(rèn)為彈體在經(jīng)過線陣像機(jī)視平面的過程中是保持固定姿態(tài)勻速運(yùn)動(dòng)的，下面推導(dǎo)該過程中彈體表面點(diǎn)的動(dòng)態(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)動(dòng)速度矢量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)的動(dòng)態(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)的對(duì)應(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ì)上是利用若干對(duì)控制點(diǎn)，求解一個(gè)含有8個(gè)未知數(shù)的非線性方程組.由動(dòng)態(tài)成像方程式（5）或式（7）可知，一對(duì)控制點(diǎn)可以建立2個(gè)方程.因此方程組可解的條件要求至少存在4對(duì)控制點(diǎn).由于彈丸目標(biāo)在飛行狀態(tài)時(shí)高速旋轉(zhuǎn)，其表面的特征點(diǎn)在線排列圖像中通常難以正確識(shí)別.而比較穩(wěn)定的特征點(diǎn)是彈尖和彈尾等特殊點(diǎn)，如圖2所示.

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

對(duì)于彈尖上的控制點(diǎn)，無論在圖像中還是在目標(biāo)中都是容易確定的；而對(duì)于彈尾上來說，圖像上容易識(shí)別的上下兩個(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對(duì)控制點(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í)，即 為 方 程 組 的唯一解：

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

2.2 初值搜索

對(duì)于式（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的線性方程.利用這一特性和解的唯一性，本文通過對(duì)α、β進(jìn)行搜索，并對(duì)每一給定值的α（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í)對(duì)目標(biāo)函數(shù)進(jìn)行對(duì)數(shù)變換，這并不影響函數(shù)的單調(diào)性.

圖3 目標(biāo)函數(shù)值的對(duì)數(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）視為觀測值，對(duì)于每一個(gè)特征點(diǎn)通過將式（5）或式（7）進(jìn)行一階泰勒展開，列各點(diǎn)的誤差方程式為

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

式中，

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

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

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

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

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

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

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

首先利用3dMax創(chuàng)建了一個(gè)彈長550mm，彈徑70mm的彈丸目標(biāo)三維表面模型；然后利用第1節(jié)線陣像機(jī)的成像原理對(duì)彈丸飛過像機(jī)視場時(shí)進(jìn)行8次仿真成像，每次成像中彈丸運(yùn)動(dòng)與姿態(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é)的測量方法，分別對(duì)仿真得到的圖像進(jìn)行測量.對(duì)合速度參數(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é)果對(duì)比表1可知，合速度的相對(duì)測量誤差小于2%，姿態(tài)的測量最大誤差小于0.4°，平均偏差為0.20°.

4 結(jié)論與展望

在本文的測量方法中，彈丸目標(biāo)的初始時(shí)刻在測量坐標(biāo)系中的位置Xo、Yo、Zo，并不是待求參數(shù)，但是為了能夠從嚴(yán)格投影方程的角度來分析彈丸的運(yùn)動(dòng)和姿態(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í)利用原方程對(duì)線性參數(shù)進(jìn)行求解，得到了方程組的一組近似解；然后進(jìn)一步使用Gauss－Newton迭代法對(duì)待求參數(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ù)先矯正，因此本文在對(duì)像機(jī)成像模型的討論中并未涉及畸變模型.

在具體的靶場實(shí)踐中采用本文測量方法，可能為測量引入較大的誤差的因素主要有：①目標(biāo)在經(jīng)過線陣像機(jī)視平面過程中的勻速特性；②像機(jī)成像中的積分時(shí)間帶來的圖像點(diǎn)模糊；③人工拾取特征點(diǎn)的不確定性.量化研究并盡量消除上述因素對(duì)測量精度造成的影響是下一步需要研究的重點(diǎn)內(nèi)容.此外，本文測量方法只針對(duì)具有旋轉(zhuǎn)體結(jié)構(gòu)的彈丸目標(biāo)而設(shè)計(jì)，而對(duì)于一般結(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)動(dòng)參數(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)動(dòng)狀態(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.