遙感統計推斷理論與應用初探

朱渭寧

浙江大學海洋學院, 東海實驗室, 浙江 舟山 316021

引 言

傳統的遙感圖像分析方法有兩種, 即分類和反演。 分類是一種定性方法, 它通常根據光譜信息, 獲取地物的類別信息, 如土地利用類別、 水體營養化類別[1-2]; 而反演是一種定量方法, 它通常根據光譜信息, 獲取地物參數的數量信息, 如地表溫度、 水體中的黃色物質濃度[3-5]。 這兩種方法在遙感中都有大量的理論與應用研究[6-7]。

本文提出一類介于分類和反演之間, 既含有定性分析, 也含有定量分析的遙感圖像分析新方法, 稱之為“遙感統計推斷”(以下有時簡稱“遙感推斷”或“推斷”)。 遙感推斷的基本過程是根據地物對象的光學屬性(如反射率)的統計信息(如均值、 方差、 概率分布函數), 去推斷地物對象的非光學屬性(如農田濕度、 水體渾濁度)的統計信息(如均值、 方差、 概率分布函數)。

在更深入地介紹遙感推斷之前, 先看一個例子。 如圖1所示, 通過遙感分類, 可以將一個湖泊或者其部分水域劃分為“寡營養水”或“渾濁水”; 通過遙感反演, 可以計算每個像素所對應的水體中含有的CDOM(黃色物質)的濃度大小; 而通過遙感推斷, 則可以獲取黃色物質濃度在湖體區域的統計和分布特征, 如均值、 方差以及對數正態分布。 由此可以看出, 遙感推斷與傳統反演和分類方法的最大區別在于, 它是基于研究地物對象的光譜分布, 去計算該對象其他屬性的分布, 因此是一種從“分布”到“分布”的過程; 而反演和分類則是基于地物對象的光譜本身, 去計算該對象其他屬性的類別或數值信息, 因此是一種從“光譜”到“類別/數量”的過程。 反演和分類是基于光譜的, 因此它們是一種基于像素(pixel-based)的方法, 也就是說, 即使對于一個像素, 也可以運用相應的分類或反演方法。 與之相反, 推斷是一種基于對象(object-based)的方法, 對于一個像素(樣本), 其分布或統計特征是沒有意義的, 因此也無法使用遙感推斷方法。

圖1 遙感統計推斷與遙感傳統分類和反演的區別(以湖泊水質研究為例)

從這三種方法的學科和理論基礎上看, 遙感分類其實就是傳統信息科學中的模式識別技術在遙感上的應用, 因此它涉及了較多統計學和數據科學的內容, 而相關的光學理論則涉及較少, 很多分類算法甚至不需要對地物光譜進行大氣校正。 而遙感反演則相反, 它是以輻射傳輸和幾何光學為基礎的, 雖然有些反演算法也是經驗統計性的(如簡單線性回歸或復雜神經網絡), 但這些經驗算法以及其他更多的半分析或分析算法在圖像處理、 波段選擇、 模型設計等方面都會涉及較多的諸如大氣校正、 地物吸收與散射特征等輻射傳輸和幾何光學的內容。 遙感推斷所對應的物理和光學過程也是以輻射傳輸為基礎的, 但是它更多涉及了統計光學的一些理論與方法[8], 即研究一個光場或一群光子參數的統計分布特性是如何形成和變化的。 如果感興趣的某種地物參數(如水體中的葉綠素濃度)呈現出某種統計分布, 那么這種分布通過輻射傳輸后產生的光學變量(如遙感反射率)也會呈現出一種相應的分布。 舉一個理想條件下的簡單例子, 如果某個湖泊在所有像素點的水質參數都是常數, 那么在任一波段上, 在所有像素點觀測到的遙感反射率也應該是一個常數。 可以想象, 從一個注滿自來水的露天游泳池中提取的每個水體像素的光譜曲線都應該一模一樣。 反之, 如果觀測到遙感反射率呈現出某種分布, 那么說明一定是地物或者參與輻射傳輸的某些變量的屬性具有某種分布, 而遙感統計推斷就是研究這些分布是怎么在地物對象的非光學屬性和光學屬性之間進行傳輸和改變的, 其最終目的就是獲知能夠引起遙感光學變量產生分布的地物屬性最有可能的分布特征。

既然已經有了分類和反演方法, 為什么還要研究和使用遙感統計推斷呢? 這是需求決定的。 通過分類, 我們可以了解對象的一些定性性質, 獲取其基本的屬性信息, 如地表是土壤還是水體, 是紅土還是黃土, 是清水還是渾水, 對于一些較為宏觀的用戶, 例如土地資源、 水資源的管理者、 決策者, 這些信息可能就足夠了。 但是如果需要進一步深入地了解地物屬性變化的機理機制, 就需要定量地分析其局部和細節特征與變化, 此時就需要運用遙感反演方法, 獲取這些屬性的數值信息, 提供給專業科研人員進一步分析和使用。 但是反演方法常會遇到兩類難題： (1)反演準度問題, 即由于地物過于復雜、 病態反演等原因, 反演模型和算法自身存在一些缺陷, 反演準確性較低, 不確定性較高[7]; (2)反演效率問題, 即使反演算法較為準確, 但是如果算法自身過于復雜, 研究對象的數據量巨大且有實時性需求時, 反演效率較低[9]。 例如, 一個湖體在高分辨的影像上可能有上百萬個像素, 對每個像素都運行一遍反演算法, 成本較高。 事實上, 對于很多諸如環境監測的用戶來說, 不需要精確地知道在每個像素位置上水體或土壤的溫度值, 他們更感興趣的往往是對象的一些統計特征, 例如一個水庫葉綠素分布的最大值、 最小值、 均值、 方差等等。 所以, 用常規的反演方法, 算出每個像素的葉綠素濃度值, 再統計出整個水庫葉綠素的均值、 方差, 可能既不準又費時。 如果傳統的反演方法是“先反演, 再統計”, 那么遙感統計推斷則是反其道： “先統計、 再推斷”。 相比于復雜的反演, 做統計是相對簡單和快速的過程, 與其對水體像素做幾百萬次反演, 再做一次統計, 不如直接先對光譜做一次統計, 再做一次推斷。

