基于高光譜數(shù)據(jù)的水體溶解氧反演研究

劉元藝 田正航 唐曉燕

（南陽理工學(xué)院信息工程學(xué)院，河南 南陽 473004）

0 引言

隨著社會(huì)的不斷發(fā)展，我國(guó)工業(yè)和城市用水量也越來越大，水資源供需日趨緊張。另外，由于很多城鎮(zhèn)生活污水和工業(yè)廢水是排放到海洋中的，因此會(huì)對(duì)海洋生態(tài)造成污染，這已經(jīng)成為工業(yè)發(fā)展乃至社會(huì)發(fā)展的障礙。

水質(zhì)參數(shù)中溶解氧（DO）的含量是海洋各類生物生存的必要物質(zhì)條件，其可以作為一個(gè)判斷污染程度的重要衡量標(biāo)準(zhǔn)。目前，國(guó)內(nèi)外已經(jīng)逐漸開始使用高光譜遙感技術(shù)進(jìn)行關(guān)于水質(zhì)監(jiān)測(cè)的研究。該技術(shù)使用的衛(wèi)星傳感器不僅具有回訪周期性，而且還具有數(shù)據(jù)獲取容易（與人工檢測(cè)相比）、水質(zhì)反演覆蓋面積較大的優(yōu)點(diǎn)，不需要消耗大量人力、財(cái)力。呂航等人采用高光譜遙感數(shù)據(jù)并結(jié)合該地區(qū)水面同步觀測(cè)數(shù)據(jù)找出了水質(zhì)參數(shù)與單波段反射率、不同波段反射率比以及不同波段反射率差值的相關(guān)性，得出溶解氧的敏感波段組合為B（636 nm）/B（669 nm）和B（636 nm）/B（674 nm）。該文選取香港近海區(qū)域的水體作為研究區(qū)，使用HJ-1A衛(wèi)星HSI數(shù)據(jù)對(duì)香港近海區(qū)域DO濃度進(jìn)行反演研究，根據(jù)高光譜反射率數(shù)據(jù)與研究區(qū)水體DO濃度的相關(guān)性分析，使用線性回歸方法建立高光譜遙感反射率與水體DO濃度的反演模型。

1 數(shù)據(jù)獲取及預(yù)處理

1.1 獲取實(shí)測(cè)數(shù)據(jù)與衛(wèi)星遙感數(shù)據(jù)

研究區(qū)地處北緯22°26′，東經(jīng)114°12′。該文所用實(shí)測(cè)DO數(shù)據(jù)來自香港特別行政區(qū)76個(gè)水質(zhì)監(jiān)測(cè)站，香港環(huán)保署每月收集海水各項(xiàng)水質(zhì)參數(shù)樣本1次。根據(jù)高光譜圖像經(jīng)、緯度范圍的限制，取其中MM區(qū)、SM區(qū)以及VM區(qū)共計(jì)42個(gè)監(jiān)測(cè)站為主要數(shù)據(jù)源。

所需遙感影像的基本要求如下：1） 圖像分辨率高。2）條帶噪聲少。3）云覆蓋量低于5%或?qū)λw區(qū)域無遮擋。4） 盡可能地保證衛(wèi)星圖像拍攝日期與水源數(shù)據(jù)日期一致，以保證地面監(jiān)測(cè)站數(shù)據(jù)和衛(wèi)星高光譜數(shù)據(jù)可以在允許的誤差內(nèi)達(dá)到同步。因此，該文使用中國(guó)資源衛(wèi)星應(yīng)用中心下載2014年1月HJ-1A衛(wèi)星HSI傳感器的2級(jí)高光譜影像產(chǎn)品數(shù)據(jù)。HJ-1A高光譜數(shù)據(jù)由HJ-1A衛(wèi)星搭載的高光譜成像光譜儀（HSI）獲取。HSI傳感器的2級(jí)產(chǎn)品是經(jīng)光譜修復(fù)、輻射校正和系統(tǒng)幾何校正的，并通過ArcGIS將校正后的圖像映射到衛(wèi)星圖像拍攝地點(diǎn)的地圖投影坐標(biāo)下的產(chǎn)品數(shù)據(jù)中。

