李 芹,鐘若飛

(首都師范大學(xué)三維信息獲取與應(yīng)用重點實驗室,北京 100048)

模擬AM SR-E數(shù)據(jù)的地表多參數(shù)反演

李 芹,鐘若飛

(首都師范大學(xué)三維信息獲取與應(yīng)用重點實驗室,北京 100048)

為了提供更適用的地表參數(shù)反演方案,對 Njoku等模擬AMSRE-E數(shù)據(jù)的多參數(shù)反演工作重新進(jìn)行正向模擬和算法改進(jìn),算法用MATLAB開發(fā)實現(xiàn)。反演結(jié)果表明,改進(jìn)后的反演方法不僅在精度上略高于原方法,而且還能反演裸露地表可實際測量的粗糙度參數(shù),并在更大的植被含水量范圍內(nèi)達(dá)到較高的反演精度。

多參數(shù)反演;Levernberg-Marquardt;迭代反演算法;微波輻射測量

0 引言

對全球環(huán)境變化進(jìn)行監(jiān)測需要獲得土壤水分、植被含水量、地表溫度及地表粗糙度等參數(shù)的時空分布信息。從遙感獲取數(shù)據(jù)的角度來看,即使在同一個傳感器平臺上,不同的觀測角度、頻率和極化方式以及地表參數(shù)的變化對亮溫的貢獻(xiàn)程度是不同的。研究適用的多參數(shù)反演方法,不但能減小對地面輔助數(shù)據(jù)的依賴程度,還可以更充分地利用微波輻射計提供的低成本、大范圍、高更新率的觀測數(shù)據(jù)。以往對于多參數(shù)反演算法的研究大部分基于最小二乘法原理[1-4],其中 Njoku[1]提出的針對 AM SR-E數(shù)據(jù)的土壤水分反演迭代算法——Levernberg-Marquardt算法,是一種經(jīng)典的多參數(shù)反演算法,但此算法的適用條件是植被含水量低于 1.5 kg/m2,并用模擬數(shù)據(jù)給出了 3個反演參數(shù) (即土壤水分、植被含水量、地表溫度)的均方根誤差[2]。最近,Shi等[5]提出了更適合AM SR-E這種大入射角、高頻率傳感器的 Qp參數(shù)化模型來代替?zhèn)鹘y(tǒng)的Q/H模型;新發(fā)射的傳感器如 SMOS等數(shù)據(jù)的應(yīng)用也需要研究反演算法。根據(jù)這些研究進(jìn)展,本研究對 N joku等模擬 AMSR-E數(shù)據(jù)的多參數(shù)反演工作重新進(jìn)行正向模擬和算法改進(jìn),以期得到更好的反演效果,為今后類似的工作提供參考依據(jù)。

1 模擬數(shù)據(jù)

本研究采用模擬數(shù)據(jù)來驗證地表多參數(shù)的反演算法。按照表 1和表 2所列的參數(shù),采用正向模型,建立裸露地表和植被覆蓋地表兩種情況下的模擬數(shù)據(jù)庫。

表 1 大氣和植被的主要經(jīng)驗參數(shù)及輔助數(shù)據(jù)①來自 2003年 GAM E-Tibe實驗數(shù)據(jù),這里僅作為模擬數(shù)據(jù)的參考經(jīng)驗值,對反演無影響。Tab.1 Mainempirical parameters and ancillary data about atmosphere and vegetation

表 2 模擬裸露地表發(fā)射率數(shù)據(jù)庫參數(shù)Tab.2 Parameters needed for bare soil’s simulated database

2 反演方法

2.1 Qp模型的引入

地表參數(shù)遙感反演的前提就是建立正向模型,即根據(jù)物理基礎(chǔ)建立起地表參數(shù)與衛(wèi)星獲得的亮溫之間的函數(shù)關(guān)系:該模型通過輸入地表參數(shù),如土壤水分、植被含水量和地表溫度,根據(jù)亮溫與地表和大氣相關(guān)變量間的關(guān)系,得到相應(yīng)的亮溫數(shù)據(jù)。一般的正向模型將衛(wèi)星與地表間的介質(zhì)分為大氣、植被和土壤3層,植被層可看作是一個位于粗糙土壤之上的單次散射層,因而需要從地表、植被和大氣 3方面考慮。地表土壤的發(fā)射率通過介電常數(shù)和觀測角與土壤水分相關(guān)。除此之外,植被覆蓋度、地表粗糙度、土壤溫度以及地形等因素都對土壤發(fā)射率有影響。

