基于Monte Carlo法的光電跟蹤測(cè)量系統(tǒng)精度分析和布站優(yōu)化方法＊

張振鐸，張 彬

（中國(guó)科學(xué)院 長(zhǎng)春光學(xué)精密儀器與物理研究所，吉林 長(zhǎng)春 130033）

引 言

靶場(chǎng)光測(cè)設(shè)備是對(duì)彈道導(dǎo)彈、飛行器、衛(wèi)星等武器設(shè)備進(jìn)行精確觀測(cè)的有效工具，它通過(guò)角度測(cè)量和交誤差分析實(shí)施對(duì)空間目標(biāo)精確定位以完成外彈道的觀測(cè)。測(cè)量中的各種誤差源直接決定定位精度［1，2］。但誤差并不是越小越好還要受到研制成本等多種因素的制約，同時(shí)各種誤差因素對(duì)光測(cè)系統(tǒng)而言，其重要性和影響方式有所不同，對(duì)各種誤差源進(jìn)行精確分析對(duì)靶場(chǎng)光測(cè)設(shè)備有重要的實(shí)際意義。一方面在設(shè)計(jì)的初期研制過(guò)程中可以對(duì)設(shè)備測(cè)量不確定度有影響的各種因素進(jìn)行嚴(yán)格的過(guò)程控制，同時(shí)對(duì)優(yōu)化材料選用、設(shè)備制造等多種環(huán)節(jié)進(jìn)行控制，以更低的研制成本更好地滿足用戶的使用要求［3－5］；另一方面通過(guò)對(duì)光測(cè)設(shè)備進(jìn)行合理的布站以提高光測(cè)設(shè)備測(cè)量精度［6，7］。光測(cè)設(shè)備的精度模型和光測(cè)設(shè)備的布站是兩個(gè)相互聯(lián)系的問(wèn)題，對(duì)光測(cè)設(shè)備的布站優(yōu)化必須以精確的誤差模型為基礎(chǔ)，同時(shí)對(duì)光測(cè)設(shè)備的誤差分析也需要給出在一定布站方式下對(duì)特定目標(biāo)軌跡的精度分析。

文中以光電經(jīng)緯儀為例，對(duì)光電經(jīng)緯儀建立了精確的Verilog－A模型，并在此基礎(chǔ)上使用 Monte Carlo法對(duì)光電經(jīng)緯儀的各種誤差源的影響進(jìn)行了詳細(xì)的統(tǒng)計(jì)學(xué)分析和研究，并針對(duì)彈道導(dǎo)彈不同的布站方式進(jìn)行數(shù)值分布，給出能適應(yīng)不同靶場(chǎng)地形的通用的布站優(yōu)化方法。

1 光電經(jīng)緯儀的誤差模型

文中提出的誤差分析方法具有一定的通用性，使用通用的模擬系統(tǒng)建模語(yǔ)言Verilog－A以對(duì)光電經(jīng)緯儀建立數(shù)學(xué)模型，并在數(shù)學(xué)模型中對(duì)測(cè)量精度有影響的各種誤差源進(jìn)行建模。Verilog－A提供了層次化的模擬系統(tǒng)模型構(gòu)架，可以一定的數(shù)學(xué)表達(dá)式在抽象的層次上對(duì)系統(tǒng)建模，以方便系統(tǒng)的分析設(shè)計(jì)和驗(yàn)證。Verilog－A支持直流、交流、瞬態(tài)、Monte Carlo法等多種分析方法。提供給用戶一定的設(shè)計(jì)參數(shù)，用戶使用Verilog－A模型可以對(duì)設(shè)備的性能進(jìn)行仿真分析進(jìn)一步指導(dǎo)光測(cè)設(shè)備的應(yīng)用。對(duì)各種光測(cè)設(shè)備建立通用的Verilog－A模型具有非常重要的實(shí)際意義。光電經(jīng)緯儀的Verilog－A模型分為基本模型和誤差模型，如圖1所示。

圖1 評(píng)估光電經(jīng)緯儀測(cè)量誤差的VerilogA模型框圖Fig.1 The verilog A model for theodolite measure error evaluation

其中基本模型實(shí)現(xiàn)光電經(jīng)緯儀的基本功能，即實(shí)現(xiàn)空間目標(biāo)3維坐標(biāo)（x，y，z）到俯仰角E和方位角A的轉(zhuǎn)換：

