王茹月,周航宇
(湘潭大學計算機學院網絡空間安全學院,湖南 湘潭 411105)
勘測地表水體并劃定其空間分布,對了解水文過程和管理水資源具有重要意義[1-2]。利用遙感多光譜數據,使用水體指數方法,人們已經研究得到廣泛使用的自動水體提取方法。遙感水體指數方法是提取全球和大區域水體信息的最有效方法[3-4],其關鍵點有兩處:一是歸一化使用的水體指數,二是自適應閾值方法[5-7]。使用水體指數方法進行水體信息提取的最終效果在很大程度上依賴于分割閾值的選取結果, 閾值的好壞會直接影響到水體信息提取精度[8-9]。Mcfeeters在提出NDWI指數的同時,將地表水的閾值看作0,即認為歸一化后正值為水,負值為非水,這是靜態閾值方法。隨后研究人員發現并非所有的圖像都是完美適用的,使用NDWI方法閾值為0時并不能準確地將建筑物類型與水體區分開。文獻[10]創建了自動提取水指數的AWEI方法,該方法以0為閾值就可以自動提取水體指數;但從現有情況來看,以0為閾值無法獲得最佳效果,仍需設定最佳閾值[4,11]。因此一些自適應閾值方法不斷被研究使用,如迭代法、大津閾值法、直方圖雙峰法等。然而,每種方法均有各自的優缺點和不適用情況。本文將對直方圖雙峰法展開研究,以期提高遙感影像水體信息提取的精度。
本文使用Landsat 8/9遙感數據,所有數據均從USGS的門戶網站Eatrh Explore上獲取。選取的L2級產品已經官方處理,無須進行幾何校正、輻射定標及大氣校正過程[12],只需轉換為地表反射率即可參與計算。試驗區選取中國東部5個區域的影像,每景影像內均有湖泊及河流。針對水體的研究均具有季節性和時間性趨勢[13],考慮季節、溫度、水生植物等對結果的影響,每個試驗區均選取冬夏兩個季節的數據,共10景影像,詳細數據參數見表1。考慮含云量對最終結果的影響較大,每景影像都已盡量選取含云量較低的數據。

表1 試驗數據參數
直方圖雙峰法最早是在1996年由Prewitt提出,全稱為灰度直方圖雙峰法。該方法是將圖像的灰度值進行分割,若灰度直方圖呈明顯的雙峰狀,則選取兩峰之間的谷底所對應的灰度級作為閾值。即對于一幅灰度值范圍在Z1到Zk的圖像,使用直方圖進行統計,x軸為[0,255]的灰度值,y軸為對應灰度值的像素點數。若直方圖形成兩個山峰形狀,則將兩峰之間y值最小時對應的x值設為Zt,Zt為該圖像的閾值。
在實際遙感影像水體提取時,可直接將歸一化后的遙感數據作為雙峰法的數據統計對象,無須再處理成0~255的灰度值。歸一化后數據范圍為[-1,1],將其分成若干個區間,統計每個區間包含的數據個數,即每個區間包含的像素點個數,再以各區間中點為x軸,各區間包含的像素點個數為y軸形成雙峰圖,取兩峰之間的最低谷所對應的x值作為最終的自適應閾值。
直方圖雙峰法的優點在于算法易于理解且容易實現,計算量小,但也存在明顯不足。使用雙峰法求遙感影像的閾值,必須有一定的圖像先驗知識,水域與非水域分布面積應較為均等。當一幅圖大面積為水體或非水體,體現在圖像結果上是單峰數據或兩峰之間差距過大;若數據統計過程不準確或噪音點過多,則會出現多峰數據情況。在雙峰法使用過程中若出現以上情況則獲取的閾值結果是無意義的。
針對數據單峰及雙峰差距過大的情況,作出如下改進:由于一景圖像過大無法直接進行處理或處理耗時較長,需將一景圖像裁剪成網格,再對每一幅裁剪后的圖像分別進行處理。圖1展示部分網格圖像歸一化后的結果,明顯存在一幅圖出現大面積水域或非水域的情況。在此前提下,本文選擇在影像裁剪成網格之后先挑選水和非水分布較為均衡的圖像進行閾值求解,取均衡圖像閾值的均值作為這一景圖像的閾值,可以有效防止在雙峰法使用過程中出現單峰及雙峰差距過大的情況。如圖1所示,mndwi_1、mndwi_35、mndwi_48、mndwi_200樣式為單一水體或非水體,需要剔除;mndwi_51、mndwi_85、mndwi_104為分布較均衡圖像,可作為選擇。