1.2 衛(wèi)星遙感數(shù)據(jù)處理

衛(wèi)星數(shù)據(jù)的前期處理主要包括以下3個(gè)步驟：1） 光譜曲線輻射亮度值轉(zhuǎn)換。該步驟主要實(shí)現(xiàn)了從遙感影像像元亮度值（值）到絕對(duì)輻射亮度的轉(zhuǎn)化，這一過程對(duì)高光譜數(shù)據(jù)處理來說是比較簡(jiǎn)單的。由于衛(wèi)星拍攝時(shí)通常與地面的距離至少為500 km，拍攝時(shí)得到的真實(shí)地物輻射亮度值通常很低，因此衛(wèi)星數(shù)據(jù)處理中心為了方便對(duì)輻射亮度數(shù)據(jù)進(jìn)行處理、保存，在預(yù)處理前使用放大因子系數(shù)來放大地物輻射率的數(shù)值。該步驟將光譜曲線所顯示的值除以衛(wèi)星處理時(shí)所放大的系數(shù)就可以得到絕對(duì)輻射亮度值，也就是真實(shí)地物輻射亮度值。2） 高光譜數(shù)據(jù)的存儲(chǔ)格式轉(zhuǎn)換。HJ-1A衛(wèi)星原存儲(chǔ)格式為H5文件，經(jīng)ENVI打開后，系統(tǒng)默認(rèn)為BSQ格式存儲(chǔ)，而在進(jìn)行FLAASH大氣校正過程中要求數(shù)據(jù)以BIL或BIP格式存儲(chǔ)，因此原始高光譜數(shù)據(jù)的存儲(chǔ)格式應(yīng)使用Convert Interleave工具進(jìn)行轉(zhuǎn)換（由BSQ格式轉(zhuǎn)換為BIL格式）。3） 波段去除。高光譜影像在衛(wèi)星拍攝時(shí)，由于大氣吸收等因素在光譜圖上會(huì)出現(xiàn)條帶噪聲，從而導(dǎo)致影像在某一波段的質(zhì)量較差，因此該文根據(jù)高光譜圖像各個(gè)波段的圖像質(zhì)量，去除了HJ-1A高光譜影像中的前30個(gè)和后5個(gè)包括條帶噪聲的波段，剩余80個(gè)波段為數(shù)據(jù)保存較為良好的波段。

由于HJ-1A光譜數(shù)據(jù)只是經(jīng)過了初步的幾何校正，為了使采樣點(diǎn)經(jīng)緯度與高光譜衛(wèi)星圖像完全匹配，并提高后續(xù)建模的精確度和準(zhǔn)確性，因此需要對(duì)圖像進(jìn)行幾何精校正，該文采用在地理空間數(shù)據(jù)云（https：//www.gscloud.cn/）網(wǎng)站上下載的同時(shí)期、同地區(qū)Landsat-8的衛(wèi)星影像，由于Landsat-8影像產(chǎn)品級(jí)別較高，且已經(jīng)經(jīng)過幾何精校正，因此可以通過Landsat-8來定標(biāo)校正HJ-1A所拍攝的高光譜影像的經(jīng)、緯度。通過ENVI國(guó)產(chǎn)衛(wèi)星匯總工具包讀取數(shù)據(jù)，使用Image Registration工具均勻選取校正點(diǎn)，通過系統(tǒng)自動(dòng)添加的矯正點(diǎn)去除誤差較大的點(diǎn)，使誤差控制在0.77即可完成幾何精校正工作。

