面向InSAR的空氣擾動影響機翼撓曲變形建模

朱莊生，張萌

（1.北京航空航天大學 儀器科學與光電工程學院，北京100083； 2.慣性技術重點實驗室，北京100083；3.新型慣性儀表與導航系統技術國防重點學科實驗室，北京100083）

干涉合成孔徑雷達（Interferometry Synthetic Aperture Radar，InSAR）是在合成孔徑雷達（Synthetic Aperture Radar，SAR）二維成像基礎上發展并逐漸成熟起來的三維遙感成像技術。InSAR技術中的相位干涉基于相位連續性假設，即要求任意相鄰像素之間的絕對相位差小于π［1］，因此雙節點（即單基線）InSAR技術適用于地形平坦區域。為解決實際測繪中高度突變區域的相位干涉問題，多節點（即多基線）InSAR技術應運而生［2］，垂直航向上基于多載荷節點構成的長短基線聯合測量，可有效克服高度突變、噪聲干擾等影響，現已成為航空遙感領域的主流發展方向。

受外部擾動的影響，InSAR載機成非理想運動，機翼也會產生變形和振動，進而導致InSAR系統各節點載荷處形成復雜運動，即引起干涉基線長度的變化。從相位干涉技術可知，干涉基線長度測量誤差直接決定了InSAR的高程分辨率［3］，目前主要借助傳遞對準技術并輔以機翼撓曲變形誤差建模補償技術［4-5］提高基線測量精度，機翼撓曲變形誤差建模是其核心問題之一。

機翼撓曲變形誤差建模方法主要有3類：

第1類是基于隨機理論建模。該方法將機翼撓曲變形當作Markov模型［6］或AR模型［7］，并將其增廣為卡爾曼濾波狀態量［8］或噪聲量［9］進行估計。如Kain和Cloutier［8］設計的42維卡爾曼濾波器中有18個撓曲運動狀態量，公式復雜，計算效率低。該方法局限性在于模型參數設定無實際依據，也未考慮外界干擾變化對機翼變形的影響分析。

第2類是直接測量法。如攝影測量法，其要求攝影機對測點通視或直視，主要用于航天領域，如美國的SRTM［10］（Shuttle Radar Topography Mission）。在航空領域，由于機翼變形過大或天氣影響，攝影機對測點難以保證直視或通視。如基于光纖光柵（Fiber Bragg Grating，FBG）傳感器的機翼變形測量，借助FBG感知應變信息，計算機翼變形，多應用在機翼健康監測［11］方面。

第3類是基于有限元仿真的變形建模。有限元仿真是基于計算流體力學（Computational Fluid Dynamics，CFD）及計算結構力學（Computational Structural Mechanics，CSM）的流固耦合理論發展起來的，其湍流模型成熟、可控性強，主要應用于計算飛機氣彈性力和結構響應［12］及安全設計［13］。如陳志敏等［14］利用有限元仿真獲得機翼受氣動載荷作用的彎曲和扭轉變形情況；張華等［15］利用Nastran軟件分析大展弦比復合材料前掠機翼和后掠機翼在氣動載荷作用下的靜氣彈變形情況。

機翼撓曲變形主要源于內部擾動（如發動機振動、機翼結構及材料等）和外部擾動（如空氣擾動等）的影響，內部擾動隨著飛行平臺的確定可認為是一些已知量，而外部擾動是一個復雜的過程。本文利用有限元軟件針對空氣擾動對機翼結構的氣動力影響與結構響應進行仿真求解，建立起空氣擾動影響下的機翼撓曲變形模型，并借助應變估計位移的模態疊加原理對模型進行驗證計算。

1 空氣擾動類型及影響因素分析

1.1 空氣擾動的復雜多源性

對于機載多基線InSAR系統，具體結構如圖1所示，DLR研制的E-SAR［16］工作高度為3～5 km，NASA研制的AIRSAR［17］和TOPSAR［18］分別工作在8 km和9 km高度；中國科學院電子學研究所研制的機載三基線InSAR［19］工作在3 km航高上，該所研制的MiNiSAR［20］工作在0.5～3 km高度。綜上，本文考慮機載InSAR系統成像高度范圍為0.5～9 km，在此工作區域內的空氣擾動主要為垂直風切變、陣風及大氣湍流［21］。

圖1 機載分布式InSAR系統Fig.1 Airborne distributed InSAR system

1.1.1 垂直風切變

垂直風切變指的是垂直氣層任意兩觀測點之間風向和風速的突然變化，計算公式為

式中：S2、S1分別為空間A、B兩點的風速值。

