基于統計學方法的衛星在軌熱變形快速分析

張也弛，張龍，羅文波，潘騰

北京空間飛行器總體設計部，北京 100094

基于統計學方法的衛星在軌熱變形快速分析

張也弛*，張龍，羅文波，潘騰

北京空間飛行器總體設計部，北京 100094

為使用在軌實測溫度數據輔助完成成像載荷視軸指向的在軌標定，首次將應用于地統計學中的普通克里格法用于衛星在軌溫度場的無偏估計，提出了一種基于統計學方法的在軌熱變形快速分析方法。應用熱試驗過程中星上少量溫度實測數據進行衛星溫度場的無偏估計，并將溫度計算結果賦值到有限元模型，進而完成了熱變形分析。應用此種方法，相比通過熱物理方法進行溫度場分析，可將熱變形分析的工期縮短數日，且載荷變形的分析結果與試驗實測結果的誤差在1′左右，可較準確地得到載荷視軸在熱變形作用下的指向變化。文章可為衛星在軌熱變形的快速分析提供參考。

溫度實測點；快速熱變形分析；統計方法；普通克里格法；變差函數；有限元；衛星

某些成像載荷與光學成像載荷不同，例如X射線望遠鏡，難以通過對天空或地面已知目標成像來進行自身指向的標定，而是始終依賴衛星平臺提供的載荷光軸指向，利用載荷觀測數據進行圖像重建，最終確定X射線源的空間位置。因此，需要選擇適當的方法進行載荷視軸指向的標定，進而保證載荷的點源定位精度。指向精度的偏差由多種因素引發，其中影響最大、也最難進行標定的就是在軌溫度場變化引發的熱變形。通常，溫度場借助熱物理方法進行分析[1-2]，為了得到載荷在軌由于熱變形引起的實時指向變化，本文擬通過在軌溫度遙測數據對熱變形進行分析。

熱變形分析首先需要溫度場作為輸入條件，但相比有限元分析所需的溫度場數據，星上所能提供的在軌溫度遙測點非常少，要依據在軌溫度實測值進行熱變形的計算，只能通過插值等方法對整個載荷的溫度場數據進行估計。另外，由于溫度場在有限元模型中的賦值需要很大的工作量，導致工期較長，文獻[3]實現了有限元模型中溫度場的自動批量導入，避免了手工賦值耗時過長的問題。然而，熱變形分析過程中還有另一個工作量很大的環節，即分析需要有限元模型上所有節點的溫度值作為輸入條件，該溫度場由熱分析計算給出。在進行熱分析前，在軌工況的選擇需要對星內設備工作情況，外熱流情況等多種因素進行統計，如果要提供衛星在大量不同時刻的溫度分布情況，需耗費大量時間。文獻[4]應用插值方法對衛星紅外遙感器輻射定標光機系統的熱變形進行了分析，證明了該系統熱變形可滿足光學準確度要求，可為結構快速熱變形分析提供借鑒，但該系統體積規模較小，也決定了其溫度分布并不十分復雜，所以該文未對所用插值方法及估計結果的準確性開展討論。

地統計學中，對于大氣溫度測算等復雜問題通常借助一些已知測點來完成，即通過選取合適的無偏估計算法，完成對未知位置的溫度估計，十分方便快捷。普通克里格法[5]作為一種無偏估計算法，在地統計學中已得到成熟應用[6]，在空間溫度估計方面的應用也得到廣泛研究[7]。由于普通克里格法得到學者和工程師的廣泛肯定，對該算法的研究也不斷深入[8-9]，并且依據空間特性的不同發展出了簡單克里格法、協同克里格法等多種推廣型算法[10-11]。

