基于SR-FSDAF多源影像融合的內陸水體NH3-N反演研究

馬玉榮，柯美如

(1.鄭州大學 水利與交通學院，河南 鄭州 450001；2.鄭州大學 圖書館，河南 鄭州 450001)

0 引言

生活與工業廢水排放尤其是NH3-N污染物激增,導致了水體發黑發臭和富營養化,嚴重破壞了河道的生態系統,影響城市居民的生活和健康。傳統的NH3-N人工監測方法工作量大、費時費力,難以實現大規模、多時段檢測,而遙感技術的快速發展,為河流大范圍水質監測提供了技術支持。利用遙感影像監測水體NH3-N濃度的基本原理是基于對水體中不同濃度NH3-N的光譜特征進行分析,選取高相關性波段建立反演模型。龔紹綺等[1]通過實驗室NH3-N溶液光譜分析得出氮在404、477 nm波長處相關性最高,對應波長的單波段反演模型的擬合度確定系數(R2)分別為0.976和0.947。段洪濤等[2]利用光譜儀對南湖總氮進行實驗,發現最大相關性出現在443 nm處,反演模型R2為0.75。徐良將等[3]利用太湖區域的高光譜數據,用波段比值對水中氨氮濃度進行反演,結果得到波段比值1 015 nm/528 nm與總氮濃度相關性最高為0.803,R2為0.839,這些相關研究為利用遙感影像監測水中氨氮提供了一定依據。

目前,遙感水質監測的衛星數據有MODIS、MERIS[4]、TM/ETM+、Landsat-8 OLI[5-6]、HJ-1 CCD、GF-1[7]、Sentinel-2[8]等,內陸水體NH3-N監測常用的遙感數據有Landsat-7 ETM+、Lamdsat-8[9-10]、SPOT 5[11]等。Landsat-8影像質量優于之前Landsat系列,空間分辨率為30 m,回訪周期為16 d,是目前水質監測的重要數據源。而MODIS衛星回訪周期為1 d,具有高時間分辨率,但空間分辨率較低,250 m (1～2波段)、500 m (3～7波段)、1 000 m (8～36波段)。Landsat和MODIS影像獲取便利,且歷史影像豐富,因此Landsat系列被廣泛用于內陸水體水質監測,MODIS則多用于海洋監測。

由于大氣條件和重訪周期的影響,滿足定量反演云量要求的中高分辨率衛星時間間隔過大,難以有效獲取河道水體污染物濃度的動態變化趨勢;MODIS等多光譜衛星時間分辨率雖高,但其空間分辨率較低,難以識別較窄的內陸河道。如何能既準確獲取內陸河道信息,又提高反演的時間頻率、加強水質污染物的動態監測,是亟待解決的問題。因此,開展多光譜和高空間分辨率影像的融合模型研究,可提高影像的空間分辨率和時間分辨率,為進一步準確、動態反演內陸河道水體的NH3-N濃度提供了保障。

現有時空融合算法根據計算原理的不同分為5類:基于解混的方法、基于權重函數的方法、基于貝葉斯原理的方法、基于學習的方法和混合型方法。基于解混的方法主要有STDFA(Spatial Temporal Data Fusion Approach)[12],其解混誤差較大,且缺乏地物的類內變化。基于權重函數的方法如STARFM(Spatial and Temporal Adaptive Reflectance Fusion Model)[13],利用圖像信息進行權重分配來估計高分辨像素值,對異質變化無效,且權重函數基于經驗獲得,缺乏遷移性。基于貝葉斯估計理論的方法[14]自定義輸入圖像和融合圖像間的關系,但函數確立過程復雜且異質表現一般。基于學習的方法使用機器學習建模高分辨率和低分辨率影像之間的映射關系,從而預測融合圖像,例如Huang等[15]提出的SPSTFM (Sparse Representation Based Spatio-Temporal Reflectance Fusion Model),其融合結果雖有提升,但學習成本較高、遷移性差,且缺乏光譜原理性支撐。混合型方法整合前4類中2種或多種方法,如王宇恒等[16]將PanSharpening和時空融合進行結合,以獲取時、空、譜互補信息。Zhu等[17]提出的FSDAF(Flexible Spatio-Temporal Data Fusion),將解混、加權函數和空間插值方法結合,減少了圖像的輸入量,對異質變化的預測有所增強,但表現仍不理想,不能滿足于更高精度的變化監測。基于此,提出了改進的時空數據融合(SR-FSDAF)模型。利用SR-FSDAF提高監測的頻次和中低分辨率圖像的利用率,生成高頻融合圖像,并應用于信陽市水質NH3-N濃度監測,對淮河流域信陽段內的內陸河道進行NH3-N濃度的反演與時空分布變化分析。

