秦 雁，吉紅香，黃本勝，陳亮雄，周 宇，楊靜學(xué)

(1.廣東省水利水電科學(xué)研究院，廣州 510635； 2. 河口水利技術(shù)國家地方聯(lián)合工程實驗室，廣州 510635； 3.廣東省粵港澳大灣區(qū)水安全保障工程技術(shù)研究中心，廣州 510635)

1 研究背景

采砂是一項劇烈干預(yù)河道自然演變的經(jīng)濟(jì)活動[1]。適當(dāng)有序的采砂不僅有利于河道行洪及航運通航，且能提供寶貴的砂石資源[1-2]。近年來，我國城市建設(shè)與基礎(chǔ)設(shè)施建設(shè)等經(jīng)濟(jì)活動對砂石需求量持續(xù)增長，高利益驅(qū)使下各種非法采砂活動屢禁不止。盲目、超量和無序采砂對防洪、涉水工程、供水、航運、漁業(yè)、水生態(tài)安全以及社會穩(wěn)定產(chǎn)生惡劣的負(fù)面影響[1-5]。然而，目前采砂監(jiān)管工作主要是以執(zhí)法船只或車輛巡邏的方式進(jìn)行的，面臨著監(jiān)管區(qū)域廣、專職人員少、證據(jù)收集難等困難，給河湖采砂動態(tài)監(jiān)管帶來了挑戰(zhàn)，急需新的監(jiān)管手段[4]。

遙感技術(shù)具有速度快、成本低、范圍廣的特點，近年來廣泛應(yīng)用于水體信息提取[6]、水質(zhì)監(jiān)測[7]、采砂監(jiān)管[8-13]等水資源調(diào)查、監(jiān)察、保護(hù)領(lǐng)域。目前，遙感技術(shù)應(yīng)用于采砂監(jiān)管主要包括采砂船識別和高渾濁水體識別2個方面。應(yīng)用高分辨率衛(wèi)星[3]或無人機遙感圖像[4]，采用人工解譯[3-4]或深度學(xué)習(xí)圖像自動識別[8]方法可識別采砂船、采砂場。另外，采砂擾動之后，水體濁度、透明度、總懸浮物濃度、懸浮泥沙濃度等相關(guān)指標(biāo)異常升高，能夠通過水色遙感技術(shù)判斷采砂活動信息[10]。鄔國鋒等[10-11]建立了鄱陽湖水體透明度(Secchi Disk Depth,SDD)與 LANDSAT TM5 圖像紅、藍(lán)波段線性經(jīng)驗關(guān)系，分析采砂對水體透明度、生物多樣性等方面的影響。張偉等[12]利用藍(lán)光、紅光、近紅外波段計算總懸浮物指數(shù)(Total Suspended Matter Index,TSMI)，并建立其與總懸浮物濃度的對數(shù)經(jīng)驗關(guān)系，分析采砂區(qū)時空變化特征。楊靜學(xué)等[13]應(yīng)用TSMI分析發(fā)現(xiàn)鶴地水庫懸浮物濃度高值區(qū)存在采砂情況。目前，常規(guī)采砂遙感監(jiān)測需經(jīng)過大氣校正、懸浮物濃度反演兩級過程，通過設(shè)定總懸浮物濃度經(jīng)驗閾值識別采砂嫌疑區(qū)。由于大氣及水質(zhì)組分時空差異大，受限于大氣糾正、總懸浮物濃度反演精度，分割閾值常需人工干預(yù)，難以實現(xiàn)多時相采砂區(qū)自動化識別。此外，總懸浮物濃度反演模型常利用近紅外波段泥沙高反射特性，不適用于僅包含紅、綠、藍(lán)3個波段的常規(guī)無人機自然圖像。

水色變化與水中光學(xué)組分濃度變化直接相關(guān)，通過水色參數(shù)可監(jiān)測水體水質(zhì)的時空變化，如黑臭水體[14]、水體富營養(yǎng)化[15]、水體透明度[16-17]、水色異常[18]等。主波長水色參數(shù)不需近紅外波段，僅采用紅、綠、藍(lán)3個波段計算，且對大氣糾正抗擾性強[17]，在監(jiān)測水體濁度、河道采砂方面極具應(yīng)用潛力，而目前這方面的研究鮮見。

本文以楓樹壩水庫、鶴地水庫水域為研究區(qū)域，應(yīng)用Landsat8、Sentinel-2衛(wèi)星遙感圖像開展采砂區(qū)主波長水色參數(shù)分析，并采用自然斷點法自動設(shè)定劃分閾值，旨在建立一種適用于紅、綠、藍(lán)3波段圖像，無需大氣糾正，自動識別采砂嫌疑區(qū)的方法，為常態(tài)化采砂監(jiān)管提供動態(tài)監(jiān)測信息。

