劉 冉，徐 鵬，劉心怡，黃 瑞，朱孟蘭

(中國電建集團貴陽勘測設(shè)計研究院有限公司，貴州 貴陽 550081)

1 引言

當(dāng)前，大部分內(nèi)陸水體受到不同程度的污染，特別是在經(jīng)濟較發(fā)達的地區(qū)水體污染更為嚴重，這不僅影響了居民的居住環(huán)境，而且制約了城市進一步發(fā)展。因此，有效的水質(zhì)監(jiān)測和治理成為了亟待解決的問題，其中監(jiān)測是基礎(chǔ)。傳統(tǒng)的監(jiān)測方法是定期實地取水樣進行分析，但是獲取的只是離散的數(shù)據(jù)，時效性也不高。遙感技術(shù)自20世紀70年代開始應(yīng)用于水體，其具有監(jiān)測范圍廣、速度快、 成本低和便于進行長期動態(tài)監(jiān)測的優(yōu)點，為實時高效監(jiān)測水質(zhì)提供了一種行之有效的方法，促進了水污染的治理和生態(tài)環(huán)境的保護，對城市、社會的可持續(xù)發(fā)展具有重要的意義。

2 研究現(xiàn)狀

隨著衛(wèi)星遙感技術(shù)的不斷發(fā)展，水質(zhì)遙感反演由海洋向內(nèi)陸水體、由定性向定量發(fā)展，可反演的水質(zhì)參數(shù)主要包括葉綠素a濃度、懸浮物濃度、CODM和濁度等，反演的模型主要有經(jīng)驗?zāi)Ｐ?、半?jīng)驗?zāi)Ｐ秃臀锢砟Ｐ蚚1]。

2.1 經(jīng)驗?zāi)Ｐ?/h3>

經(jīng)驗?zāi)Ｐ椭饕腔谛l(wèi)星多光譜波段或波段組合和水質(zhì)參數(shù)之間的相關(guān)性分析，然后選擇與水質(zhì)參數(shù)相關(guān)性最大的波段或波段組合建立經(jīng)驗?zāi)Ｐ?，廣泛應(yīng)用于內(nèi)陸水質(zhì)遙感監(jiān)測。如：Steven M.Kloiber等[2]根據(jù)TM影像和地面同步觀測數(shù)據(jù)，利用TM1和TM3的比率建立一個三參回歸模型來預(yù)測SDT(Secchi disk transparency)。余豐寧等[3]利用TM各波段數(shù)據(jù)和太湖同步觀測數(shù)據(jù)進行相關(guān)性分析，選擇最佳波段建立了葉綠素含量的遙感定量模型。

經(jīng)驗?zāi)Ｐ碗m然在特定條件、區(qū)域能取得比較理想的結(jié)果，但必須以大量實驗為基礎(chǔ)，而且物理學(xué)概念、機理不清，模型適應(yīng)性差。

2.2 半經(jīng)驗?zāi)Ｐ?/h3>

半經(jīng)驗?zāi)Ｐ蚚4]是針對高光譜數(shù)據(jù)源的水質(zhì)遙感監(jiān)測方法，對實測水體光譜或者衛(wèi)星高光譜進行相關(guān)處理，選取和水體水質(zhì)參數(shù)相關(guān)性最好的特征波段或變量，然后和水質(zhì)參數(shù)建立反演模型，該模型具有一定的物理意義。

Koponen等[5]通過對航空高光譜影像和MERIS影像的多個波段進行組合或比值后和葉綠素a濃度建立各種模型，根據(jù)反演結(jié)果和測點數(shù)據(jù)計算RMSE，選出葉綠素a濃度最佳的反演模型；疏小舟等[6]在太湖地區(qū)采用高分辨率地物光譜儀實地測量湖水的反射光譜曲線，并同時采樣分析葉綠素等水質(zhì)參數(shù)，研究發(fā)現(xiàn)內(nèi)陸水體葉綠素a濃度>5 μg/L時，光譜比R705nm/R675nm、R700nm附近反射峰位置適合用于葉綠素a濃度的遙感監(jiān)測。