1 時空融合模型SR-FSDAF

FSDAF由Zhu等[17]于2016年提出,該方法對地物進行分類,利用不同類型地物在MODIS上體現出的變化,粗略得到每類地物的時間變化值,對預測日期的MODIS影像進行薄板樣條采樣(TPS),得到粗略的地物空間變化值。TPS基于空間上相關性對數據進行差值重采樣,能保留圖像的局部變化信息。將2種變化值利用鄰域信息賦予不同權重,減少預測結果的偏差,具有更強的穩定性和空間連續性。但FSDAF對于空間突變的預測表現相對無突變情況仍然下降較多。基于FSDAF的思想,對空間突變預測的主要信息源MODIS影像,使用高效亞像素卷積神經網絡(Efficient Sub-Pixel Convolutional Neural Network, ESPCN)替代TPS,保留更多的紋理信息,從而提高算法對于空間異質變化的預測精度。

1.1 時空融合模型SR-FSDAF構建

基于超分辨率和靈活時空數據融合模型的SR-FSDAF分為6個步驟:基于無監督分類的MODIS像元純度計算、像元時間變化粗估計、時間變化殘差計算、基于ESPCN的超分辨率空間變化預測圖像重構、殘差分布計算以及基于鄰域信息的增強和融合。為方便描述,將高空間、低時間分辨率的Landsat圖像表述為“細分辨率圖像”,將高時間、低空間分辨率圖像MODIS表述為“粗分辨率圖像”,使用t1時刻的Landsat和MODIS影像及t2時刻的MODIS影像預測t2時刻的Landsat影像。算法流程如圖1所示,模型的變量及定義如下:

m——一個MODIS像元內對應的Landsat-8影像的像元數;

(xi,yi)——第i個MODIS像元的下標;

i——MODIS的像元下標;

j——在一個MODIS像元中的Landsat-8像元的下標(j=1,2,…,m);

M1(xi,yi,b),M2(xi,yi,b)——在t1和t2時刻MODIS(xi,yi)位置上第b波段的值;

L1(xij,yij,b),L2(xij,yij,b)——在t1和t2時刻MODIS(xi,yi)位置上第j個Landsat-8像元在第b波段的值;

Pc(xi,yi)——第i個MODIS像元中第c類Landsat-8像元所占的比例;

ΔM(xi,yi,b)——第i個MODIS像元在第b波段上t1—t2時刻的變化;

ΔF(c,b)——Landsat-8影像c類地物第b個波段t1-t2時刻的變化,c為水體、農田、建筑物和林地其中之一。

圖1 SR-FSDAF模型算法流程Fig.1 Flowchart of SR-FSDAF algorithm

(1)基于無監督分類MODIS像元純度計算

MODIS影像的空間分辨率為500 m,Landsat影像空間分辨率為30 m。經過圖像配準,一個MODIS像元大約與16個Landsat像元覆蓋面積相同,基于此對應關系,對t1時刻的Landsat圖像進行K-means地物分類,將分類數設置為水體、農田、建筑物和林地4類。通過Landsat像元分屬哪幾類及各類占比Pc(xi,yi),c=1,2,3,4表示分類之一。

(2)像元時間變化粗估計