2 數(shù)據(jù)與方法

2.1 研究區(qū)介紹

楓樹壩水庫、鶴地水庫均為廣東大(I)型水庫，因位于省際交界區(qū)域，交通不便，管理薄弱，造成非法采砂屢禁不止，難以監(jiān)管。楓樹壩水庫集雨面積5 150 km2，總庫容19.4億 m3，位于東江上游的龍川縣楓樹壩鎮(zhèn)，是集防洪、發(fā)電、供水功能于一體的綜合性水庫。近年來，楓樹壩水庫入庫河流尋鄔水常出現(xiàn)水體異常渾濁的情況，疑似存在非法采砂現(xiàn)象。鶴地水庫集雨面積1 495 km2，總庫容11.44億 m3，庫區(qū)地跨廣東廉江、化州和廣西陸川、博白4個縣(市)。鶴地水庫是一座以灌溉為主，兼顧防洪、發(fā)電、旅游的水利樞紐骨干工程，承擔(dān)著雷州半島97 333 hm2農(nóng)田灌溉和近400萬人口生活用水供水任務(wù)。2012—2018年，鶴地水庫交界水域非法采砂行為時有發(fā)生，僅依靠人工巡查，難以掌握采砂的時空分布信息。引入遙感監(jiān)測方法后，通過加大跨區(qū)執(zhí)法力度，并經(jīng)過長效治理，自2019年之后，鶴地水庫非法采砂現(xiàn)象得到了有效控制。

2.2 數(shù)據(jù)源

采用Sentinel-2衛(wèi)星T50RLN標(biāo)準(zhǔn)分幅數(shù)據(jù)，包括2019-11-08、2019-11-23、2019-12-08、2020-02-21、2020-04-16共5個時相圖像，于https:∥scihub.copernicus.eu/dhus/#/home網(wǎng)站下載獲得，用于楓樹壩水庫采砂監(jiān)測。Sentinel-2衛(wèi)星圖像為L1C級產(chǎn)品，是經(jīng)過幾何糾正的表觀反射率數(shù)據(jù)。同時采用Sentinel-2衛(wèi)星T49QDE標(biāo)準(zhǔn)分幅數(shù)據(jù)，以及相近時相行列號為124045的Landsat8衛(wèi)星數(shù)據(jù)監(jiān)測鶴地水庫，檢驗方法對不同衛(wèi)星數(shù)據(jù)的適用性。Landsat8衛(wèi)星圖像于https:∥earthexplorer.usgs.gov/網(wǎng)站下載，為L1級數(shù)據(jù)，需進(jìn)行幾何糾正和輻射糾正處理獲得表觀反射率數(shù)據(jù)。

2.3 研究方法

2.3.1 CIE-XYZ顏色系統(tǒng)介紹

水體中非藻類懸浮物、浮游植物和有色可溶性有機物，這3種光活性物質(zhì)及水分子共同決定了水體顏色[17]。由于水中物質(zhì)組分的物理化學(xué)性質(zhì)不同、濃度不同，水體顯示出不同的顏色，因此水體顏色具有表征水體水質(zhì)的潛力。

為了建立水色與水質(zhì)的關(guān)系，需要對顏色進(jìn)行數(shù)字化表示。國際照明委員會(CIE)于1931年建立了一套標(biāo)準(zhǔn)顏色系統(tǒng)CIE-XYZ顏色系統(tǒng)[16]。該系統(tǒng)將CIE-RGB系統(tǒng)中的3個原色R、G、B轉(zhuǎn)換為3個理想的原色X、Y、Z，關(guān)系式[19]為

(1)

CIE-XYZ顏色系統(tǒng)中光譜的三刺激值X、Y、Z全為正值，并且X、Z兩原色只代表色度，沒有亮度，光亮度只與三刺激值Y成比例[16]。為了便于目視三刺激值定義的顏色，國際照明委員會規(guī)定了CIE-xy二維色度圖，定義x+y+z=1，對z進(jìn)行降維處理，使每種顏色都對應(yīng)一個色度坐標(biāo)。根據(jù)CIE-XYZ系統(tǒng)中的三刺激值X、Y、Z計算色度圖上二維坐標(biāo)(x,y)，計算公式為

(2)

圖1 CIE-xy色度圖Fig.1 CIE-xy chromaticity diagram

