周毅 劉瑤 田淑芳

(1 中國(guó)地質(zhì)大學(xué)(北京) 地球科學(xué)與資源學(xué)院，北京 100083)(2 自然資源部國(guó)土衛(wèi)星遙感應(yīng)用中心，北京 100048)

透明度作為反映水體清澈/渾濁程度的水質(zhì)參數(shù)，被廣泛應(yīng)用于湖泊學(xué)和海洋學(xué)研究中，對(duì)評(píng)估水體內(nèi)水生植物多樣性、生產(chǎn)力和水體營(yíng)養(yǎng)程度等方面具有重要意義[1]。過(guò)去幾十年來(lái)，受人口增長(zhǎng)、經(jīng)濟(jì)發(fā)展中不合理的開發(fā)利用及污染等因素影響，我國(guó)許多內(nèi)陸水體(湖泊、河流和水庫(kù))呈富營(yíng)養(yǎng)化狀態(tài)，水體質(zhì)量急劇下降，帶來(lái)水體生態(tài)系統(tǒng)的結(jié)構(gòu)破壞與功能退化，對(duì)人民生活及國(guó)家經(jīng)濟(jì)造成巨大損失和不良影響[2]。近年來(lái)，國(guó)家為改善水生態(tài)環(huán)境，開展了重點(diǎn)流域保護(hù)修復(fù)、重點(diǎn)工程水質(zhì)保障等一系列工作。有必要定期開展內(nèi)陸水體透明度等水質(zhì)參數(shù)的監(jiān)測(cè)工作，服務(wù)于生態(tài)修復(fù)措施的水環(huán)境改善效果評(píng)價(jià)，以及為后續(xù)水資源的保護(hù)和管理提供決策數(shù)據(jù)支持。透明度通常定義為塞氏盤(Secchi Disk)在水中垂直下放至剛剛在視野中消失的深度。由于塞氏盤獲取成本低且操作簡(jiǎn)便，使用塞氏盤進(jìn)行水體透明度測(cè)量已有上百年的歷史。然而，通過(guò)布設(shè)地面監(jiān)測(cè)點(diǎn)或者船載測(cè)量的方式難以實(shí)現(xiàn)對(duì)大區(qū)域和多湖泊的透明度重復(fù)觀測(cè)。相比而言，能夠周期性穩(wěn)定獲取和實(shí)現(xiàn)空間覆蓋的遙感衛(wèi)星，成為了開展長(zhǎng)期大規(guī)模水體透明度反演的重要數(shù)據(jù)源。

國(guó)內(nèi)外已有眾多學(xué)者基于遙感數(shù)據(jù)開展水體透明度反演研究，常見方法有經(jīng)驗(yàn)/半經(jīng)驗(yàn)?zāi)Ｐ秃桶敕治瞿Ｐ?。文獻(xiàn)[3]基于陸地衛(wèi)星-5(Landsat-5)專題成像儀(TM)數(shù)據(jù)的藍(lán)波段、紅藍(lán)波段反射率比值，使用相關(guān)系數(shù)矩陣、多元逐步回歸等方法完成了1973～1998年間美國(guó)大型湖泊的透明度評(píng)價(jià)研究；文獻(xiàn)[4]使用Landsat-5和陸地衛(wèi)星-7(Landsat-7)增強(qiáng)型專題成像儀(ETM+)數(shù)據(jù)的藍(lán)紅波段，通過(guò)建立經(jīng)驗(yàn)?zāi)Ｐ蛯?duì)1990～2010年間美國(guó)緬因州湖泊開展透明度反演研究；文獻(xiàn)[5-6]建立了一種基于半分析算法(Quasi Analytical Algorithm，QAA)的水體透明度反演模型，可應(yīng)用在包括中分辨率成像光譜儀(MODIS)、中等分辨率成像頻譜儀(MERIS)、陸地衛(wèi)星-8(Landsat-8)等衛(wèi)星的多種傳感器上實(shí)現(xiàn)水體透明度的反演。國(guó)內(nèi)方面，文獻(xiàn)[7]使用Landsat-5中紅藍(lán)波段建立自然對(duì)數(shù)與透明度自然對(duì)數(shù)的線性模型，對(duì)2004年6個(gè)時(shí)期鄱陽(yáng)湖的水體開展透明度反演，取得了較好精度；文獻(xiàn)[8]應(yīng)用Landsat-8 陸地成像儀(OLI)數(shù)據(jù)的紅、綠波段建立對(duì)數(shù)模型，反演了2013年的東平湖水體透明度，取得了19.77%的總體相對(duì)誤差；文獻(xiàn)[9]基于透明度和水體漫衰減系數(shù)關(guān)系建立半分析模型，并通過(guò)MODIS、近海高光譜成像儀(HICO)數(shù)據(jù)對(duì)近岸水體透明度進(jìn)行反演。