圖1 部分歸一化結果圖輪廓展示
針對數據多峰情況,在試驗過程中發現多峰情況出現的嚴重程度與直方圖區間個數的多少有較大關系,在逐步試驗過程中最終確定可通過改變劃分區間個數調節最大峰的像素點所占比重,進而可以有效緩解多峰情況的出現。為研究得出最大峰占比多少時閾值為最佳結果,進行對比試驗,試驗流程如圖2所示。
本文試驗研究過程如下:
(1)下載獲得Landsat 8/9圖像,選出b2、b3、b4、b6波段分別裁剪成128×128大小的網格圖片,從序號1開始按照順序命名圖片。
(2)對b3、b6波段使用MNDWI方法進行歸一化,歸一化后數據范圍為[-1,1]。
(3)將歸一化后圖像中水和非水分布不均勻的圖像剔除,保存剩余圖像。
(4)對每幅圖像選取最大峰占比值從0.7到0.2每隔0.05遞減的數據進行試驗,即分別取最大峰像素點數占總點數比重為0.2、0.25、0.3、0.35、0.4、0.45、0.5、0.55、0.6、0.65、0.7的條件作對比試驗。對圖像數據不斷循環調整區間份數以改變最大峰占比,達到要求的條件后取最大峰與次大峰之間最低處x值作為該圖像閾值。取所有圖像的均值作為這一景圖像的最終閾值。
(5)根據最終閾值,將大于閾值的像素點面積求出即為水體的面積。
(6)精度判定。
將b2、b3、b4波段融合成RGB圖像,在裁剪后的網格圖像中隨機挑選10幅圖進行人工面積標注,與自動求得的水域面積進行精度對比。基于不同最大峰占比值的雙峰法精度對比結果見表2。為宏觀、清晰地反映整體變化趨勢,圖3和圖4為10景圖像的精度變化折線圖,圖3為冬季數據結果,圖4為夏季數據結果。

圖3 冬季數據精度變化

圖4 夏季數據精度變化

表2 基于不同最大峰占比值的精度對比 (%)
為獲得最終的結論,合理的綜合數據得出結果,本文計算了精度變化的均值及方差,如圖5、圖6所示。由試驗結果可得出,最大峰像素點占比0.3時為精度最佳結果,即使用雙峰法求自適應閾值時,應先通過改變區間份數調節最大峰像素點數占總點數到0.3,再取最大峰與次大峰之間最低谷為閾值時所求的閾值結果為分類精度最高值。改進之后,雙峰法再出現單峰、雙峰等不適用情況的比例由0.775降到0.07,即參與最終閾值計算的數據已基本不再出現不適用雙峰法的情況。

圖5 精度均值變化

圖6 精度方差變化
為展示改進的雙峰法與其他自適應閾值方法的效果對比,另外使用迭代閾值法進行水體提取,其他試驗步驟相同,只替換了自適應閾值方法,數據結果見表3,在相同試驗面積的情況下,改進的雙峰法精度總要高于迭代法。計算各研究區的Kappa系數、總體分類精度參數[14],平均結果分別為0.858、0.931,說明改進的雙峰法的水體提取結果處于高度一致性范圍內。

表3 改進雙峰法與迭代法精度比較 (%)
本文應用改進的直方圖雙峰法,使用2022年4月21日Landsat 8遙感影像數據計算得到鄱陽湖總水域面積為2 001.5 km2,水域結果如圖7所示。其中,黑色部分為水體,與人工標注面積相比精度高達99.30%,證明該方法的可用性。

圖7 鄱陽湖水體提取過程展示
本文詳細研究了圖像閾值分割方法之一的直方圖雙峰法,其優點在于算法簡單,計算量小,但若使用雙峰法需有一定的圖像先驗知識,如圖像中水、非水分布面積應較為均等,數據直方圖應形成兩個山峰狀數據才能較準確提取閾值,因此雙峰法的使用有一定的局限性。通過本文試驗研究,針對雙峰法的局限性做出如下改進。
(1)對一景圖像裁剪后挑選水、非水分布均勻的圖像參與閾值計算,取其閾值均值作為這一景圖像的閾值,可以有效防止單峰及雙峰差距過大情況的出現。
(2)通過改變區間份數調節最大峰像素點數占總點數比值到0.3時,既可以有效防止出現多峰數據,也能在此條件下取得最好的精度結果。
改進后的雙峰法最高精度可達99.14%,高于原雙峰法精度,且有效減少雙峰法使用時依賴水域分布以及圖像先驗知識的情況,大大提高了雙峰法的適用性。