根據GJB 5601—2006［22］給出的垂直風強度，計算H＝30、50、70、100m工作高度突變所引起的風切變大小，如圖2所示。

圖2 不同高度變化引起的垂直風切變Fig.2 Vertical wind shear caused by different height changes

由圖2可知，在10 km的工作高度內，100 m的高度變化最大可引起0.35m／s的垂直風切變；航空中經常取30m的高度變化判斷其對飛行器的影響［23］，本文計算可得垂直風切變大小為0.1m／s；根據國際民航組織規定垂直風切變小于1m／s的情況下與風切變強度小于0.033／s的風切變等級為微弱影響［23］。除此之外，由于垂直風切變影響載機升力最終改變載機工作高度，需控制載機產生相反的工作高度變化，如垂直風切變引起工作高度下降，則需提高當前工作速度以產生更大升力來提高載機的工作高度；而起飛和降落階段的工作速度小，機動余量小，改變工作高度又不夠，通常認為垂直風切變主要影響起飛和降落階段［24］。

1.2.1 仿真條件設置

建立如圖4所示的機翼蒙皮骨架結構。蒙皮厚度為2mm，骨架結構由翼梁、翼肋和桁條組成，圖4中從左端0.1m處起每隔0.2m分別為測點1～測點15；采用BOEING103翼型，翼弦為0.6m，總長為3m，材料為鋁合金7075，材料屬性如表1所示。

基于有限元軟件ANSYSWorkbench分析載機工作速度、高度引起的大氣湍流變化對機翼變形量的影響，主要利用ICEM CFD軟件、Fluent模塊及Static Structural模塊進行流固耦合求解。應用ICEM CFD軟件建立如圖5所示的機翼飛行流場網格區域，分為壁面區域、湍流入口區域、湍流出口區域及機翼結構。在Fluent模塊中計算機翼不同工作狀態下所受的氣動載荷，湍流模型采用單方程Spalart-All maras模型；氣體設置為constant，不同工作高度上的氣體密度根據表2［22］確定；邊界條件采用速度入口、壓力出口條件，壓力、溫度根據表2［22］確定。在Static Structural模塊中將Fluent計算載荷施加在機翼上，計算機翼結構應變及變形。

圖4 機翼模型Fig.4 W ing model

表1 鋁合金7075材料屬性Table 1 A luminum alloy 7075 m aterial p roperties

圖5 機翼飛行流場區域Fig.5 W ing flight flow field area

以載機工作高度3 km、速度200 m／s為例，Fluent計算蒙皮氣動載荷及湍流速度變化，如圖6所示，Static Structural計算骨架結構Y向（橫向）位移，如圖7所示，圖7中左端至右端的變化表示機翼變形量的增大。

表2 不同工作高度上的大氣物理參數［22］Table 2 Atmospheric physical parameters at different working altitudes［22］

圖6 機翼所受氣動載荷及湍流速度變化Fig.6 Changes in aerodynamic loads and turbulent velocity on wing

圖7 機翼結構Y向位移Fig.7 Y-direction displacement of wing structure

1.2.2 大氣湍流影響機翼撓曲變形仿真分析

載機工作高度一定時，即大氣密度、壓力、溫度均未變化時，僅僅是工作速度變化影響大氣湍流，仿真3 km航高上150 m／s起每隔5 m／s直至250m／s的工作速度下大氣湍流對機翼撓曲變形的影響，結果如圖8所示。

為分析航空遙感系統工作高度對機翼變形量的影響，采用載機從2.5 km 起每隔0.5 km 到6 km的工作高度變化，以及從125 m／s起每隔25m／s到300m／s的工作速度變化來獲取仿真數據，具體的湍流參數按照表2設置。在此僅列出測點15的位移，如圖9所示。

圖8 不同工作速度時的機翼撓曲變形Fig.8 Deflection of wing at different working speeds

圖9 不同工況下測點15的機翼撓曲變形Fig.9 Wing deflection deformation atmeasuring point 15 under different conditions

可以看出，隨著工作高度遞增（壓力減小），位移量逐漸減小且近似一階線性關系；綜合圖8可看出，位移量隨工作速度變化近似二階關系。下文將借助機翼受力變形模型并結合實驗數據來分析大氣湍流對機翼撓曲變形的影響。

2 模型建立及參數辨識

2.1 大氣湍流影響機翼撓曲變形模型建立