目前，水體透明度反演研究多以多光譜數(shù)據(jù)為數(shù)據(jù)源，采用星載高光譜數(shù)據(jù)開展研究相對(duì)較少。相較于多光譜數(shù)據(jù)，高光譜數(shù)據(jù)具有更高的光譜分辨率，在用于水體透明度反演時(shí)，能夠提供更豐富的波段選擇，且具備更精準(zhǔn)的光譜特征，因此在水質(zhì)參數(shù)反演精度的提升方面具有優(yōu)勢(shì)和潛力。

我國(guó)于2018年5月9日成功發(fā)射高分5號(hào)(GF-5)衛(wèi)星，其搭載新一代高光譜相機(jī)(Advanced Hyper Spectral Imager，AHSI)的光譜覆蓋范圍為0.40～2.50 μm，空間分辨率為30 m，可見光近紅外(VNIR)與短波紅外(SWIR)分別有150和180個(gè)波段，光譜分辨率分別為5 nm和10 nm。此外，資源一號(hào)02D衛(wèi)星(又稱為5米光學(xué)業(yè)務(wù)衛(wèi)星)也于2019年9月12日成功升空。該星搭載的AHSI，具有與GF-5衛(wèi)星相同的光譜覆蓋范圍和空間分辨率，可見光近紅外波段有76個(gè)通道(光譜分辨率為10 nm)，短波紅外波段包括90個(gè)通道(光譜分辨率為20 nm)。我國(guó)后續(xù)還將發(fā)射搭載高光譜相機(jī)的衛(wèi)星，通過(guò)多星組網(wǎng)，未來(lái)有望能夠基于星載高光譜數(shù)據(jù)開展內(nèi)陸水體水質(zhì)監(jiān)測(cè)。

本文以北京市官?gòu)d水庫(kù)為研究區(qū)，以資源一號(hào)02D衛(wèi)星高光譜影像為數(shù)據(jù)源，應(yīng)用半經(jīng)驗(yàn)和半分析模型方法反演研究區(qū)的水體透明度。通過(guò)對(duì)結(jié)果精度進(jìn)行評(píng)價(jià)，分析資源一號(hào)02D衛(wèi)星高光譜數(shù)據(jù)在反演水體透明度方面的表現(xiàn)與應(yīng)用潛力。

1 材料與方法

1.1 研究區(qū)概況

官?gòu)d水庫(kù)位于北京西北部永定河上游(見圖1)，水庫(kù)面積約130 km2。1970年以來(lái)，受入庫(kù)流量下降、污染負(fù)荷增加等影響，官?gòu)d水庫(kù)水質(zhì)不斷惡化，于1997年退出北京飲用水水源地之列；后經(jīng)治理，水庫(kù)在2007年再度成為北京備用飲用水水源地。近年來(lái)，隨永定河綜合治理和生態(tài)修復(fù)推進(jìn)，官?gòu)d水庫(kù)水質(zhì)逐年好轉(zhuǎn)[10]。

圖1 研究區(qū)示意圖Fig.1 Map of study area

1.2 數(shù)據(jù)及處理

2020年8月13日在官?gòu)d水庫(kù)開展了與資源一號(hào)02D衛(wèi)星同步的地表試驗(yàn)，獲取了地表實(shí)地的光譜數(shù)據(jù)和水質(zhì)參數(shù)(葉綠素a、懸浮物濃度、透明度等)。高光譜遙感數(shù)據(jù)的成像時(shí)間為2020年8月13日上午11時(shí)21分，地面試驗(yàn)則在當(dāng)天10～14時(shí)內(nèi)完成，以保證星地?cái)?shù)據(jù)獲取的一致性。本次試驗(yàn)在研究區(qū)范圍內(nèi)布設(shè)了18個(gè)均勻分布的點(diǎn)位(見圖2)，其中實(shí)測(cè)水體透明度范圍為0.4～0.7 m，均值為0.58 m。

圖2 水體透明度實(shí)測(cè)位置圖Fig.2 Location map of measured water transparency

對(duì)獲取的星載高光譜數(shù)據(jù)，進(jìn)行了輻射定標(biāo)、大氣校正、水體范圍提取和遙感反射率(Rrs)計(jì)算等一系列預(yù)處理工作，具體的處理過(guò)程如圖3所示。

