星光成像的大氣影響研究（Ⅲ）：大氣折射

陶志煒，戴聰明，武鵬飛，任益充，梅海平，馮云松，饒瑞中，魏合理

（1 中國科學院合肥物質科學研究院 安徽光學精密機械研究所 中國科學院大氣光學重點實驗室， 合肥 230031）（2 國防科學技術大學 電子對抗學院 紅外與低溫等離子體安徽省重點實驗室， 合肥 230037）（3 先進激光技術安徽省實驗室， 合肥 230037）（4 合肥工業大學 物理學院， 合肥 230601）

0 引言

人類利用恒星來進行導航，最早可追溯到古代人們通過北極星來確定方位。直到20 世紀50年代，星敏感器的橫空出世，大大提升了恒星導航的精度。星敏感器是一種高精度的姿態敏感測量儀器，它通過成像系統對星空成像，測量恒星矢量在星敏感器坐標系中的分量，利用已知的恒星精確位置來確定載體相對于慣性坐標系的三軸姿態［1］。星敏感器導航技術的高精度、抗干擾性強、可不依賴其它系統進行獨立導航等優點在近地空間各類機載、艦載、車載平臺上有著廣泛的應用［2］。

早期，星敏感器應用于衛星平臺，承擔了衛星姿態測量的任務，是衛星平臺不可或缺的測量設備。由于大氣層外可近似看成是真空環境，因此，星載星敏感器在工作時幾乎不受大氣的影響，其測姿精度通?？蛇_到角秒級甚至亞角秒量級。隨著空天觀測平臺的發展以及大氣層內星敏感器觀測高度的降低，大氣對恒星探測的影響愈發顯著，星敏感器在成像時不可避免的會受到地球大氣背景輻射、湍流、折射等因素的影響。白天大氣分子與氣溶膠粒子會對太陽光產生散射，使得天空產生復雜的背景光，對星敏感器在白天觀星時產生極大的干擾，嚴重降低成像的信噪比，使得恒星目標湮沒在背景中，無法識別［3］。另一方面，大氣分子與氣溶膠粒子的存在會使得星光傳輸產生衰減，再加上大氣本身存在折射率起伏，還會影響光波的相位，并對振幅產生調制，引發星光的閃爍和抖動［4-10］。由于大氣折射率隨高度存在一定變化，因此當星光經過大氣傳輸時，大氣折射會延長星光的傳輸路徑，改變星光的傳輸方向，嚴重影響星光導航的精度。

為了定量評估大氣對星光成像的影響，本系列文章Ⅰ、Ⅱ著重研究了背景輻射導致星光“看不見”以及大氣湍流導致星光“看不準”這兩個重要問題。為了進一步精確獲取恒星的位置，采用星敏感器進行定位和導航，本文針對導致星光“看不準”中另一重要因素大氣折射展開深入的研究。關于星光大氣折射影響的研究最早可追溯到大氣折射計算模型［11］，該模型中大氣折射使天體的表面位置向天頂移動。隨后，國內外相關學者提出了眾多模型用于星光大氣折射的計算［12-18］，并劃分了幾個使用不同星光折射計算方法的區域［19］：當觀測天頂角小于70°時，計算僅需要溫度，壓力和濕度等參數便可獲得較為不錯的精度。當觀測天頂角大于70°且在地平線的20°和5°之間時，溫度梯度成為主導因素，計算需要使用文獻［20］的方法，并使用標準大氣的溫度梯度和觀測者的測量條件進行數值積分。當觀測角度在更接近地平線的地方，需要在數值積分中使用局部溫度梯度隨高度變化的實際測量值進行計算。

本文是星光成像大氣影響三部分研究（背景輻射、湍流及折射）的第三部分，目標在于選取最佳的星光折射計算數理模型，并針對我國幅員遼闊特點，選取典型地區，不斷獲取測試數據，優化完善我國典型區域星光折射的數據模型。首先，通過美國標準大氣的溫壓數據，結合Ciddor 折射率模型，從計算精度、迭代次數以及計算速度三個方面對比分析了不同折射計算模型的優缺點，選取了最佳折射計算模型。其次，采用Ciddor 折射率模型計算得到的我國典型地區的折射率廓線數據，結合最佳折射計算模型，研究了星光折射引起的折射角、色散角、橫向位移以及路徑延長因子和延時等物理量的變化規律，并深入分析了輸入參數不確定性對大氣折射角計算的影響。本文的研究成果有望為實現星敏感器的準確觀測，星光導航定位的精確校正提供理論和數據支撐。

1 理論模型

1.1 折射率模型

為了準確描述大氣引起的折射和色散現象，建立準確的折射率模型是必要的。2014年，BARRELL H等比較了最常見的折射率模型，以評估光學設計軟件Zemax OpticStudio 中的大氣表面模型的準確性［21］。與此同時，作者還證實了OpticStudio 使用了Barrell&Sears 折射率模型［22］中的一個過時的方程式。與CIDOR P E 的最新工作相比［23］，在I 波段（800～934 nm）和天頂角為30°的情況下，色散角可能存在0.8 mas 的差異。另一方面，若將Ciddor 模型與Birch&Downs 模型［24］和Bonsch&Potulski 模型［25］相對比，此時誤差可能會縮小到0.2 mas 左右。本文選取Ciddor 模型作為計算折射率計算中最為精確的模型，因為其被國際大地測量學協會認為是計算大氣折射率最為標準的方程［26］。除此之外，該模型對于本文的計算波段范圍內的溫度、濕度、壓強都是準確的。其基本原理是假設大氣折射率n與大氣密度ρ存在Lorentz-Lorenz 關系，即［23，27］，通過將n和ρ換算成一組參考條件，計算出感興趣的大氣條件的折射率。關于Ciddor 模型介紹讀者請詳見文獻［23］，這里不做過多的介紹。若考慮更長波段，例如中紅外波段，或者OH和H2O 吸收線開始影響折射率變化的情況，這時采用Mathar 模型［28］可能是不錯的選擇，關于Mathar 模型的實際檢驗，可見SKEMER A J 等在中紅外天空中的測試［29］。