在粗糙地表情況下,傳統(tǒng)模擬裸露地表微波輻射特性的方法是采用Q/H模型以及一系列由此衍生出來的模型。本研究引入 Shi等在 2006年研究并提出的Qp參數(shù)化模型。該模型以A IEM模型提供的模擬數(shù)據(jù)庫為參照,對比傳統(tǒng)的Q/H模型和 Qp參數(shù)化模型在 H、V兩個極化和極化比值下的模擬精度,發(fā)現(xiàn)Qp參數(shù)化模型更適合AM SR-E等大入射角輻射計,并且給出了從 6～37 GHz的Qp計算公式。該公式建立了 Qp參數(shù)與粗糙度參數(shù) s(均方根高度)和 l(相關(guān)長度)的商 (即 s/l)之間的函數(shù)關(guān)系[5]。

在有植被覆蓋的情況下,選擇τ-ω模型來模擬植被頂層的亮溫 Tbp[7]。該模型適用于低矮植被覆蓋地表(如農(nóng)田),且對低頻下植被覆蓋地表的微波輻射描述逼近真實情況。植被的覆蓋雖然減弱了來自土壤的輻射,但同時增加了植被本身的輻射,使得土壤水分變化對植被頂層亮溫的變化貢獻(xiàn)減小,即在有植被覆蓋的時候,土壤水分的敏感性降低而植被含水量的敏感性升高,從而土壤水分反演變得困難。此時的土壤水分反演思路主要有兩種:一是估算并去除植被的影響,轉(zhuǎn)化為在裸露地表下反演土壤水分;二是將土壤水分和植被含水量同時反演。本研究采用第二種思路。

衛(wèi)星觀測的亮溫實際上是大氣頂層的亮溫,在微波輻射信號從地表到達(dá)傳感器的傳輸過程中,還受到大氣的影響,一般采用經(jīng)驗關(guān)系來消除[7]。

2.2 參數(shù)敏感性分析

參數(shù)敏感性分析是為了輔助反演時選擇合適的通道組合以及分析反演結(jié)果,通常表示為在給定的參數(shù)基準(zhǔn)值 x0和變化范圍 Xj下,歸一化的參數(shù)敏感性系數(shù)值[1,2]。由于引入了 Qp模型代替 Q/H模型作為正向模擬的模型,同時增加了一個待反演參數(shù)——粗糙度,并且擴(kuò)大了植被含水量的范圍,因此,以下將要重新分析所有通道 (除了 89 GHz以外)對各個待反演參數(shù)的敏感性,并將其與 Njoku等在 6～18 GHz 3個通道的工作進(jìn)行對照。

表 3列出了裸露地表情況下AM SR-E 10個通道 (除 89 GHz以外)的亮溫對地表粗糙度 (m)、土壤含水量(sw)、地表溫度 (ts)的敏感性。當(dāng)實際地表為裸露土壤或接近裸露土壤 (植被覆蓋率很低)時,相比植被含水量,更受關(guān)注的地表參數(shù)是 m、sw和 ts。為了簡化正向模型,這里固定植被含水量vwc=0 kg/m2;為了方便比較不同頻率和極化的地面發(fā)射率對這 3個參數(shù)的敏感性,表中的亮溫是尚未考慮大氣影響的近地面亮溫,與表 4和表 5中針對植被覆蓋地表的大氣頂層的亮溫是不同的。這樣的前提并不影響對反演算法的檢驗。

表 3 裸露地表情況下的各參數(shù)敏感性Tab.3 Parameters’sensitivities in case of bare soil surface