與經(jīng)驗?zāi)Ｐ拖啾?，半?jīng)驗?zāi)Ｐ途哂幸欢ǖ奈锢硪饬x，模型精度有所提高。

2.3 分析模型

分析模型是基于太陽光在水體中的傳輸過程以及水體中的葉綠素、懸浮物、污染物和水分子等組分的吸收、散射特性對傳輸過程的影響，構(gòu)建水體反射光譜與污染濃度之間函數(shù)關(guān)系的物理模型。

Hoogenboom等[7]利用水體各組分的相互作用，構(gòu)建一種快速的、可行的反演水體葉綠素濃度的矩陣模型；李四海等[8]利用SeaWIFS數(shù)據(jù)，通過現(xiàn)場光譜測試及準同步采樣，建立了遙感參數(shù)與含沙量的定量反演模式；鄧孺孺等[9]在研究污染水體反射光譜成像的物理過程的基礎(chǔ)上, 建立了葉綠素不敏感波段遙感數(shù)據(jù)像元反射率與污染物濃度的函數(shù)關(guān)系, 用多波段數(shù)據(jù)進行求解，并采用TM數(shù)據(jù)運用該方法對珠江口海域水污染進行定量分析，取得了很好的效果。

分析模型具有清晰的物理概念、機理，不需要像經(jīng)驗、半經(jīng)驗?zāi)Ｐ鸵蕾嚧罅康耐綄崪y數(shù)據(jù)，適用性強，反演精度高，是目前定量遙感常采用的方法。

3 模型原理及推理

3.1 輻射傳輸過程分析

太陽光與水體相互作用的輻射傳輸過程如圖1所示。

圖1 太陽光在水中輻射傳輸過程

從圖1中可以看出，到達水面的入射光—包括太陽直射光和天空散射光(天空光)在水面發(fā)生鏡面反射，一部分反射光向上傳輸進入傳感器；其余的入射光發(fā)生折射、透射進入到水中，其中一部分與水體中的組分如：水分子、懸浮泥沙、CODM等發(fā)生吸收散射，部分散射光經(jīng)過水體二次衰減后，穿出水面進入傳感器；在水深不大的前提下，另一部分入射光繼續(xù)向下傳輸?shù)竭_底部，被水底反射，反射光同樣經(jīng)過水體二次衰減后穿過水體，進入傳感器，這部分反射光與水體的散射光共同組成水中光或叫離水反射輻射[10]。

3.2 水體輻射傳輸方程

基于上面的分析可以得知水域之上衛(wèi)星傳感器接收到的行星反射率主要由鏡面反射比Rm和水體出水反射率Rwo兩部分組成，即：

Rw=Rwo+Rm

(1)

水面鏡面反射帶有少量的水體本身的信息，它的強度主要與表面粗糙度有關(guān)，但是這部分反射光幾乎不攜帶任何水質(zhì)信息，對于水質(zhì)遙感來說是干擾信息，應(yīng)該予以剔除。鏡面反射比Rm與波長無關(guān)，水在紅外波段為強吸收，Rwo近似等于零，所以Landsat8數(shù)據(jù)的紅外波段OLI6(1.560～1.660μm)的反射率約等于鏡面反射比Rm，每個波段減去OLI6即可消除鏡面反射比Rm即得到水體出水反射率Rwo：

Rwo=Rw-Rm=Rw-ROLI6

(2)

而水體出水反射率Rwo又由水中散射率Rws和水底對水體反射率的貢獻Rg組成，即：

Rwo=Rws+Rg

(3)

又由具有朗伯反射地物的反射率與輻亮度的關(guān)系可得：

(4)

式(1)～(2)中，Lwo為水體的離水輻亮度；Eo為大氣層外太陽平均輻照度；θ為太陽天頂角；

而水體離水輻亮度Lwo由兩部分構(gòu)成，即水中散射光向上的部分Lws和水底反射光Lg；

如果只考慮一次散射，并且假設(shè)水體中的組分只有水分子、懸浮泥沙和污染物，則水深h處的薄層水的散射對離水輻亮度Lwo的貢獻為：