1.2 折射計算模型

星光大氣折射計算的精度不僅取決于折射率計算模型的精度，同時也依賴于折射計算模型的準確性。本節將構建星光成像大氣折射影響的折射角和色散模型，大氣折射角引起的橫向位移模型以及路徑延長模型等等。

1.2.1 折射角和色散

折射角和色散的計算依賴于模型幾何和大氣參數的選取。針對不同的模型幾何，可將地球大氣分為平面平行大氣和球形大氣。如圖1 所示，平面平行大氣通常是假設天頂角z0≈η。該模型僅對于z0→0 的情況下成立，隨著z0的增大，折射角計算引起的誤差也越大。平面平行大氣條件下折射角和色散角的計算公式表示為［30］

圖1 各向同性球形整層大氣折射示意圖Fig.1 Diagram illustrating refraction in a homogeneous spherical whole-layer atmosphere

式中，ζ代表光線折射前的角度，n0代表對應觀測高度的大氣折射率，λ代表波長。值得注意的是，λ對于折射角和色散角的影響由折射率模型決定。針對球形大氣，國內外學者提出眾多模型和計算方法，例如Cassini模型［31］，Corbard 誤差函數模型［32-34］，Mathar 氣壓指數模型［35］，Oriani 定理［36］以及基于SLALIB 包的計算方法［37］等等，下面將主要介紹Cassini 模型（原因見模型精度小節中算法速度分析），其余模型詳見附錄。

當z0≈η不成立時，各向同性球形大氣條件下折射角計算公式為

式中，r⊕代表折算地球半徑，r⊕的大小由緯度和海拔高度共同決定（具體計算方法可參考文獻［38，39］），hr則代表折算高度。在各向同性大氣中，僅需要觀測高度所對應的大氣密度和壓強數據便可計算出hr，表示［30］為

式中，p0和ρ0分別代表觀測高度處對應的大氣壓強和密度，g代表重力加速度。若考慮等溫大氣，hr通常由式（4）決定［30］。

式中，kb代表Boltzmann 常數，T0代表觀測高度對應處的溫度，m代表空氣分子的平均質量。

值得注意的是，上述介紹的折射角和色散計算模型無論對于哪種模型，都是將大氣看成整層來進行處理，實際上，這種處理方法一定程度上會喪失計算精度，因此，下面我們介紹兩種折射角計算方法，他們通常是將大氣劃分成若干層來進行計算的，其在計算速度上稍有欠缺（更多關于算法精度和速度的比較見下一小節）。一般而言，大氣分層的個數是不確定的，一般以前后兩次折射角的差值小于0.01″作為終止迭代的標準［40］?？紤]到低層大氣對于折射計算貢獻較大，因此每層大氣所對應高度一般采用式（5）進行計算［41］。

式中，ri代表每一層的厚度，t=i/N，i=1，2，…，N代表節點數，N代表總層數，r′⊕=r⊕+h0，h0代表觀測高度，A=r′⊕sinz0。若假設劃分為若干層后的大氣每一層的折射率以及折射率梯度（等同于光線曲率）為常數，將計算方法劃分為等折射率光線追跡法和等曲率光線追跡法［40］。對于等折射率光線追跡法而言，通常是假設光線是在每一層的邊界處發生偏折，通過Snell 折射定律，特別地，對于第i層和第i+1 層，滿足ni(hi+r⊕)sinθi=ni+1(hi+1+r⊕)sinθi+1，其中hi+1=hi+ri，如圖2（a）所示，建立相鄰兩層的折射角度變化關系。當上述關系遍歷所有層后，得到折射角計算公式［42］

圖2 球形多層大氣折射示意圖Fig.2 Diagram illustrating refraction in a spherical multiple-layer atmosphere

對于等曲率光線追跡法，通常是假設在每一層內折射率梯度為常數，即曲率半徑在每一層的邊界處發生變化，特別地，如圖2（b）所示，對于第i層而言，曲率半徑滿足1/ki=?sinφidni/nidhi，其中ki代表第i層的曲率半徑，φi代表光線沿第i層傳播方向與該層半徑矢量方向的夾角。一般地，曲率半徑還可表示為微分方程［40，43］

除此之外，通過圖2（b）中的幾何關系，還可得到微分方程組［43］

式中，?和d?的定義詳見圖2（b）。結合等式（7）和等式（8），通過圖2（b）中的幾何關系設置初值，使用四階Runge-Kutta 法便求解上述一階非線性微分方程組，得到光線發生折射前入射角度z∞。至此，等曲率光線追跡法計算得到的折射角可表述為

1.2.2 橫向位移

上節我們通過求解微分方程組獲得了光線經過大氣折射后所產生的折射角和色散角。一般地，若要通過星敏感器確定恒星位置，通?？梢酝ㄟ^恒星視在位置和真實位置所產生的折射角進行彌補。進一步，若要通過星敏感器確定大氣中其他目標時，特別地當目標位置處于地球表面和低仰角時，若此時按原有方法計算其位置便會產生一定的誤差，使得定位的精度產生偏差，因此有必要計算該誤差，以彌補和精準確定目標所在位置［43］。以光線入射到大氣層前的方向為基準，當目標位于一定高度時，真實光線方向與折射前光線方向存在一定的距離偏差，我們將此位移稱為橫向位移b，如圖2（b）所示。通過圖中的幾何關系可得到微分方程［43］

