楊 志，董坤明，王 珊，陳 槐，王乃茹，范紅霞

(1. 安徽省現(xiàn)代交通設(shè)計(jì)研究院有限責(zé)任公司，安徽 合肥 230011； 2. 水發(fā)規(guī)劃設(shè)計(jì)有限公司，山東 濟(jì)南 250100； 3. 江寧區(qū)水務(wù)局，江蘇 南京 210000； 4. 南京水利科學(xué)研究院 水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210029)

隨著南京沿江經(jīng)濟(jì)的發(fā)展，對(duì)長江岸線、水資源開發(fā)利用的要求愈發(fā)迫切，對(duì)長江河勢(shì)穩(wěn)定、防洪安全的要求也越來越高。探討南京河段的河床與洲灘演變特征，可為南京段河勢(shì)穩(wěn)定、岸線資源的開發(fā)利用等提供科學(xué)依據(jù)，具有重要的研究意義[1]。長江南京段是典型的分汊型河道，相關(guān)的河床演變研究較多，但受水文測(cè)站布點(diǎn)、地形資料等制約，研究往往只局限于局部河段[2]。臧英平等[3]根據(jù)監(jiān)測(cè)資料分析了長江南京段近20 年5 個(gè)時(shí)期的河勢(shì)及岸段變化，表明河勢(shì)基本穩(wěn)定，深泓線橫向擺動(dòng)較小，深槽區(qū)平面位置基本穩(wěn)定，岸線沖淤變化不大；朱宇馳[4]利用長江南京段近50 年4 個(gè)時(shí)期的航道圖，結(jié)合GIS、RS 和數(shù)字高程模型（DEM），分析了南京段的岸線變遷及水下地形沖淤演變，總體表現(xiàn)出淤積-沖刷-穩(wěn)定的特征；徐韋等[5]利用1998 年和2013 年南京段的歷史水下地形數(shù)據(jù)，結(jié)合2015—2016 年多波束測(cè)深、流速與河床沉積物數(shù)據(jù)，探討了南京段河槽演變對(duì)三峽水庫等人類活動(dòng)的響應(yīng)規(guī)律。

傳統(tǒng)的河床演變分析很大程度上受制于水文站點(diǎn)和野外實(shí)測(cè)資料，僅利用傳統(tǒng)觀測(cè)手段很難全面有效地掌握河流基本信息。隨著科學(xué)技術(shù)的發(fā)展，衛(wèi)星遙感可獲取多時(shí)相觀測(cè)影像數(shù)據(jù)，能為河床演變提供直接有效的方法，越來越多的研究人員采用衛(wèi)星遙感影像研究河床平面形態(tài)的變化[6-8]。河床平面形態(tài)的變化與河床演變關(guān)系密切，河床平面形態(tài)通過平面位置遷移、河寬變化、沙洲的形態(tài)變化及其產(chǎn)生和消亡過程、汊道增減等，來反映河床演變過程[9]。已有關(guān)于長江的衛(wèi)星遙感研究成果極多，但極少針對(duì)南京河段的詳細(xì)研究。嚴(yán)夏青等[10]提取了2016—2020 年月度長江干流上海至宜賓段水域面積，并分析其年際、年內(nèi)變化規(guī)律；姜玲玲等[11]研究了長江干流江心洲在1984—2019 年的變化過程，發(fā)現(xiàn)江心洲在過去35 年平面形態(tài)總體上以擴(kuò)張為主，但存在一定的區(qū)域性差異。

鑒于目前南京段長時(shí)間序列演變分析的研究成果較少，且遙感影像演變方面研究不足的現(xiàn)狀，通過廣泛收集南京河段長時(shí)間遙感影像序列（共約800 幅），并引入本征正交分解法分析河床邊界和洲灘的時(shí)空精細(xì)演變特征，以期為河流的合理開發(fā)利用、科學(xué)保護(hù)和有效管理提供支撐。

1 研究區(qū)域及數(shù)據(jù)分析方法

長江南京段位于江蘇省境內(nèi)，上起貓子山，下至三江口，全長93 km，是長江中下游干流14 個(gè)重點(diǎn)河段之一（圖1）。長江兩岸干堤總長191.05 km（不含洲堤、以山代堤段和通江河港堤），其中江南97.72 km，江北93.33 km。徑流下行受潮水上行頂托形成兩種水勢(shì)交互，在南京段形成了新生洲、新濟(jì)洲、梅子洲及八卦等諸多沙島和淺灘，使長江在南京段多次分汊，主、支汊交替變化。