基于式(1)計算MODISt1～t2影像反射率變化。

ΔM(xi,yi,b)=Mt2(xi,yi,b)-Mt1(xi,yi,b)。

(1)

根據光譜解混原理[18],ΔM可被表示為式(2),用最小二乘反向求解式(2)可得到ΔF(c,b)。

(2)

(3)時間變化殘差計算

(3)

如果2個時刻之間地物發生了光譜信息的變化,所造成的偏差定義為R(xi,yi,b),計算見式(4),其中n=16。

R(xi,yi,b)=ΔM(xi,yi,b)-

(4)

(4)基于ESPCN超分辨率空間變化預測圖像重構

ESPCN是對低分辨率圖像卷積的高效超分算法,其通過學習低分辨率影像和高分辨率影像之間的特征映射,恢復空間特征信息。ESPCN的輸入為裁剪并重采樣至與高分辨率圖像大小一致的低分辨率圖像和高分辨率圖像的圖像對數據集。ESPCN將卷積神經網絡直接應用于低分辨率圖像,利用亞像素卷積將低分辨率特征映射放大。3次卷積后,得到通道數為r2的與輸入圖像大小一樣的特征圖,再將特征圖像重排列為rH×rW×1的高分辨率圖像,具體網絡訓練參數設置參照文獻[19]。

(5)

(5)殘差分布計算

將殘差R(xi,yi,b)和殘差ESR(xij,yij,b)結合,使用4×4移動窗口計算與中心細分辨率像元地物類別一致的像元數和窗口內像元總數的比值I(xij,yij)作為同質性指標。計算見式(6),其中(xij,yij)為中心像元,當窗口內像元類別一致時,Ik的值為1,否則為0。

(6)

計算2個預測結果之間的差值,為殘差的差異值ESR-TP(xij,yij,b):

(7)

依據同質性指標I(xij,yij)對預測的細分辨率影像的殘差分布進行權重分配,得到權重分布Ew(xij,yij,b),見式(8)。將Ew進行歸一化得到歸一化殘差分布W(xij,yij,b),見式(9)。t2時刻的預測細分辨率像元的殘差EL(xij,yij,b),見式(10)。

Ew(xij,yij,b)=ESR-TP(xij,yij,b)×I(xij,yij)+

R(xi,yi,b)×[1-I(xij,yij)],

(8)

(9)

EL(xij,yij,b)=n×R(xi,yi,b)×W(xij,yij,b)。

(10)

t1～t2的細分辨率像元在波段b的像素變化值ΔL(xij,yij,b)計算如下:

ΔL(xij,yij,b)=EL(xij,yij,b)+ΔF(c,b)。

(11)

(6)基于鄰域信息的增強和融合

使用鄰域信息增加預測穩定性,減少移動窗口的區塊效應。對于在t1時刻細分辨率圖像像元(xij,yij),選取n個同類且與像元在其鄰域內光譜差最小的細像元。第k個細像元與像元的光譜差為Sk,計算如下:

(12)

相似像元對中心像元的貢獻權重Dk與距離呈負相關,距離越遠,貢獻值越少,計算方法見式(13)。歸一化處理后,權重wk的計算方法見式(14):

(13)

(14)

(15)

1.2 模型評價指標

為了比較融合效果的精度和優劣,使用STARFM、FSDAF和SRCNN嵌入FSDAF三種算法作為對比融合模型,分別定性和定量對結果進行分析。定性評價方法通過視覺比較將t2時刻模型融合結果和真實Landsat圖像進行對比,得到局部細節的模擬真實程度和整體的色彩差異性是否過大。定量評價方法使用3類評價指標,分別從融合圖像的結構整體相似性、反射率還原度和光譜保真度進行綜合評價。

結構整體相似性指標為結構相似度(SSIM),被廣泛用于評估2張相似圖像的線性關系強度,計算如下:

(16)

式中:μx、μy是t2時刻Landsat真實圖像和融合的細分辨率圖像的平均值,σx、σy是二者的圖像方差,C1、C2是用于保證指標計算結果的有理性的非零常數。2幅圖像的整體結構越相似,SSIM的值就越接近1。