結合式（7）、式（8）和式（9），通過圖2（b）中的幾何關系設置初值，其中?∞，φ∞代表根據上小節設置的初值求解得到的大氣最外層的?，φ，使用四階Runge-Kutta 法便可求解上述一階非線性微分方程組，得到任意觀測高度的橫向位移b(∞)(h0)。

1.2.3 路徑延長因子和路徑延時

如圖2 所示，星光經過大氣傳輸后光線會因大氣折射現象發生彎曲，延長了星光的傳輸路徑，造成一定程度的路徑延時［44］。本小節將通過等折射率光線追跡法構建星光傳輸的幾何模型，計算星光經大氣傳輸后的實際路徑長度，給出大氣折射引起的路徑延長因子的計算公式。如圖3 所示，大氣按照式（5）被劃分為若干層，以星光真實位置點S為起點，觀測者位置O為終點連線，兩者之間的距離代表星光的斜程距離，可通過式（11）得到。

圖3 大氣折射引起的星光路徑延長示意圖Fig.3 Diagram illustrating the elongation of the optical ray trajectory caused by atmospheric refraction

根據圖3 中給出的幾何關系并結合正弦定律、Snell 折射定律，可得到

式中，Ci和Cmax表示為

δi和δN表示為

?i代表第i?1 層大氣到第i層之間的折射角，根據直接積分法，?i可近似解析表示為［45］

結合式（12）、（15），對圖3 中的三角形△CPi?1Pi，△CPN S使用余弦定律可得大氣每一層的光線路徑長度

根據式（16）、（17）可計算出光線經大氣傳輸后的實際路徑長度為，因此，星光經過大氣折射后所產生路徑延長因子和路徑延時分別可表示為

式中，c代表光速，c=299 792 458.0 m/s。

1.3 模型精度

1.3.1 美國標準大氣

地球的大氣層可以看作是一個球狀分層的介質，在不同高度層上有特定的折射率值分布。為了獲得Ciddor 模型下的折射率分布并通過該分布對比不同折射計算模型的精度，首先使用美國標準大氣模型進行計算。一般來說，美國標準大氣是一個理想化的穩定狀態的代表，它給出了大氣壓強、溫度和其他參數隨高度的變化情況，其中最大高度可達1 000 km［44］。美國標準大氣假設溫度隨高度是滿足線性分布的，而壓強隨高度的變化可以通過氣體定律和流體靜力學方程中獲得［46］，圖4（a）、（b）給出了美國標準大氣0～86 km 情況下的溫度廓線和壓強廓線，根據Ciddor 模型，可計算得到美國標準大氣折射率（圖4（c））及折射率梯度（圖4（d））隨高度的變化曲線。根據上節的大氣折射計算模型，計算得到美國標準大氣下不同模型的折射角隨天頂角的變化情況，從圖4（e）可以看到折射角總體是隨著天頂角的增大而不斷變大，當觀測仰角較小時，此時大氣折射的影響最為顯著。需要指出的是，不同模型計算的折射角差異在圖中并不顯著，不同模型的計算精度需放大對比。最后，圖4（f）和（g）還計算了美國標準大氣下星光經大氣產生的橫向位移、路徑延長因子和路徑延時，這都為計算我國典型地區不同時間段所對應的大氣折射參數提供了參考。

圖4 美國標準大氣下的大氣參數廓線和折射計算結果Fig.4 Atmospheric parameter profiles and refraction calculations in the U.S.standard atmosphere

1.3.2 計算精度

圖5 是不同折射計算模型計算的折射角（包括等折射率光線追跡法、等曲率光線追跡法、Cassini 模型、Corbard 誤差函數模型、Mathar 氣壓指數模型、Oriani 定理、平面平行大氣模型以及基于SLALIB 包計算結果）與直接積分法［31，45，47］計算結果偏差隨天頂角的變化情況。這里以直接積分法作為基準進行比較是因為在美國標準大氣下使用直接積分法獲得折射角可近似替代Pulkovo 折射表［48］給出的折射角的標準值［49］。當越小，說明此時采用的折射計算模型與標準值越接近。如圖5 所示，平面平行大氣模型計算得到的折射角與標準值相差最大，說明使用該模型計算折射角和色散所產生的誤差也最大。除此之外，還可以看出當天頂角較小時，使用Oriani 三次方定理和五次方定理計算得到折射角與標準值相差最小，此時選取Oriani 定理計算星光經大氣所產生的折射角最為準確。另外，相比于Oriani 定理， 天頂角較小時Cassini 模型和等折射率光線追跡法計算結果也表現出較為精確的結果。當天頂角位于45°附近時，使用Cassini 模型計算得到的結果最為精確，此時，與Oriani三次方定理相比，五次方定理表現出明顯的優勢。然而，當天頂角不斷增大時，這時使用Oriani 定理計算大氣折射角將較Cassini 模型和等折射率光線追跡法產生較大的誤差。與Cassini模型相比，雖然在大天頂角時使用等曲率光線追跡法、Corbard 誤差函數模型以及SLALIB 包計算將取得更為精準的計算結果，但其上述三種模型在天頂角較小將引入較大的計算誤差。綜上所述，若選取等折射率光線追跡法和Cassini模型來計算星光大氣折射所產生的折射角和色散角，可獲得較為不錯的計算精度。