1.1 數(shù)據(jù)來源

本文收集的數(shù)據(jù)主要來自Landsat（地球資源技術(shù)衛(wèi)星），包括1984—2020 年間南京河段所屬區(qū)域（Landsat 軌道號(hào)（120，38））約800 幅衛(wèi)星遙感圖片（衛(wèi)星重現(xiàn)期為15 d，每月2 幅），僅利用其中兩個(gè)波段（Green 綠、NIR 近紅外）的數(shù)據(jù)，空間分辨率為30 m。南京潮位資料來源于歷年的《長江流域水文資料》第6 卷第6 冊(cè)。

1.2 計(jì)算方法

采用基于綠波段與近紅外波段的歸一化差分水體指數(shù)（NDWI）來凸顯影像中的水邊線信息[12]。以2020-02-06 的南京段遙感影像為例，經(jīng)過歸一化差分水體指數(shù)計(jì)算后，采用otsu 閾值分割算法進(jìn)行二值化（圖2(a)）?？梢?，NDWI 指數(shù)可清晰凸顯出水邊線信息。圖2 中黑色區(qū)域取值為0，白色區(qū)域取值為1。

圖2 長江南京河段河床邊界與洲灘提?。?020-02-06 遙感影像）Fig. 2 The extracted shorelines and sandbars in the Nanjing reach of the Yangtze River

當(dāng)水邊線信息凸顯后，采用形態(tài)學(xué)方法提取河床邊界和洲灘，主要利用開操作與閉操作。利用形態(tài)學(xué)方法提取長江南京段河床邊界和5 個(gè)大洲灘（圖2(b)）。為清晰顯示，圖2(b)中河床邊界和洲灘的比例尺不同。由于Landsat 衛(wèi)星影像分辨率為30 m，在提取河床邊界和洲灘面積時(shí)會(huì)引發(fā)誤差，本文按最不利情況考慮，即所有的水邊線提取誤差都為30 m。

本征正交分解法（Proper Orthogonal Decomposition，POD）主要對(duì)高維度特征數(shù)據(jù)進(jìn)行降維，僅保留高維度數(shù)據(jù)最重要的特征，去除噪聲和不重要的特征，從而簡化數(shù)據(jù)分析。本文主要用于河床邊界或洲灘形態(tài)序列中含能模態(tài)的提取。POD 方法將信號(hào)序列投影到最優(yōu)函數(shù)空間，使信號(hào)在該最優(yōu)函數(shù)空間投影的平均能量最大，則該最優(yōu)函數(shù)空間的解為信號(hào)序列的空間相關(guān)系數(shù)矩陣的特征向量空間，可通過特征分解進(jìn)行求解，具體步驟見文獻(xiàn)[13]。

2 計(jì)算結(jié)果

2.1 遙感影像篩選

采用目視解譯法人為挑選了289 張不受云層影響的遙感影像，分析可用遙感影像的時(shí)域分布特征。由圖3 可見，夏季6—9 月降雨較多，云層覆蓋頻繁，導(dǎo)致可用圖片數(shù)量較少；可用圖片的年際分布也不均。水邊線提取需要在同一水位（擬定為85 高程5 m 水位，接近南京河段造床流量45 000 m3/s 對(duì)應(yīng)水位5.31m）下，南京河段岸坡坡比一般為1∶5，按Landsat 衛(wèi)星圖片30 m 的分辨率，小于6 m 的水位差反映在圖片上不超過1 個(gè)像素；同理，五大洲灘的坡比一般為1∶13，只要水位差小于2.3 m，則誤差不超過1 個(gè)像素。以2015 年為例，圖4 給出了南京潮位站的日均水位，從2015 年可用遙感影像中篩選出4 月和10 月各1 張可用照片，其拍攝日期的水位與擬定5 m 水位差為±1 m，滿足河道及洲灘的水邊線提取要求。同理，從每年圖片中選取2 幅組成新的序列，共74 幅，其中最大的水位差為2 m，滿足計(jì)算精度要求。

圖3 可用遙感影像的時(shí)域分布特征（總數(shù)289）Fig. 3 Temporal distribution characteristics of available remote sensing images (total 289)

圖4 2015 年南京潮水位站日均水位（以水位5 m±1 m 篩選遙感影像）Fig. 4 Daily average water level at Nanjing tidal level station in 2015 (filtering at water level 5 m±1 m)

圖5 為最終篩選出的歷年遙感影像拍攝時(shí)間。1984—2020 年，篩選的影像時(shí)間基本分布在每年的3—4 月和10—11 月，只有極個(gè)別年份的可用影像拍攝時(shí)間在2、7 和12 月。