在完成數(shù)據(jù)前期處理工作后，由于HJ-1A高光譜衛(wèi)星所攜帶的傳感器不可避免地受大氣的吸收和散射作用，地表反射數(shù)據(jù)和實(shí)際測(cè)量值之間有一定的差異，因此使用高光譜圖像反演地物的反射光譜信息，需要對(duì)其進(jìn)行大氣校正。該文選用FLAASH大氣校正方法，該方法對(duì)校正參數(shù)的要求比較嚴(yán)格，稍有差錯(cuò)都會(huì)對(duì)最終的結(jié)果造成很大的影響，在進(jìn)行大氣校正前應(yīng)該對(duì)衛(wèi)星數(shù)據(jù)進(jìn)行處理，使其符合以下要求：經(jīng)過輻射定標(biāo)之后的輻射率數(shù)據(jù)的頭文件中包括中心波長(zhǎng)（）值和波段寬度（）值；要求是ENVI標(biāo)準(zhǔn)格式文件；以格式BIL存儲(chǔ)，以便FLAASH讀取處理。在ENVI的FLAASH大氣校正模塊中輸入從頭文件中讀取的各項(xiàng)定標(biāo)參數(shù)， ENVI就可以根據(jù)各項(xiàng)信息完成大氣校正處理工作。由于衛(wèi)星拍攝地點(diǎn)與時(shí)間的差別，因此氣溶膠模型的大小以及環(huán)境模型的類別都會(huì)對(duì)最終大氣校正的結(jié)果造成影響，應(yīng)在中國(guó)資源衛(wèi)星中心仔細(xì)查詢環(huán)境一號(hào)衛(wèi)星的不同模型參數(shù)，以便縮小后續(xù)處理的誤差。大氣校正前后光譜曲線如圖1所示。

圖1 大氣校正前后光譜曲線

經(jīng)過大氣校正后的光譜圖像數(shù)據(jù)不能準(zhǔn)確反映水體物質(zhì)的本質(zhì)特征，尤其對(duì)高光譜數(shù)據(jù)源來說，其波段較多，且光譜曲線的波動(dòng)范圍具有相似性，會(huì)使后續(xù)分析出現(xiàn)一定的誤差，且不便于直接從中提取易于計(jì)算的光譜曲線特征，因此需要在大氣校正之后通過包絡(luò)線去除操作對(duì)光譜曲線進(jìn)行再次處理，以突出光譜曲線的特征。

該文通過使用ENVI中的Continuum Removal模塊工具讀取經(jīng)大氣校正操作后的高光譜影像頭文件并進(jìn)行包絡(luò)線去除操作，得到去除包絡(luò)線后的高光譜影像，最終得到歸一化的水體反射率。

2 半經(jīng)驗(yàn)線性回歸建模

2.1 皮爾遜相關(guān)性分析

采用衛(wèi)星高光譜傳感器進(jìn)行水域水質(zhì)參數(shù)的光譜測(cè)量，并使用皮爾遜相關(guān)性分析法尋找水質(zhì)參數(shù)的光譜特征，選擇具有最佳相關(guān)性的波段或波段組合，使用線性回歸等方法對(duì)該地區(qū)高光譜遙感數(shù)據(jù)和水質(zhì)參數(shù)進(jìn)行反演研究的方法被稱為半經(jīng)驗(yàn)法。該文使用半經(jīng)驗(yàn)線性回歸方法，通過水體表面反射率與水質(zhì)參數(shù)濃度之間的關(guān)系構(gòu)建水質(zhì)參數(shù)反演模型。當(dāng)建立水質(zhì)參數(shù)反演模型時(shí)，先從高光譜圖像的光譜曲線中尋找水質(zhì)參數(shù)的光譜特性，并對(duì)水質(zhì)參數(shù)濃度和水體反射率進(jìn)行皮爾遜相關(guān)性分析，尋找其中相關(guān)性絕對(duì)值最大的波段，并借此建立相應(yīng)的反演模型。對(duì)2個(gè)變量、來說，它們之間的相關(guān)性可以用相關(guān)系數(shù)模型來計(jì)算。該文正是基于此關(guān)系，使用皮爾遜相關(guān)系數(shù)計(jì)算光譜波段的大小與歸一化反射率之間的相關(guān)性，皮爾遜相關(guān)性系數(shù)是用來衡量2組數(shù)據(jù)與之間線性關(guān)系的強(qiáng)度。且當(dāng)皮爾遜相關(guān)性的絕對(duì)值在不同區(qū)間時(shí)，所代表的相關(guān)性也各不相同，應(yīng)使相關(guān)性絕對(duì)值的大小盡可能地接近1，計(jì)算方法如下。