圖5 美國標準大氣下不同折射模型計算的折射角與直接積分法計算結果的偏差隨天頂角的變化曲線Fig.5 The deviation of the refraction angle calculated by different refraction models and the direct integration as a function of zenith angle in the U.S.standard atmosphere

1.3.3 迭代次數

為了進一步從迭代次數角度對比多層大氣折射計算模型（等折射率光線追跡法和等曲率光線追跡法）之間的區別，圖6 給出了在天頂角為85°情況下使用的等折射率、等曲率光線追跡法計算的折射角隨大氣層數的變化情況（這里選取85°進行計算，是因為在美國標準大氣下根據這兩種大氣折射計算模型得到折射角與標準值誤差基本一致，因此排除了計算誤差的影響）。需要說明的是，圖6 的計算是在假設hmax= 86 km 的基礎上進行的，由式（7）和（8）可以看出，若要使用等曲率光線追跡法求解折射角和色散角，通常需要折射率和折射率梯度曲線隨高度變化的函數關系的。為此，盡管本文的驗證計算只使用了86 km 以下的大氣作為折射影響的主要因素，但在等曲率光線追跡法的計算當中仍使用美國標準大氣下1 000 km 以下的數據作為生成n(h)和dn( )h/dh函數的來源。從圖6（a）中可以看出，隨著大氣層數的不斷增大，使用兩種多層大氣折射計算得到折射角逐漸趨于常數，相比于等折射率光線追跡法，等曲率光線追跡法收斂的速度較快。與等折射率光線追跡法不同的是，當大氣層數較小時，等折射率計算結果隨著層數的增大單調遞增，而等曲率法計算結果卻出現明顯的振蕩，這可能是由于數值求解微分方程組帶來的計算誤差導致的。若考慮兩者的迭代誤差，如圖6（b）所示，等曲率法較等折射率法而言，在大氣層數較小時，其迭代誤差就基本趨向于0，說明等曲率光線追跡法在迭代次數方面較等折射率法表現出明顯的優勢。然而，值得說明的是，等曲率光線追跡法依賴折射率梯度的精度，需要超過最大觀測高度以上的折射率和折射率梯度的數據（通常采用指數擬合，但會影響折射角計算的精度，最好有實測數據），另一方面，等曲率光線追跡法在處理實測數據時，需要對觀測高度以下的折射率梯度數據進行分段平均，這時采用的平均方法（分段平均，滑動平均）也將很大程度決定折射計算的精度，但好處是相對于等折射率光線追跡法而言，大氣分層個數較少，需要的迭代次數也較少?？傮w來看，如圖5 所示，兩種算法計算精度基本接近，盡管等曲率光線追跡法需要較少的大氣分層便能獲得不錯的結果，但由于其不確定性因素較多且計算較為復雜，因此采用等折射率光線追跡法進行計算依然是最佳的選擇。

圖6 等折射率、等曲率光線追跡法計算的折射角隨大氣層數的變化曲線Fig.6 The variation curve of refraction angle as a function of the number of atmospheric layers， which is calculated by the method of the equivalent refractive index ray tracing and the equivalent curvature one

1.3.4 算法速度

第一小節從計算精度的角度對比了不同整層大氣折射模型和多層大氣折射模型的計算精度，然而面對多經緯度、多時段、多個觀測仰角情況下折射和色散計算的應用場景時，選取一個速度較快且精度較高的計算方法對于實現實時星光導航尤為重要。因此，圖7 左圖對比了不同折射計算模型的運算時間隨計算天頂角個數的變化情況?？梢钥吹?，隨著計算天頂角個數的不斷增加，運算時間不斷增加，其中多層大氣折射計算模型運算時間較其他模型較長，這可能是由于多層大氣折射模型的多次迭代計算造成的。對比等折射率法，可以發現等曲率法運算時間更長，這說明使用4 階Runge-Kutta 的計算程序仍有待優化，后續可考慮使用Pytorch 庫［50］在GPU 上進行運算來提升求解方程組的計算速度。與多層大氣折射模型相比，基于SLALIB包計算速度相比其他解析模型最慢，這可能因為SLALIB 包計算需要從python 中不斷調用Fortran 代碼導致的［37］。圖7 右圖給出了不同解析折射計算模型在計算91 個天頂角后的運算時間對比，可以看到Cassini 模型和平面平行大氣模型相比其他計算模型速度最快，而Mathar 氣壓指數模型相比其他模型速度最慢，Oriani三次方定理和五次方定理計算速度基本一致。因此，綜合計算精度、計算速度和迭代次數，本文選取等折射率光線追跡法和Cassini 模型作為折射計算的模型。

圖7 不同折射計算模型的運算時間隨天頂角個數的變化和91 個天頂角情況下不同理論折射計算模型的時間直方圖Fig.7 The variation curves of the calculation time of different refraction calculation models as a function of the number of zenith angles and the time histograms of different theoretical refraction calculation models in the regime of calculating 91 zenith angles

2 我國典型區域折射計算結果