圖3 高光譜數(shù)據(jù)處理過(guò)程Fig.3 Hyperspectral data process flow

1)輻射定標(biāo)

通過(guò)ENVI5.3軟件，將原始影像的VNIR數(shù)據(jù)(76個(gè)波段)和SWIR數(shù)據(jù)(90個(gè)波段)合并為一個(gè)文件，應(yīng)用Apply Gain and Offset工具導(dǎo)入絕對(duì)定標(biāo)參數(shù)后進(jìn)行輻射定標(biāo)。

2)大氣校正

FLAASH大氣校正模型是目前在多/高光譜影像中使用效果最好的大氣校正工具之一，對(duì)于地物波譜信息可以高保真地恢復(fù)，簡(jiǎn)潔快速地獲取比較準(zhǔn)確的物理模型參數(shù)[11]。因此本文在ENVI5.3中應(yīng)用FLAASH模塊對(duì)影像數(shù)據(jù)進(jìn)行大氣校正。綜合考慮數(shù)據(jù)特點(diǎn)、成像時(shí)間及研究區(qū)位置，本文在FLAASH模塊中選擇中緯度夏季(Mid-Latitude Summer)大氣模型，水體反演波段設(shè)置為940 nm，并選用鄉(xiāng)村(Rural)氣溶膠模型，以及2-波段(考夫曼-泰倫)(2-Band(K-T))氣溶膠反演方法。

3)水體范圍提取

本文選擇歸一化水體指數(shù)(Normalized Difference Water Index，NDWI)方法進(jìn)行水體范圍提取[12]，其公式為

(1)

式中：ρGreen、ρNIR分別代表綠波段和近紅外波段的地表反射率。

式(1)最初以Landsat TM數(shù)據(jù)進(jìn)行計(jì)算，而資源一號(hào)02D衛(wèi)星高光譜數(shù)據(jù)信噪比相對(duì)較低，選擇單一波段計(jì)算容易受到噪聲的影響[13]。因此，本文參照Landsat TM數(shù)據(jù)的綠光(0.52～0.60 μm)和近紅外波段范圍(0.76～0.90 μm)，通過(guò)波段合并，分別計(jì)算資源一號(hào)02D衛(wèi)星數(shù)據(jù)在Landsat TM這兩個(gè)波段范圍內(nèi)的地表反射率平均值，用于歸一化水體指數(shù)的計(jì)算。在此基礎(chǔ)上，運(yùn)用直方圖雙峰法(Histogram Bimodal Method，HBM)分離水體和其他類型地物，并將小于1000像元(約0.9 km2)的微小水體消除，只保留主體水域，從而實(shí)現(xiàn)影像中水體范圍的提取。

4)遙感反射率計(jì)算

遙感反射率計(jì)算公式為

(2)

式中：λ為波長(zhǎng)，ρ(λ)為經(jīng)過(guò)FLAASH大氣校正后的地表反射率，ρsky(λ)為天空光反射率，Lu(λ)為水面上行輻亮度，Ld(λ)為水面下行輻亮度，Lsky(λ)為天空光輻亮度。

即使是反射率相對(duì)較高的渾濁水體，在短波紅外的反射率也趨于0，因此通常認(rèn)為圖像中短波紅外的反射率是來(lái)自天空光的影響，將短波紅外的反射率用于天空光的去除[14]。但是，由于AHSI的輻射定標(biāo)系數(shù)是基于地表亮目標(biāo)計(jì)算得到，而水體反射率相對(duì)較低，造成水體像元在短波紅外波段噪聲相對(duì)較大，不適宜用于天空光的糾正。因此，本文采用忽略天空光的影響，以ρ(λ)/π作為遙感反射率的近似值。盡管采用這樣近似計(jì)算的結(jié)果可能略大于實(shí)際的遙感反射率，但由于通常認(rèn)為天空光的影響與波長(zhǎng)無(wú)關(guān)，因此，近似計(jì)算的遙感反射率與真實(shí)反射率具有一致的光譜形狀?；诖?，本文采用波段比值的經(jīng)驗(yàn)?zāi)Ｐ?，和半分析模型進(jìn)行透明度反演。這兩類方法對(duì)于透明度的反演，主要取決于光譜的波形特征，受遙感反射率的量級(jí)影響較小[15]。

1.3 水體透明度反演方法

本文分別運(yùn)用半經(jīng)驗(yàn)?zāi)Ｐ秃桶敕治瞿Ｐ头椒ㄟM(jìn)行官?gòu)d水庫(kù)水體透明度反演。