本文首次探討了普通克里格法應用于衛星熱變形分析的可行性，提出一種基于在軌溫度遙測數據進行衛星熱變形快速分析的方法，用以跳過繁瑣的熱分析工況選擇過程。文章使用熱平衡試驗過程中衛星結構上的少量溫度實測數據，利用普通克里格法進行衛星溫度場估計，再進行溫度場賦值與變形計算，最后將分析結果與熱試驗過程中的實測變形結果進行比較，說明這種結構熱變形分析方法的正確性和可行性。

1 熱變形快速分析方法

傳統熱變形分析采用熱分析結果作為溫度場輸入，分析的工期往往需要數日乃至更長的時間，本文提出的方法與其最大的不同，在于本方法跳過熱分析步驟，直接使用統計學理論，基于少量溫度實測點完成結構的溫度場分析，用于估計衛星在軌任一時刻的變形情況，為衛星在軌熱變形分析提供一種快速有效的方法。

1.1 普通克里格法

普通克里格法依據溫度實測點即可完成溫度場估計，因此該方法可依據星上的溫度遙測數據，提高熱變形分析工作的效率，在航天任務中十分有利于在軌應用。該方法分為三步，如圖1所示：1)根據已知測點數據，對空間場進行結構分析，選取并提出變差函數模型；2)通過優化算法(本文采用拉格朗日乘數法)求解克里格方差的極小值，若克里格方差滿足要求，則進一步計算插值權值，若克里格方差不滿足要求，則回到步驟1，重新選取變差函數；3)得到計算插值權值后，進行克里格插值計算。

圖1 普通克里格法計算示意Fig.1 Calculation process of ordinary Kriging

普通克里格法一般公式為

待定點值的期望誤差應為0，即有

變差函數根據實驗變差函數擬合得到，實驗變差函數由采樣點及其坐標計算，算式如下：

(5)

得到實驗變差函數后，可用其擬合理論變差函數，擬合方法見文獻[6]。對變差函數的擬合，已經有成熟的專用軟件及Matlab工具箱，用于根據采樣點及其坐標擬合變差函數。

在式(2)的約束條件下，通過優化方法求解式(4)中克里格方差的極小值，也就得到了相應變差函數模型下的權重λi，即可進一步開展未知點的溫度計算。

進行二維估計時，對某個點進行數據估計可采用四分法，如圖2所示：一般對平面某點進行估計時至少存在4個采樣點，且應保證4個象限均有采樣點存在(邊界點除外)。進行三維數據估計時，宜采用八分法，即最好保證空間8個象限均有采樣點存在。

圖2 采樣點分布Fig.2 Sampling point distributions

對某個區域內進行溫度場估計，就是根據鄰近數據和變差函數進行數據平滑的過程。如果溫度場存在分布不均勻甚至叢聚現象，且采樣點數據較少，甚至在溫度突變位置無采樣數據，則會對估計結果造成較大偏差。因此在結構上存在熱源，或者在隔熱層的兩側存在溫度突變的區域都應布置采樣點。

1.2 應用普通克里格法的熱變形快速分析

由以上方法替代傳統熱變形計算中的熱分析步驟，完成溫度場估計后，即可進行熱變形計算，本文的熱變形計算在Nastran中完成。Nastran在進行有限元計算時，所使用的bdf文件是Patran根據有限元模型及工況條件生成的，這其中就包含了熱變形分析所需的溫度場數據。要在Nastran中進行熱變形分析，需要將在外部生成的溫度場數據導入bdf文件。

本文使用統計學方法推算有限元模型上所有節點的溫度值后，采用VC++編寫接口文件，將計算得到的溫度場寫入bdf文件，實現有限元模型溫度場的賦值，最后將熱變形計算結果與熱平衡試驗中的實測變形結果進行比較，進而說明方法的正確性。整個熱變形運算方法的流程如圖3所示。

圖3 快速熱變形計算方法流程Fig.3 Process of the rapid thermal deformation analysis

2 分析結果和試驗驗證

2.1 結構形式