誤差模型使用坐標(biāo)變換的方法，光電跟蹤經(jīng)過(guò)了基座坐標(biāo)系、水平軸坐標(biāo)系、照準(zhǔn)部坐標(biāo)系、相面坐標(biāo)系，分別引入豎軸誤差、橫軸誤差、照準(zhǔn)軸誤差和傳感器誤差，同時(shí)編碼器也會(huì)引入編碼器誤差，在坐標(biāo)系X1（x，y，z）經(jīng)過(guò)旋轉(zhuǎn)變換轉(zhuǎn)到坐標(biāo)系X2（x′，y′，z′），旋轉(zhuǎn)θ角的旋轉(zhuǎn)變換矩陣為M（θ），則同一點(diǎn)在X1中坐標(biāo)與X2中坐標(biāo)的關(guān)系為：

由于旋轉(zhuǎn)變換矩陣的逆矩陣為：

得到：

光電跟蹤依次經(jīng)過(guò)了基座坐標(biāo)系、水平軸坐標(biāo)系、照準(zhǔn)部坐標(biāo)系、相面坐標(biāo)系等旋轉(zhuǎn)變換。各誤差源所引入的旋轉(zhuǎn)變換矩陣定義為Mi（i＝1，2…n），則在相面坐標(biāo)系中單位光軸方向矢量（1，0，0）T在基座坐標(biāo)系中的坐標(biāo)可以表達(dá)為：考慮到方位角和俯仰角旋轉(zhuǎn)方向的定義，各坐標(biāo)系間旋轉(zhuǎn)矩陣如表1所示

其中rot x，rot z，rot y分別為繞x，z，y軸的旋轉(zhuǎn)變換矩陣，在旋轉(zhuǎn)角很小時(shí)可以表示為：

表1 光電經(jīng)緯儀誤差模型中旋轉(zhuǎn)變換Tab.1 The rotation transformation used for theodolite intercross measure model

坐標(biāo)變換法是一種對(duì)光測(cè)設(shè)備誤差源進(jìn)行表達(dá)的通用的方法，使用此方法可以方便地對(duì)光電經(jīng)緯儀建立包含各種誤差源的準(zhǔn)確的數(shù)學(xué)模型，并進(jìn)一步建立光電經(jīng)緯儀雙站交匯測(cè)量模型如圖2所示，輸出dx，dz，dy分別為x，z，y三個(gè)方向上的交匯測(cè)量誤差。其中交匯測(cè)量模塊在異側(cè)布站時(shí)使用K公式，在同側(cè)布站時(shí)使用L公式，在Verilog－A模型中需要對(duì)除法運(yùn)算分母做限幅處理以避免在平衡點(diǎn)處出現(xiàn)收斂問(wèn)題。

圖2 光電經(jīng)緯儀交匯測(cè)量模塊Fig.2 Theodolite intercross measure model

2 基于Verilog－A模型的誤差分析法

針對(duì)上述光電經(jīng)緯儀的數(shù)學(xué)模型進(jìn)行統(tǒng)計(jì)分析，基本模型參數(shù)如表2所示。使用兩種基本的誤差分析方法，即最壞情況法分析和Monte Carlo法分析。計(jì)算了各種單項(xiàng)誤差對(duì)系統(tǒng)性能的影響?，F(xiàn)使用的分析方法具有一定的通用性，可用于各種光測(cè)儀器的誤差計(jì)算中。

2.1 最壞情況法分析法

在最壞情況下假定所有的誤差源都取最大誤差，由于在實(shí)際情況下所有的誤差源不是同時(shí)取得最大誤差數(shù)值，因此最壞情況分析將得到較為悲觀的預(yù)測(cè)數(shù)值，但這種分析方法可以快速估計(jì)各單項(xiàng)誤差對(duì)系統(tǒng)誤差影響程度和影響方式，同時(shí)給出理論上最大誤差。使用此模型仿真了各單項(xiàng)誤差都取最大數(shù)值的最壞情況，誤差主要影響方位角A，而對(duì)俯仰角E影響較小。其中編碼器誤差直接加入到方位角和俯仰角誤差中，而傳感器誤差、照準(zhǔn)差、橫軸差對(duì)系統(tǒng)誤差的影響相對(duì)較為復(fù)雜，可以通過(guò)掃描相關(guān)誤差參數(shù)計(jì)算出各誤差參數(shù)對(duì)系統(tǒng)性能的影響。

表2 光電經(jīng)緯儀模型參數(shù)選取Tab.2 Parameter selection for intercross measure model