從表 3可以看出,m和 sw在 H極化時的敏感性明顯強(qiáng)于 V極化時的敏感性;ts在兩種極化時的敏感性相當(dāng)。m和 ts在各個頻率的敏感性相當(dāng),sw在低頻的敏感性明顯強(qiáng)于在高頻的敏感性。所以,同樣的通道數(shù)目,低頻組合比高頻組合更容易獲得較好的反演結(jié)果。此外,由于正向模型包括了 Qp模型,則在 H極化的粗糙地表反射率中V極化的光滑地表反射率的貢獻(xiàn)不為 0,同時使用雙極化數(shù)據(jù)可以消除極化差,較之單一的 H極化的頻率組合更容易獲得較好的反演結(jié)果。

在有植被覆蓋的地表,隨著植被覆蓋度的增大,vwc會增加,同時冠層對地表輻射的衰減效應(yīng)增強(qiáng)。表 4和表 5分別列出了 vwc≤1.5 kg/cm2和vwc>1.5 kg/cm2兩種情形下大氣頂層亮溫對 sw、vwc及 ts的敏感性。

表 4 低植被覆蓋情況下的各參數(shù)敏感性Tab.4 Parameters’sensitivities in case of lower vegetation cover surface

表 5 高植被覆蓋情況下的各參數(shù)敏感性Tab.5 Parameters’sensitivities in case of higher vegetation cover surface

從表 4、5可以看出,sw的敏感性在 H極化時顯著強(qiáng)于在 V極化時的;vwc和 ts的敏感性在兩種極化時基本相當(dāng);敏感系數(shù)最小的是 ts,預(yù)測 3個參數(shù)反演結(jié)果誤差最大的是 ts。隨著 vwc值的變大,sw和 vwc的敏感性反而大大變?nèi)?ts的敏感性也變?nèi)?反演的難度增加;sw敏感性強(qiáng)弱變化趨勢為:6GHz>10GHz>18GHz,而 vwc敏感性強(qiáng)弱的變化趨勢為:6GHz<10GHz<18GHz,ts的敏感性隨頻率改變的變化幅度很小。

2.3 算法改進(jìn)和實現(xiàn)

為了與 Njoku算法進(jìn)行對比,本文所采用的迭代算法仍然是基于 Levernberg-Marquardt原理[2]。

這里的反演問題歸結(jié)為一個結(jié)合了微波遙感定量反演理論的非線性系統(tǒng)的最優(yōu)化問題,涉及到一些數(shù)學(xué)原理,故考慮使用 MATLAB來輔助實現(xiàn),以便于調(diào)用其工具箱的矩陣運(yùn)算和求解方程組等函數(shù)。MATLAB的最優(yōu)化工具箱提供的 lsqnonlin函數(shù)雖然可以選擇Levernberg-Marquardt方式來實現(xiàn)優(yōu)化,但層層調(diào)用,且封裝了其他的優(yōu)化方式,不利于修改和調(diào)用。本研究在MATLAB幫助菜單下的web source中找到Miroslav Balda提供的 LM Fnlsq函數(shù)[6],參考并自定義算法。

算法的輸入?yún)?shù)為最大迭代次數(shù)、阻尼系數(shù)初始值、待反演參數(shù)值的邊界范圍、代價函數(shù)的容差及待反演參數(shù)的容差,具體實現(xiàn)步驟如下:

(1)定義代價函數(shù) Fun和雅克比矩陣 J。Fun即各個通道的亮溫觀測值 (用模擬值加上高斯隨機(jī)噪聲表示)與用估計值 xk(第 k次迭代后的 x)計算的亮溫估計值之差的加權(quán)平方和,用各個通道觀測噪聲標(biāo)準(zhǔn)差平方的倒數(shù)作為權(quán)重;J可以用解析函數(shù)的一階泰勒導(dǎo)數(shù)或者有限數(shù)值差分求導(dǎo)。由于模型公式的復(fù)雜性,若直接計算一階偏導(dǎo)數(shù)作為雅克比矩陣,不但計算復(fù)雜且可能會導(dǎo)致無法收斂,所以這里采用后一種方式,但其缺陷是計算量大,速度慢。