圖5 篩選出的遙感影像（總數(shù)74，每年2 幅）Fig. 5 Filtered remote-sensing images (total 74, two images per year)

2.2 河床平面演變特征

圖6 給出了近40 年南京段河道水面面積隨時(shí)間變化情況，河道面積提取相對(duì)誤差均值為2.0%。圖中的水面面積指南京河段河床邊界包圍的區(qū)域，包括了洲灘面積。近40 年來，南京河段的水面面積緩慢下降，從最大的340 km2下降至320 km2，減幅約為6%。說明人類活動(dòng)不斷侵蝕河流，使河床邊界向內(nèi)收縮。

圖6 南京河段水域面積隨時(shí)間變化（相對(duì)誤差均值2.0%）Fig. 6 Temporal variation of mainstream area in the Nanjing reach in recent 40 years

圖7 為南京河段河床邊界的均值和方差。近40 年來，貓子山對(duì)岸、南京青奧體育公園、太子山公園及長江村附近的河床邊界出現(xiàn)了較大改變，尤以長江村的河床邊界變化最為劇烈，減小面積約10.35 km2，占近40 年總變化面積的52%，主要由水產(chǎn)養(yǎng)殖和圍河造田導(dǎo)致。

圖7 近40 年南京河段河床邊界統(tǒng)計(jì)參數(shù)Fig. 7 Statistical parameters of the river boundary sequence of the Nanjing reach in last 40 years

利用本征正交分解法提取了南京河段河床邊界序列的主要含能模態(tài)。由圖8 可知，一階模態(tài)含能比例很高（97%），貢獻(xiàn)了絕大多數(shù)能量；二階及更高階模態(tài)占比很少。圖9 僅給出了前2 階模態(tài)，可見一階模態(tài)反映了河段河床邊界的時(shí)均形態(tài)，二階模態(tài)反映南京段河床邊界在貓子山對(duì)岸和長江村附近出現(xiàn)了較大變化。

圖8 南京河段河床邊界序列各階模態(tài)含能Fig. 8 Energy content of POD modes of the river boundary sequence of the Nanjing reach

圖9 南京河段河床邊界序列的主要含能模態(tài)（0 為非河床，1 為河床）Fig. 9 The most energetic POD modes (1st-2nd) in river boundary sequence of the Nanjing reach

圖10 為南京段河床邊界序列的前兩階模態(tài)投影系數(shù)隨時(shí)間的變化?？芍浑A模態(tài)的投影系數(shù)隨時(shí)間基本不變，數(shù)值約為600，表明一階模態(tài)反映恒定不變的空間特征，對(duì)應(yīng)河床邊界的時(shí)均形態(tài)。二階模態(tài)的投影系數(shù)隨時(shí)間遞減，表明在1984—2005 年間長江村附近的河床邊界不斷向內(nèi)收縮，并在2005 年后基本保持不變。投影系數(shù)反映了長江村附近河床邊界變化的時(shí)間演變特征。

圖10 一階及二階模態(tài)投影系數(shù)隨時(shí)間的變化Fig. 10 Variation of projection coefficients of 1st and 2nd modes with time

2.3 洲灘平面演變特征

圖11 為南京河段近40 年的洲灘面積隨時(shí)間變化情況，八卦洲、梅子洲、新潛洲、新濟(jì)洲及新生洲的面積提取相對(duì)誤差分別為1.5%、4.4%、6.8%、4.7%和4.3%。八卦洲、梅子洲及新濟(jì)洲較為穩(wěn)定、變化較小；新生洲變化較大，表現(xiàn)為洲頭消退和洲尾淤長，但兩者程度相近，導(dǎo)致面積變化很小；新潛洲洲尾淤長嚴(yán)重，面積增幅200%，但此處幾乎無人類活動(dòng)，洲尾淤長主要由水沙變化所致。

圖11 近40 年南京河段洲灘面積隨時(shí)間變化Fig. 11 Temporal variation of sandbar areas in the Nanjing reach in recent 40 years

由于洲灘數(shù)量較多，為避免文章冗余，僅以變化較大的新潛洲為例進(jìn)行分析。由圖12 可知，近40年內(nèi)，新潛洲的洲頭變化較小，洲尾出現(xiàn)劇烈變化，洲尾增長甚至超過了原洲島長度的1.5 倍。

圖12 近40 年新潛洲形態(tài)統(tǒng)計(jì)參數(shù)Fig. 12 Statistical parameters of the Xinqian Sandbar shape sequence in last 40 years