如果2組數(shù)據(jù)： {，，， …，X} 和： {，，，…，Y} 具有線性關(guān)系，那么皮爾遜相關(guān)系數(shù)如公式（1）、公式（2） 所示。

式中：P、P分別為數(shù)據(jù)集、的標(biāo)準(zhǔn)差；P為標(biāo)準(zhǔn)差的數(shù)學(xué)符號(hào)；Сov為協(xié)方差的數(shù)學(xué)符號(hào)；為、中數(shù)據(jù)的總個(gè)數(shù)；X、Y為數(shù)據(jù)集、中的數(shù)據(jù)；、分別為數(shù)據(jù)集、的平均值。

皮爾遜相關(guān)系數(shù)即2個(gè)數(shù)據(jù)集協(xié)方差與標(biāo)準(zhǔn)差的比值，取值為-1～1，其絕對(duì)值越接近1，說明相關(guān)性越顯著。當(dāng)建立水質(zhì)參數(shù)反演模型時(shí)，先從高光譜圖像的光譜曲線中尋找水質(zhì)參數(shù)的光譜特性，并對(duì)水質(zhì)參數(shù)濃度和水體反射率進(jìn)行皮爾遜相關(guān)性分析，尋找其中相關(guān)性絕對(duì)值最大的波段，并建立相應(yīng)的反演模型。如果效果不好，則對(duì)光譜曲線進(jìn)行一階導(dǎo)數(shù)和二階導(dǎo)數(shù)或包絡(luò)線去除等處理，尋找具有明顯的特征波段，并對(duì)其進(jìn)行相關(guān)性分析，該文最終選擇了包絡(luò)線去除的方法。該文使用SPSS 25軟件分析了HJ-1A遙感衛(wèi)星的單波段反射率和DO濃度之間的相關(guān)關(guān)系，在SPSS軟件中輸入DO濃度作為一個(gè)變量，將HJ-1A高光譜80個(gè)波段的反射率作為另一個(gè)變量，建立雙變量皮爾遜相關(guān)關(guān)系，再建立實(shí)測(cè)數(shù)據(jù)與遙感反射率的相關(guān)關(guān)系，去除異常值之后的皮爾遜相關(guān)系數(shù)如圖2所示。

對(duì)去除包絡(luò)線后的實(shí)測(cè)溶解氧數(shù)據(jù)與HJ-1A影像光譜曲線建立的相關(guān)性圖像進(jìn)行分析發(fā)現(xiàn)，在598 nm和656 nm處出現(xiàn)明顯的反射峰，在609 nm和669 nm處出現(xiàn)明顯的吸收谷。對(duì)這兩處的反射峰波段和吸收谷波段進(jìn)行多次組合測(cè)試，尋找相關(guān)系數(shù)高的波段組合。分別對(duì)609 nm、669 nm、598 nm和656 nm進(jìn)行反射峰與吸收谷的相加減、相乘除或其他的組合方式，最終獲得了以下的相關(guān)關(guān)系，波段組合相關(guān)系數(shù)如圖3所示。

圖3 波段組合皮爾遜的相關(guān)性系數(shù)

2.2 DO反演模型建立