(2)初始值設(shè)定。初始值包括待反演參數(shù) x0(可以任意選取,但為了盡快收斂,一般取它的中間值)、迭代次數(shù)最大值 Max Iter、阻尼系數(shù) l(缺省值是0。接近 0時,切換到 GN方法來搜索解;接近無窮大時,切換到 SD方法來搜索解。如果代價函數(shù)的值逐漸變小,1/scale變小,加快收斂速度,反之,scale變大,增加解的搜索范圍,scale的取值一般在 2～10之間)、代價函數(shù)容差 FunTol及待反演參數(shù)容差XTo l。

(3)調(diào)用函數(shù)。在MATLAB下自定義。

(4)限制待反演參數(shù)值的邊界范圍。邊界范圍參見表 2。如果迭代反演的結(jié)果超出邊界范圍,則調(diào)整算法的輸入?yún)?shù),重復(fù)步驟 (3)。

(5)滿足結(jié)束條件時 (達(dá)到最大迭代次數(shù)、小于代價函數(shù)容差 FunTol及小于待反演參數(shù)容差XTo l),結(jié)束迭代。

3 驗證與分析

給模擬亮溫數(shù)據(jù)加上高斯白噪聲來代表衛(wèi)星實際觀測亮溫。在對所有數(shù)據(jù)進(jìn)行反演前,先抽樣反演,對算法的輸入?yún)?shù)進(jìn)行訓(xùn)練,得到對應(yīng)最高反演精度輸入?yún)?shù)的經(jīng)驗配置。

3.1 裸露地表參數(shù)的反演

在裸露地表條件下,待反演參數(shù)有綜合粗糙度(m)、土壤含水量 (sw)及地表溫度 (ts)。此時在反演誤差的貢獻(xiàn)中少了植被層的影響,土壤水分反演的精度較高。實驗發(fā)現(xiàn),在 6GHz和 10GHz的雙極化觀測組合下,地面溫度的反演誤差較大。這里假設(shè)已經(jīng)獲得了 6GHz、10GHz、18GHz、23GHz及37GHz雙極化的大氣參數(shù)。大氣校正后,由大氣頂層的亮溫得到近地面的亮溫,然后再進(jìn)行反演。實驗發(fā)現(xiàn),加入 37GHz雙極化觀測后,地表溫度反演的精度明顯提高,這與前面的參數(shù)敏感性分析結(jié)果是一致的。4通道算法的最佳頻率組合是 6GHz和37GHz的雙極化,共 4個通道的近地面亮溫,并且,10個通道的觀測組合得到的反演精度最高,圖 1顯示了分別添加 0 K和 0.3 K的高斯白噪聲后 3個地表參數(shù)的反演結(jié)果。

圖 1 裸露地表情況下 m、s w、t s的反演結(jié)果Fig.1 Retrieval results of m、s w and t s in case of bare soil surface

從圖 1可見,高斯隨機(jī)噪聲從 0 K增加到 0.3 K時,反演精度下降,說明該算法對噪聲比較敏感。表 6列出了不同噪聲條件下,3個地表參數(shù)的反演精度。

表 6 不同高斯噪聲條件下的地表參數(shù)反演精度Tab.6 Retrieval RM SE of m、s w and s t under various guassnoise conditions

任意改變初始值,反演結(jié)果基本不變。將反演結(jié)果輸入到正向模型中,各個通道亮溫的方差也均低于添加的噪聲,表明多通道反演具有降噪效果。

3.2 植被覆蓋地表的參數(shù)反演

在植被覆蓋條件下,待反演參數(shù)有土壤含水量(sw)、植被含水量 (vwc)、地表溫度 (ts)。N joku等在算法的模擬研究和實測數(shù)據(jù)驗證中發(fā)現(xiàn),當(dāng) vwc<1.5 kg/m2時,數(shù)據(jù)反演結(jié)果誤差較小,并且在模擬研究中,在給模擬數(shù)據(jù)庫添加了 0.3 K的高斯白噪聲后,3個參數(shù)的反演誤差分別小于 0.06 g/cm3、0.1 kg/m2和 1.5℃。但該算法不適合反演 vwc>1.5 kg/m2的數(shù)據(jù)。筆者認(rèn)為其原因為:①反演算法的實現(xiàn)僅從數(shù)學(xué)原理出發(fā),沒有考慮遙感模擬數(shù)據(jù)庫的參數(shù)范圍以及遙感模型的適用范圍,迭代可能得到了無意義解,從而在調(diào)用遙感模型代碼時提前跳出程序;②當(dāng)植被含水量大于 1.5 kg/m2時,輻射亮溫對土壤含水量的敏感性大大減小,而代價函數(shù)在逼近最小值的過程中存在多解,在滿足某一個迭代結(jié)束條件之前,反演得到的結(jié)果隨機(jī)性較大。