前一節我們介紹了不同的折射率模型以及折射計算模型，并分析和對比了不同折射計算模型的精度和速度，本節使用上一節選取的折射率模型以及折射計算模型來計算我國典型地區不同時刻、不同觀測高度、不同觀測天頂角以及不同波長下大氣折射引起的折射角、色散角、橫向位移、路徑延長因子以及路徑延時等物理量，為了解我國典型地區不同時刻的大氣折射規律，實現星敏感器的準確觀測、星光導航定位的精確校正提供數據支持。圖8 是我國西北地區（以大柴旦為典型代表，90°10′ E～96°22′ E， 37°35′ E～39°12′ E，海拔2 829～5 655 m）和東南沿海（以臺州為典型代表，124°34′ E， 28°50′ N，海拔162 m）某一天早晨和夜晚的實時觀測數據繪制的大氣溫度、壓強以及濕度隨著高度的變化曲線，可以看到西北地區和東南沿海大氣溫度隨高度整體呈現先降低后升高，且溫度轉折點通常位于對流層頂。另一方面，還可以發現早晨溫度廓線明顯小于夜晚溫度廓線。除此之外，從圖8 可以看到兩地大氣壓強隨高度整體呈現單調指數遞減趨勢，且早晚大氣壓強無明顯差別，但相比于西北地區，東南沿海大氣壓強要高于西北地區。最后，從圖8 最右側一列，可以看出東南沿海大氣濕度隨高度整體呈現不斷變小趨勢，然而大柴旦早晚的濕度廓線并無明顯規律，這可能是由于當日實測過程中天氣變化所致，導致大氣濕度在對流層出現反常增大。圖9 是根據Ciddor 模型計算得到的折射率和折射率梯度的廓線，相比于美國標準大氣，實測數據計算得到的梯度曲線抖動更加劇烈，為此我們對兩地折射率梯度廓線進行分段平均處理。

圖8 西北地區、東南沿海某天早晨和夜晚的大氣溫度廓線、壓強廓線、相對濕度廓線Fig.8 Atmospheric temperature profile， pressure profile and relative humidity profile in the morning and night of a day of the northwest area and southeast coast of China

圖9 根據Ciddor 模型計算的西北地區、東南沿海某天早晨和夜晚的折射率廓線，折射率梯度廓線Fig.9 Refractive index profiles， refractive index gradient profiles calculated from the Ciddor model in the morning and night of a day of the northwest area and southeast coast of China

2.1 折射角

2.1.1 不同觀測高度

圖10 是使用Cassini 模型和等折射率光線追跡法計算的西北地區、東南沿海早晨和夜晚不同觀測高度情況下（即0 km，1.5 km，3 km，4.5 km，6 km ，7.5 km，9 km 和10 km，這里指的是相對于海平面以上的觀測高度）折射角隨天頂角的變化情況，本節計算在λ=0.9 μm 的基礎上進行的?？梢钥吹讲煌Ｐ?、不同地區和不同時間段折射角整體隨天頂角的增加而不斷變大，且Cassini 模型和等折射率光線追跡法的計算結果基本一致，僅僅在西北地區7.5 km～9 km 處出現幾個微角秒的偏差，這可能是由于計算誤差造成的。對比不同地區折射角的計算結果，可以看出東南沿海折射角明顯大于西北地區，但整體變化趨勢基本一致。進一步，對比兩地不同時刻的折射角差異，還可以看到隨著天頂角的增大，早晨和夜晚折射角出現明顯的差異，這種差異的變化范圍大約在5～10 個微角秒的數量級上變化。最后，我們還可看出西北地區早晚折射角差異要略大于東南沿海，但這種差異隨著觀測高度的增大逐漸消失，隨著觀測高度的增加，早晚折射角以及折射角差異隨天頂角的變化也將逐漸減小，這是因為隨著觀測高度，光線因大氣折射所走過的路徑長度和彎曲程度不斷變小，從而導致折射角逐漸減小導致的。

圖10 使用Cassini 模型和等折射率光線追跡法計算的西北地區、東南沿海早晨和夜晚不同觀測高度情況下折射角隨天頂角的變化曲線Fig.10 The variation curve of refraction angle as a function of zenith angle for different observation altitudes in the morning and night of a day of the northwest area and southeast coast of China using the Cassini model and the equivalent refractive index ray tracing method

2.1.2 不同波長

圖11 是使用Cassini 模型和等折射率光線追跡法計算的西北地區、東南沿海早晨和夜晚不同波長情況下（即0.9 μm，1 μm，1.1 μm，1.2 μm，1.3 μm，1.4 μm，1.5 μm 和1.65 μm）折射角隨天頂角的變化情況，本節計算的觀測高度為3 km?？梢钥吹讲煌ㄩL情況下使用Cassini 模型和等折射率光線追跡法計算的折射角基本無偏差。除此之外，還可以看到相比于西北地區，東南沿海的折射角要略大，不同波長下早晚折射角差異隨天頂角的增大不斷增大，且西北地區早晚折射角之差略大于東南沿海，這與上小節的規律基本一致。對比不同觀測波長折射角隨天頂角的變化曲線，可以看出折射角隨觀測波長在角秒量級無明顯變化，需進一步調整變化單位進行比較。圖12 是使用Cassini 模型和等折射率光線追跡法計算的西北地區、東南沿海早晨和夜晚不同觀測波長下兩波長計算的折射角偏差隨天頂角的在微角秒量級上變化情況（其中λ1=1 μm，1.1 μm，1.2 μm，1.3 μm，1.4 μm，1.5 μm 和1.65 μm，λ2=0.9 μm），可以看到不同波長下兩波長折射角偏差隨天頂角的增加不斷增大，且相比西北地區而言，東南沿海在不同波長下的折射角偏差較西北地區偏大。除此之外，從圖12 中可以明顯的看到隨著波長的增大，不同波長計算得到折射角差異也逐漸變大，但這種變化是在微角秒量級，這可能是因為波長對折射角的影響是通過折射率間接影響而導致的。

圖11 使用Cassini模型和等折射率光線追跡法計算的西北地區、東南沿海早晨和夜晚不同波長下折射角隨天頂角的變化曲線Fig.11 The variation curve of refraction angle as a function of zenith angle for different values of wavelength in the morning and night of a day of the northwest area and southeast coast of China using the Cassini model and the equivalent refractive index ray tracing method