遙感推斷的另一個益處是為反演提供輔助信息, 這就如同分類給反演提供輔助信息。 由于地物和輻射傳輸的復雜性, 很多反演算法的局限性都較大。 在某種時空背景下針對某種地物特征的較為準確的算法, 當地物特征和時空背景轉變時, 可能會變得不準確, 因此算法的參數、 模型、 甚至算法自身都需要做相應的調整。 而此時分類就可以發揮作用, 例如二類水體黃色物質的反演一直是一個難題, 近年來提出了先分類再反演的算法[10], 如先將水體分為清水和渾水兩類, 對這兩類水體分別運用有針對性的不同算法, 從而提高整個水體的反演準度[11]。 遙感推斷扮演的也是類似的角色, 針對特別復雜的對象時, 可以采取“先分類、 再推斷、 后反演”的策略。 遙感反演一般遵從“實地采樣”、 “建模驗證”和“圖像應用”三個步驟。 實地采樣獲取的數據, 如在某水庫采集的30個水樣, 往往是對象(水庫)的一個“樣本集”而不是“總體”, 第二步反演建模也是基于這個樣本集的, 而將建立好的模型應用于遙感影像時, 輸入的數據卻是對象的“總體”, 即整個水庫的所有像素, 因此當“樣本集”與“總體”的統計和分布特征存在偏差時, 反演模型的結果也會產生相應的偏差。 而遙感推斷正是對地物對象“總體”的一種描述, 我們根據推斷結果校正和調整反演算法的模型和參數, 使得將算法應用于總體時, 能夠得到合理的、 與分布推斷相吻合的、 沒有統計偏差的準確結果。

遙感推斷的第三個用途是對地物對象進行定性分析。 從分類和定性的角度看, 推斷方法給出的是一個地物對象屬性的統計和分布信息, 而這些信息可以為傳統遙感分類分析提供一些新的視角, 也就是說, 不僅僅可以依據對象光譜自身的簽名特征對像素分類, 也可以依據地物對象自身屬性的分布特征對地物整體進行分類。 顯然, 根據地物屬性特征分類是更有指向性、 更為精確的分類方法, 因為根據輻射傳輸, 遙感所觀測到的光譜簽名在某種意義上是一種“二手”信息, 它是由作為“第一手”信息的地物屬性特征所決定的。 根據葉綠素、 黃色物質、 懸浮泥沙的濃度分布特征直接給一個水體分類, 顯然要比間接使用其光譜簽名特征分類更為準確。 由于地物非光學屬性的分布特征是通過其光學屬性的分布特征推斷而來, 因此有時甚至可以省略具體的推斷過程, 直接將地物光學屬性的分布特征與其類別相關聯, 即在遙感分類時, 輸入變量不是像素的光譜, 而是地物的光譜分布。

除了傳統基于像素光譜的遙感分類, 面向對象遙感分類也是一種較常見的分類方法[12-14]。 需要說明的是, 面向對象遙感分類中的對象和遙感推斷中的對象的語義往往不同。 面向對象遙感分類中的對象往往指具有一定相似性的像素集合, 生成這些集合的過程實質上是一種數字圖像分割, 而基于圖像像素亮度的統計分布特征(閾值)正是圖像分割的一種常見方法[15], 如大津法[16]。 與之類似, 遙感統計推斷也是基于像素亮度統計分布特征, 但是其目標不是分類和分割, 而是在已經確定地物對象或類別的前提下(即在一定的空間范圍內), 利用像素亮度的統計分布特征去推斷引起某種亮度統計分布的地物屬性的分布特征。

遙感統計推斷與傳統統計推斷也有所不同, 傳統統計推斷往往是根據樣本集的統計特征去推斷總體的統計特征[17-18], 例如去一塊農田或一個水庫采集了30個土壤或水體樣本, 根據這30個樣本去推斷整個農田或水庫的狀況, 這就是傳統統計推斷。 傳統統計推斷在成本和效率上顯然更低, 而通過遙感統計推斷, 則可以節省成本和提高效率。

表1總結了上面陳述的關于遙感分類、 反演和推斷的特點和異同之處。 從分類, 到推斷, 再到反演, 是一個由“質”到“量”逐漸變化的過程, 推斷作為一種介于定性和定量之間的方法, 有其自身的特點和優勢, 相比于分類, 推斷在定性上的優勢是更“準”、 更“?！? 相比于反演, 推斷在定量上的優勢是更“快”、 更“精”, 因此可以作為傳統遙感分類與反演方法的一種有益補充。

表1 遙感分類、 反演與推斷的特征與比較