本文所分析的星上載荷，其主體結構由下板、中板、上板以及三層板間的支撐筒組成。其中支撐筒內部采用主動溫控，上板采用散熱涂層進行被動控制。在有限元分析中，對上板、中板和下板的建模采用了二維網格單元，對支撐筒的建模采用了二維網格和三維HEX8單元。上板、中板和下板采用二維數據估計；支撐筒形狀較復雜，采用三維估計。

在諸多因素中，有限元模型對熱變形的影響占據了很大的比重。本文所用星上載荷的有限元模型為經過多次試驗修正的模型，包括正弦振動試驗修正與熱變形修正，這些修正對于最后分析精度的保證都起到了重要作用。

圖4 結構形式示意Fig.4 General view of the structure

2.2 試驗結果

在該載荷熱平衡試驗過程中，星上安裝了溫度遙測點，這些測點溫度即為普通克里格法進行溫度場估計的輸入條件；同時，為了解科學成像載荷在熱變形作用下的視軸指向變化，在星上安裝了傾角傳感器。溫度及傾角測點位置見圖5～圖6，CD1～CD26為溫度測點，測點溫度見表1。AT01～AT05為傾角傳感器。5個傾角傳感器(AT01～AT05)實測的Y、Z方向(水平方向)傾角變化結果見圖7，圖中第47～67 h為低溫工況時段。所測位置最終傾角可由Y、Z方向傾角進行計算：

式中：α、β分別為Y、Z方向傾角測量值；γ為總傾角變化。

圖5 上板溫度及傾角測點示意Fig.5 Temperature and dip angle sensors on the upper plate

圖6 中板溫度及傾角測點示意Fig.6 Temperature and dip angle sensors on the middle plate

2.3 分析結果及討論

根據已知點數據，需要建立變差函數模型來進行圖1中后續步驟的計算。對于溫度分布來說，距離測點較近的點，其溫度值與測點溫度值關聯性更大；反之，距離測點較遠的點，其溫度值與測點溫度值關聯性則較小，根據這種特點，本文中的變差函數初步確定采用高斯型：

式中：θi為比例參數，該參數依據已知測點數據，由最大似然估計確定。根據實測溫度數據擬合的變差函數見圖8。

圖8 根據采樣點擬合的變差函數Fig.8 Variogram fitted by the sampling points

使用變差函數完成插值后，采用交叉驗證法驗證結果，即假設部分測點溫度未知，用其他采樣點對這些點進行估計，再與實測值進行比較。變差函數選取的越符合實際，那么克里格估值結果也就越符合實際結果。分別對CD1～CD4、CD5～CD8進行估計：在估計CD1～CD4溫度時，采用CD5～CD17作為已知溫度點；在估計CD5～CD8溫度時，采用CD1～CD4、CD9～CD17作為已知溫度點，將CD1～CD4、CD5～CD8實測溫度和采用普通克里格法得到的溫度值進行比對，結果見表2，可見估計值與實測值偏差很小，均在0.6℃以內。說明變差函數選用高斯型能夠較準確地還原溫度場，可以進行下一步計算。

使用所有試驗實測數據進行估計，整個載荷上板的溫度估計結果見圖9。

圖9 上板溫度估計結果Fig.9 Estimate results of temperature on the upper plate

項目測點實測溫度/℃克里格插值估計溫度/℃第一組估計CD1-32 6-32 9CD2-33 5-33 7CD3-32 0-32 1CD4-34 0-34 4第二組估計CD5-33 3-33 5CD6-34 6-34 8CD7-32 4-32 97CD8-34 2-34 4

在溫度穩定后，圖4中載荷溫度變化較大的結構為上板和中板。而對于板間的同一支撐筒來說，在溫度穩定后，支撐筒上部與下部溫差僅在1.5℃左右，無論采用克里格插值還是更簡單的泰森多邊形插值，對分析結果的影響僅反映在一階小量。鑒于支撐筒的溫度分布特點，在同一支撐筒上最多只布置了兩個測點。以其中的一個支撐筒為例說明其三維溫度場建模步驟，如圖10所示：首先建立一個包絡該支撐筒的簡單有限元網格，將支撐筒頂部與底部測點溫度分別賦值于該有限元網格的頂部節點和底部節點；之后使用Patran中自帶的Field生成方法(這里采用泰森多邊形法[12])生產空間場；最后使用該Field將溫度場賦值給支撐筒。