該文選擇相關(guān)系數(shù)最高的A+A/B波段，并以該組合波段與溶解氧的關(guān)系建立線性回歸方程，為了保證模型的準(zhǔn)確性，盡可能地減少隨機(jī)誤差，該文隨機(jī)抽取了29 個(gè)監(jiān)測(cè)點(diǎn)數(shù)據(jù)，其中20個(gè)監(jiān)測(cè)點(diǎn)數(shù)據(jù)為訓(xùn)練集，其余9個(gè)監(jiān)測(cè)點(diǎn)數(shù)據(jù)為預(yù)測(cè)集，取平均結(jié)果建立散點(diǎn)圖，并擬合回歸方程，最終得到的反演模型為= 39.68- 290.83+787.52-931.03+ 409.6，并通過計(jì)算得出擬合系數(shù)=0.639，反演效果較好，建立的模型圖如圖4所示。

圖4 線性回歸模型

利用模型預(yù)測(cè)剩下的9個(gè)監(jiān)測(cè)站的溶解氧濃度對(duì)溶解氧模型的精度進(jìn)行驗(yàn)證，由實(shí)測(cè)值與預(yù)測(cè)值之間的擬合曲線可以看出模型的適用性，DO濃度的實(shí)測(cè)值和預(yù)測(cè)值的擬合曲線如圖5所示。根據(jù)實(shí)測(cè)值與預(yù)測(cè)值的線性回歸曲線分析可以看出模型的擬合程度較好，預(yù)測(cè)擬合系數(shù)指數(shù)也達(dá)到了0.669 8，說明運(yùn)用數(shù)學(xué)模型建立溶解氧含量與遙感圖之間的關(guān)系的效果較好。

圖5 DO 實(shí)測(cè)值與預(yù)測(cè)值擬合曲線

3 結(jié)語

綜上所述，基于遙感的水質(zhì)監(jiān)測(cè)已經(jīng)隨著水體污染的問題變得越來越重要，而基于常規(guī)水質(zhì)監(jiān)測(cè)的方法有許多缺點(diǎn)，會(huì)受各種因素的制約和限制，很難滿足較大水域范圍、實(shí)時(shí)空間動(dòng)態(tài)以及便利快速的監(jiān)測(cè)要求，該文通過查閱國(guó)內(nèi)外基于遙感進(jìn)行水質(zhì)監(jiān)測(cè)的理論、方法和技術(shù)，充分發(fā)揮了高光譜遙感的優(yōu)勢(shì)，以香港近海區(qū)域?yàn)檠芯繉?duì)象，利用實(shí)地監(jiān)測(cè)數(shù)據(jù)、環(huán)境與災(zāi)害監(jiān)測(cè)預(yù)報(bào)HJ-1A衛(wèi)星HSI高光譜數(shù)據(jù)影像，通過建立線性回歸方程的方法對(duì)香港近海區(qū)域溶解氧濃度進(jìn)行監(jiān)測(cè)。

此外，該文還總結(jié)了國(guó)內(nèi)外有關(guān)溶解氧反演技術(shù)的最新進(jìn)展，水質(zhì)反演理論由分析法向半經(jīng)驗(yàn)方法發(fā)展，對(duì)遙感數(shù)據(jù)也提出了更高的要求。對(duì)HJ-1A高光譜遙感數(shù)據(jù)進(jìn)行數(shù)據(jù)處理、大氣校正、幾何精校正及包絡(luò)線去除等預(yù)處理。以香港近海區(qū)域?yàn)檠芯浚Y(jié)合環(huán)境保護(hù)署的實(shí)測(cè)數(shù)據(jù)、HJ-1A影像數(shù)據(jù)對(duì)DO進(jìn)行遙感分析。使用線性回歸方法對(duì)DO實(shí)測(cè)值和光譜數(shù)據(jù)進(jìn)行建模，指數(shù)達(dá)到0.669 8。