遙感統計推斷是一個在遙感理論和應用方面都有較多擴展空間的新方向。 本文旨在提出這一方法, 并做一些理論和應用上的初探工作, 引入與之相關的一些概念、 原理和技術, 展現一些簡單示例, 以期拋磚引玉, 推動該領域進一步地深入研究和發展。

1 基礎概念和理論

1.1 基本原理

對于一個給定地物(研究對象), 假設其空間范圍是A,A在圖像上包含的像素數量是n, 傳感器觀測到的原始DN范圍是R=[0, 255]。 我們有兩種從概率統計的角度看待圖像的方式： (1)將該地物所有可能生成圖像的集合看作是一個樣本空間S, 即對地物的觀測是一種隨機實驗, 將某一時刻獲取的一副圖像作為這個集合的一個元素, 即該樣本空間的一個樣本點s, 對地物的一次觀測是一次隨機事件。 可知樣本空間中的元素數量是一個極大的數值256n, 而每幅圖像出現的概率基本都是1/256n, 即我們幾乎不可能從兩副圖像上觀測到每個像素的DN值都完全一樣的地物A; (2)將所有可能的DN值看作一個樣本空間, 即該空間中含有256種元素, 對地物的一次觀測(成像)相當于對樣本空間做了n次抽樣(一種隨機實驗)。 將每次抽樣當作一個隨機事件, 該事件的結果即抽出的DN值, 因此n次抽樣的結果為x1,x2, …,xn∈R, 這些結果組成了一個隨機變量X=(x1,x2, …,xn), 也就是說, 在觀測對象的空間范圍A內所有像素的DN值就是一個隨機變量, 其值域是[0, 255]的離散整數, 且值域中每個值被抽中的次數(概率)是不一樣的, 例如對于DN=0的元素, 其被抽中(在A范圍被觀測到)的次數是m0, 則其對應的概率就是P(0)=m0/n, 這樣256個元素的概率就形成了一種概率密度為p(X)=(P(0),P(1), …,P(255))的離散型分布。

遙感推斷的理論基礎涉及較多光學(統計光學)和統計(推斷方法)方面的內容。 統計光學是現代光學的一個分支, 其主要是應用概率統計來分析光學問題, 研究光場的統計性質(Goodman, 2015)。 傳統統計光學中研究的光場主要包括自發輻射和受激輻射產生的熱光和激光等, 其光場本身含有較為明顯或細微的統計特征。 遙感成像本質上也是一個光學過程, 不過一般可以認為或假定, 在一定范圍內(即感興趣的地物對象區域), 大氣頂端的太陽入射光場是一種均勻分布, 當入射光經與大氣和地物的吸收和反射后, 在傳感器上成像所形成的二維(如單波段圖像)或三維(如多光譜或高光譜數據立方體)圖像光場的分布會因此改變, 從而具有某些統計性質。 在給定時刻t, 傳感器對空間范圍A內的地物成像, 則圖像像元ai接收的信號Li為

Li=f(λi,Θi,pi,Ψai,Ψsi)

(1)

式(1)中,f是輻射傳輸方程,λ是入射信號光譜,Θ是入射和觀測幾何,p是極化,Ψa和Ψs分別是大氣和地表特性的參數集(梁順林等, 2019)。 假定在給定的時空范圍(A,t)內,Θi和pi不隨ai變化, 并且如果進一步假設, 入射光(大氣頂層)和大氣參數(在A范圍內的水平方向)也是均勻分布, 則入射輻射在近地表的入射值λsurf也是一個不隨ai變化的均勻分布, 則有

Li=f(λsurf,Θ,p,Ψa,Ψsi)

(2)

所以, 在A范圍內的n個Li則形成一個變量L=(L0,L1, …,Ln)。 將Li歸一化到[0, 255]的范圍內, 則L即等于上面從抽樣角度描述過的地物隨機變量X。

如果將地表特性參數Ψs也看作是一個隨機變量, 則從式(2)可以看出,L是Ψs的函數, 且以λsurf,Θ,p和Ψa為函數參數。 由隨機變量函數的分布特性可知, 假設Ψs的概率密度為gS(s), 且輻射傳輸函數f是單調函數, 則L的概率密度函數gL(l)為

(3)

若f是非單調函數, 則有

(4)

式(4)中,s1,s2, …,sn是方程l=f(sk)的所有n個解。

事實上, 求解式(3)或式(4)極為困難, 一是因為輻射傳輸過程f是一個極其復雜的迭代過程, 我們無法獲知其反函數f-1的解析表達, 二是因為參與輻射傳輸的地表特性參數Ψs可能不止一個,Ψs是一個多元隨機變量, 每個參數的概率分布都會對最終L的分布產生影響。

在遙感推斷中, 我們已知的是Li的分布gL(l), 需要獲知的是Ψs的分布gS(s), 因此假設式(2)有反函數Ψs=h(Li)=f-1(λsurf,Θ,p,Ψa,Li), 則式(3)和式(4)分別變為

(5)

(6)

其中l1,l2, …,ln是方程Ψs=h(lk)的所有n個解。

由于函數h一般沒有解析形式且L受多變量Ψs決定, 式(5)和式(6)也較難求解。 但是如果針對某個單變量地物參數Ψs, 如果能獲取其與Li關系近似的解析表達, 則式(5)和式(6)可以相應地解出。 事實上,Ψs與Li之間的近似關系h往往可以表達為反演所獲取的算法或模型, 然而, 正向我們前面所指出的, 反演算法通常是基于幾何光學的, 其在建模過程中, 受實地采樣數據數量、 質量以及地物參數復雜性的影響, 具有一定的不確定性和誤差; 再者, 很多表達近似關系h的反演模型也不是解析可導的, 因此根據式(5)或式(6)計算地物的概率分布函數gS往往有較大的不確定性或無法求解。

