星載微波輻射計的線性與非線性反演算法比較

王兆徽 劉宇昕，2 宋清濤,2 吳彬鋒

(1 國家衛星海洋應用中心, 北京 100081) (2 國家海洋局空間海洋遙感與應用研究重點實驗室, 北京 100081)

?

星載微波輻射計的線性與非線性反演算法比較

王兆徽1劉宇昕1，2宋清濤1,2吳彬鋒1

(1 國家衛星海洋應用中心, 北京 100081) (2 國家海洋局空間海洋遙感與應用研究重點實驗室, 北京 100081)

針對星載微波輻射計反演海洋大氣參數的過程，分析了反演通道的選擇，比較了線性與非線性反演算法的區別。非線性反演算法逐步舍去低分辨率的低頻通道，能夠直接有效地提高反演產品的空間分辨率。對幾種反演算法的分析和實驗表明，非線性反演算法中的基于Nelder-Mead的最優化反演算法相對于線性反演算法，其精度較高，具有一定的開發價值和應用潛力。

星載微波輻射計； 反演算法； Nelder-Mead算法

1 引言

星載微波輻射計是測量地表或大氣通過熱輻射輻射出的電磁波的儀器，其應用始于1962年美國發射的水手2號(Mariner-2)衛星進行金星飛行探測。從1972年開始，微波輻射計應用于海洋方面的研究。美國1972年發射的雨云5號(Nimbus-5)衛星和1975年發射的雨云6號(Nimbus-6)衛星搭載了單通道電磁掃描微波輻射計(Electrically Scanning Microwave Radiometer，ESMR)，驗證了被動微波遙感反演海洋大氣參數的潛力。1978年發射的雨云7號(Nimbus-7)衛星搭載了采用扇形機械掃描方式的掃描多通道微波輻射計(Scanning Multichannel Microwave Radiometer，SMMR)。1987年之后國防氣象衛星計劃(Defense Meteorological Satellite Program，DMSP)搭載了專用遙感器微波成像儀(Special Sensor Microwave/Image，SSM/I)。SSM/I改進了SMMR的許多儀器硬件問題，并為后續的“熱帶降雨測量任務”(Tropical Rainfall Measuring Mission，TRMM)衛星搭載的熱帶微波成像儀(Tropical Microwave Imager,TMI)和先進微波掃描輻射計(Advanced Microwave Scanning Radiometer,AMSR)的設計提供了基礎。AMSR[1]是日本宇宙航空研究開發機構(JAXA)和美國航空航天局(NASA)聯合設計發射的星載微波掃描輻射計。包括搭載于地球觀測-水(EOS-Aqua)衛星的AMSR-E，搭載于日本先進地球觀測衛星-2(ADEOS-II)的AMSR，以及搭載于全球環境變化觀測-W1(GCOM-W1)衛星的AMSR-2。

中國的“風云”系列衛星和“海洋”系列衛星均搭載了微波輻射計。2000年6月25日發射的風云二號(FY-2B)是靜止軌道氣象衛星，搭載了3通道的微波輻射計。風云二號的后續衛星和風云三號都[2]增設了微波輻射計的通道數目。中國于2011年8月發射了首顆海洋動力環境衛星——海洋二號(HY-2A)衛星，首次實現集主被動傳感器于一身，其中掃描微波輻射計是該星的4個微波載荷之一，采用圓錐掃描方式[3]，主要用于測量全球海洋表面溫度、海面風速、海洋上空水汽和降雨等參數。

微波輻射計反演算法的研究為微波輻射計的應用提供了技術支持。使用微波輻射計的觀測數據和可靠的反演算法能夠獲得海面溫度、海面風速、大氣中的水汽含量，云液態水含量等特定信息，這些信息對全球水循環以及氣候變化等問題的研究提供了幫助。除此之外，對反演算法的深入研究還能夠為儀器的設計(例如觀測頻率、觀測天線位置等)提供合理的需求分析。在軌運行的微波輻射計系統中，AMSR系列由于其連續穩定的工作狀態和相對廣泛的應用基礎，是目前研究最深入的星載微波輻射計。因此，基于AMSR-E開展的反演算法研究對于我國海洋動力環境衛星的發展具有一定的參考價值。