1.3.1 半經(jīng)驗(yàn)?zāi)Ｐ?/p>

根據(jù)文獻(xiàn)[16]研究成果，半經(jīng)驗(yàn)?zāi)Ｐ头椒ㄖ饕ㄟ^(guò)對(duì)地面同步點(diǎn)位的影像光譜和地表測(cè)量的透明度進(jìn)行相關(guān)性分析，選擇合適的波段或波段組合進(jìn)行建模。本文通過(guò)分析資源一號(hào)02D衛(wèi)星圖像和地表采集的水體光譜特征，選擇中心波長(zhǎng)為653 nm的波段，采用對(duì)數(shù)回歸模型(見圖4)得到水體透明度反演公式為

Zsd=-0.438×ln (Rrs(653))-1.335 8

(3)

式中：Rrs(653)為中心波長(zhǎng)為653 nm波段處的遙感反射率。

圖4 對(duì)數(shù)回歸模型Fig.4 Logarithmic regression model

1.3.2 半分析模型

通過(guò)半分析模型反演水體透明度主要步驟為：①估算固有光學(xué)量(IOPs)，即總吸收系數(shù)a和后向散射系數(shù)bb；②基于a和bb反演漫衰減系數(shù)(Kd，m-1)；③用Kd和Rrs確定Zsd。

本文分別運(yùn)用第5版及第6版半分析算法(QAAv5、QAAv6)[17]計(jì)算總吸收系數(shù)a和后向散射系數(shù)bb；并采用更新后的半分析模型[18]對(duì)漫衰減系數(shù)Kd進(jìn)行反演

m1×(1-m2×e-m3×a(λ))bb(λ)

(4)

式中：θs為太陽(yáng)天頂角，bbw(λ)為純水的后向散射系數(shù)，m0-3與γ分別取0.005、4.26、0.52、10.8和0.265的固定值，a(λ)為總吸收系數(shù)，bb(λ)為后向散射系數(shù)。

對(duì)于水體透明度的反演，本文通過(guò)公式[19]進(jìn)行計(jì)算

(5)

式中：Kd(tr)表示水體在可見光譜范圍(410～665 nm)上的漫衰減系數(shù)，Rrs(tr)表示該波長(zhǎng)下相應(yīng)的遙感反射率。

1.4 評(píng)價(jià)方法

本文通過(guò)建立線性回歸模型來(lái)分析反演結(jié)果與實(shí)測(cè)值之間的相關(guān)性，并計(jì)算絕對(duì)誤差(AE)、相對(duì)誤差(RE)以及均方根誤差(RMSE)用于對(duì)結(jié)果的精度評(píng)價(jià)為

IAE=|dm-dp|

(6)

(7)

(8)

式中：dm和dp分別表示透明度的實(shí)測(cè)值和反演值；n為樣本數(shù)。

2 結(jié)果與分析

2.1 不同反演結(jié)果與實(shí)測(cè)值對(duì)比

本文分別運(yùn)用1.3節(jié)中的方法采用資源一號(hào)02D衛(wèi)星數(shù)據(jù)對(duì)研究區(qū)水體透明度進(jìn)行計(jì)算，并提取了與實(shí)測(cè)點(diǎn)位匹配的相應(yīng)反演值。應(yīng)用線性回歸模型對(duì)反演結(jié)果與實(shí)測(cè)值進(jìn)行了相關(guān)性分析并得到相應(yīng)的相關(guān)系數(shù)值(見圖5)。

圖5 線性回歸分析圖Fig.5 Linear regression analysis results

此外，運(yùn)用1.4節(jié)中各誤差分析方法計(jì)算得到各反演結(jié)果的誤差分析情況(見表1)。

表1 反演結(jié)果誤差分析表Table 1 Error analysis of inversion results