2.2 輸入參數不確定性對折射角計算的影響

由前文的分析的可知，折射計算的準確性和精度不僅取決于折射率模型和折射計算模型的準確性，同樣取決于輸入數據的精度和準確度，然而在實際大氣中進行探測測量時，所得的大氣參數廓線將不可避免地受到隨機噪聲的影響，從而影響最終折射計算的準確度。為了定量分析大氣參數廓線隨機噪聲誤差對實際折射計算的影響，精確計算折射角計算誤差和隨機噪聲之間的關系，本小節分別溫度、壓強、相對濕度、CO2濃度以及波長等角度計算了輸入參數不確定對折射角計算的影響。圖13～16 分別給出了使用Cassini模型和等折射率光線追跡法計算的西北地區、東南沿海早晨和夜晚不同參數對應的不同輸入參數誤差情況下折射角誤差隨天頂角的變化情況。需要說明的是，本節計算假設溫度、壓強和相對濕度廓線所帶來的誤差滿足高斯隨機的正態分布，這也滿足實際測量中大氣參數廓線所附加的誤差類型。

對比圖13～16 中Cassini 模型和等折射率光線追跡法的計算結果，可以看出，使用Cassini 模型和使用等折射率光線追跡法計算得到輸入參數不確定性對折射角的影響規律基本一致，只有在計算不同壓強誤差對折射角的影響表現出異常，可以看到使用Cassini 模型計算得到折射角誤差明顯高于等折射率光線追跡法，且誤差是后者的好幾倍。除此之外，對比圖13 和圖14 西北地區早晚折射角的誤差數據，可以看出在較大的天頂角時，輸入溫度參數不確定性對早晨折射角計算影響偏大，輸入相對濕度參數不確定性對夜晚折射角計算的影響偏大。另一方面，對比圖15 和圖16 東南沿海早晚折射角的誤差數據，可以看出輸入相對濕度參數不確定性對早晨折射角計算的影響偏大?？傮w而言，輸入參數誤差對西北地區折射計算的影響要明顯小于東南沿海。最后，從圖13～16 整體折射角誤差的變化情況可以看到，隨著觀測天頂角和輸入參數誤差（包括溫度、壓強、相對濕度、CO2濃度以及波長）的增加，折射角計算所帶來的誤差也逐漸增加，這一點是不言而喻的。

圖14 使用Cassini 模型和等折射率光線追跡法計算的西北地區夜晚不同參數對應的不同輸入參數誤差情況下折射角誤差隨天頂角的變化曲線Fig.14 The variation curve of refraction angle error as a function of zenith angle under different input parameter errors corresponding to different parameters in the evening of the northwest area of China calculated using the Cassini model and the equivalent refractive index ray tracing method

圖15 使用Cassini 模型和等折射率光線追跡法計算的東南沿海早晨不同參數對應的不同輸入參數誤差情況下折射角誤差隨天頂角的變化曲線Fig.15 The variation curve of refraction angle error as a function of zenith angle under different input parameter errors corresponding to different parameters in the morning of the southwest coast of China calculated using the Cassini model and the equivalent refractive index ray tracing method

圖16 使用Cassini 模型和等折射率光線追跡法計算的東南沿海夜晚不同參數對應的不同輸入參數誤差情況下折射角誤差隨天頂角的變化曲線Fig.16 The variation curve of refraction angle error as a function of zenith angle under different input parameter errors corresponding to different parameters in the evening of the southwest coast of China calculated using the Cassini model and the equivalent refractive index ray tracing method

2.3 色散

圖17 是使用Cassini 模型和等折射率光線追跡法計算的西北地區、東南沿海早晨和夜晚不同觀測高度情況下（即0 km，1.5 km，3 km，4.5 km，6 km，7.5 km，9 km 和10 km，這里指的是相對于海平面以上的觀測高度）色散角隨天頂角的變化情況。需要說明的是，本節計算的0.9 μm 和1.65 μm 之間的色散現象。從圖中可以看出東南沿海不同觀測高度層早晚的色散情況要明顯大于西北地區，另一方面，隨著觀測高度的增加，兩地不同時刻的色散現象明顯減弱，這與上一小節折射角的變化規律是一致的。除此之外，從圖中還可以看出，不同高度層西北地區的早晚色散角差異要大于東南沿海，且隨著觀測高度的增加，這種趨勢逐漸喪失。最后，與不同地區早晚折射角差異隨觀測高度變化規律一致的是，觀測高度越低，兩地早晚色散差異越大，這可能是由于色散差異隨著傳輸路徑增大不斷累積導致的。

圖17 使用Cassini 模型和等折射率光線追跡法計算的西北地區、東南沿海早晨和夜晚不同觀測高度情況下色散角隨天頂角的變化曲線Fig.17 The variation curve of dispersion angle as a function of zenith angle for different observation altitudes in the morning and night of a day of the northwest area and southeast coast of China using the Cassini model and the equivalent refractive index ray tracing method

2.4 橫向位移

2.4.1 不同觀測高度