本文參考AMSR-E的儀器參數和說明文件(Algorithm Theoretical Basis Document，ATBD)，對比了幾種線性反演算法和非線性反演算法。實驗結果表明，星載微波輻射計的非線性反演算法雖然計算效率低于線性反演算法，但其物理意義明確且準確度更高，具有應用的價值。

2 反演算法概述

目前國際最常用的星載微波輻射計反演算法由“遙感系統”(Remote Sensing System，RSS)開發完成，Wentz參考Goodberlet，Swift和Wilkerson等人為SSM/I開發的GSW算法[4]，為AMSR開發了多元線性回歸[5](Multiple Linear Regression，MLR)算法(該算法結構簡單、計算快速)。RSS在2007年公布了與之類似的新的回歸模型[6]，并增加了插值技術的使用，由于新模型中的二次項可以視為單獨的一次項進行回歸分析，因此這種新模型可以視為升級版多元線性回歸模型。除了線性模型的反演方法，國內外的學者還研究了星載微波輻射計的非線性反演算法[7-8]。中國國家衛星海洋應用中心(National Satellite Ocean Application Services，NSOAS)的空間海洋遙感與應用研究重點實驗室(Key Laboratory of Space Ocean Remote Sensing and Application，LORA)開發了基于Nelder-Mead的最優化反演算法[9]。

反演算法的理論基礎是輻射傳輸模型(Radiative Transfer Model，RTM)。輻射傳輸模型建立了海面溫度(Ts，單位K)，海面風速(W，單位m/s)，大氣垂直積分水汽含量(V，單位mm，以下簡稱水汽含量)，大氣垂直積分液水含量(L，單位mm，以下簡稱液水含量)等海洋大氣參數與遙感衛星接收到的輻射亮溫之間的關系，是反演計算的理論依據。基本關系可以表達為

TB=TBU+τ[ETs+TBΩ]

(1)

式中：TB為星載傳感器接收到的總的輻射亮溫，TBU為上行有效溫度，τ為大氣透過率，E為海面發射率，TBΩ為大氣輻射的海面散射。此外，RTM還包括了觀測頻率、入射角等可以作為常數計算的儀器參數。輻射傳輸模型應用于不同微波輻射計的反演工作還需要考慮儀器參數的差別。兩種典型的星載微波輻射計的主要參數如表1所示。

表1 在軌運行的微波輻射計相關參數的比較

3 線性模型反演海洋大氣參數

3.1 多元線性回歸反演海洋大氣參數

多元線性回歸算法可以表述為如下形式：

(2)

式中:Pj為反演所求的不同的海洋大氣參數，R、I為線性方程。下標i為輻射計的通道(1=6.9V，2=6.9H，…)，下標j為反演的海洋大氣參數(1=TS，2=W，3=V，4=L)，c為反演系數。反演系數可以采用模擬計算獲得，也可以采用實際觀測的微波輻射計亮溫和現場觀測的海洋大氣參數配對所得的數據獲得。

由于兩個低頻波段(6.9 GHz，10.7 GHz)與三個高頻波段(18.7 GHz,23.8 GHz，36.5 GHz)對空氣的吸收不同，因此回歸方程采用如下形式：

I(TB)=TB,v=6.9,10.7 GHz

I(TB)=-ln(290-TB),

v=18.7,23.8,36.5 GHz

(3)

R(X)=X

(4)

在反演計算之前，需要對觀測的輻射亮溫進行關于入射角、海面風方向和海面鹽度的修正。入射角和海面鹽度可以簡單地使用線性模型進行修正，而海面風方向在多元線性回歸模型中只能采用180°的模糊修正。

3.2 升級版多元線性回歸反演海洋大氣參數

升級版多元線性回歸算法的第一步是獲得反演初步結果。反演系數的獲取與多元線性回歸算法類似：