使用該模型仿真了因俯仰角變化而引起的傳感器誤差、軸系照準(zhǔn)差和水平軸誤差A(yù)分量的變化曲線，如圖3所示。從圖3中可以看出，隨著俯仰角的增加，三種誤差源快速增加，文獻(xiàn)表明照準(zhǔn)差和水平軸誤差的δE分別與sec（E）和tan（E）成正比，而水平軸誤差與tan（E）與仿真結(jié)果相一致。同時(shí)可以看出在俯仰角較大時(shí)照準(zhǔn)差對(duì)系統(tǒng)影響約為橫軸誤差的2倍。其中圖3（a）表示了不同誤差像素?cái)?shù)目對(duì)系統(tǒng)誤差的影響，傳感器誤差將隨著誤差像素?cái)?shù)目的提高而增加，誤差像素主要來(lái)自量化誤差、細(xì)分誤差和拖尾誤差。

圖3 俯仰角變化對(duì)傳感器誤差、軸系照準(zhǔn)差和橫軸誤差的影響Fig.3 The influence of collimate error，vertical axes error，horizontal axes error

圖4表示了豎軸誤差隨方位角和俯仰角變化曲線，其中左圖表示豎軸方位誤差隨方位角和俯仰角的變化，可以看出豎軸差與方位角呈現(xiàn)正弦變化規(guī)律，且在豎軸傾斜方向上取得最大值，而隨俯仰角E以tan（E）的規(guī)律快速增加；右圖豎軸俯仰誤差與俯仰角E無(wú)關(guān)且隨方位角呈現(xiàn)正弦規(guī)律變化。此仿真結(jié)果與文獻(xiàn)［1］論述相一致。圖5表示總誤差源隨著方位角和俯仰角變化曲線，其中總方位誤差隨方位角A以正弦形式變化，其原因主要來(lái)自豎軸誤差的影響，總方位誤差隨俯仰角E快速增加，主要受到橫軸差、照準(zhǔn)差以及傳感器誤差的影響。

2.2 Monte Carlo法分析法

在實(shí)際情況下各個(gè)單項(xiàng)是以一定概率分布的形式出現(xiàn)的，誤差合成不是簡(jiǎn)單的線性疊加關(guān)系，傳統(tǒng)的方法是基于統(tǒng)計(jì)學(xué)的t分布和χ2分布不確定度分析法，這種方法在處理測(cè)量誤差傳遞時(shí)是基于線性化近似模型，同時(shí)假設(shè)各種誤差源間是相互獨(dú)立的。由于光測(cè)設(shè)備本身是一個(gè)復(fù)雜的非線性系統(tǒng)，同時(shí)誤差源間也不是完全獨(dú)立的，因此決定了這種方法具有一定的局限性。Monte Carlo法是一種通用的誤差分析工具，它將設(shè)備的各項(xiàng)誤差源表達(dá)為一定概率分布函數(shù)的形式，以相同概率分布產(chǎn)生隨機(jī)數(shù)進(jìn)行仿真計(jì)算，從而得到各種誤差對(duì)系統(tǒng)的影響。系統(tǒng)的輸出呈現(xiàn)出一定統(tǒng)計(jì)分布，通過(guò)使用MatLab統(tǒng)計(jì)學(xué)工具箱擬合的方法，可以得到系統(tǒng)輸出的準(zhǔn)確的統(tǒng)計(jì)分布函數(shù)，并以一定的數(shù)字特征如期望和方差的形式來(lái)表達(dá)。Monte Carlo算法的準(zhǔn)確度主要取決于采樣點(diǎn)的數(shù)目，可以通過(guò)合理的選用采樣點(diǎn)的數(shù)目以達(dá)到所需的計(jì)算精度。

圖4 豎軸誤差隨方位角和俯仰角的變化曲線Fig.4 The relationship of vertical axes error with azimuth and pitch angle

圖5 最壞情況下經(jīng)緯儀總誤差隨方位角和俯仰角變化曲線Fig.5 The total error with azimuth and pitch angle in the worst case

圖6 Monte Carlo誤差分析的一般步驟Fig.6 The general step for Monte Carlo error analysis

圖7 Monte Carlo分析結(jié)果的正態(tài)分布擬合Fig.7 The normal distribution fitting result

圖6表示了對(duì)系統(tǒng)進(jìn)行Monte Carlo法分析的一般步驟。在Monte Carlo法分析中假設(shè)各誤差源符合正態(tài)分布，精度范圍為±δ，其中δ為正態(tài)分布標(biāo)準(zhǔn)差，各誤差精度范圍如表1所示。圖7為方位角135°，俯仰角為45°時(shí)使用Monte Carlo法分析得到統(tǒng)計(jì)直方圖以及使用MatLab統(tǒng)計(jì)工具箱擬合的結(jié)果，擬合得到正態(tài)分布δ＝8.53。使用此方法進(jìn)一步分析了各種誤差源對(duì)系統(tǒng)性能的影響如表3所示，是在考慮各種單項(xiàng)誤差的情況下計(jì)算得到的系統(tǒng)的方位角誤差和俯仰角誤差。從分析結(jié)果看，照準(zhǔn)差和傳感器誤差對(duì)總誤差的影響最大。