反射率還原程度的評價指標為均方根誤差(RMSE),其反映像素值的還原程度和細節信息的模擬融合結果,計算如下:

(17)

式中:x(i,j)為真實圖像,y(i,j)為融合圖像。

光譜保真程度使用的評價指標為光譜角(SAM),其將單像素的光譜視為高維向量,計算2幅圖像同一位置像素的光譜向量的向量夾角,夾角越小,像素之間的光譜越相似。具體的夾角計算如下:

(18)

2 NH3-N遙感反演模型

2.1 反演模型

特征波段的準確選擇是獲得高反演精度的基礎。利用水質參數NH3-N的實測數據計算和波段組合的相關系數,并選擇相關性高的波段或波段組合作為建模參數。

傳統的統計回歸模型(TSR)已廣泛用于水質參數反演,采用一次擬合、二次擬合、對數擬合和指數擬合4類回歸函數來構建反演模型,通過對比分析,選擇精度最高的一種模型作為后續反演的模型。

2.2 評價指標

反演模型的準確性通過R2、平均絕對百分比誤差(MAPE)和RMSE進行評估。MAPE計算如下:

(19)

式中:x為實測數據,y為反演模型的估計值,n為檢驗的數據數量。

3 研究區域與數據

3.1 研究區域

信陽市(113°45′E～115°55′E,30°23′N～32°27′N)地處河南省南部,面積18 916 km2,地形南高北低,丘陵平原相間,屬于亞熱帶和溫帶季風氣候過渡帶以及濕潤和半濕潤地區過渡帶,四季分明,雨熱同期。信陽市位于淮河流域,區域內主要支流有浉河、潢河、灌河和史河等。

3.2 實驗數據

(1)遙感影像數據

Landsat-8空間分辨率為30 m,回訪周期為16 d。MODIS空間分辨率為500 m,回訪周期為1 d。MODIS和Landsat具有明顯的時間分辨率與空間分辨率差異,且部分光譜波段重合度高,是時空融合算法常用的融合數據源,使用MODIS和Landsat OLI影像進行實驗。

時空融合模型的實驗數據為信陽市2017年11月8日及11月24日的Landsat圖像與MODIS圖像對,期間圖像沒有發生異質突變。NH3-N反演的遙感數據選取云量少于10%且已做過大氣校正的Landsat-8 OLI和MODIS日地表反射率數據,影像波段等具體信息如表1和表2所示。

表1 Landsat影像信息

表2 MODIS影像信息

對所選取的Landsat-8 OLI和MODIS數據進行預處理,并將MODIS重投影和重采樣至480 m,以便于與空間分辨率為30 m的Landsat影像像素進行匹配和計算。

(2)水質監測數據

實測數據使用2016年9—12月以及2017年2月的信陽區域內信陽市水利局設置的監測站點的實測數據資料,共計773組數據。選取43個監測站點的數據進行實驗,站點覆蓋信陽市境內淮河主要干支流。

4 實驗結果

4.1 不同時空融合模型融合及其對比

將同一時刻的Landsat圖像和MODIS圖像裁剪至16∶1的大小,將Landsat的B2、B3、B4、B5、B6、B7和MODIS的B3、B4、B1、B2、B6、B7波段按照順序合成新的多光譜圖像,并輸入SR-FSDAF模型,得到預測日期的融合圖像,融合結果如圖2所示。可以看出,SR-FSDAF的局部區域的細節地物捕捉和保留能力更好。

圖2 圖像融合結果對比Fig.2 Comparison of image fusion results

SR-FSDAF、STRAFM、SRCNN嵌入模型和FSDAF四種融合方法的RMSE、SSIM、SAM的具體值如表3所示,每個波段的最好融合結果值都采用加粗下劃線進行突出。分析可知:

(1)光譜保真程度評價指標SAM為圖像像素的向量夾角,光譜角越小,像素之間的光譜越相似,圖像保真度越高。SR-FSDAF的SAM得到了最優值3.417,說明有著相對最高的光譜保真度。

(2)在Band1上,FSDAF模型的反射率還原程度指標RMSE最小為0.004 5,SSIM也最小為0.982,獲得了較好的融合效果。而在Band2～4、Band7上,SR-FADAF模型的RMSE和SSIM獲得相對較好的值,和最優值接近。

綜合考慮3個指標,SR-FSDAF模型融合影像精度最高,局部區域細節地物捕捉和保留能力更好,計算結果更優。

表3 圖像融合精度對比

4.2 基于融合圖像的反演

利用信陽市水文局2016年12月5日與2017年1月1日對信陽市重點水功能區布設的43個采樣斷面采集到的86組水樣數據,隨機抽取其中62組作為模型訓練集用于模型的構建,結合2016年12月7日及2016年12月30日的Landsat與MODIS影像生成的融合影像,分析不同波段和波段組合運算值與NH3-N濃度值之間的Pearson相關系數。共對28個波段的計算值進行了相關性分析,并對相關波段進行了不同組合,具體如表4所示。

表4 融合影像和Landsat-8 OLI影像波段組合相關性對比Tab.4 Correlation of different band combinations of fused image and Landsat-8 OLI image

28個波段模型中波段組合(R-G)/(R+G)的r值最高為0.805,且為p<0.01,為顯著相關。因此最終使用波段組合(R-G)/(R+G)作為NH3-N濃度反演模型的x值。

本文以常見的一次函數、二次函數、指數函數和對數函數4種統計模型作為擬合函數進行反演模型構建,擬合結果如圖3～圖6所示。

圖3 一次函數擬合結果Fig.3 Linear function fitting results

圖4 二次函數擬合結果Fig.4 Quadratic function fitting results

圖5 對數函數擬合結果Fig.5 Logarithmic function fitting results

圖6 指數函數擬合結果Fig.6 Exponential function fitting results

從圖中可以看出,二次函數模型不論在濃度高值區還是低值區實際值和函數值擬合度高,整體效果最好,將其作為最終的NH3-N的反演模型,R2為0.664:

(20)

式中:y為NH3-N的估計濃度值。

為驗證反演模型在信陽淮河流域的適用性,使用未參與擬合的另外24組實測站點的數據與反演模型的結果進行對比驗證,利用MAPE與RMSE進行定量精度評定,總體二次回歸的NH3-N的濃度反演模型的平均相對誤差為12.37%,RMSE計算結果為0.065 6,滿足反演需要的精度要求,以此為基礎可展開較為可靠的進一步研究。實測數據和模型預測數據具體如表5所示。

表5 實測數據與模型預測結果對比Tab.5 Comparison between measured data and model prediction results

4.3 不同時期NH3-N濃度整體水平

根據相關調查資料[20],淮河水系每年的7—8月為豐水期,每年12月—次年2月為枯水期,其他月份為平水期。

不同時期(由于豐水期影響,將3—6月和9—11月看作平水期的2個不同階段進行分析,分別表述為平水期Ⅰ和平水期Ⅱ)的水體NH3-N濃度類別占比變化如圖7所示。可以看出,整體上2017年NH3-N濃度呈現略微下降趨勢,并具有明顯的季節特征。枯水期、平水期Ⅰ、豐水期和平水期Ⅱ的NH3-N濃度平均值為0.753、0.706、0.533、0.544 mg/L。整體趨勢為:豐水期<平水期<枯水期。

圖7 不同水期NH3-N濃度水體面積比例Fig.7 Proportion of water area with NH3-N concentration in different water periods

豐水期時期NH3-N濃度對應標定的6類水的水面積平均占比分別為:13.89%、54.26%、24.6%、9.56%、3.73%、0.31%,水體主要成分為Ⅱ類NH3-N濃度水,水質情況相對較好。豐水期的NH3-N濃度相對最低,變化趨勢平緩。