(5)

ti=TBi-150,v=6.9,10.7,18.7,36.5 GHz

ti=-ln(290-TBi),v=23.8 GHz

(6)

第二步，RSS以海面溫度和海面風速為例進行雙線性插值：

(7)

理論上對于4個待求的反演結果，可以推廣為希爾伯特空間的“四線性插值”，但這樣的方法具有極大的計算規模，并不一定具有實際應用價值。式(7)的下標k表示3～34℃，下標l表示0～37m/s，以上插值范圍按照實際需求做等分處理，一個方便的方法是將插值節點間隔取為1。據此，式(5)反演的初步結果將位于由等分間隔的海面溫度和海面風速構造的插值格網。最終的反演結果可以寫為這樣的插值形式：

(8)

式中：權重系數w的計算與式(5)反演的初步結果距插值節點的距離相關。

考慮到海面鹽度、入射角和海面風向的影響，最終反演的回歸系數需要進行修正。鹽度校正可以簡單地使用線性模型進行修正，與多元線性回歸模型相同。入射角和海面風向的校正采用雙線性插值方法直接修正反演系數。入射角插值節點選為55°±1°，間隔0.5°；海面風向插值節點選為0°～180°，間隔20°。這種方法需要逐一計算每個反演時的溫度和風速的插值節點處(38×38=1444個節點)的不同入射角和海面風向構成的另外一個維度上的插值節點(5×10=50個節點)，共需要計算72 200(1444×50)組反演系數。

這種升級的反演算法，其實質是使用局部的線性關系描述整體空間內的函數關系。式(7)和式(8)聯立，可以看出對局部反演結果的插值同樣等同于對反演系數進行插值。這與反演系數關于入射角和海面風向的校正的本質是一致的。值得一提的是，當這種方法在插值間隔上趨于無窮小(格網無窮密)，維度拓展到包含所有RTM的參量(一般情況，包括海面溫度、海面風速、水汽含量、液水含量、海面鹽度、入射角、海面風向等)的希爾伯特空間時，就與LORA開發的最優化方法反演海洋大氣參數的算法核心相一致了。

4 非線性方法反演海洋大氣參數

輻射計觀測的輻射亮溫可以表達為海面溫度(Ts)、海面風速(W)、水汽含量(V)和液水含量(L)的函數：

TBi=fi(Ts,W,V,L)+ei

(9)

ei=eobsi+emi

(10)

式中：ei為觀測亮溫和RTM模擬亮溫之間的誤差，該誤差的來源包括儀器產生的觀測誤差eobsi和模型誤差emi。模型誤差emi包括RTM與真實環境的誤差ertmi和反演模型與RTM之間的誤差erei。

4.1 牛頓法反演海洋大氣參數

牛頓法是最常用的求解方程及方程組的方法。

牛頓方法中，式(9)被寫為多元一階泰勒展開的形式：

(11)

在不考慮觀測亮溫和RTM模擬亮溫之間的誤差(AMSR-E的ATBD所述，按照最小二乘法，這一項忽略與否并不關鍵)以及忽略高階無窮小量時，式(11)可以寫成如下形式：

AΔP+ΔT=0

(12)

式中：A為雅克比矩陣，有

(13)

(14)

已知x*是Ax=b的最小二乘解的充要條件為[10]，x*是ATAx=ATb的解。因此式(12)的最小二乘解為

ΔP=-(ATA)-1ATΔT

(15)

即

Pj=Pj(0)-(A(0)TA(0))-1A(0)TΔT(0)

(16)

因此有

Pj(k+1)=Pj(k)-(A(k)TA(k))-1A(k)TΔT(k)

(17)

式(17)就是牛頓法反演的迭代公式。迭代計算中，計算量最大的步驟在于雅克比矩陣的計算。這主要是由于fi的偏導數不方便獲得解析式，只能使用數值微分來計算。這是十分不易進行的，因為無法準確地在一定的計算精度下找到合適的差商步長[11]：步長過大結果必然不準確；步長過小會導致切線和割線無法分別，造成計算錯誤。一個有用的方法是變步長差商。但由于計算量的原因，這顯然不實用。除了這種情況，雅克比矩陣和(14)這兩者都導致需要大量的計算fi。