遙感統計推斷需要做的就是在不建立反演模型的前提下, 直接建立gS與gL之間的關系, 以求達到準確且快捷地獲取地物參數的統計分布特征。

1.2 遙感推斷的適用對象

原則上遙感推斷可以應用于任意空間范圍A, 但是在實際運用時, 對A有一定的要求和限制：

(1)A必須包含具有統計意義的像素數量, 如前所述, 如果A僅僅包含一個或幾個像素, 則無法對A做出統計描述或統計意義不顯著;

(2)A應當對應于圖像中有共同大類屬性的地物對象, 如土壤、 水體、 植被等, 顯然不能將這些對象的光譜混為一體進行推斷, 就如同一種反演算法往往是估算同一類地物的參數;

(3)在范圍A內, 研究對象的屬性往往具有某些統計分布特征, 如水體渾濁度、 土壤濕度、 農作物健康度等等。 有些人工地物, 如瀝青、 水泥或石料鋪設的路面、 場地等等, 由于面積較小、 材質均勻, 其屬性一般情況下不會呈現較明顯的統計分布特征;

(4)A一般需要有明確的、 有意義的邊界, 即A對應于一個對象的整體, 而不是其局部。 例如A是一片森林或一個湖泊, 如果A僅僅是森林或湖泊的局部, 則對其推斷出的結果缺乏指示研究對象整體分布特征的意義;

(5)A一般應是連通的區域, 即通常不將兩個不連通的區域當作一個區域進行推斷, 例如一副影像上有兩個湖泊, 即使名稱相同, 但如果水體不連通, 其某種屬性的分布特征可能有較大差異, 將其混合后推斷出一個單一的結果對兩個湖泊都沒有指示意義。 相反, 有時候兩個湖泊可能名稱不同, 但是如果存在寬闊的連通水道, 水體交換頻繁, 此時反而須將其作為一個整體進行推斷。

1.3 地物對象光譜的統計描述

對地物光譜做統計描述, 首先要面臨的問題是對何種級別的光譜信息進行統計描述。 這里的光譜信息可以是原始DN、 經過標定后TOA幅亮度、 經過大氣校正后的地表反射率, 對于水體來說, 還可以是除去水面反射后的離水幅亮度或遙感反射率Rrs(remote sensing reflectance)。 事實上, 這些光譜信息在不同級別上進行轉換時, 其對應的概率統計特征也會相應地改變, 因此最適合于遙感推斷的光譜信息自然應該是剛離開地物且未受大氣及其他干擾要素(如水面反射)影響, 如地表反射率或水體的遙感反射率。 但是, 這并不意味著其他級別的光譜信息不能用于遙感推斷, 正如遙感分類可以基于地表反射率, 也可以基于未經大氣校正的TOA反射率。 如前所述, 大氣或其他干擾要素對區域A內地物光譜概率分布的影響可能是各向均勻的, 因此這些影響往往可以被統計消除。 例如, 如果區域A內DN呈現一個正態分布(μ,σ2), 則DN×gain+offset標定后的變量仍然是一個正態分布, 只不過其均值和方差變為gain×μ+offset和(gain×σ)2。

在地物對象的空間范圍A內, 并在給定的波長λ上, 可以對其所包含的n個像素的光學屬性L(如DN、 地表反射率R、 離水輻亮度Lw等)做出各種統計描述, 常見的參數有最大值、 最小值、 中值、 均值、 方差等等。 更為可視化的描述是根據L的數值繪制統計直方圖, 從而可以觀察和獲取更多的分布特征與細節。 這一直方圖可以看作是隨機變量L的概率分布函數的近似, 我們稱之為對象的“光譜概率分布”(spectral possibility distribution, SPD)。 基于SPD, 可以進一步計算諸如偏度(skewness)、 峰度(kurtosis)等統計變量。

當用統計直方圖描述SPD特征時, 還需要有幾點注意事項： (1)在計算SPD之前, 有時候需要從區域A對應的n個像素中減去部分(例如1%～4%)有可能引入不確定性的像素。 例如在水質遙感中,A中一些極亮或極暗的反射率是由水面的白帽或深色船舶引起, 此時可以將n中對應反射率最大和最小2%的像素排除后, 對剩余96%的像素進行統計描述; (2)當SPD用于定性/分類分析時, 有時還需要將其歸一化, 以便比較不同研究對象之間的相對差異。 例如兩個湖體在某個波段上的SPD看似不同的正態分布, 但是歸一化后(即分布函數經某種放大、 縮小和平移后), 兩者可能是幾乎相同的; (3)除了描述單波段的SPD外, 還可以根據所研究地物參數的光學特征, 描述并計算不同波段的組合信息, 如NDVI已經可以直接指示研究區域的植被分布信息; 研究水體中的葉綠素時, 可以借鑒其兩波段反演算法[19], 描述水體紅和近紅外波段比的統計分布特征。

1.4 入射分布的正向傳輸

入射分布的正向傳輸研究是指已知地表參數在區域A內的統計分布特征, 估算由此引起的地表反射光譜的分布特征。 這其中須假設A范圍內的大氣在水平方向上是各向同性的, 因此入射光在與地物交互之前是一個常數(均勻分布), 而需要估算的則是剛離地的、 未受大氣影響的反射光譜的分布, 如離水幅亮度或遙感反射率, 即通過大氣校正后從圖像中獲取的地表的光譜信息。