InSAR載機的工作速度、高度變化將引起大氣湍流變化，由于湍流在飛行中測量難度較大，所以本文將InSAR載機的工作速度、高度作為輸入量，大氣湍流作為未知中間量，機翼變形量作為輸出量，建立起載機不同工作速度、高度影響下的機翼撓曲變形模型。

變形建模主要有機理建模法和測試分析法。機理建模法根據對象特性分析其因果關系尋找內部機理，所建模型具有物理數學意義。如詹斌［31］通過分析石墨烯納米金屬基復合材料的變形機理建立起其塑性變形模型。該方法一般基于大量簡化及假設前提，難度較大，不具有普適性。測試分析法又稱參數辨識法［32］，是將研究對象特性視為“黑箱”，只能利用大量表征系統規律和狀態的實驗數據，運用統計方法建立與原始數據擬合度最好的模型，所建模型在線校正能力強；該方法無需了解系統機理，但需設計合理實驗以獲得系統最大信息量。所以，通常采用機理分析結合參數辨識法［33］建立模型。如關永亮［34］基于復合材料變形機理綜合參數辨識法建立起復合材料機翼的非線性動力學變形模型。

對于本文建立機翼撓曲變形模型而言，首先需要明確大氣湍流與其引起的機翼所受載荷之間的關系，其次需要明確機翼載荷與機翼撓曲變形之間的關系，從而初步建立起機翼撓曲變形機理模型，最后通過獲取的仿真實驗數據對模型進行參數辨識并建立模型。

圖10 蒙皮表面氣動載荷分布Fig.10 Aerodynamic load distribution on skin surface

對于圖4所示的機翼結構，蒙皮感受如圖10所示的沿翼展分布載荷qS、翼弦分布載荷和qC，為飛機提供升力并將外載荷傳遞給骨架，并不是主要的承力結構。對于骨架結構而言，桁條對蒙皮起支撐作用，并承受由蒙皮傳遞的氣動載荷；翼肋保持機翼的氣動外形并承受一定剪力；翼梁在根部與機身相連，主要承受機翼的彎矩和剪力［35］，并且多節點InSAR成像載荷主要分布在翼梁上。綜上，考慮骨架整體結構承受機翼的彎矩和剪力，即InSAR基線長度誤差將直接取決于機翼骨架結構的撓曲變形建模誤差［3］。

實際情況下，翼梁結構與機身相連，形成一個懸臂梁的狀態，本文對于骨架整體結構也是基于懸臂梁假設考慮的，受力變形情況如圖11所示。

2009年，美籍華裔科學家高錕獲得諾貝爾物理學獎，然而在得到自己摘取“科學王冠上的明珠”的消息時，高錕的阿爾茨海默病病情已經非常嚴重：他不能理解自己為何被稱為“光纖之父”，也不知道諾貝爾獎是一個科學家所能夠獲得的最高榮譽。

圖11 機翼骨架結構變形示意圖Fig.11 Schematic diagram of wing skeleton structure deformation

文獻［36］給出了懸臂梁在載荷q作用下與其產生的橫向彎曲位移（撓度）W之間的關系：

式中：pw（z）＝［（l－z）4＋4l3z－l4］／（24EI），l為懸臂梁長度；z為距固定端距離；E為彈性模量；I為截面慣性矩。

根據梁的彎曲理論，距中性層距離d處的軸向位移U與該截面處撓度W的關系為：U（z，h）＝將式（4）代入計算可得設軸向位移為

已知機翼蒙皮感受載荷qL、qD，以及載荷作用下機翼變形量W、U，并根據式（1）可確定載機工作速度影響湍流速度。綜上，初步建立起載機工作速度引起大氣湍流的變化所影響的機翼撓曲變形模型為

式中：AY、AZ項隨機翼位置變化，與機翼材料、尺寸及氣動性能有關；vB體現骨架感受氣動載荷的大小，B為一常值，不隨機翼位置變化，表示載機工作速度、湍流速度對載荷的影響。由機翼受力變形機理建立模型式（6）、式（7），借助實驗數據辨識模型參數可建立起載機工作速度影響機翼撓曲變形模型；而工作高度影響湍流尺度、強度機理復雜，只能借助參數辨識法建立。

2.2 基于回歸分析的參數辨識

基于參數辨識法建立機翼撓曲變形模型的方法主要有時序分析法和回歸分析法。時序分析法是對實驗數據基于ARMA模型、AR模型、MA模型等研究分析、識別參數。如解春明等［7］將機翼撓曲變形當作AR模型，并利用卡爾曼濾波對其進行在線補償。該方法的局限性在于驗證數據的平穩有序性，重點是模型階次n的選擇，即考慮前n時刻對當前時刻的影響。對本文建立載機工作高度與速度影響機翼撓曲變形的模型而言，并不適用。