載荷溫度場最終映射結果見圖11。

圖10 支撐筒溫度場賦值Fig.10 Temperature assignment of the supporting column

圖11 溫度場映射結果Fig.11 Mapping results of the temperature field

完成溫度場映射后，即可使用Nastran進行有限元分析，上板變形見圖12，可見變形結果為上板四周向上翹曲，距離中心基準越遠處變形值越大。

利用變形分析結果，對圖5～圖6中的傾角傳感器AT01～AT05所在位置進行熱變形計算，分析得到的傾角變化與試驗實測傾角變化的比較見表3。

圖12 載荷上板變形Fig.12 Thermal deformation of the upper plate of the payload

由表3中比較可見，AT01、AT03、AT05的分析值與實測值之差均在1′以內，AT02、AT04分析值與實測值之差在1′～2′之間，分析結果與試驗結果吻合良好。對于科學成像載荷熱變形的計算已經達到了較高的精度，可滿足該載荷的要求，因此文中所用分析過程可應用于該載荷視軸在熱變形作用下的指向變化估計。

表3 傾角變化的試驗值與分析值的比較

同時，本方法使用在軌溫度遙測數據估計結構溫度場，因此分析結果也對應著溫度遙測數據采集時刻的變形情況。所以，需要知道在軌某一時刻的結構變形情況時，只需提供其附近時刻的溫度遙測數據，就可對該時刻的變形情況進行分析。

3 結束語

文章首次應用一種統計學方法——普通克里格法，提出了衛星在軌熱變形計算的快速方法。該方法省略了繁雜的熱分析過程，將溫度場估計工期由數天甚至數周減少為幾小時，可用于衛星的實時在軌變形估計。基于上述方法，本文對衛星載荷熱變形進行了分析，結果表明，分析值與實測值差別在1′左右，滿足該載荷的指向精度要求。

基于本文的工作，可在后續研究中進一步探索如何應用變差函數的結構套合、層次套合及各向異性套合等方法，針對結構溫度場分布更復雜的情況進行估計。達到快速性計算并預估復雜溫度場和熱變形的目的，最終以能實施溫度或變形的主動控制為目標開展下一步工作。

References)

[1] 游思梁，陳桂林，王淦泉. 星載輻射計掃描鏡太陽輻射熱-結構建模與仿真[J].中國空間科學技術，2011，31(1)：62-69.

YOU S L，CHEN G L，WANG G Q. Thermal-structural modeling and simulation for scan mirror on imager in solar radiation[J].Chinese Space Science and Technology，2011，31(1)：62-69(in Chinese).

[2] 張鏡洋，常海萍，王立國. 小衛星瞬態熱分析模型修正方法[J]. 中國空間科學技術，2013，33(4)：24-30.

ZHANG J Y，CHANG H P，WANG L G. Correction method for transient thermal analysis model of small satellite[J]. Chinese Space Science and Technology，2013，33(4)：24-30(in Chinese).

[3] 張思斯. 掃描輻射計及冷凝管道的熱分析[D]. 上海：上海交通大學，2008：11-21.

ZHANG S S. The thermal analysis of radiometer and HDPE pipe[D]. Shanghai：Shanghai Jiaotong University，2008:11-21(in Chinese).

[4] 肖慶生，楊林華，趙壽根.衛星紅外遙感器輻射定標光機系統熱-結構耦合變形分析[J]. 航天器環境工程，2011，28(1)：46-51.