研究入射分布正向傳輸的主要方法是蒙特卡特模擬或基于輻射傳輸的數值模擬, 如涉及水體輻射傳輸的Hydrolight模型。 首先模擬各種地物參數可能的分布特征, 例如, 葉綠素呈某種正態分布, 黃色物質和非藻類顆粒物呈某種指數分布, 然后將這些參數輸入模擬模型, 模型運算后輸出與之對應的光譜分布, 從而有助于了解如果觀測到某種光譜分布, 則其理論上可能對應于何種地物參數的統計分布形式。 注意, 把這種對應應用于推斷時需謹慎, 因為它們可能并不是一一對應, 即不同的地物參數分布特征可能產生相同的光譜分布。 類似情況在反演中也常有出現, 特別是在水色遙感領域, 水色組分(葉綠素、 黃色物質、 非藻類顆粒物)濃度的不同組合可能會產生相同的水面光譜。 同理, 水色組分濃度分布的不同組合也可能會產生相同的水面光譜分布。

我們不但要研究地物參數由于概率分布函數形式(例如正態分布、 對數正態分布、 指數分布)的不同所引起的光譜分布形式的不同, 也要研究分布形式相同, 但分布參數不同所引起的光譜分布形式和參數的差異, 而這些分布參數正是我們感興趣的, 是具有定量指示作用的關鍵參數。 例如, 當水體中懸浮顆粒呈現正態分布, 但分布的均值和方差發生變化時, 光譜是否也是正態分布? 如果是正態分布, 分布的均值和方差如何相應地變化。

當光譜與地物參數之間存在輻射傳輸的簡化形式時, 如對于水體有水面反射率R∝bb/a, 其中a和bb分別是水體的總吸收與總后向散射系數, 則可以利用式(5)或式(6), 理論上計算與分析入射分布的正向傳輸過程, 當理論計算過于復雜時, 也可以采用蒙特卡洛模擬法。

1.5 地物分布的反向推斷過程

地物分布的反向推斷過程是遙感統計推斷的常規過程, 即通過觀測到的光譜信息的統計特征和分布去推斷地物參數的統計特征和分布, 其具體包含兩類不同的推斷情形： (1)統計特征推斷： 僅僅推斷地物參數的若干統計特征參數, 如均值、 方差、 最大最小值等; (2)統計分布推斷： 推斷地物參數在研究范圍內完整的概率分布形式。 可以看出, 上面的情形(2)的結果是包含情形(1)的, 即如果獲知了某地物參數的概率密度分布, 易知其相應的統計特征參數; 而反之則不行, 因為不同的概率分布形式可能具有相同的統計特征參數。

通過光譜的概率分布直接推斷地物參數的概率分布, 即上面描述的情形(2), 是較為困難的過程, 其框架性、 系統性的方法有待于后續的深入研究。 情形(1)的推斷則相對簡單, 本文介紹一種簡易方法, 即根據實測數據樣本集的重采樣法(bootstrap method)。 假設實地采樣的樣本集為S, 從S中抽取一個子集s1, 然后分別計算s1中光譜和地物參數的某個感興趣的統計參數(如均值), 例如計算水體遙感反射率在670 nm的均值Rrs1和顆粒物濃度C1; 重復這一過程, 抽取子集s2, 計算Rrs2和C2; 將此過程重復n次, 形成一個光譜均值集合{Rrs1,Rrs2, …,Rrsn}以及與之對應的濃度均值集合{C1,C2, …,Cn}, 然后根據這兩個集合擬合出遙感反射率均值和顆粒物濃度均值之間的推斷模型。 地物參數的其他統計特征也可以用相同方法進行推斷建模。 因此, 通過bootstrap方法建立的是一個推斷算法集, 推斷不同的統計參數需要使用算法集中的不同算法。 將推斷算法應用于遙感影像時, 只需要計算地物空間范圍A內的光譜統計參數(如均值), 然后代入其相對應的推斷模型, 即可獲得地物參數的某種統計特征(如顆粒物濃度均值)。

利用bootstrap推斷還有一些注意事項： (1)從樣本集S中抽取子集s時, 應適當選擇s的大小, 使得其包含的光譜和地物參數的統計特征既有一定的統計意義, 也有一定的差異性。 抽樣方式可以采用隨機重采樣, 使子集S和s的地物和光譜參數具有類似的統計分布特征; 如果S集包含的樣本數較多, 也可以進行“重分布”采樣, 使得S和s中的地物和光譜參數具有不同的統計分布特征。 例如樣本集S中地物和光譜參數呈現均勻分布, 我們從中抽取具有不同特征的正態、 對數正態以及指數分布的子集s, 從而可以研究當地物參數呈此類分布時, 相對應的光譜分布的特征, 使得生成的推斷模型更為準確, 適用性更廣; (2)遙感推斷依據的光譜統計分布信息不僅僅是單波段的, 還可以是多個波段的, 或者使用不同的波段組合或模型, 即根據地物參數的光學性質, 利用以往反演建模的先驗知識, 確定合適的光譜分布信息。 例如, 葉綠素反演的兩波段模型是基于紅綠波段的比值, 則可以依據紅綠波段比的分布去推斷葉綠素的分布; (3)很多常見的統計概率密度分布形式是由若干分布參數(如均值和方差)所確定的, 如果假設地物參數是已知的某種常見分布函數, 當通過bootstrap方法推斷出該分布的參數后, 實際上該地物的完整分布就確定了, 即我們完成的是一個參數推斷。 即使該地物的具體分布函數未知, 獲取其均值、 方差等關鍵統計信息, 也可以大致描述其完整分布的主要信息。