根據(jù)線性回歸分析結(jié)果可以看出，運(yùn)用半經(jīng)驗(yàn)?zāi)Ｐ头椒ǚ囱莸玫降乃w透明度結(jié)果與實(shí)測(cè)結(jié)果相關(guān)性最高，R2為0.473 3；而基于QAAv5與QAAv6的半分析模型方法反演結(jié)果相關(guān)性較低，其中基于QAAv5計(jì)算的反演結(jié)果與實(shí)測(cè)值相關(guān)性最低，R2為0.419 3。圖5(b)和(c)中，由于QAAv5與QAAv6算法主要是針對(duì)一類大洋水體提出，是以550～560 nm作為標(biāo)準(zhǔn)波長(zhǎng)；將其應(yīng)用于渾濁水體，會(huì)造成對(duì)總吸收系數(shù)和后向散射系數(shù)的低估，進(jìn)而低估漫衰減系數(shù)，導(dǎo)致反演的水體透明度偏高[20-21]，因此擬合線低于1∶1線。通過(guò)表1可知，半經(jīng)驗(yàn)?zāi)Ｐ腿〉玫姆囱萁Y(jié)果，相較2種半分析模型方法結(jié)果，整體精度較高，絕對(duì)誤差平均值為0.05 m，平均相對(duì)誤差為9.5%，而RMSE僅為0.07 m；兩種半分析模型得到的反演結(jié)果精度相對(duì)較低，其中，基于QAAv5計(jì)算的反演結(jié)果精度較差，最大相對(duì)誤差達(dá)到了174.0%，其平均相對(duì)誤差為94.7%，RMSE為0.55 m。而基于QAAv6反演結(jié)果平均相對(duì)誤差及RMSE分別為31.4%、0.19 m，相較QAAv5結(jié)果精度稍高，主要原因是QAAv6算法對(duì)渾濁水體的情況進(jìn)行了適當(dāng)?shù)膬?yōu)化，對(duì)渾濁的內(nèi)陸水體的適用性稍有提升。

上述3種方法反演的透明度結(jié)果與實(shí)測(cè)結(jié)果的R2都比較低，這主要是由于本研究區(qū)的同步點(diǎn)位數(shù)量有限；另一方面，反演結(jié)果的精度分析可以看到，表明該方法在應(yīng)用于官?gòu)d水庫(kù)透明度反演中具有較好的效果。

2.2 官?gòu)d水庫(kù)水體透明度分析

通過(guò)上述分析，運(yùn)用半經(jīng)驗(yàn)?zāi)Ｐ头椒ǚ囱莸玫降乃w透明度具有較高精度，因此本文基于該方法得到官?gòu)d水庫(kù)水體透明度(Zsd)狀態(tài)圖(見圖6)。

圖6 官?gòu)d水庫(kù)Zsd狀態(tài)圖Fig.6 Zsd state of Guanting reservoir

從圖6可以看出，官?gòu)d水庫(kù)水體透明度范圍在0～1 m之間，相比東北部水體，西南部水體透明度更高，說(shuō)明該區(qū)域水體質(zhì)量更好。該結(jié)果與實(shí)地測(cè)量得到的水質(zhì)狀態(tài)基本一致。

3 結(jié)論

本文以北京市官?gòu)d水庫(kù)為研究區(qū)，采用資源一號(hào)02D衛(wèi)星高光譜數(shù)據(jù)，應(yīng)用半經(jīng)驗(yàn)?zāi)Ｐ?、半分析模型方法進(jìn)行研究區(qū)的水體透明度反演。通過(guò)反演結(jié)果與實(shí)地測(cè)量結(jié)果的相關(guān)性分析和精度評(píng)價(jià)，得到主要結(jié)論如下：

(1)本文分別使用半經(jīng)驗(yàn)?zāi)Ｐ?、基于QAAv5與QAAv6的半分析模型反演得到水體透明度，通過(guò)結(jié)果與實(shí)測(cè)數(shù)據(jù)相關(guān)性及精度評(píng)價(jià)，半經(jīng)驗(yàn)?zāi)Ｐ头椒ǖ玫降姆囱葜蹬c實(shí)測(cè)值相關(guān)性最高，相關(guān)系數(shù)為0.473 3，且結(jié)果精度也最高，平均絕對(duì)誤差、平均相對(duì)誤差及RMSE分別為0.05 m、9.5%、0.07 m；其次為基于QAAv6半分析模型方法的反演結(jié)果，平均絕對(duì)誤差、平均相對(duì)誤差及RMSE分別為0.17 m、31.4%、0.19 m；而基于QAAv5的半分析模型反演結(jié)果精度較差。

(2)基于資源一號(hào)02D衛(wèi)星的高光譜數(shù)據(jù)，可實(shí)現(xiàn)水體范圍的提取并用于內(nèi)陸渾濁水體的透明度反演。相比較常規(guī)的水色衛(wèi)星，資源一號(hào)02D衛(wèi)星具有更高的空間分辨率，可應(yīng)用于中小型湖庫(kù)的水質(zhì)監(jiān)測(cè)。下一步的研究是驗(yàn)證該數(shù)據(jù)應(yīng)用于更多不同光學(xué)特性水體中的透明度反演能力。