CIE-xy二維色度圖如圖1所示。圖1中S稱作等能白光點，其色度坐標(biāo)為(1/3，1/3)，表示3種原色等量混合。C點是任一顏色色度坐標(biāo)點，坐標(biāo)為(x,y)，從S指向C的直線與光譜軌跡的交點D所對應(yīng)的波長，即該顏色的主波長λd。補色波長λc為從C指向S的直線與光譜軌跡的交點M所對應(yīng)的波長。從等能白光點到色度圖外圍線的每一條射線上的點都具有相同的顏色主波長。顏色主波長λd是顏色量化的重要指標(biāo)，它將可見光顏色從380～700 nm以1 nm間隔來標(biāo)識，能夠以具體的波長形式表示一種顏色的色調(diào)。

2.3.2 主波長水色參數(shù)采砂識別算法

水庫水域在采砂擾動之后，水體濁度升高，水體光譜反射峰向紅光通道偏移，水體表現(xiàn)為黃棕色。水體逐漸渾濁，水色定量表征體現(xiàn)為主波長水色參數(shù)由藍(lán)綠波段向黃紅波段逐漸遞增，即水體濁度升高，主波長增加。

原始遙感圖像經(jīng)過預(yù)處理后獲得表觀反射率圖像，再依據(jù)改進(jìn)歸一化差異水體指數(shù)(Modified Normalized Difference Water Index,MNDWI)進(jìn)行水陸分離[20]，獲得水體掩模。將水體掩模范圍內(nèi)紅、綠、藍(lán)3波段表觀反射率代入式(1)，計算X、Y、Z三刺激值。根據(jù)式(2)，將XYZ轉(zhuǎn)為色度圖坐標(biāo)(x,y)。采用張晨朝[21]提出的分區(qū)法，由色度坐標(biāo)(x,y)計算獲得水體主波長水色參數(shù)圖像。

進(jìn)而采用自然斷點分類法Natural Breaks(Jenks)[22]對主波長圖像進(jìn)行分類測試，分析確定高濁度水體的劃分閾值，可實現(xiàn)采砂嫌疑區(qū)識別。自然斷點分類法的思想是在確定分類個數(shù)的情況下，計算各種閾值組合下分類的方差和，其值最小的就是最優(yōu)的分類結(jié)果，即可自動確定最優(yōu)分類閾值。根據(jù)已知試驗區(qū)的情況，主波長圖像分為高、中、低3個等級，數(shù)值最高的等級可識別為采砂嫌疑區(qū)。主波長水色參數(shù)計算及自然斷點分類程序算法均采用python語言，調(diào)用gdal開源庫編程實現(xiàn)。

3 結(jié)果與分析

3.1 鶴地水庫

采用本方法對Landsat8衛(wèi)星、Sentinel-2衛(wèi)星多時相遙感圖像進(jìn)行處理，測試不同傳感器、不同大氣狀況下，自動識別鶴地水庫采砂嫌疑區(qū)的效果，結(jié)果如圖2、圖3所示。

圖3 鶴地水庫Sentinel-2衛(wèi)星遙感采砂嫌疑區(qū)識別Fig.3 Identification of suspected sand mining area by Sentinel-2 satellite remote sensing in Hedi Reservoir

Landsat8衛(wèi)星不同時相主波長水色參數(shù)圖像顯示了水體濁度差異，九洲江水域及化州庫灣異常渾濁區(qū)域主波長值域范圍為525～590 nm，渾濁程度較小的庫心區(qū)域主波長值域范圍為480～490 nm。根據(jù)自然斷點法分別以最高類閾值>525.45、525.50、525.47、525.47 nm將2015年1月1日、2016年11月3日、2017年12月2日、2017年12月24日共4個不同時相主波長高值區(qū)識別為采砂嫌疑區(qū)。由4個不同時相的結(jié)果可知，采砂區(qū)劃分閾值十分穩(wěn)定，即使2017年1月22日及2017年12月24日存在薄云，閾值差異仍然非常微小，說明方法抗擾性強，無需對圖像進(jìn)行大氣糾正預(yù)處理，易于使用。

Sentinel-2衛(wèi)星與相近時相Landsat8衛(wèi)星的主波長水色參數(shù)圖像反映的水體濁度一致，2016年10月30日、2017年1月28日、2017年10月30日、2017年12月19日共4個不同時相采砂嫌疑區(qū)識別閾值分別是522.55、522.64、522.47、525.72 nm，差異也非常小。由此說明主波長水色參數(shù)對于不同傳感器具有很強的適應(yīng)性和穩(wěn)定性。