第2 節介紹了星光經大氣折射后產生橫向位移的概念，表明了確定橫向位移的大小有利于更精確的使用星敏感器定位大氣層內的目標位置。本文將基于我國典型地區的一次探空數據，計算不同觀測高度情況下星光折射橫向位移的大小。圖18 給出了西北地區和東南沿海早晨和晚上不同觀測高度下橫向位移b(∞)(h0)隨天頂角的變化情況。在給出計算結果之前，需要說明的是，由式（7）、（8）和（10）可以看出，使用微分方程組求解橫向位移時，需要折射率和折射率梯度曲線隨高度變化的函數關系。為此，盡管橫向位移的計算只使用西北地區和東南沿海30 km 以下的實測數據作為折射影響的主要因素，如圖9 所示，但在使用四階Runge-Kutta 法計算橫向位移時通常需要高度超過30km 的數據作為生成n(h)和dn(h)/dh函數的來源。因此，在本節計算過程中，當h<30 km 時，我們使用實測氣象探空數據根據Ciddor 模型計算得到的折射率和分段平均后折射率梯度曲線的線性插值函數作為生成n(h)和dn(h)/dh函數的來源（不同平均方法影響計算精度，經計算分段平均計算得到的折射角較滑動平均更接近直接積分法得到的結果）；當h>30 km 時，由于折射率和折射率梯度隨高度變化通常滿足指數衰減規律，因此本文采用h<30 km 時折射率和分段平均后的折射率梯度曲線的指數擬合函數作為生成n(h)和dn(h)/dh函數的來源（在美國標準大氣的驗證計算中，我們使用了1 000 km 以下數據作為生成梯度函數的來源，這時計算的折射角更為接近標準值，這也是使用微分方程組進行折射計算的弊端，它通常需要更多的實測數據才能獲得更高的精度）。除此之外，由圖9（d）可以看出，東南沿海夜晚的折射率梯度曲線在2 km 高度以下存在數值躍變，因此在對東南沿海夜晚折射率梯度數據進行擬合時，采用了2 km 以上的數據進行指數擬合，此時計算得到折射角與Cassini 模型給出的結果基本接近，驗證該方法的可行性。

圖18 西北地區和東南沿海早晨和晚上不同觀測高度下橫向位移隨天頂角的變化曲線Fig.18 The variation curves of lateral shift as a function of zenith angle under different observation heights in the morning and evening of the northwest area and southeast coast of China

從圖18 中可以看出西北地區和東南沿海早晚不同觀測高度下的橫向位移整體與天頂角呈現正相關關系，且隨著天頂角的增大，橫向位移增加的速率也不斷變大，尤其當天頂角為85°時，橫向位移呈現量級上的增長。與此同時，通過對比西北和東南沿海兩地橫向位移的大小，可以發現西北地區折射引起橫向位移的大小較東南沿海地區偏大。另一方面，對比不同高度層（即0 km，1.5 km，3 km，4.5 km，6 km，7.5 km，9 km、10 km 和15 km，這里指的是相對于海平面以上的觀測高度）的觀測結果，可以發現兩地觀測的橫向位移的大小隨著觀測高度的增加不斷減小，這也是合乎情理的，因為觀測高度越高，星光因折射而產生的彎曲效應也越小。值得注意的是，當天頂角較小時，觀測高度從0 km～15 km 變化時橫向位移基本減小了一個數量級。除此之外，對比兩地不同高度下早晚橫向位移隨天頂角的變化規律可以看出，夜晚相比早晨折射產生的橫向位移偏大，隨著觀測高度的增大，早晚兩時橫向位移產生的偏差逐漸減小，且當天頂角較小時，早晨的橫向位移會逐漸大于夜晚的結果。

2.4.2 不同波長

圖19 是西北地區、東南沿海早晨和夜晚不同觀測波長下兩波長計算的橫向位移偏差隨天頂角的在微米量級上變化情況（其中λ1=1 μm，1.1 μm，1.2 μm，1.3 μm，1.4 μm，1.5 μm 和1.65 μm，λ2=0.9 μm，且本節計算在觀測高度為3 km 的基礎上進行的），可以看到不同波長下西北、東南兩地早晨和夜晚的兩波長橫向位移偏差整體隨天頂角的增加而不斷變大。除此之外，對比兩地早晨和夜晚的兩波長位移偏差數據可以看出當天頂角較大時，夜晚較早晨而言兩波長橫向位移的偏差更大。另一方面，對比西北地區和東南沿海的結果，可以看到西北地區波長引起的橫向位移偏差要大于東南沿海，說明西北地區波長對橫向位移的影響更大。最后，通過對比不同波長下兩波長產生的橫向位移偏差，可以看到隨著波長的增大，不同波長計算得到橫向位移差異也逐漸變大，且當天頂角較小時，這種差異較為明顯。值得注意的是，這種變化是在厘米量級進行的，這可能也是由于波長對折射角的影響是通過折射率間接影響導致的。

圖19 西北地區、東南沿海早晨和夜晚不同觀測波長下兩波長計算的橫向位移偏差隨天頂角的變化曲線Fig.19 The deviation of lateral shift calculated at two wavelengths as a function of zenith angle under different values of wavelength in the morning and night of the northwest area and southeast coast of China

2.5 路徑延長因子和路徑延時

2.5.1 不同觀測高度

大氣折射使得光線經過大氣傳輸產生彎曲現象，延長了光線傳輸的路徑，增加了光束信號的延時，是星光導航定位技術的主要誤差來源之一，為此，本小節針對我國典型地區不同時段，不同觀測高度計算了大氣折射引起的路徑延長因子以及路徑延時，并針對不同區域、時間和高度下路徑延長因子給出了擬合得到的多項式函數。圖20 給出了西北地區、東南沿海早晚不同觀測高度層下路徑延長因子和路徑延時隨天頂角的變化情況，通過計算得到的數據，給出了兩地不同高度層下早晚路徑延長因子隨天頂角變化的多項式擬合函數（以西北地區早晨和東南沿海地區夜晚觀測高度為3 km 時的數據為例）。

圖20 西北地區、東南沿海早晚不同觀測高度層下路徑延長因子和路徑延時隨天頂角的變化曲線Fig.20 The variation curves of path elongation ratio and path delay as a function of zenith angle under different observation heights in the morning and evening of the northwest area and southeast coast of China