本文建立載機工作速度與高度自變量與機翼撓曲變形因變量之間的數學關系，與利用回歸分析辨識參數的理論思想相契合，所以基于回歸分析法借助有限元仿真實驗數據對機翼撓曲變形機理模型辨識參數。

回歸分析中，常用線性回歸模型：LR（xm）＝Am1xm＋Am2；非線性回歸模型：指數函數等。利用回歸分析可對同一組數據擬合出多種不同模型，常用相關指數R2及殘差平方和SSE選擇最優擬合模型，SSE越小及R2越接近1，說明實際數據與預測值越相近，擬合程度就越好。

根據上述分析，對于載機工作速度影響下的機翼撓曲變形模型，直接利用非線性模型中的辨識參數建立模型；對于工作高度變化所引起的大氣湍流影響下的機翼撓曲變形模型，利用上述回歸模型進行擬合處理，選擇與原始實驗數據擬合程度最好的模型進行下一步的計算及預測。

2.3 大氣湍流影響載機撓曲變形模型建立

對于3 km工作高度上載機工作速度變化的實驗數據，借助冪函數u＝AvB辨識參數建立起載機工作速度影響機翼撓曲變形模型，在此列出機翼骨架結構上測點2、測點7、測點12、測點15的Y向位移數據依據u＝AvB的參數辨識結果，如表3所示。

對比表3中3個不同擬合模型，根據SSE越小和R2越接近1的判斷準則選擇最優擬合模型，確定InSAR載機在3 km工作高度上工作速度從150m／s變化至250m／s的情況下，Y向、Z向機翼變形量滿足：

式中：AY、AZ項與機翼測點位置、機翼結構的材料、尺寸及氣動性能有關；v2.015項體現載機工作速度和湍流速度與機翼承受載荷間的關系。

對于建立工作高度影響機翼撓曲變形模型而言，高度僅僅體現壓力、溫度和密度的變化，所以該模型轉換為容易測量獲取的不同工作高度上所對應的壓力值與機翼變形量之間的關系。對同一工作高度2.5～6 km上載機工作速度變化引起的機翼變形數據進行處理，發現均滿足式（8）和式（9），僅系數AY、AZ發生變化，變化趨勢如圖12所示?？梢钥闯?，AY、AZ隨高度變化趨勢保持一致且近似一階關系，與距固定端距離z關系近似三階關系。

選擇測點2、測點7、測點12、測點15的Y向位移進行參數辨識，統計結果如表4所示。

表3 載機工作速度影響的參數辨識結果Table 3 Parameter iden tification results of carrier speed

圖12 高度變化時的AY、AZ 系數Fig.12 AY and AZ at different heights

表4 載機工作高度影響的參數辨識結果Table 4 Parameter identification results of car rier’s work ing height

對比表4中模型，根據SSE越小及R2越接近1的判斷準則，確定工作高度變化時系數AY、AZ滿足：

式中：aY、aZ隨機翼結構及測點位置z變化；P表示不同工作高度上的壓力。

綜上，載機工作高度以及工作速度影響機翼撓曲變形模型可表示為

式中：P、v項為與機翼載荷大小相關量。

3 仿真實驗

3.1 基于建立模型的機翼撓曲變形誤差仿真分析

利用式（14）、式（15），對于不工作高度、不同工作速度下變形數據進行預測計算，表5和表6列出了測點7、測點12、測點15的變形數據。

表5 本文模型預測Y向變形Table 5 Y-direction deform ation pred icted by p roposed model

表6 本文模型預測Z向變形Table 6 Z-direction deform ation pred icted by p roposed model

3.2 基于模態疊加原理的機翼撓曲變形誤差仿真分析

對于本文所建模型，適用于已知InSAR載機工作速度、高度的機翼撓曲變形量的計算，本文借助模態疊加原理（又稱模態法）對其適用性及精度進行驗證。模態疊加原理［38］基于載荷作用下結構應變響應估計位移，存在：

式中：D為待估計位移值；S為測點應變值；DST為應變、位移轉換矩陣，軸向位移矩陣DST u＝［фu］（［ψu］T［ψu］）－1［ψu］T，橫向位移矩陣DSTω＝［фω］（［ψu］T［ψu］）－1［ψu］T，［фu］、［ψu］和［фω］分別為軸向應變、位移和橫向位移矩陣。

圖13 骨架結構的應變、位移模態Fig.13 Skeleton structure strain and displacementmode