3 基于Verilog－A模型的光電經(jīng)緯儀布站優(yōu)化

表3 A＝135°，E＝45°情況下 Monte Carlo法分析結(jié)果Tab.3 The Monte Carlo analytical result when A＝135°and E＝45°

基于Monte Carlo法的精度分析可以進(jìn)一步應(yīng)用于光電經(jīng)緯儀的布站設(shè)計(jì)中，可以針對(duì)特定的地形特點(diǎn)和彈道軌跡使用Monte Carlo法仿真光電經(jīng)緯儀的交匯測(cè)量過(guò)程，彈道軌跡公式如下：

現(xiàn)選擇彈道發(fā)射后2s到7s進(jìn)行理論彈道軌跡優(yōu)化，彈道發(fā)射方向?yàn)閄方向，彈道方程參數(shù)為v0＝4.5km／s，θ＝0.13，采用異側(cè)布站方式，兩站點(diǎn)坐標(biāo)為（10km，20km，0）和（10km，－20km，0），使用K公式進(jìn)行交匯測(cè)量，在只考慮經(jīng)緯儀的測(cè)量誤差的情況下，測(cè)量得到的x，z，y三個(gè)方向的誤差曲線隨時(shí)間變化曲線如圖8所示。對(duì)其在仿真時(shí)間內(nèi)的平均誤差繪制統(tǒng)計(jì)直方圖，并使用MatLab統(tǒng)計(jì)工具箱進(jìn)行正態(tài)分布擬合得到如圖9的計(jì)算結(jié)果，x，z，y三個(gè)方向的正態(tài)分布方差分別為1.168m、2.756m和0.515m。計(jì)算表明在垂直于飛行方向的z方向上誤差最大，而在高度y方向下誤差最小。

圖8 仿真x，z，y方向測(cè)量誤差隨時(shí)間變化曲線Fig.8 The time variable measure error in the direction of x，z，y

圖9 x，z，y方向平均誤差直方圖Fig.9 The average measure error histogram in the direction of x，z，y

在異側(cè)布站情形下，設(shè)置站點(diǎn)z坐標(biāo)分別為＋20km和－20km，使用Monte Carlo法仿真得到站點(diǎn)x坐標(biāo)變化對(duì)交匯測(cè)量誤差的影響，如圖10所示，圖中4條曲線分別對(duì)應(yīng)于x，z方向站點(diǎn)定位誤差d分別為0m，1m，2m和3m時(shí)系統(tǒng)三個(gè)方向上的測(cè)量誤差，從計(jì)算結(jié)果可知站點(diǎn)定位誤差在垂直于飛行方向（z方向）上影響較大，且在站點(diǎn)定位誤差小于1m時(shí)，經(jīng)緯儀本身誤差起主要作用；而當(dāng)站點(diǎn)定位誤并差大于1m時(shí)，站點(diǎn)定位誤差將成為主要的誤差源。同時(shí)注意到當(dāng)站點(diǎn)x坐標(biāo)為20km時(shí)三個(gè)方向交匯測(cè)量誤差均呈現(xiàn)最小值，為最優(yōu)布站位置。仿真結(jié)果表明在飛行段x方向中部位置布站可以獲得最小測(cè)量誤差。

圖10 異側(cè)布站時(shí)不同站點(diǎn)位置精度下站點(diǎn)x方向距離與測(cè)量誤差關(guān)系Fig.10 The relationship between measure error xdistance for different side station position optimization

設(shè)置站點(diǎn)x坐標(biāo)位＋20km，－20km時(shí)，使用Monte Carlo法仿真了站點(diǎn)z坐標(biāo)變化對(duì)交匯測(cè)量誤差的影響，如圖11所示，其中z坐標(biāo)變化范圍設(shè)置為1～20km，從仿真結(jié)果可以看出：

（1）x，z方向的誤差隨z方向站點(diǎn)距離減小而減小。

（2）而y方向誤差隨z方向站點(diǎn)距離呈現(xiàn)拋物線形態(tài)，且受布站位置誤差的影響較大。在不存在站點(diǎn)定位誤差時(shí)，站點(diǎn)z坐標(biāo)為2.25km處，y方向誤差最??；在布站定位誤差為1m時(shí)，站點(diǎn)z坐標(biāo)為5m處，y方向誤差最小；站點(diǎn)定位誤差為3m時(shí)，站點(diǎn)z坐標(biāo)為5m處，y方向誤差最小。定義y方向誤差最小的點(diǎn)為最優(yōu)布站點(diǎn)。