實際的遙感推斷還要面對一個問題： 如何結合實地采樣獲取的信息進行推斷? 面向推斷的采樣和面向反演與分類的采樣有一些異同, 并且針對不同的地物特征也要采取不同的策略。 如果將衛星觀測到的地物參數在空間范圍A內的光譜信息看作是一個總體(population), 則實地采樣獲取的僅僅是從總體中抽樣的一個樣本集。 這個樣本集如果足夠大且沒有偏差, 例如采樣點密集且在整個區域A內隨機分布, 則可以利用傳統統計推斷方法, 根據地物參數及光譜在實測樣本集中的分布推斷出其在總體A中的各自分布。 由此可以看出, 對于面向推斷的采樣, 往往要求樣本集對于總體是無偏的, 而對于面向遙感分類和反演的采樣, 往往要求光譜特征的差異與多樣性, 采樣前須預估地物參數(類別)和光譜可能存在明顯差異的子區域, 并在該區域內采集若干代表性樣本, 但是這些代表性樣本的數量往往不和其對應子區域的面積成正比, 因此根據樣本集去推斷總體的統計特征與分布往往會引起一定的偏差。

另外, 由于遙感影像是對地物的瞬間觀測, 而實地采樣不可能在短時間內采集較大空間范圍內的數據, 因此地物參數的時變特征以及空間范圍A也會對樣本集的代表性有一定影響。 一些較為固定, 短時內(例如三五天內)變化不明顯的地物(如農田作物)參數, 每天采集一些子區域, 三五天內的數據則可以覆蓋整個研究區域, 然后將各個子區域的數據合并成整個區域的無偏樣本集; 但是對于水體等流動性、 時變性較快的地物參數, 按照同樣方式采集的樣本集則可能存在一定偏差, 將其用于推斷時可能會有一定的不確定性和偏差。

2 應用初探

2.1 國內湖泊、 海灣水體的SPD特征

湖泊有明確的邊界與范圍, 是遙感統計推斷理想的研究對象, 其水質參數, 如溫度、 濁度、 葉綠素和黃色物質濃度等, 也往往具有一定的分布特征。 以往的一些觀測表明, 一個生態系統中生物種群數量或個體年齡(生長或活動狀態)有時呈正態或對數正態分布, 而葉綠素正是水體中浮游植物數量和生長狀態的指示參數。 這些水質參數如果直接或間接影響水體內在光學屬性的分布, 就會相應地改變遙感觀測到的表觀光學屬性的分布。

我們利用Landsat-8衛星影像, 提取分析了國內130多個主要水體的SPD特征。 這些水體主要以湖泊為主, 也包含一些相對封閉的海岸帶水體, 如膠州灣、 三都澳等。 Landsat-8原始數據從美國地質調查局USGS網站獲取, 經由大氣校正模塊ACOLITE處理后, 可以直接提取水體區域在Landsat-8五個可見光/近紅外波段(443、 483、 551、 655和865 nm)的遙感反射率Rrs, 除去所有像素中反射率最大及最小2%的像素, 然后計算剩余96%像素在每個波段的統計直方圖。 為了便于在相同范圍內比較不同水體的SPD特征, 設置bin數量為20, 并將Rrs歸一化到0-1, 即

(7)

式(7)中,Rrs_norm是歸一化后的遙感反射率,Rrs是像素處的觀測值,Rrs_max和Rrs_min是水體區域中遙感反射率的最大與最小值, 其歸一化后的Rrs_norm值分別對應1和0。

圖2顯示了我國一些主要水體基于Landsat-8影像的SPD特征。

圖2 國內主要水體在Landsat-8五個可見光/近紅外波段上的光譜概率分布

(1)單一波長SPD的統計分布特征。 單一波長SPD不僅包含一些常見的概率分布形式, 如正態分布[圖2(b)賽里木湖865 nm、 圖2(f)洱海443 nm]、 對數正態分布[圖2(d)扎陵湖655 nm、 圖2(b)賽里木湖463 nm]、 指數分布[圖2(c)移山湖865 nm、 圖2(o)冬給措納湖655 nm]、 均勻分布[圖2(i)巢湖483 nm、 圖2(h)查干湖655 nm]等, 而且包含很多有一定特征, 但無法歸結為已知的常見概率分布形式, 如星云湖551 nm的平臺狀分布[圖2(p)]、 日居錯443、 483、 551 nm的雙峰分布[圖2(k)]、 巢湖865 nm的階梯狀分布[圖2(i)]、 微山湖483 nm的鞍狀分布[圖2(u)]、 鄱陽湖551 nm的不規則分布等[圖2(a)]。 有些SPD雖然具有相近的概率分布形式, 但是其統計參數卻存在一定差異, 如都是正態或對數正態分布, 其參數均值μ和方差σ2不同(如賽里木湖、 洱海、 冬給措納湖443 nm); 同是指數分布, 其參數斜率λ不同(扎布耶錯和青海湖865 nm); 一些分布具有一定對稱性(如三都澳655 nm), 另一些分布則具有不同偏度(如興凱湖655 nm); 一些分布較為平緩(如膠州灣483 nm), 另一些分布則具有較大峰度(如巢湖483 nm)。