4.2 最優化方法反演海洋大氣參數

已知在式(9)中，TBi表示輻射計的第i個通道觀測亮溫，fi表示第i個通道的RTM。n個通道的亮溫數據就可以構建一個n元非線性超定方程組。于是，得到一個關于觀測亮溫和均方誤差(Mean Square Error，MSE)E的表達式：

TBi=fi(Ts,W,V,L)+ei

(18)

(19)

式中：ki是權重系數，表示各個觀測通道的貢獻程度。絕大多數情況，令ki為1或0是比較方便且合適的。

顯然，當E取得最小值時，對應的(Ts，W，V，L)為這個非線性超定方程組的最優解。固定步長牛頓法是解決這類問題的一種常用的、基礎的方法：

對于最優化問題minf(X)，k次迭代后得到Xk。將X=Xk在f(X)處二階泰勒展開，有

(20)

令

Q(X)=f(Xk)+2f(Xk)(X-Xk)=0

(21)

解(21)可得

Xk+1=Xk-[2f(Xk)]-1f(Xk)

(22)

然而，這種方法有顯而易見的缺陷。首先，由RTM和觀測亮溫構成的方程組的黑塞矩陣的局部正定在數學上是不容易嚴格證明的，通常只能夠從物理意義上說明它的正定性。其次，與式(13)類似，式(19)的二階偏微分方便的計算方法仍然是使用數值微分技術，這種技術的缺點就是計算量大且依賴于差商計算中兩點距離趨于無窮小的程度。

需要提出的是：式(19)中ki的設置是十分必要的。最優化方法與牛頓法，即傳統的最小二乘法，最大的區別，就是最優化方法區分了觀測誤差的貢獻程度。顯然，最小二乘法可以認為是ki=1的最優化方法的解法。另一個實際的例子是，當某一通道的數據被認為是錯誤時，式(13)的雅克比矩陣顯然會因此降低規模。與之對應，最優化方法中這一通道的k=0。因此，如果各個通道觀測數據的測量誤差分布不同以及對反演某一參數的貢獻不同時，采用合理的方法推算ki的值將會是一件很有意義的工作。

基于以上分析，LORA的研究團隊提出了使用Nelder-Mead(NM)算法的最優化方法反演大氣海洋參數的方法，力求解決上述各種方法在反演時的缺陷。NM算法的最優化反演能夠方便地求得式(19)中的最小值。NM算法[12-13]是一種直接搜索方法，借用了幾何學中單純形的概念。單純形是一種N維空間中具有N+1個頂點的特殊多面體。常見的單純形有1維空間中的直線，2維空間中三角形，3維空間中的四面體等。NM算法通過使用反射(reflection)，擴展(expansion)，收縮(contraction)和緊縮(shrink)四種運算技巧，不斷調整N維空間的單純形頂點使之逼近最優解。

5 仿真實驗

基于AMSR的ATBD提供的參數化的RTM，建立了仿真數據庫。通過這些仿真數據，使用參數化的RTM計算出模擬亮溫。仿真數據庫包含400 000組仿真數據，每組數據包含的物理量及分布范圍如下：海面風向(0°～180°)，觀測角(54.7°～55.3°)，海面鹽度(32‰～37‰)，海面溫度(0～30 ℃)，海面風速(0～20 m/s)，水汽含量(0～75 mm)，液水含量(0～0.3 mm)。數據分布采用平均分布的分布形式。

相對于多元線性回歸反演每個參數使用全部10個通道的數據，基于NM算法的最優化方法反演不同的海洋大氣參數時使用了不同的通道組合[14-15]。這樣的處理方法在提高反演結果分辨率的同時，還能有效地減少計算量、提高計算效率。表2記錄了反演不同參數的反演通道組合。

表2 非線性反演的通道選擇

5.1 線性方法和非線性方法的反演結果比較