Landsat8衛(wèi)星2015年和2016年2個時相、Sentinel-2衛(wèi)星2016年一個時相主波長水色參數(shù)遙感結(jié)果顯示，廣西、廣東交界九洲江水域異常渾濁，2017年Landsat8、Sentinel-2衛(wèi)星不同時相的遙感結(jié)果均顯示，除了九洲江水域，茂名化州與湛江廉江交界庫灣水域也異常渾濁。主波長水色參數(shù)反映的鶴地水庫水域渾濁程度時空分布特征與楊靜學(xué)等[13]的研究結(jié)果相似，均體現(xiàn)了九洲江及化州庫灣存在采砂嫌疑。但由于主波長水色參數(shù)模擬人眼對水色的感官，將紅、綠、藍(lán)3個波段信息壓縮，受大氣狀況影響小，可通過自然斷點分類法自動劃分高值區(qū)識并將其別采砂嫌疑區(qū)，具有人工干預(yù)少、抗擾性強的優(yōu)勢，更適用于多時相遙感反演結(jié)果對比分析。

圖4 楓樹壩水庫Sentinel-2衛(wèi)星遙感采砂嫌疑區(qū)識別Fig.4 Identification of suspected sand mining area by Sentinel-2 satellite remote sensing in Fengshuba Reservoir

3.2 楓樹壩水庫

對2019年11月8日至2020年4月16日共5個時相Sentinel-2衛(wèi)星遙感圖像計算主波長水色參數(shù)，結(jié)果如圖4所示。采用自然斷點法劃分主波長水色參數(shù)為3個等級，經(jīng)測試各時相最高值分類閾值非常穩(wěn)定，為525.4±0.1 nm。以主波長高值區(qū)識別為采砂嫌疑區(qū)，由圖4可見，方法自動性強、穩(wěn)定可靠，既能識別水體較為渾濁時的情況，也未對2020年2月21日及2019年11月8日水體較為清澈時的情況產(chǎn)生誤判。

多時相主波長水色參數(shù)遙感結(jié)果顯示，廣東、廣西交界尋鄔水河段經(jīng)常出現(xiàn)異常渾濁現(xiàn)象，與周圍支流水色形成明顯對比，極可能存在采砂現(xiàn)象。根據(jù)2019年11月6日以及2020年12月17日高分辨率衛(wèi)星遙感圖像顯示，該河段在圖4(a)中A、B點位存在采砂船。如圖5所示，在B點位清晰可見采砂船周圍水體渾濁，正在進(jìn)行采砂活動。由高分辨率衛(wèi)星遙感圖像采砂船解譯結(jié)果，證實了主波長水色參數(shù)采砂區(qū)識別結(jié)論。

圖5 楓樹壩水庫采砂船識別驗證Fig.5 Verification and identification of sand dredgerin Fengshuba reservoir

4 結(jié) 論

經(jīng)楓樹壩水庫、鶴地水庫不同地區(qū)，Landsat8、Sentinel-2衛(wèi)星不同遙感圖像試驗發(fā)現(xiàn)，主波長能夠反映水色，其值隨水體濁度升高而增大。采用自然斷點法劃分主波長為高、中、低3個等級，高值區(qū)分類閾值非常穩(wěn)定。采砂區(qū)高渾濁水體可采用自然斷點法自動選擇閾值判定，或者依據(jù)本文試驗采用主波長>525 nm作為判斷條件。試驗證明主波長水色參數(shù)在識別采砂方面具有方便快捷、自動識別、適用多種衛(wèi)星圖像及無人機自然圖像的優(yōu)勢。方法可為采砂動態(tài)監(jiān)管提供時空分布信息，輔助管理人員有針對性地開展重點巡查，有效提升常態(tài)化采砂監(jiān)管的能力。

由于缺少實測數(shù)據(jù)，本文尚未建立濁度與主波長水色參數(shù)的定量關(guān)系，后續(xù)需進(jìn)行光譜采樣試驗，選取水色參數(shù)和光譜指數(shù)建立深度學(xué)習(xí)水體濁度遙感反演模型。此外，可進(jìn)一步提升數(shù)據(jù)處理效率，應(yīng)用Google Earth Engine遙感云平臺高頻率計算Landsat8、Sentinel-2等衛(wèi)星圖像主波長水色參數(shù)，統(tǒng)計多時相采砂嫌疑區(qū)識別頻度圖像，快速評估經(jīng)常性渾濁水體，以此排除因雨洪影響而渾濁的水域，識別異常情況，提供采砂監(jiān)管輔助信息。