針對上述情況,本研究對反演方法進(jìn)行了改進(jìn),具體措施為:

(1)對地表參數(shù)進(jìn)行抽樣,試用反演算法,獲得最高反演精度所對應(yīng)的一些參數(shù)配置,并將這些參數(shù)配置作為經(jīng)驗值;

(2)在反演算法中增加對反演參數(shù)值邊界條件的限制(迭代算法的缺陷就在于對超定方程組求最優(yōu)解時會產(chǎn)生多解,為了保證最優(yōu)解落在地表參數(shù)的物理范圍或者實際范圍內(nèi),給出了一個數(shù)學(xué)角度上的修正方案,即一旦超過范圍就取上下限),假如迭代點的土壤水分含量小于 0,則令它為 0.01,如果大于 0.45,則令它為 0.44;

(3)計算有限差分雅克比矩陣代替。如果將正演的植被模型看作一個關(guān)于待反演的 3個地表參數(shù)的函數(shù),這個函數(shù)是非常復(fù)雜的非線性函數(shù),不方便計算每次迭代后 3個地表參數(shù)的偏導(dǎo)數(shù),所以在計算雅克比矩陣時,用有限差分雅克比矩陣代替,這樣設(shè)計有利于保證反演結(jié)果的收斂和算法的穩(wěn)定性,但是算法的時間和空間復(fù)雜度變大。

從 2.2節(jié)參數(shù)敏感性分析中得到用 3個低頻、雙極化、共 6個通道數(shù)據(jù)一起參加反演能獲得最好的反演結(jié)果的推斷,在實驗中證實了這一判斷的正確性。當(dāng)用 6個通道同時反演 3個參數(shù)時,圖2顯示了分別在噪聲 0 K和 0.3 K的條件下的反演結(jié)果。

圖 2 植被覆蓋地表情況下 s w、t s、v wc的反演結(jié)果Fig.2 Retrieval results of s w t s and v wc in case of vegetation cover surface

分析認(rèn)為,出現(xiàn)少數(shù)反演誤差比較大的點,主要是因為迭代反演的不穩(wěn)定性。vwc和 sw的敏感性強(qiáng)弱相當(dāng),且明顯強(qiáng)于 ts,所以反演精度 RM SE(sw)≈RM SE(vwc)>RM SE(st)。隨著 vwc的增大,vwc和 sw的敏感性都變?nèi)?反演誤差相對較大的點集中在 vwc值靠近邊界最大值的一端,這一現(xiàn)象在噪聲增大為0.3 K時更明顯,同時也說明算法對噪聲比較敏感。

表 7 不同 v wc變化范圍內(nèi) 3個地表參數(shù)的反演精度 (高斯白噪聲為 0.3 K)Tab.7 Retrieval RM SE of s w、v wc and t s under various v wc value range(Noise=0.3 K)

任意改變初始值,反演結(jié)果基本不變。將反演結(jié)果輸入到正向模型中,各個通道亮溫的方差也均低于添加的噪聲,表明多通道的反演方法取得了降噪效果。

4 結(jié)論

(1)本文對 Njoku等的 Levernberg-Marquard t迭代方法進(jìn)行了改進(jìn),改進(jìn)后的反演方法在一定程度上提高了反演精度。

(2)在實際反演中,衛(wèi)星的通道噪聲不是簡單的高斯白噪聲。本文用模擬數(shù)據(jù)加上高斯噪聲來代替實際數(shù)據(jù),目的是便于與 Njoku等的研究工作進(jìn)行對比,從理論意義上為相關(guān)研究工作提供參考。