平水期的NH3-N濃度擾動較大,3—6月的平水期前期受枯水期的影響,NH3-N濃度相對較高,9—11月的平水期后期受到豐水期的影響,NH3-N濃度相對較低。

在枯水期,NH3-N濃度最高,高于豐水期和平水期,且變化平緩。枯水期間NH3-N相對濃度最低為12月,主要為Ⅱ類水和Ⅲ類水,占比達到82.77%,從12—次年2月NH3-N濃度逐漸上升,在1月劣Ⅴ類NH3-N濃度水的比例有所增加,出現在竹竿河的部分區域;2月Ⅴ類NH3-N濃度水面積增長較多,主要出現在小潢河和浉河以及灌河下游。枯水期Ⅰ類NH3-N濃度水占總水面平均面積為6.69%,Ⅱ類NH3-N濃度水占整體面積的35.41%,Ⅲ類NH3-N濃度水的平均占比為44.44%,Ⅳ類NH3-N濃度水的平均面積占比為10.15%,Ⅴ類為2.76%,劣Ⅴ類為0.54%主要出現在個別小支流。

淮河流域信陽段的NH3-N濃度在春季濃度最高,主要是冬季枯水期河流水量少且流速緩慢,NH3-N污染物易產生累積造成的,其濃度均值從0.77 mg/L上升至0.81 mg/L。隨著水量增加,NH3-N濃度在5—6月開始下降,并在7—9月的夏季時期因雨季NH3-N濃度全年濃度最低,濃度均值下降至0.52 mg/L。10—12月進入冬季枯水期,NH3-N濃度又緩慢攀升,12月濃度均值為0.69 mg/L。

4.4 重點區域NH3-N濃度分布

依據局部NH3-N濃度差異程度的不同,NH3-N濃度較高的局部小區域主要出現在潢河和浉河的下游。

浉河NH3-N濃度局部長期在1.6～2.3 mg/L,8月時最大值達到了4.5 mg/L,一年內NH3-N濃度變化情況如圖8所示。

圖8 浉河局部區域NH3-N濃度變化Fig.8 Variation of NH3-N concentration in local area of the Shihe River

潢河NH3-N濃度為0.1～1.9 mg/L,小潢河濃度為0.3～2.5 mg/L,且濃度在1.5 mg/L區域基本集中在中下游段,與季節變化趨勢相關弱,一年內NH3-N濃度變化情況如圖9所示。

圖9 潢河局部區域NH3-N濃度變化Fig.9 Variation of NH3-N concentration in local area of the Huanghe River

根據調查,信陽市主要河流周邊大部分為農田,市內主要排污口分布在浉河和潢河的工業生產區和生活區,其面源污染主要受農業農田影響,局部V類及以上的嚴重污染主要受人類工業活動影響。

在豐水期,由于被雨水稀釋,NH3-N整體濃度更低;枯水期NH3-N濃度最高,濃度變化平緩;除去排污的影響,在3—9月的農業生產期,農田肥料等污染會造成水體的NH3-N濃度上升和局部變化。對于浉河、潢河等污染物濃度異于季節規律且常有高濃度NH3-N分布的情況,則主要是由工業生產引起的。

5 結束語

提出了改進的SR-FSDAF模型,通過不同融合模型融合結果對比,發現SR-FSDAF模型比STARFM、SRCNN嵌入模型和FSDAF具有更優的視覺效果和指標結果,圖像RMSE、SSIM和SAM分別為0.009 9(均值)、0.975 (均值)和3.417。基于SR-FSDAF方法,利用Landsat-8 OLI和MODIS進行時空影像融合,并進行了NH3-N濃度的反演及其時空分析,在一定程度上彌補了傳統水質監測方法獲得NH3-N污染物時空變化特征難的問題。

提出的SR-FSDAF方法在提高NH3-N濃度反演頻率上具有一定參考價值,但仍需隨季節性差異和區域性差異進行調整和改進,可為小區域水質監測提供模型與參考。