利用本征正交分解法提取新潛洲時(shí)間序列的主要含能模態(tài)，一階模態(tài)含能比例較高為83%，二階至四階模態(tài)含能比例在1%~10%，更高階模態(tài)含能比例均在1%以下。圖13 僅給出前四階模態(tài)，可見模態(tài)分析法可以很好地提取不同時(shí)段新潛洲的淤長特征，其中一階模態(tài)反映了洲灘的初始形態(tài)，三階、二階及四階模態(tài)分別反映了洲尾的初期、中期及后期增長規(guī)律。

圖13 新潛洲形態(tài)序列的主要含能模態(tài)（0 為非洲灘，1 為洲灘）Fig. 13 The most energetic POD modes (1st-4th) in shape sequence of the Xinqian Sanbar

圖14 為新潛洲形態(tài)序列的前四階模態(tài)的投影系數(shù)隨時(shí)間的變化情況?？芍浑A模態(tài)的投影系數(shù)隨時(shí)間變化較小，表明一階模態(tài)反映1 個(gè)較為恒定的空間特征，對(duì)應(yīng)時(shí)均形態(tài)。三階模態(tài)的投影系數(shù)在1984—2000 年呈遞增趨勢(shì)，此后基本不變，反映出洲尾初期增長的時(shí)間演變特征。二階模態(tài)的投影系數(shù)在2000—2005 年呈遞增趨勢(shì)，此后改變很小，反映出洲尾中期的增長時(shí)間特征。四階模態(tài)的投影系數(shù)在2005—2020 年呈遞增趨勢(shì)，反映出洲尾后期增長的時(shí)間演變特征。綜合可知，新潛洲在近40 年內(nèi)始終處于洲頭穩(wěn)定、洲尾淤長的狀態(tài)。

圖14 一階至四階模態(tài)投影系數(shù)隨時(shí)間的變化Fig. 14 Variation of projection coefficients of 1st-4th modes with time

3 結(jié) 語

本文利用收集的1984—2020 年間Landsat 衛(wèi)星的南京段衛(wèi)星遙感圖片，通過目視解譯及水位（拍攝日期）篩選，挑選出74 幅圖片（每年2 張）。以85 高程5 m 水位為水邊線提取基準(zhǔn)，采用歸一化差分水體指數(shù)來凸顯河段水邊線信息，采用形態(tài)學(xué)方法提取河床邊界和內(nèi)河洲灘，并利用本征正交分解法分析兩者的空間與時(shí)間演變特征，其中含能模態(tài)反映了河道的空間演變特征，投影系數(shù)反映了時(shí)間演變特征，主要結(jié)論如下：

（1）近40 年來，長江南京段水域面積從最大的340 km2下降至320 km2，減幅6%，河床邊界改變最大的區(qū)域位于長江村附近，此處邊界向內(nèi)收縮范圍占近40 年總變化面積的52%，主要由水產(chǎn)養(yǎng)殖和圍河造田導(dǎo)致。本征正交分解一階模態(tài)反映了河床邊界的時(shí)均形態(tài)，二階模態(tài)反映了貓子山對(duì)岸和長江村附近河床邊界的收縮變化，投影系數(shù)表明該變化主要出現(xiàn)在1984—2005 年間。

（2）長江南京段的八卦洲、梅子洲及新濟(jì)洲的洲灘面積變化較小，均值分別為67、18、10 km2；新生洲洲頭消退和洲尾淤長的面積相近，導(dǎo)致面積變化較小，為11 km2；而新潛洲自洲頭向下游逐年淤長，增加的長度甚至超過了原洲島長度的1.5 倍，面積為6 km2，總增幅200%。一階模態(tài)反映了洲灘的初始形態(tài)，三階、二階及四階模態(tài)分別反映了洲尾的初期、中期及后期增長規(guī)律。模態(tài)投影系數(shù)可知洲尾初期、中期及后期增長分別發(fā)生在1984—2000 年、2000—2005 年及2005—2020 年。

（3）本文應(yīng)用的本征正交分解方法可以從看似無序的河床演變時(shí)間序列中提取其主要空間演變特征，并刻畫這些主要空間演變特征隨時(shí)間的演化規(guī)律。鑒于目前中等分辨率衛(wèi)星圖像的易獲取性，對(duì)于其他河流，特別是偏遠(yuǎn)地區(qū)缺少實(shí)測(cè)資料的河流，應(yīng)用本文方法可以很好地研究其時(shí)空演變特征。