(3)本次研究發(fā)現(xiàn)用 Qp參數(shù)化模型代替 Q/H模型,使得待反演的粗糙度和作為已知信息的粗糙度與實際的粗糙度測量聯(lián)系起來,而不是 N joku等使用的 Q/H模型中僅僅有粗糙度意義的 H參數(shù)。粗糙度信息可以由 DEM數(shù)據(jù)庫或者主動傳感器數(shù)據(jù)反演提供,為主被動傳感器聯(lián)合反演搭建了橋梁。由于 AM SR-E對粗糙度的敏感性不如主動傳感器數(shù)據(jù),限制了地表參數(shù)反演精度的進(jìn)一步提高,如何結(jié)合主動數(shù)據(jù),如高分辨率雷達(dá)在粗糙度估算方面的優(yōu)勢,來提高土壤水分反演精度,也是當(dāng)前 SMAP任務(wù)的研究熱點。

(4)從數(shù)學(xué)角度來看,引入遙感反演迭代算法的缺陷在于初始值和迭代結(jié)束條件等 (對應(yīng)程序中函數(shù)的輸入?yún)?shù))對最優(yōu)解的影響較大。本次研究僅僅是在非線性迭代反演算法的實現(xiàn)上做了一點小的改進(jìn),如何采取措施進(jìn)行更理想的改進(jìn)是當(dāng)前非線性迭代反演的一個難題,今后仍然有不少研究工作需要開展。

[1] Njoku E G,L L i.Retrieval of Land Surface Parameters Using Passive Microwave Measurements at6～18 GHz[J].IEEE Transactions on Geoscience and Remote Sensing,1999,7(3):79-93.

[2] Njoku E G.AM SR Land Surface Parameters[M/OL].USA:NASA Jet Propulsion Laboratory,1999.http://eospso.gsfc.nasa.gov/ftp_ATBD/REV IEW/AMSR/atbd-amsr-land.pdf.

[3] Paloscia S,Macelloni G,Santi E,et al.A Multifrequency Algorithm for the Retrieval of Soil Moisture on a Large Scale Using Microwave Data from SMMR and SSM/I Satellites[J].IEEE Transactions on Geoscience and Remote Sensing,2001(39)8:1655-1661.

[4] Aurelio Cano.The SMOS Mediterranean Ecosystem L-Band Characterisation Experiment(MELBEX-I)over Natural Shrubs[J].Remote Sensing of Environment,2010,(114)4:844-853.

[5] Shi JC,Jiang L M,Zhang L X,et al.A Parameterized Multifrequency-Polarization Surface Emission Model[J].IEEE Transactions on Geoscience and Remote Sensing,2005,(43)12:2831-2641.

[6] MATLABW eb Source[EB/OL].[2009-01-27].http://www.mathworks.cn/matlabcentral/.

[7] Richard A M de Jeu.Retrieval of Land Surface Parameters Using Passing Microwave Remote Sensing[D].Thesis VR IJE University Amsterdam,2003.

(責(zé)任編輯:刁淑娟)

Multiple Surface Parameters Retrieval of Simulated AM SR-E Data

LI Qin,ZHONG Ruo-fei
(Capital Normal University Key Lab of3D Information Acquisition and Application,Beijing 100048,China)

In order to provide more suitable surface retrieval methods,the authors reconducted forward modeling and made algorithm improvement for multiple parameters retrieval of AM SRE-E data simulated by Njoku developed Levernberg-Marquardt iteration algorithm,and realized projectcode in MATLAB.The result show s that the imp roved retrieval method not only has somewhat higher precision than the original method but also can make retrieval of measurable roughness parameters of the exposed surface. In addition,it can attain relatively high precision in a wider water content range of vegetation.

Multiple parameter retrieval;Levernberg-Marquardt;Iterational inversion algorithm;Microwave radiometer measurement

李 芹 (1986-),女,碩士研究生,主要從事微波遙感定量反演地表參數(shù)、G IS應(yīng)用與開發(fā)工作。

TP 79

A

1001-070X(2011)01-0042-06

2010-06-19;

2010-06-26

國家自然科學(xué)基金項目 (編號:40701127)。