XIAO Q S，YANG L H，ZHAO S G. Thermal-structure coupled deformation in an optical-mechanical system for radiometric calibration of satellite IR remote sensor[J]. Spacecraft Environment Engineering，2011，28(1)：46-51(in Chinese).

[5] NOEL C. Spatial prediction and ordinary Kriging[J]. Mathematical Geology，1988，20(4)：405-421.

[6] 王家華，余元華，高海余. 克里金地質繪圖技術[M]. 北京：石油工業出版社，1999： 39-110.

WANG J H，YU Y H，GAO H Y. Beijing:Petroleum Industry Press[M]. Beijing：Petroleum Industry Press，1999：39-110(in Chinese).

[7] 彭思嶺. 氣息要素時空插值方法研究[D]. 長沙：中南大學，2010：5-16.

PENG S L. Developments of spatio-temporal interpolation methods for meteorological elements[D]. Changsha：Central South University，2010：5-16(in Chinese).

[8] MARCOTTE D. Cokriging with MATLAB[J]. Computers & Geosciences，1991，17(9)：1265-1280.

[9] MARTIN H T. Matlab Recipes for Earth Sciences[M]. German：Springer，2007：28-39.

[10] SIMPSON T W，BOOKER A J，GHOSH D，et al.Approximation methods in multidisciplinary analysis and optimization:a panel discussion[J]. Struct. Multidisc Optim.，2004，27(5)：302-313.

[11] DANIEL B，CHRIS L F，ARMIN I. Hierarchical nonlinear approximation for experimental design and statistical data fitting[J]. SIAM Journal on Scientific Computing，2007，29(1)：49-69.

[12] 李少華，劉遠剛，王延忠. 泰森多邊形在地質數據去叢聚中的應用[J]. 物探與化探，2011，35(4):562-564.

LI S H，LIU Y G，WANG Y Z. The application of Thiessen polygen to the declustering effect of geological data[J]. Geophysical & Geochemical Exploration，2011，35(4)：562-564(in Chinese).

(編輯：車曉玲)

Rapid thermal deformation analysis of on-orbit satellites on the basis of statistic method

ZHANG Yechi*，ZHANG Long，LUO Wenbo，PAN Teng

Beijing Institute of Spacecraft System Engineering，Beijing 100094，China

To calibrate the orientation position of the satellite orbit by the temperature sensors，the ordinary Kriging method widely used in geo-statistics was applied into the temperature estimation process of the satellite. A rapid thermal deformation analysis method was developed based on the statistic methods. Based on a small amount of temperature sensors，estimation of the temperature field was carried out. Then，the temperature field was used for the finite element analysis to accomplish the thermal deformation analysis. Differences between the analytical and the experimental results of the deformation of the payload are around 1′，which confirms that the method can be used to obtain the variation of optical axis of the payload. Furthermore，compared to the thermo physical method，the work period can be reduced by several days. This work provides guidance for rapid thermal deformation analysis of on-orbit satellites.

temperature measuring points；rapid thermal deformation analysis；statistic method; ordinary Kriging；variogram；finite element；satellite

10.16708/j.cnki.1000-758X.2016.0062

2016-03-29；

2016-06-17；錄用日期：2016-08-22；

時間：2016-12-16 11:29:12

http:∥www.cnki.net/kcms/detail/11.1859.V.20161216.1129.007.html

戰略性先導科技專項

張也弛，張龍，羅文波，等. 基于統計學方法的衛星在軌熱變形快速分析[J].中國空間科學技術，2016，36(6)：

55-61.ZHANGYC，ZHANGL，LUOWB，etal.Rapidthermaldeformationanalysisofon-orbitsatellitesbasedontemperaturemeasuringpoints[J].ChineseSpaceScienceandTechnology，2016，36(6)：55-61(inChinese).

V414.3

A

http:∥zgkj.cast.cn

*通訊作者：張也弛(1983-)，男，博士，高級工程師，zhangyechi83@163.com，研究方向為航天器總體設計