利用ANSYSWorkbench模態分析可獲得骨架結構的應變、位移模態，如圖13所示。基于仿真分析最終計算出DST矩陣。將仿真獲取的應變值利用式（16）便可估計Y向、Z向位移。表7和表8列出了不同飛行條件下測點7、測點12、測點15利用模態疊加原理的計算結果。

表7 Y向變形模態疊加原理計算結果Table 7 Y-d irection deform ation calcu lated by modal superposition principle

表8 Z向變形模態疊加原理計算結果Table 8 Z-d irection deform ation calcu lated by modal superposition p rinciple

3.3 仿真結果及分析

從表5～表8可知，利用本文模型及模態疊加原理計算Y向、Z向位移與ANAYS仿真有一定誤差。表9列出了2種方法的計算誤差。

從表9可以看出，利用模態疊加原理及參數辨識建立模型計算橫向最大誤差均在0.6mm以內，軸向最大誤差均在0.015mm以內，即利用模態疊加原理與回歸預測2種方法均能很好地估計出Y向、Z向位移量，既驗證了本文方法的可行性，又對比了2種方法的計算精度。

表9 本文模型與模態疊加原理計算誤差值對比Table 9 Com parison of error values calculated by p roposed model and modal superposition p rinciple

4 實測實驗

4.1 實驗系統及步驟

為驗證本文方法，搭建起如圖14所示的機翼結構分布式光纖光柵測量系統，將InSAR載機在不同大氣湍流情況下所承受氣動載荷轉換為對機翼結構施加不同加載量。以0.5 kg每隔0.5～5 kg的加載量代表3 km工作高度上150m／s起每隔10m／s直至240m／s工作速度的機翼所受載荷；仿真實驗中不同湍流擾動下的機翼變形量由ANSYS直接計算獲取，本節實驗利用光纖光柵傳感器實時獲取其變形量；由此建立起湍流擾動模擬加載量與機翼變形量之間的映射關系，從而進行預測計算。

圖14 實驗室機翼結構Fig.14 W ing structure in laboratory

在機翼結構上進行加載實驗，如圖15所示，在實驗過程中從0.5 kg開始，每次以0.5 kg的重量遞增直到5 kg，每次加載保持靜止狀態5min以上，利用光纖光柵測量加載過程中的應變變化及雙全站儀測量系統測量2個靶標測點的位移量。

圖15 機翼加載實驗Fig.15 Loading experiment in wing

4.2 實驗結果及分析

對懸臂梁進行10次加載，全站儀測量2個靶標點的位移量如表10所示。

利用光纖光柵測量應變量并借助模態疊加原理計算可得加載過程當中的位移變化，并利用前7次加載基于模態值進行非線性回歸分析建立模型，擬合值如圖16所示。預測計算在機翼上施加39.2、44.1、49 N的位移值，與全站儀測量值進行對比計算，結果如表11所示。

表10 全站儀測量值Table 10 Total station measurement values

從表11可以看出，利用模態值與全站儀測量值對比，Y向、Z向誤差均在0.35 mm以內，既驗證了本文方法的可行性，又驗證了光纖光柵測量機翼撓曲變形的精度。

圖16 模態值與擬合模型值對比Fig.16 Comparison of modal and fitted values

表11 預測值與全站儀測量值對比Table 11 Com parison of p redicted and total station measurement values

5 結 論

本文通過仿真分析了大氣湍流對多節點In-SAR系統的載機機翼撓曲變形的影響，提出了一種基于機理建模綜合參數辨識方法建立湍流影響機翼撓曲變形模型的方法，利用基于模態疊加原理的結構變形估計方法對其進行驗證計算，并在實驗室搭建起機翼撓曲變形測量系統對本文方法進行實驗驗證。

1）對不同工作高度、速度下的大氣湍流仿真實驗，構建了影響大氣湍流因素與機翼結構撓曲變形之間的數學關系，從而實現將湍流對機翼撓曲變形的影響量化并建模分析以及預測。

2）與傳統的大氣湍流對于機翼撓曲變形影響的隨機理論建模分析相比，本文方法借助大量仿真實驗數據作為計算基礎，提高了所建模型的可靠性。

本文的主要目的是基于工作高度、速度量測量來構建大氣湍流與機翼撓曲變形量間的映射關系，其研究思路和研究方法對所有機翼結構具有普適性和推廣性。但是，針對不同機翼結構，其升力、阻力等相關系數不同，以及在同等外力條件下，機翼變形量也不同；所以，針對不同的機翼結構本文所建模型系數應有所調整，這也是后續的重點研究工作。