（3）最優(yōu)布站點(diǎn)z方向距離定義為z＊，當(dāng)站點(diǎn)z向距離z＜z＊時(shí)，y方向誤差隨z減少快速增加；當(dāng)站點(diǎn)z向距離z＞z＊時(shí)，x方向和z方向誤差隨z增加快速增加。

圖11 異側(cè)布站時(shí)不同站點(diǎn)位置精度下站點(diǎn)z方向距離與測(cè)量誤差關(guān)系Fig.11 The relationship between measure error z distance for different side station position optimization

文獻(xiàn)［1］表明對(duì)定點(diǎn)測(cè)量時(shí)，當(dāng)站點(diǎn)交匯角在90°左右交匯測(cè)量的精度最高，使用Monte Carlo法計(jì)算表明針對(duì)特定的彈道軌跡，在考慮經(jīng)緯儀的各項(xiàng)誤差源以及布站位置誤差時(shí)，交匯測(cè)量精度將呈現(xiàn)更為復(fù)雜的關(guān)系，Monte Carlo法可以為光電經(jīng)緯儀的布站優(yōu)化提供指導(dǎo)。

使用Monte Carlo法進(jìn)一步繪制了站點(diǎn)位置（x，z）變化時(shí)異側(cè)對(duì)稱布站情形下三維誤差分布曲面如圖12所示，使用此方法可以快速計(jì)算得到最優(yōu)布站位置為（5km，20km，0）附近。如果給定一定的靶場(chǎng)地形信息，使用此方法可以針對(duì)特定靶場(chǎng)地形給出最優(yōu)布站方式。

圖12 站點(diǎn)位置精度1m時(shí)，x，z，y方向測(cè)量誤差分布Fig.12 The x，z，y measure error distribution with 1mthe station position precision

4 結(jié) 論

對(duì)光電經(jīng)緯儀建立了包含各種誤差源的準(zhǔn)確的數(shù)學(xué)模型，使用此模型進(jìn)行最壞情況分析和Monte Carlo法分析，分析了各種誤差源影響系統(tǒng)性能的程度和方式。進(jìn)一步使用此模型對(duì)雙站異側(cè)布站情形進(jìn)行了優(yōu)化，分析表明針對(duì)沿x方向發(fā)射的彈道軌跡，飛行段x方向中部位置布站可以獲得最小測(cè)量誤差；布站z方向存在最優(yōu)位置z＊，此處y方向交匯測(cè)量誤差最小，當(dāng)布站z向距離z＜z＊時(shí)，y方向誤差隨z減少快速增加；當(dāng)布站z向距離z＞z＊時(shí)，x方向和z方向誤差隨z增加快速增加。提出的基于Monte Carlo法的布站優(yōu)化可以進(jìn)一步推廣到特定靶場(chǎng)地形的情況，對(duì)經(jīng)緯儀總體設(shè)計(jì)及布站方式的選擇具有一定的理論指導(dǎo)意義。

［1］趙學(xué)顏，李迎春.靶場(chǎng)光學(xué)測(cè)量［D］.北京：裝備指揮技術(shù)學(xué)院，2001：35－38.

［2］王家騏.光學(xué)儀器總體設(shè)計(jì)［D］.長(zhǎng)春：長(zhǎng)春光機(jī)所研究生部，2003：27－32.

［3］杜俊峰.光電經(jīng)緯儀測(cè)量精度指標(biāo)的確定［J］.應(yīng)用光學(xué)，2006，5（6）：506－509.

［4］張 寧，沈湘衡.應(yīng)用跟蹤誤差等效模型評(píng)價(jià)光電經(jīng)緯儀跟蹤性能［J］.光學(xué)精密工程，2003，18（3）：677－684.

［5］李 慧，沈湘衡.光電經(jīng)緯儀數(shù)字化模型的建立方法及其應(yīng)用［J］.光電子技術(shù)，2007，27（2）：97－100.

［6］侯宏錄，周德云.光電經(jīng)緯儀異面交會(huì)測(cè)量及組網(wǎng)布站優(yōu)化設(shè)計(jì)［J］.光子學(xué)報(bào)，2008，37（5）：1023－1028.

［7］郭麗華，朱元昌.基于遺傳算法的光電經(jīng)緯儀布站優(yōu)化設(shè)計(jì)［J］.儀器儀表學(xué)報(bào)，2010，5（4）：741－746.