(2)單一水體內不同波長SPD之間的組合特征。 單一水體內不同波長SPD之間的關系具有相似性和多樣性。 在有些水體中, 不同波長的SPD曲線在形狀(概率分布函數)和位置(分布函數的統計參數)都比較接近, 如在移山湖[圖2(c)]和扎布耶錯[圖2(e)]中, 五條曲線都呈指數分布, 且形狀和位置基本重合; 而在有些水體中, SPD的形狀較為相似, 但位置存在差異, 例如在賽里木湖中[圖2(b)], 各個波段的SPD都呈正態或對數正態分布, 但是從443到483 nm再到551 nm, 隨著波長增加, SPD曲線的均值左移, 峰度增加, 但是從551到655 nm再到865 nm, SPD曲線的均值轉而右移, 峰度降低, 形成一個以551、 483、 443、 655、 865 nm為順序, 峰度依次降低, 均值逐漸右移的對數正態SPD組合, 在峽山水庫、 陽宗海、 駱馬湖等水體中也觀測到類似的SPD組合。 還有一種常見的特征是五條SPD曲線形成不同的組, 組內SPD的形狀和位置都相同, 而組外則不同。 例如在興凱湖、 洪澤湖和三都澳中, SPD在443、 483、 551和655 nm的曲線較為相近, 形成一組, 偏度偏右, 而近紅外865 nm的SPD則單獨形成另一組, 形狀雖與其他波長的SPD相近, 但是偏度偏左。 這種組合在洱海中也類似, 兩組偏度都是偏左, 但是形狀卻不同, 865 nm呈指數分布, 而另一組中的4個SPD呈正態分布。 其他一些組合特征還包括： 扎陵湖中, 藍波段(443和483 nm)一組, 綠波段(551 nm)一組, 而紅和近紅外波段(655和865 nm)一組; 巢湖和納木錯中, 藍綠一組, 紅/近紅外一組; 冬給措納湖中, 綠和近紅外形成一組, 而紅和兩個藍波段各自形成一組; 青海湖中, 紅波段單獨形成一組, 其余四個波段形成一組; 在日居錯和星云湖中, 各個波段的SPD都存在較大差異;

(3)同一水體在不同時間的SPD特征。 同一水體在不同時間或季節, 也具有不同的SPD特征, 例如太湖在冬季[2017年12月5日, 圖2(r)]和春季[2014年3月16日, 圖2(s)], 微山湖在冬季[2013年11月29日, 圖2(u)]和春季[2014年3月14日, 圖2(t)]呈現的SPD都不盡相同。 太湖在冬季藍綠形成一組, 峰度較高, 右偏, 紅和近紅外分別形成一組。 在春季, 綠色從藍綠組中分離出來, 峰度有所降低, 單獨形成一組; 而藍紅三個波段的峰度降低十分明顯, 整個曲線較為平展, 而近紅外波段的變化則不明顯。 微山湖從2013年冬季到2014年春季, 各個波段的SPD特征都發生了較大變化, 冬季分布大致呈現一種兩端突出, 中間平緩的鞍狀曲線, 到了春季后, 突出的兩端向中間移動, 且峰度降低。 太湖春季3月中旬SPD綠波段的分離可能是由于氣溫變暖后, 浮游植物和藻類活動增加, 濃度增高, 引起水體IOP發生變化, 從而造成SPD相應改變; 而微山湖未觀測到類似變化, 可能是因為微山湖緯度較高, 3月中旬水溫仍較低, 水生生物還處于非活躍期。

從上面結果可以看出, 不同水體或同一水體在不同時間單一波段SPD的分布特征以及不同SPD的組合特征不僅具有一定的相似性, 也具有一定的多樣性。 我們推測, 這些相似性與多樣性與水體及其周邊的物理、 化學、 生物、 環境、 生態等屬性或因素有一定關聯, 如水體溫度、 鹽度、 深度、 集水區所處的氣候帶、 海拔、 植被、 土地利用與覆蓋及人類活動強度等。 SPD特征可用于潛在地指示這些水質相關屬性和因素, 無須推斷出這些屬性和因素統計分布的量化特征, 即前述所說的基于SPD的遙感推斷可以像基于光譜的遙感分類一樣對研究對象進行分類和定性分析。 當然, 如果進一步推斷出這些屬性的量化統計分布特征, 則可以對研究對象進行更深入更有針對性的定量分析。

本節僅僅是針對Landsat-8觀測到的水體SPD做一些直觀的描述性的初步分析, 以顯示遙感統計推斷潛在的定性分析功能。 我們建議將來收集更多數據, 對SPD的統計特征做一些定量描述, 量化其一些基本參數, 如均值、 方差、 偏度、 峰度等等, 使得基于這些參數的分類和定性分析結果更為準確。

2.2 西湖懸浮顆粒物分布的初步推斷

我們利用在西湖采集的數據, 建立了一個簡單的統計推斷模型, 用于推斷西湖懸浮顆粒物在412 nm的吸收系數aNAP(412)的主要統計參數——均值。 2017年—2019年期間, 我們在西湖各湖區采樣十余次, 共采集可用樣本108個, 在每個樣本點測量了光譜、 透明度、 懸浮顆粒物和CDOM吸收系數等水質參數[20]。