(5)

假設(shè)水深為H，對式(5)從0到H積分，得

(6)

將式(6)代入式(4)，得到水中的散射率Rw s：

(7)

同理，根據(jù)水體反射的物理機制，可以求出水底反射光對水體離水輻亮度的貢獻為：

(8)

式(8)中Rb為水底底質(zhì)的反射率；將式(8)代入(4)，可得到水底對水體反射率的貢獻Rg：

Rg=Rbe-μ(α+β)H

(9)

將式(7)和式(9)代入(3)，得到體的出水反射率R_wo：

(10)

3.3 污染物和懸浮泥沙提取

由于河口水體比較渾濁，污染也比較嚴重，一般水域(水深>2m)到達水底的反射光極其微弱，可以忽略，其效果相當(dāng)于水深極大，即H→∞，則式(10)可簡化為[9]：

(11)

式(11)為本文所使用的最終水質(zhì)反演模型。

4 模型在珠江口水質(zhì)參數(shù)定量提取中的應(yīng)用

本文所采用的遙感數(shù)據(jù)是2016年2月7日廣州地區(qū)Landsat8的OLI傳感器的多光譜數(shù)據(jù)，處理流程如圖2所示。

圖2 處理流程

4.1 輻射定標

(12)

式(12)中，Qcal為影像以16位量化的亮度值；Mρ為波段λ的反射率調(diào)整因子；Aρ為波段λ的反射率調(diào)整參數(shù)，Mρ和Aρ可從頭文件中讀取，其數(shù)值如表1所示。

表1 各波段定標系數(shù)

(13)

式(13)中θ為太陽天頂角，與太陽高度角互余，太陽高度角可從頭文件中讀取。

4.2 大氣校正

大氣輻射校正是進行定量遙感的先決條件, 其目的是將衛(wèi)星遙感數(shù)據(jù)轉(zhuǎn)換為地物真實信息的表觀反射率，對后續(xù)的定量反演有著非常重要的影響。一般可通過三種途徑進行大氣散射校正，即輻射傳遞方程式計算法、野外波譜測度回歸分析法及多波段圖像的對比分析法[12]。

本文選擇基于大氣傳輸物理機制的暗像元法進行大氣校正。該方法以山區(qū)陰影部分的植被作為黑體，先假設(shè)陰影區(qū)植被的反射率為0，估算出大氣散射對程輻射的貢獻，再利用迭代法對程輻射進行校正[13]。

傳感器接收到的像元反射率R(λ)為：

(14)

式中，T'為入射方向大氣透過率；ω為大氣散射系數(shù)；Rg為地面反射率；T為垂向上大氣透過率；Ρ(Θ)為散射相函數(shù)，Θ為散射角；

由式(14)可得陰影區(qū)植被行星反射率Rv為：

(15)

根據(jù)式(15)對所選取的一個暗像元進行計算，所求得的校正參數(shù)如表2。

表2 第一次計算后的校正參數(shù)

表3 第二次計算后的校正參數(shù)

由于第一次計算時，我們先假設(shè)在陰影區(qū)植被的反射率為0，把在陰影區(qū)植被處所讀得的數(shù)值直接作為程輻射反射率處理，沒有考慮植被對下行的天空光的反射，會存在一定誤差，因此需要進行一步訂正[13]。利用第一次計算的參數(shù)進行第二次迭代，計算得到的校正參數(shù)如表3所示。

為了檢驗大氣校正的效果，從影像上選取了植被這一典型地物樣本，讀取植被某點的反射率值，并繪制大氣校正前后反射率曲線的對比圖，如圖3所示。

波段

從圖中可以看出，校正前植被反射率曲線與標準植被反射率曲線有很大偏差，特別是在反射峰附近偏差更大，而大氣校正后植被反射率有了明顯的改正，反射率曲線則與標準植被反射率曲線相似，這也說明了暗像元法是一種比較有效的大氣校正方法。

4.3 污染物和懸浮泥沙濃度的提取