使用400 000組仿真數據、2000×243組真實數據和3種反演算法進行了實驗，分析了多種反演算法的反演精確度及計算效率。表3是采用NM算法、多元線性回歸算法、升級版多元線性回歸算法分別對仿真數據進行反演實驗的結果。可見，升級版的多元線性回歸算法的反演精度優于多元線性回歸算法，NM算法的反演精度優于升級版的多元線性回歸算法。這幾種反演算法的反演結果與目標數據之間的偏差(Bias)都很小，在1×10-4的數量級左右。NM方法反演的2010年5月1日升軌數據的海面溫度結果如圖1所示，反演結果的分布合理。和線性回歸模型相比，非線性反演沒有引入新的誤差，因此,這種高精度的反演方法除了能提高反演產品的質量，還能夠為在軌儀器定標提供準確的數據。

表3 3種反演算法的反演結果與仿真數據的比較

圖1 NM算法反演AMSR-E升軌數據反演海面溫度Fig.1 AMSR-E ascending SST by Nelder-Mead algorithm

5.2 線性方法和非線性方法的計算效率比較

圖2 NM的最優化算法反演的反演效率Fig.2 Retrieval time of the non-linear Nelder-Mead algorithm

多元線性回歸由于采用簡單的線性計算方式，在業務化應用中具有計算效率高的優點。與之相比，非線性迭代方法由于算法本身的復雜性，其計算效率遠遠低于多元線性回歸。統計了使用基于NM算法的最優化反演400 000組數據時每一步的迭代次數以及反演AMSR-E的半軌2000×243組真實數據的反演耗時。對于4個不同的海洋大氣參數，反演的迭代次數按照海面溫度、海面風速、水汽含量和液水含量的順序遞減。使用單臺HP-820Z圖形工作站，基于NM的最優化方法反演2000×243組真實數據[16](見圖2)的耗時約為90 min；多元線性回歸反演這組數據的計算時間遠小于非線性方法，約為0.5 min。因此，基于NM算法的最優化方法應用于業務化工作需要解決計算耗時的問題。

6 結束語

從線性算法和非線性算法兩方面對反演方法進行了歸納討論，具體包括：多元線性回歸算法與升級版的多元線性回歸算法；牛頓迭代算法與最優化算法。實驗表明，LORA開發的基于NM的最優化反演算法與傳統的星載微波輻射計的線性反演算法相比，物理意義明確、精度高1個數量級。在解決了計算耗時問題之后，該方法具有一定的研究價值和應用潛力。未來，隨著星載微波輻射計系統的完善與發展，我們計劃將計算快速的線性反演算法和具有明確物理意義的非線性反演算法相結合。預計新的聯合反演方法能夠為我國海洋系列衛星產品的綜合應用提供更好的技術保障，同時為在軌定標等工作提供可靠的數據支持。

References)

[1]Kawanishi T, Sezai T, Ito Y, et al. The advanced microwave scanning radiometer for the earth observing system (AMSR-E), NASDA’s contribution to the EOS for global energy and water cycle studies[J]. IEEE Transactions on Geoscience and Remote Sensing, 2003, 41 (2)：184-194

[2]張升偉, 李靖, 姜景山, 等.風云3號衛星微波濕度計的系統設計與研制[J]. 遙感學報, 2008, 12 (2): 199-207

Zhang Shengwei, Li Jing, Jiang Jingshan,et al. Design and development of microwave humidity sounder for FY-3 meteorological satellite[J]. Journal of Remote Sensing, 2008, 12(2):199-207 (in Chinese)

[3]蔣興偉,林明森,宋清濤. 海洋二號衛星主被動微波遙感探測技術研究[J]. 中國工程科學, 2013, 15 (7):4-11

Jiang Xingwei, Lin Mingsen, Song Qingtao. Active and passive microwave remote sensing technology of the HY-2A ocean satellite mission[J]. Engineering Sciences, 2013, 15 (7): 4-11 (in Chinese)