樣本點基本覆蓋西湖的整個湖區, 實測數據大致呈指數衰減分布[見圖3(a)], 推斷建模采用前面介紹的Bootstrap方法, 即(1)從108個樣本中隨機抽取50個樣本, 形成一個重采樣樣本集s, 計算s中aNAP(412)、 5個Landsat-8單波段以及這個5個單波段兩兩形成的10個波段比的均值; (2)重復過程(1)500次, 形成一個含有500條記錄的建模集M; (3)基于建模集M, 在aNAP(412)均值和15個光譜均值之間分別做簡單的線性回歸, 形成15個推斷模型; (4)根據推斷模型的統計參數(如相關系數R2), 從15個模型選擇一個最佳模型作為最終的推斷模型。

圖3 西湖懸浮顆粒物吸收系數aNAP(412)均值的推斷建模與驗證

基于實測數據獲取的西湖懸浮顆粒物吸收系數aNAP(412)均值的Landsat-8推斷模型為

(8)

式(8)中, 下標mean表示該變量在整個湖區像素的均值, 模型的R2=0.44, 見圖3(b)。 我們仍然用實測數據進行驗證： 從108個樣本中隨機抽取k(30≤k≤80)個樣, 由這k個樣直接計算aNAP(412)的均值, 即aNAP(412)均值的真實值, 再根據k個樣的光譜均值和推斷模型式(8)計算aNAP(412)的推斷值, 重復這一過程500次, 則可根據這500對aNAP(412)的真實值和推斷值之間的誤差獲知推斷模型的準確性。 驗證結果顯示, 推斷模型即式(8)的相對誤差僅7.1%, 見圖3(c)。

建模結果顯示, 波段比B2/B4是用于推斷aNAP(412)均值的最佳變量, 其次是B1/B4(R2=0.36), 而其他14個單一波段或波段比變量所擬合的模型的推斷準確性較差。 另外我們也測試了基于Bootstrap法, 用單波段或波段比的方差去推斷aNAP(412)的方差, 但是發現建模效果很差,R2的最大值僅有0.027(以B4作為輸入變量), 說明方差這一統計參數無法通過簡單的單波段或波段比方差進行線性擬合推斷。 由于我們只測試了單波段和波段比這兩類變量, 在推斷模型的擬合上也只用了最簡單的線性回歸, 所以推斷結果并不是特別理想, 將來如果能借鑒相關的反演模型, 使用更為復雜的變量組合或擬合方式也許能進一步提高均值和方差的推斷準確度。

影響推斷準確性的還包括用于建模的實測數據的質量。 在本研究中, 實測數據的空間范圍雖然覆蓋了整個湖區, 但是卻是在不同時間采集的。 另外, 西湖整個湖區包含中間大湖、 北面的北里湖, 西面的西里湖和茅家埠, 西南的小南湖等幾個分湖區, 每個分湖區的水體雖然相通, 但是橋下連通的水道較為狹窄, 加之各個湖區上游船的活動強度也不同, 因此各個分湖區的水體有一定的獨立性, 如果將其當作一個整體進行推斷, 也可能會降低推斷的準確性。

在對實測數據進行bootstrap重采樣時, 由于實測數據數量有限, 我們也沒有進行如1.5節所描述的帶分布特征的重采樣, 而只是進行隨機采樣, 使得每個重采樣后的子集中數據的分布都與原始數據集的分布相類似, 這樣可能會造成得到的推斷模型準確度降低且具有一定的局限性。

3 結 論

遙感統計推斷的一個特點是避免大量反演計算又能獲取地物參數關鍵的定量信息, 即先快速獲取地物的光譜信息的統計特征, 再輸入推斷模型中, 獲取地物其他參數的統計特征。 注意, 我們不能將推斷模型和反演模型相混淆, 雖然有時候兩者在形式上有些類似。 對象的光譜統計信息不能輸入到反演模型中, 而像素的光譜簽名信息也不能輸入到推斷模型中。 反演模型往往是基于像素光譜, 如果我們根據區域A內所有像素生成一個“均值”或“中值”光譜, 然后將該光譜輸入反演模型運行一次, 雖然速度也相對較快, 但這樣只能獲取地物參數的一個“均值”或“中值”, 無法獲知此參數在研究區域內的分布狀況(如方差), 且反演模型對所謂的“均值”或“中值”光譜是否能正確響應, 獲取的均值或中值與地物參數的實際統計值是否吻合, 或者其反演的準確度是否比推斷的準確度高, 也是待進一步研究的問題。

本文提出了遙感統計推斷的基本概念和原理, 即一種基于輻射傳輸、 統計光學和概率分布變換的, 兼顧定性與定量分析功能的遙感數據分析技術與方法。 與傳統的遙感分類和反演相比, 遙感推斷具有“快”、 “準”、 “精”、 “?！钡葍瀯?。 基于Landsat-8提取的水體SPD特征表明, 研究對象的光譜概率分布是其非光學屬性和環境因子的一種指示與表征, 很多水體的SPD既有多樣性, 也有相似性。 基于西湖實測數據的簡單推斷模型則表明, 遙感統計推斷能夠較為快速地、 準確地估算出研究對象的一些關鍵的概要性的定量信息。

通過初步的理論探索和應用測試表明, 基于對象光譜統計分布特征的定性分析與定量推斷是遙感科學領域一個有意義的新研究課題, 仍有很多基礎理論工作以及一些具體的實用技術與應用實踐還未涉及, 希望感興趣的讀者能深入研究, 在將來對這一課題做進一步的拓展與完善。