可見光波段能夠很好地反映水體信息，對水質(zhì)信息比較敏感。OLI2波段為藍光波段，對水的穿透力最大，包含更多的水下細節(jié)和水深信息，不利于水質(zhì)信息的提?。欢鳲LI3波段為綠光波段，對水體具有一定的穿透能力，可反映水體渾濁度，對水體污染的研究效果好，OLI4波段為紅光波段，是信息量最豐富的一個波段，對水中的懸浮泥沙反應(yīng)敏感，因此本文選用OLI3和OLI4這兩個波段進行計算。

本文研究的對象是水域的水質(zhì)信息，所以對遙感圖像進行水陸分離，剔除各波段的陸地信息部分，提取水體部分，可加快運算速度。由于水體在近紅外波段強烈吸收，為低反射，可通過觀察圖像在近紅外波段水陸交界處的實際反射率值，選取一個合適的值作為水陸分界的閾值，再根據(jù)植被在近紅外波段的反射率高于紅光波段這一特征將植被陰影區(qū)域分離出來[14]。

所使用的水體各組分的散射和吸收參數(shù)如表4所示。

表4 水質(zhì)參數(shù)

將表4水體各組分的參數(shù)代入式(11)，建立一個二元一次方程組，可求解出污染物濃度和懸浮泥沙濃度的解析解，再利用PCI Geomatica V9.0圖像處理軟件建模對影像進行處理，得到珠江口水域的污染物和懸浮泥沙濃度分布圖，最后使用ArcMap進行分級渲染，結(jié)果如圖4和圖6所示。

圖4 珠江口污染物濃度

從圖4中可以直觀地看出珠江口西部海域污染物濃度總體比東側(cè)的高，對比圖5的實際影像中，也可以看出西岸水體呈現(xiàn)灰黃色，確實比東岸的渾濁。

但是其中東側(cè)水域中，在三角電廠、東寶河入?？?、深圳機場福永碼頭和深圳灣附近污染物濃度明顯比周圍海域的都高，從影像中也可以觀察到該區(qū)域的水體呈現(xiàn)黑黃色，污染物濃度高。

圖5 珠江口部分水域影像(真彩色合成)

圖6 珠江口懸浮泥沙濃度分布

從圖6中，可以看出珠江口西部海域懸浮泥沙濃度明顯比東側(cè)的高，特別是在南沙十九涌以南的水域，從影像中也可以清晰地看出該區(qū)域水呈泥黃色，比較渾濁；一些河流的入?？谔幩w也比較渾濁，懸浮泥沙濃度較高。

5 結(jié)論

(1)從最終的反演結(jié)果來看，珠江口污染物和懸浮泥沙濃度的分布與實際情況基本相符，西部水域這兩種污染物的濃度大體比東岸的高，但是東岸存在著典型的污染源，這與附近的地區(qū)的生產(chǎn)活動相關(guān)，這就為水污染治理提供了有效的信息，有助于治理方案的制定和實施，這也是水質(zhì)定量遙感所帶來的社會效益和經(jīng)濟效益。

(2)本文使用暗像元法進行大氣校正，假設(shè)大氣是均質(zhì)的，而實際上大氣的結(jié)構(gòu)和成分是很復(fù)雜的，在時空分布上也很不均勻，而且只選取單一暗像元對整幅影像進行校正，會產(chǎn)生一定的誤差。另外最終的結(jié)果也缺乏衛(wèi)星過境時地面同步實測數(shù)據(jù)來進行精度檢驗，明顯還存在不足之處。

(3)總的來說，該模型物理機理清晰，是在分析入射光在水體中的傳輸過程的基礎(chǔ)上建立的，模型中的參數(shù)也比較少，也不像經(jīng)驗半經(jīng)驗?zāi)Ｐ托枰龃罅康膶嶒?，而且有些參?shù)，如清潔水體的吸收、散射系數(shù)等適用性比較高，不需要花費大量人力物力進行重復(fù)測量，可以應(yīng)用于其他數(shù)據(jù)或區(qū)域的反演，模型有較好的適用性，可行性也比較高。