西北地區早晨3 km時的擬合函數為

東南沿海夜晚3 km時的擬合函數為

由圖20 可以看出路徑延長因子和路徑延時隨天頂角的增大不斷增大，當天頂角較大時，這種增大最為劇烈。除此之外，路徑延時基本在微秒量級上變化，且西北地區和東南沿海不同高度下早晨和夜晚和路徑延長因子和路徑延時基本完全重合，通過早晚路徑延長因子和延時偏差隨天頂角的變化曲線可以看出，隨著天頂角的增加，延長因子和延時的偏差整體不斷增加。另一方面，可以看到延長因子的變化主要在10?8數量級以下范圍內變化，而延時主要在納秒量級左右上下波動，并且這種變化隨著觀測高度的增加量級上將變得更小，因此，可以指出不同高度下兩地早晚的延長因子和延時的差異可以幾乎忽略不計。對比西北地區和東南沿海的計算結果，可以看出東南沿海大氣折射造成的路徑延長比和路徑延時相比與西北地區較大。最后，還可以看出隨著觀測高度的降低，不同天頂角下延長因子和延時整體逐漸減小，這與之前折射角和橫向位移得到的結論是一致的。

2.5.2 不同波長

圖21、22 是西北地區、東南沿海早晨和夜晚不同觀測波長下兩波長計算的路徑延長因子偏差和延時偏差隨天頂角的變化情況（其中λ1=1 μm，1.1 μm，1.2 μm，1.3 μm，1.4 μm，1.5 μm 和1.65 μm，λ2=0.9 μm，且本節計算在觀測高度為3 km 的基礎上進行的），可以看到不同波長下西北、東南兩地早晨和夜晚的兩波長延長因子和延時的偏差整體隨天頂角的增加而不斷變大，且不同波長下兩波長延長因子偏差主要在10?10以下范圍進行波動，延時偏差主要在納秒量級上左右波動，因此可以認為波長對路徑延長因子和路徑延時的影響可以忽略不計。除此之外，還可看出兩地早晨和夜晚的兩波長延長因子和延時的偏差基本一致，這說明不同時段波長對延長因子和延時的影響并不差別。最后，從圖20 還可以看出東南沿海地區兩波長路徑延長因子偏差較西北地區偏大，這說明東南沿海波長對延長因子的影響更大，但由于波長對延長因子和延時的影響幾乎可以忽略不計，因此這種影響并不重要。

圖21 西北地區、東南沿海早晨和夜晚不同觀測波長下兩波長計算的路徑延長因子偏差隨天頂角的變化曲線Fig.21 The deviation of the path elongation ratio calculated at two wavelengths as a function of zenith angle under different values of wavelength in the morning and night of the northwest area and southeast coast of China

圖22 西北地區、東南沿海早晨和夜晚不同觀測波長下兩波長計算的路徑延時偏差隨天頂角的變化曲線Fig.22 The deviation of the path delay calculated at two wavelengths as a function of zenith angle under different values of wavelength in the morning and night of the northwest area and southeast coast of China

3 結論

本文選取最佳的星光大氣折射模型，著重研究了大氣折射對星光成像的影響。根據美國標準大氣給出的大氣參數廓線數據，計算了平面平行大氣、整層球形大氣以及多層球形大氣模型下的折射角的變化特性，對比分析了不同折射計算模型的計算精度和計算速度，選取了精度和速度最佳的折射計算模型。與此同時，針對多層球形大氣折射模型，研究了計算結果與迭代次數（即大氣層數）之間的關系，計算了美國標準大氣情形下大氣折射產生的橫向位移、路徑延長因子和路徑延時，驗證了多層球形大氣折射模型的可靠性和普適性。結合典型地區不同時段實測的大氣參數廓線數據，計算了不同觀測條件以及不同波長情況下的折射角、色散、橫向位移、路徑延長因子以及路徑延時的變化情況，評估了不同大氣參數由于輸入參數的不確定性對折射計算的影響。研究發現：1）使用Ciddor 折射率模型并結合Cassini 折射計算模型或等折射光線追跡法計算得到的折射角最為準確；2）相比于改變觀測波長而言，提升星敏感器的觀測高度或減小星敏感器的觀測天頂角能極大程度上減輕星光成像的大氣折射影響；3）當不同輸入參數存在噪聲和不確定性時，溫度參數不確定性對折射計算的影響最大。因此，為了能一定程度上減輕折射計算的誤差，有效提高溫度的測量精度相比修正其他參數的不確定性而言更為重要。

附錄

A.1 Corbard 誤差函數模型

Corbard 誤差函數模型假設大氣密度隨高度呈現指數衰減，表示為

式中，δn0=n0?1，t=hr/r⊕，為式（3）或（4）給出的大氣折算高度與地球半徑（為海平面觀測點對應的地球的曲率半徑）的比值，

A.2 Oriani 定理

Oriani 定理也稱為Oriani 三次方定理

式中，A和B分別代表常數。若對等式（A1）中的Ψ(x)采用Laurent 級數進行展開，可得到關于大氣折射角計算的另一個等式，即Oriani 五次方定理

A.3 Mathar 氣壓指數模型

Mathar 氣壓指數模型通過折射率（磁化率）隨高度呈指數衰減的大氣模擬望遠鏡位置上方的折射率，它通常指的是將折射角的直接積分公式Taylor 展開成奇數階的tanz0的函數

式中，σi(i=1，3，5，…)可根據磁化率隨高度呈指數衰減的假設計算得到

式中，μ0=n20?1，代表磁化率。