[4]Goodberlet M A,Swif C T,Wilkerson J. Ocean surface wind speed measurements of the Special Sensor Microwave/Imager (SSM/I)[J]. IEEE Transactions on Geoscience and Remote Sensing, 1990, 28 (5): 823-828

[5]Wentz F, Meissner T. Algorithm theoretical basis document:AMSR ocean algorithm,Version 2 [R]. Santa Rosa:Remote Sensing System,2000

[6]Wentz F J, Meissner T. Supplement 1 algorithm theoretical basis document for AMSR-E ocean algorithms[R]. Santa Rosa:Remote Sensing System, 2007

[7]王振占. 用19.35 GHz星載微波輻射計(SSM/I)亮溫反演海面風速[J]. 海洋技術, 2003, 22(2): 1-6

Wang Zhennzhan. A model for inversing sea surface wind speeds by using SSM/I 19.35 GHz vertical and horizontal brightness temperatures[J]. Ocean Technology, 2003, 22(2): 1-6 (in Chinese)

[8]王振占. 海面風場全極化微波輻射測量[D]. 北京. 中國科學院空間科學與應用研究中心, 2005

Wang Zhenzhan,Sea surface wind vector measured by polarimetric microwave radiometer[D]. Beijing:National Space Science Center,CAS, 2005 (in Chinese)

[9]王兆徽, 宋清濤, 蔣興偉, 等. 星載微波輻射計的非線性反演算法[J]. 高技術通訊, 2015(5)

Wang Zhaohui, Song Qingtao, Jiang Xingwei, et al. Non-linear retrieval algorithm for passive satellite microwave radiometer[J]. Chinese High Technology Letters, 2015(5) (in Chinese)

[10]Lawson C L, Hanson R J. Solving least squares problems[M]. Philadelphia: SIAM: 1974

[11]Mathews J H, Fink K D. Numerical methods using MATLAB[M]. Upper Saddle River: Prentice Hall, 1999

[12]Nelder J A, Mead R. A simplex method for function minimization[J]. Computer Journal, 1965, 7(4): 308-313

[13]Powell M J. On search directions for minimization algorithms[J]. Mathematical Programming, 1973, 4(1): 193-201

[14]Gentemann C L, Meissner T,Wentz F J. Accuracy of satellite sea surface temperatures at 7 and 11 GHz[J]. IEEE Transactions on Geoscience and Remote Sensing, 2010, 48(3): 1009-1018

[15]Martin S, An introduction to ocean remote sensing[M]. Cambridge: Cambridge University Press: 2004

[16]National Snow & Ice Data Center. NSIDC. AMSR-E L2A data[ED/OL]. [2014-11-01]. Ftp://n5eil01u.ecs.nsidc.org/AMSA/AE_L2A.003

(編輯：張小琳)

Comparison of Linear and Nonlinear Retrieval Algorithm for Spaceborne Microwave Radiometer

WANG Zhaohui1LIU Yuxin1,2SONG Qingtao1,2WU Binfeng1

(1 National Satellite Ocean Application Service, Beijing 100081, China) (2 Key Laboratory of Space Ocean Remote Sensing and Application, SOA, Beijing 100081, China)

The paper analyzes the channel choice of ocean retrieval algorithm and compares the differences between linear and nonlinear retrieval algorithm by spaceborne microwave radiometer data process. To discard the low-frequency channel in nonlinear retrieval algorithm gradually is a direct and effective way to improve the spatial resolution of retrieval products. Several analyses and experiments for retrieval algorithm indicate that compared with the linear retrieval algorithm, the optimization algorithms based on Nelder-Mead method has higher accuracy, application and development potential.

spaceborne microwave radiometer; retrieval algorithm; Nelder-Mead algorithm

2015-06-09；

2015-07-15

國家自然科學基金(41276019)；海洋公益性行業科研專項(201305032)；中法海洋衛星預先研究項目

王兆徽，男，碩士，研究方向為微波海洋遙感。Email:wzh@mail.nsoas.org.cn。

V443

A

10.3969/j.issn.1673-8748.2015.04.021