基于DEM與河網(wǎng)密度的水系提取與應(yīng)用
——以瀾湄流域?yàn)槔?/h1>
2022-12-08 06:18:02包紅軍張雨鳳
水力發(fā)電 2022年11期
曹 爽，王 蒙，包紅軍，張雨鳳
(1.國(guó)家氣象中心，北京 100081；2.中國(guó)氣象局-河海大學(xué)水文氣象研究聯(lián)合實(shí)驗(yàn)室，北京 100081；3.江蘇省水文水資源勘測(cè)局蘇州分局，江蘇 蘇州 215011)
0 引 言
流域水系和流域邊界能有效反映流域地形地貌特征，是流域面雨量監(jiān)測(cè)預(yù)報(bào)和分布式水文模型構(gòu)建的重要基礎(chǔ)數(shù)據(jù)，是數(shù)字孿生流域建設(shè)的主要依據(jù)[1]。它們也是完成面雨量實(shí)況監(jiān)測(cè)、預(yù)報(bào)預(yù)測(cè)及預(yù)警服務(wù)的關(guān)鍵環(huán)節(jié)和前提條件，數(shù)據(jù)精度對(duì)相關(guān)計(jì)算和模擬結(jié)果有著直接影響[2]。隨著地理信息系統(tǒng)(GIS)和數(shù)字高程模型(DEM)的迅猛發(fā)展，國(guó)內(nèi)外不斷有學(xué)者對(duì)提高提取流域河網(wǎng)精度和子流域劃分進(jìn)行深入研究和改進(jìn)，以期得到最接近真實(shí)河網(wǎng)的數(shù)字化水系。一些學(xué)者[3-6]以DEM數(shù)據(jù)源為基礎(chǔ)通過優(yōu)化改進(jìn)算法提出了移動(dòng)窗口法、坡面徑流模擬法、谷線搜索法、從DEM直接提取河網(wǎng)與劃分子流域的方法等；還有一些學(xué)者[7-12]引入矢量河網(wǎng)作為DEM控制條件，通過Agree算法、Burn-in算法、Stream Burning算法等對(duì)主要河流進(jìn)行高程與流向修正，從而解決平原區(qū)和受人類活動(dòng)干擾的城市化的水系提取問題；另外，設(shè)置合理的網(wǎng)格數(shù)和閾值范圍也是控制河網(wǎng)精度的決定性因素之一。有大量研究表明[13-19]，可以通過河網(wǎng)密度法、水系分形維數(shù)法、均值變點(diǎn)分析法、河道平均坡降法、流域?qū)挾确植挤ǖ确椒ù_定最佳集水面積閾值；在子流域劃分方面，有些學(xué)者提出了利用河流分汊和拓?fù)潢P(guān)系[20]、多閾值虛擬河網(wǎng)融合技術(shù)[21]等劃分子流域，還提出了考慮高山流域[22]、湖泊水庫(kù)范圍[23]、山區(qū)平原地貌差異[24]、反映地表下墊面類型特征[25]的子流域劃分和編碼方法。
在前人研究成果的基礎(chǔ)上，本文以瀾滄江-湄公河流域(以下簡(jiǎn)稱“瀾湄流域”)為研究對(duì)象，從遙感影像獲取矢量河網(wǎng)數(shù)據(jù)，采用Agree算法、應(yīng)用Arcgis軟件疊加矢量河網(wǎng)從而修正主要河道DEM高程數(shù)據(jù)；分析河網(wǎng)密度與集水面積閾值的變化關(guān)系，以曲線割線斜率法定量確定最佳集水面積閾值，提取流域水系、邊界并劃分子流域；基于河網(wǎng)“套合差”法進(jìn)行精度評(píng)價(jià)，驗(yàn)證提取河網(wǎng)的可靠性；再以一次降雨過程為例分析流域面雨量的空間分布，將結(jié)果應(yīng)用于流域面雨量監(jiān)測(cè)、預(yù)報(bào)服務(wù)業(yè)務(wù)，為水文氣象預(yù)報(bào)預(yù)警工作提供技術(shù)支撐。
1 研究區(qū)概況與數(shù)據(jù)
本文以瀾湄流域作為研究區(qū)，流域范圍介于9°26′～34°2′ N，93°46′～108°52′ E，河長(zhǎng)4 900 km左右，流域面積約81萬(wàn)km2，是東南亞地區(qū)最重要也是最大的國(guó)際河流。該流域發(fā)源于中國(guó)青海省唐古拉山脈，流經(jīng)中國(guó)、緬甸、老撾、泰國(guó)、柬埔寨、越南等六國(guó)，以云南省南臘河口為界，中國(guó)境內(nèi)河段稱為瀾滄江，境外河段稱為湄公河，終由越南胡志明市流入南海。瀾湄流域地形地勢(shì)呈北高南低態(tài)勢(shì)，高度落差大，表現(xiàn)出很強(qiáng)的垂直地帶性，高程隨緯度減小逐漸下降，地貌由北向南依次是高山峽谷區(qū)、中低山寬谷區(qū)和沖積平原，海拔跨度從3 500～5 000 m降至1 000～3 000 m，再降至1 000 m以下。受地形地貌影響，瀾滄江-湄公河表現(xiàn)為典型的南北向狹長(zhǎng)型河流，流域形狀似“帚狀”；上游支流短小且少，呈 “樹枝狀”、“羽狀”水系；中游東岸支流水系發(fā)育較好且西岸少有大支流，呈“梳狀”水系；下游多為平原和三角洲，河網(wǎng)特別發(fā)育，呈“辮狀”、“格狀”水系特征[26-27]，見圖1。此外，在氣候變化和人類活動(dòng)的共同作用下，瀾湄流域頻發(fā)旱澇災(zāi)害[28]，坡地山區(qū)易發(fā)山洪常伴有泥石流，河谷壩區(qū)、平原和三角洲的耕地城鎮(zhèn)被洪水侵襲，瀾滄江中游常有“焚風(fēng)”效應(yīng)導(dǎo)致干旱河谷等災(zāi)害事件。
圖1 瀾湄流域概況
本文采用NASA SRTM的空間分辨率為90 m的DEM數(shù)據(jù)，以及國(guó)家基礎(chǔ)地理信息中心提供的全球1∶1 000 000包含行政區(qū)劃、湖泊、水系等信息的矢量數(shù)據(jù)，而降水?dāng)?shù)據(jù)采用中央氣象臺(tái)5 km×5 km分辨率的智能網(wǎng)格降水。
2 研究方法
2.1 河網(wǎng)提取方法
2.1.1 河網(wǎng)修正Agree算法
DEM數(shù)據(jù)是高程信息的反演產(chǎn)品，其本身就存在一定誤差。在提取水系研究河網(wǎng)分布時(shí)無(wú)法綜合考慮地形地貌、雨水沖刷、人工河道等其他因素的共同作用，會(huì)直接影響河網(wǎng)提取精度。尤其是在地勢(shì)平坦、受人類活動(dòng)影響大的地區(qū)表現(xiàn)明顯[9,11]。
瀾湄流域下游地區(qū)即為地勢(shì)落差小的沖積平原，且耕地、人工開鑿水渠、城市化等人類活動(dòng)影響大。以DEM作為單一要素提取的河網(wǎng)水系與實(shí)際分布會(huì)存在較大差異，遂引入Agree算法修正DEM。該算法原理是通過疊加矢量河網(wǎng)，降低實(shí)際河網(wǎng)所在柵格的高程，控制水流方向，增加提取河網(wǎng)柵格的“匯流”能力。
算法主要步驟包括：①通過衛(wèi)星遙感影像解譯、Google earth清繪河道中線或國(guó)家基礎(chǔ)地理信息中心提供的地理信息數(shù)據(jù)等途徑獲取所需的實(shí)際主要河網(wǎng)數(shù)字化水系的矢量文件；②疊加矢量河網(wǎng)和填洼后的DEM，以矢量線要素為中心做緩沖區(qū)，緩沖區(qū)鄰域半徑的設(shè)置應(yīng)不小于1/2個(gè)DEM分辨率；③運(yùn)用轉(zhuǎn)換工具將矢量緩沖區(qū)范圍轉(zhuǎn)為柵格文件，進(jìn)而運(yùn)算降低重疊部分DEM柵格高程，獲取修正后DEM，命名為DEM-A。
2.1.2 流域水系提取
基于DEM的水系和流域邊界提取技術(shù)路線見圖2。工作原理是基于地表徑流漫流模型，關(guān)鍵步驟包括：迭代計(jì)算洼地填平；D8單流向算法確定水流流向；沿水流流向計(jì)算每個(gè)柵格單元的上游匯流能力即匯流累積量；設(shè)置匯流閾值，匯流量不小于閾值的柵格是潛在河網(wǎng)；根據(jù)Strahler水系分級(jí)法對(duì)河流進(jìn)行分級(jí)；捕捉傾斜點(diǎn)和計(jì)算分水嶺，獲取自然全流域和自然子流域出口和流域邊界；水系河網(wǎng)和流域邊界矢量化[14]。
圖2 水系和流域邊界提取技術(shù)路線
2.2 確定最佳集水閾值方法
在流域提取過程中有一個(gè)重要的變量參數(shù)是集水面積閾值。閾值的變化控制著水系河網(wǎng)的疏密程度，目視判讀的方式存在一定的主觀性，計(jì)算結(jié)果難以達(dá)成統(tǒng)一。于是，本文采用河網(wǎng)密度法，運(yùn)用擬合指數(shù)函數(shù)曲線的割線所對(duì)應(yīng)的斜率與擬合曲線相切的數(shù)學(xué)方法來確定切點(diǎn)，該切點(diǎn)的物理意義即為河網(wǎng)密度變化趨于平穩(wěn)時(shí)的最佳集水面積閾值。
擬合集水面積閾值與河網(wǎng)密度關(guān)系曲線，構(gòu)建集水面積閾值與河網(wǎng)密度關(guān)系方程f(x)，定義域?yàn)閇x1，x2]，值域?yàn)閇f(x2)，f(x1)]，令f(x)的一階導(dǎo)數(shù)為f′(x)。在定義域內(nèi)，曲線割線的斜率k=[f(x1)-f(x2)]/(x1-x2)，以該斜率作為擬合曲線的切線，即k=f′(x0)，求解x0，切線與擬合曲線的切點(diǎn)(x0，y0)代表了河網(wǎng)密度隨集水面積閾值由劇烈變化變?yōu)槠骄徸兓霓D(zhuǎn)折點(diǎn)。即，x0為目標(biāo)閾值[16]。
2.3 河網(wǎng)精度評(píng)價(jià)
本文對(duì)河網(wǎng)精度評(píng)價(jià)分為兩個(gè)方面：一方面，檢驗(yàn)提取的水系和數(shù)字化水系主要河道河長(zhǎng)的相對(duì)誤差，有效反映偏離真值的實(shí)際大?。涣硪环矫?，引入“河網(wǎng)套合差”的概念，指的是疊加提取水系和數(shù)字化水系兩部分，兩水系由于位置偏移會(huì)形成細(xì)碎的多邊形，計(jì)算這些細(xì)碎多邊形面積占流域總面積的比值即為套合差，該比值越小代表兩條水系的位置偏差越小[14,29](見圖3)。
圖3 河網(wǎng)套合差示意
2.4 流域面雨量計(jì)算
流域面雨量是指某一流域或區(qū)域整個(gè)面上的平均降雨量。它能客觀地反映流域的降水情況，在水文模型和水情預(yù)報(bào)等工作中應(yīng)用廣泛。本文的降水資料數(shù)據(jù)采用的是中央氣象臺(tái)5 km×5 km網(wǎng)格的智能網(wǎng)格預(yù)報(bào)降水產(chǎn)品，以算術(shù)平均法估算流域面雨量預(yù)報(bào)結(jié)果適用可行，計(jì)算簡(jiǎn)單。具體計(jì)算原理及公式不再詳述[30]。
3 結(jié)果與分析
3.1 集水面積閾值與水系河網(wǎng)的變化關(guān)系
本文選擇1 000、10 000、15 000、20 000、25 000、30 000、40 000、50 000、60 000、100 000共10個(gè)河網(wǎng)柵格數(shù)為閾值樣本，在提取水系河網(wǎng)的過程中以此判斷匯流累積量柵格數(shù)量的可變參數(shù)影響河網(wǎng)疏密程度。不同集水面積閾值對(duì)水系提取的影響特征見表1。
表1 受集水面積閾值影響提取水系河網(wǎng)特征變化
由表1可看出，隨著閾值面積由8.1 km2逐漸增加至810 km2的變化過程中，河道總數(shù)從44 429條減少為433條，河源數(shù)從22 261個(gè)減少為217個(gè)，總河長(zhǎng)從198 854.57 km減至23 363.51 km，河道級(jí)別從八級(jí)河流減少至五級(jí)河流；在集水面積閾值擴(kuò)大了100倍的情況下，河道總數(shù)、河源數(shù)相應(yīng)減少為1/100，總河長(zhǎng)則減少成了近原來的1/9，流域面積雖受到一定影響有一些減少但變化不大。另外，流域河網(wǎng)分級(jí)結(jié)果也存在一定變化規(guī)律。在同一集水面積閾值下，河道數(shù)和河長(zhǎng)隨河網(wǎng)級(jí)別增加逐漸減少且減小幅度逐漸變緩，一級(jí)河道數(shù)占河道總數(shù)50%左右，二級(jí)河道數(shù)和河長(zhǎng)是一級(jí)河道數(shù)和河長(zhǎng)的1/2，一級(jí)和二級(jí)的河長(zhǎng)和河道數(shù)在總數(shù)的占比達(dá)到75%左右。
3.2 流域最佳閾值的確定與精度評(píng)價(jià)
依據(jù)表1中河網(wǎng)閾值柵格數(shù)對(duì)應(yīng)的河網(wǎng)密度繪制圖4關(guān)系曲線。由圖4可知，隨閾值增大，河網(wǎng)密度逐漸減小，閾值范圍為103～104，最大減小幅度呈斷崖式；閾值范圍為104～4×104，仍呈減小趨勢(shì)但下降幅度逐漸變緩；閾值范圍為4×104～105，變化逐漸趨于穩(wěn)定。
圖4 河網(wǎng)柵格數(shù)閾值—河網(wǎng)密度關(guān)系曲線
本文引入多種函數(shù)關(guān)系擬合河網(wǎng)柵格數(shù)閾值和河網(wǎng)密度，最終確定采用擬合效果最好的冪函數(shù)關(guān)系，得到擬合方程
y=6.165 7x-0.465(R2=0999 3)
(1)
式中，y為河網(wǎng)密度，km/km2；x為河網(wǎng)柵格數(shù)閾值；R2為相關(guān)系數(shù)，代表擬合程度，R2=0.999 3表明擬合程度高，兩者相關(guān)性強(qiáng)。
如圖5所示，連接閾值樣本首尾兩端點(diǎn)做割線，以割線斜率為擬合曲線切線斜率做切線，求得切點(diǎn)即為擬合曲線由劇烈到平緩變化的轉(zhuǎn)折點(diǎn)。對(duì)擬合方程(1)求導(dǎo)，得到擬合曲線一階導(dǎo)數(shù)
圖5 最佳閾值的確定
k=y′=-2.867x-1.465
(2)
進(jìn)一步求解方程(2)與斜率k的關(guān)系，獲得切點(diǎn)坐標(biāo)為(15 933，0.068)，此時(shí)河網(wǎng)柵格數(shù)閾值15 933為提取水系的最佳閾值，即最佳集水面積閾值約為129.1 km2，以該閾值提取水系結(jié)果如圖6所示。切點(diǎn)縱坐標(biāo)0.068其物理意義為，擬合曲線上河網(wǎng)柵格數(shù)閾值15 933對(duì)應(yīng)的河網(wǎng)密度是0.068 km/km2；利用Arcgis水文分析工具以15 933為閾值提取水系，計(jì)算得到河網(wǎng)密度為0.069 km/km2，相對(duì)誤差為1.5%，擬合效果好。
圖6 90 m分辨率下瀾湄流域最佳集水面積閾值提取水系
為驗(yàn)證提取水系精度，反映提取水系與矢量化數(shù)字水系偏移程度，對(duì)河流總長(zhǎng)的相對(duì)誤差和河網(wǎng)“套合差”兩個(gè)指標(biāo)進(jìn)行計(jì)算分析。提取水系河流總長(zhǎng)52 061.6 km，數(shù)字水系河流總長(zhǎng)48 395.5 km，相對(duì)誤差為7.5%?；贏rcGIS圖層加載1∶1 000 000矢量化數(shù)字水系圖作為標(biāo)準(zhǔn)做檢驗(yàn)，對(duì)比數(shù)字化水系和提取水系之間的位置偏移，統(tǒng)計(jì)由此偏移產(chǎn)生的兩水系之間的細(xì)碎多邊形面積(見圖7)。統(tǒng)計(jì)結(jié)果顯示，整個(gè)瀾湄流域內(nèi)的細(xì)碎多邊形面積約為20 000 km2，河網(wǎng)套合差為2.5%，小于3%，吻合程度較好。由圖7a可見，上游地區(qū)海拔較高、地勢(shì)變化明顯且無(wú)較大湖泊，兩水系疊加形成的多邊形相對(duì)較少，多邊形面積約為2 100 km2，對(duì)應(yīng)的流域面積為166 000 km2，河網(wǎng)套河差為1.3%；圖7b顯示為流域中下游地區(qū)，該范圍內(nèi)地勢(shì)逐漸變緩，水系豐富，在柬埔寨境內(nèi)有較大湖泊洞里薩湖，兩水系之間的多邊形較上游明顯增多且面積增大，計(jì)算得中下游的河網(wǎng)套河差為2.9%，上游和中下游的河網(wǎng)套河差結(jié)果均在有效范圍內(nèi)，但中下游地區(qū)河網(wǎng)偏移程度較上游大，提取水系精度略差于上游地區(qū)?？傮w而言，本文提取水系結(jié)果與實(shí)際水系較為吻合，能有效、客觀地反映瀾湄流域水系整體特征，可為后續(xù)研究提供數(shù)據(jù)支撐。
圖7 提取水系與數(shù)字化水系河網(wǎng)套河差
3.3 河網(wǎng)分級(jí)和子流域劃分
由90 m分辨率高程數(shù)據(jù)提取最優(yōu)集水面積閾值的瀾湄流域水系情況如圖6所示，河網(wǎng)分為1～6級(jí)，一級(jí)為河道最低級(jí)別的河源水系，共計(jì)1 342條，占總河道數(shù)的50%，一級(jí)河網(wǎng)發(fā)育系數(shù)約為2.9，河系不均勻系數(shù)約為1.3。
基于Arcgis自動(dòng)化提取流域出口，共劃分小流域3 259個(gè)，流域面積普遍在50～900 km2范圍。為滿足業(yè)務(wù)應(yīng)用需要，依據(jù)水系上下游匯流關(guān)系，建立河道水系與小流域以及各小流域間的拓?fù)潢P(guān)系；參考流域內(nèi)地形特征以及湖泊范圍，滿足水系連續(xù)性、分水線完整性、出水口準(zhǔn)確性等關(guān)鍵條件；對(duì)流域面積在5 km2以下的微小流域或中間匯流區(qū)，以匯流關(guān)系為主要依據(jù)合并到相鄰小流域內(nèi)[31-32]；另外，又以流域內(nèi)干流上的8個(gè)主要水文站點(diǎn)作為流域出口對(duì)其進(jìn)行流域劃分；進(jìn)而完成人工合并小流域的工作。將瀾湄流域的子流域劃分出3個(gè)等級(jí)，根據(jù)子流域面積由小到大依次定義為一級(jí)、二級(jí)、三級(jí)流域。一級(jí)流域共計(jì)9個(gè)子流域，子流域面積最小的是昌都—舊州區(qū)間，約為1.6萬(wàn)km2，最大的是穆達(dá)漢—上丁區(qū)間，約為24萬(wàn)km2。二級(jí)流域共計(jì)61個(gè)子流域，子流域面積小于5 000 km2的占比25%，在5 000～10 000 km2的占23%，在10 000～20 000 km2的占33%，大于20 000 km2的占19%。三級(jí)流域共計(jì)328個(gè)子流域，子流域面積普遍控制在500～5 000 km2，占比達(dá)到90%。其中，小于500 km2的為不適合與相鄰合并的中間匯流區(qū)，另存在一個(gè)大于5 000 km2的子流域，受洞里薩湖水域限制，包含湖域區(qū)域，子流域面積約6 700 km2。
3.4 不同子流域尺度流域面雨量預(yù)報(bào)分析
選取2021年7月24日～25日的一次降水過程預(yù)報(bào)為例，即為24日8時(shí)起報(bào)的24 h面雨量，分析不劃分子流域的瀾湄流域以及劃分出一級(jí)流域、二級(jí)流域的面雨量結(jié)果對(duì)此次過程在空間分布上的響應(yīng)。未劃分子流域的全流域經(jīng)算術(shù)平均計(jì)算面雨量為25 mm，全流域面積大僅以一個(gè)數(shù)值表達(dá)面雨量，在降水分布均勻且降雨量較小時(shí)可起到一定參考作用，但在降雨空間分布不均時(shí)無(wú)法體現(xiàn)降雨中心。圖8a中一級(jí)流域面雨量計(jì)算結(jié)果顯示，昌都—舊州和穆達(dá)漢—上丁面雨量量級(jí)為小雨，舊州—允景洪量級(jí)為中雨，允景洪—清盛、清盛—瑯勃拉邦和瑯勃拉邦—萬(wàn)象段量級(jí)為大雨，舊州—允景洪和萬(wàn)象—穆達(dá)漢量級(jí)達(dá)到暴雨。圖8b顯示，二級(jí)流域進(jìn)一步縮小子流域面積，面雨量結(jié)果與智能網(wǎng)格預(yù)報(bào)降雨結(jié)果在暴雨中心的空間體現(xiàn)上更為一致，降雨較大范圍集中于中下游地區(qū)，泰國(guó)中東部和老撾中南部交界的位置，有6個(gè)子流域面雨量超過60 mm，量級(jí)為大暴雨，最大在泰國(guó)境內(nèi)達(dá)到143 mm，與一級(jí)流域萬(wàn)象—穆達(dá)漢段暴雨級(jí)面雨量相呼應(yīng)。
圖8 面雨量空間分布
4 結(jié)論與討論
本文分析了瀾湄流域集水面積閾值與水系特征變化的數(shù)量關(guān)系，建立了集水面積閾值與河網(wǎng)密度之間的冪函數(shù)，通過曲線割線斜率法客服主觀性，確定了90 m分辨率高程下的最佳集水面積閾值為129.1km2。并以1∶1 000 000的矢量化數(shù)字水系圖做驗(yàn)證依據(jù)，總河長(zhǎng)相對(duì)誤差為7.5%，河網(wǎng)套河差為2.5%，均說明提取水系與實(shí)際情況吻合程度較高，提取水系有效，進(jìn)一步劃分了瀾湄流域1～6級(jí)的水系和3個(gè)等級(jí)的子流域?；谥悄芫W(wǎng)格預(yù)報(bào)降水，以一次降水過程展現(xiàn)面雨量在不同等級(jí)子流域上的空間分布，子流域等級(jí)越大、劃分越細(xì)致對(duì)暴雨中心的體現(xiàn)越明顯。本文的研究成果，可在水文氣象業(yè)務(wù)上為瀾湄流域的面雨量監(jiān)測(cè)、預(yù)報(bào)和水文模型的構(gòu)建提供參考，為數(shù)字孿生流域建設(shè)與流域防洪減災(zāi)奠定基礎(chǔ)。
雖然本文對(duì)DEM進(jìn)行了修正工作，但是在處理較大湖泊、入?？谌侵薜忍厥獾匦紊吓c真實(shí)情況方面仍存在差異；研究區(qū)域面積大，在高精度高程數(shù)據(jù)情況下運(yùn)行速度慢，選取的集水閾值數(shù)據(jù)量有限。下一步工作可以考慮按照地形、水系類型等因素先將整個(gè)流域劃分為幾個(gè)不同的小區(qū)域，以求在保證高程精度情況下提高運(yùn)算效率，進(jìn)而綜合模擬流域水系。
猜你喜歡
基于小世界網(wǎng)絡(luò)的海河流域河網(wǎng)結(jié)構(gòu)及功能響應(yīng)水科學(xué)進(jìn)展(2023年4期)2023-10-07 11:23:44 寧夏紅柳溝流域水沙變化及產(chǎn)沙分析農(nóng)業(yè)科學(xué)研究(2022年2期)2022-08-01 05:25:30 纖維素基超疏水材料的制備與應(yīng)用研究進(jìn)展中國(guó)造紙(2022年2期)2022-04-04 16:40:13 板栗樹下打坑集水技術(shù)河北果樹(2020年1期)2020-02-09 12:31:40 基于PSR模型的上海地區(qū)河網(wǎng)脆弱性探討華東師范大學(xué)學(xué)報(bào)(自然科學(xué)版)(2019年3期)2019-06-24 05:29:09 基于小波去噪的稱重雨量數(shù)據(jù)分析成都信息工程大學(xué)學(xué)報(bào)(2019年1期)2019-05-20 09:14:26 集水明排與深井井點(diǎn)聯(lián)合降水工程應(yīng)用中國(guó)港灣建設(shè)(2017年11期)2017-12-19 12:27:05 不同引水水源對(duì)平原河網(wǎng)影響分析水利科技與經(jīng)濟(jì)(2017年5期)2017-04-22 02:39:28 SL—1 型雨量傳感器故障分析排除和維護(hù)西藏科技(2016年5期)2016-09-26 12:16:40 MIKE11模型在城市河網(wǎng)生態(tài)調(diào)水工程中的應(yīng)用水利科技與經(jīng)濟(jì)(2016年9期)2016-04-22 01:07:22
曹 爽，王 蒙，包紅軍，張雨鳳
(1.國(guó)家氣象中心，北京 100081；2.中國(guó)氣象局-河海大學(xué)水文氣象研究聯(lián)合實(shí)驗(yàn)室，北京 100081；3.江蘇省水文水資源勘測(cè)局蘇州分局，江蘇 蘇州 215011)
0 引 言
流域水系和流域邊界能有效反映流域地形地貌特征，是流域面雨量監(jiān)測(cè)預(yù)報(bào)和分布式水文模型構(gòu)建的重要基礎(chǔ)數(shù)據(jù)，是數(shù)字孿生流域建設(shè)的主要依據(jù)[1]。它們也是完成面雨量實(shí)況監(jiān)測(cè)、預(yù)報(bào)預(yù)測(cè)及預(yù)警服務(wù)的關(guān)鍵環(huán)節(jié)和前提條件，數(shù)據(jù)精度對(duì)相關(guān)計(jì)算和模擬結(jié)果有著直接影響[2]。隨著地理信息系統(tǒng)(GIS)和數(shù)字高程模型(DEM)的迅猛發(fā)展，國(guó)內(nèi)外不斷有學(xué)者對(duì)提高提取流域河網(wǎng)精度和子流域劃分進(jìn)行深入研究和改進(jìn)，以期得到最接近真實(shí)河網(wǎng)的數(shù)字化水系。一些學(xué)者[3-6]以DEM數(shù)據(jù)源為基礎(chǔ)通過優(yōu)化改進(jìn)算法提出了移動(dòng)窗口法、坡面徑流模擬法、谷線搜索法、從DEM直接提取河網(wǎng)與劃分子流域的方法等；還有一些學(xué)者[7-12]引入矢量河網(wǎng)作為DEM控制條件，通過Agree算法、Burn-in算法、Stream Burning算法等對(duì)主要河流進(jìn)行高程與流向修正，從而解決平原區(qū)和受人類活動(dòng)干擾的城市化的水系提取問題；另外，設(shè)置合理的網(wǎng)格數(shù)和閾值范圍也是控制河網(wǎng)精度的決定性因素之一。有大量研究表明[13-19]，可以通過河網(wǎng)密度法、水系分形維數(shù)法、均值變點(diǎn)分析法、河道平均坡降法、流域?qū)挾确植挤ǖ确椒ù_定最佳集水面積閾值；在子流域劃分方面，有些學(xué)者提出了利用河流分汊和拓?fù)潢P(guān)系[20]、多閾值虛擬河網(wǎng)融合技術(shù)[21]等劃分子流域，還提出了考慮高山流域[22]、湖泊水庫(kù)范圍[23]、山區(qū)平原地貌差異[24]、反映地表下墊面類型特征[25]的子流域劃分和編碼方法。
在前人研究成果的基礎(chǔ)上，本文以瀾滄江-湄公河流域(以下簡(jiǎn)稱“瀾湄流域”)為研究對(duì)象，從遙感影像獲取矢量河網(wǎng)數(shù)據(jù)，采用Agree算法、應(yīng)用Arcgis軟件疊加矢量河網(wǎng)從而修正主要河道DEM高程數(shù)據(jù)；分析河網(wǎng)密度與集水面積閾值的變化關(guān)系，以曲線割線斜率法定量確定最佳集水面積閾值，提取流域水系、邊界并劃分子流域；基于河網(wǎng)“套合差”法進(jìn)行精度評(píng)價(jià)，驗(yàn)證提取河網(wǎng)的可靠性；再以一次降雨過程為例分析流域面雨量的空間分布，將結(jié)果應(yīng)用于流域面雨量監(jiān)測(cè)、預(yù)報(bào)服務(wù)業(yè)務(wù)，為水文氣象預(yù)報(bào)預(yù)警工作提供技術(shù)支撐。
1 研究區(qū)概況與數(shù)據(jù)
本文以瀾湄流域作為研究區(qū)，流域范圍介于9°26′～34°2′ N，93°46′～108°52′ E，河長(zhǎng)4 900 km左右，流域面積約81萬(wàn)km2，是東南亞地區(qū)最重要也是最大的國(guó)際河流。該流域發(fā)源于中國(guó)青海省唐古拉山脈，流經(jīng)中國(guó)、緬甸、老撾、泰國(guó)、柬埔寨、越南等六國(guó)，以云南省南臘河口為界，中國(guó)境內(nèi)河段稱為瀾滄江，境外河段稱為湄公河，終由越南胡志明市流入南海。瀾湄流域地形地勢(shì)呈北高南低態(tài)勢(shì)，高度落差大，表現(xiàn)出很強(qiáng)的垂直地帶性，高程隨緯度減小逐漸下降，地貌由北向南依次是高山峽谷區(qū)、中低山寬谷區(qū)和沖積平原，海拔跨度從3 500～5 000 m降至1 000～3 000 m，再降至1 000 m以下。受地形地貌影響，瀾滄江-湄公河表現(xiàn)為典型的南北向狹長(zhǎng)型河流，流域形狀似“帚狀”；上游支流短小且少，呈 “樹枝狀”、“羽狀”水系；中游東岸支流水系發(fā)育較好且西岸少有大支流，呈“梳狀”水系；下游多為平原和三角洲，河網(wǎng)特別發(fā)育，呈“辮狀”、“格狀”水系特征[26-27]，見圖1。此外，在氣候變化和人類活動(dòng)的共同作用下，瀾湄流域頻發(fā)旱澇災(zāi)害[28]，坡地山區(qū)易發(fā)山洪常伴有泥石流，河谷壩區(qū)、平原和三角洲的耕地城鎮(zhèn)被洪水侵襲，瀾滄江中游常有“焚風(fēng)”效應(yīng)導(dǎo)致干旱河谷等災(zāi)害事件。
圖1 瀾湄流域概況
本文采用NASA SRTM的空間分辨率為90 m的DEM數(shù)據(jù)，以及國(guó)家基礎(chǔ)地理信息中心提供的全球1∶1 000 000包含行政區(qū)劃、湖泊、水系等信息的矢量數(shù)據(jù)，而降水?dāng)?shù)據(jù)采用中央氣象臺(tái)5 km×5 km分辨率的智能網(wǎng)格降水。
2 研究方法
2.1 河網(wǎng)提取方法
2.1.1 河網(wǎng)修正Agree算法
DEM數(shù)據(jù)是高程信息的反演產(chǎn)品，其本身就存在一定誤差。在提取水系研究河網(wǎng)分布時(shí)無(wú)法綜合考慮地形地貌、雨水沖刷、人工河道等其他因素的共同作用，會(huì)直接影響河網(wǎng)提取精度。尤其是在地勢(shì)平坦、受人類活動(dòng)影響大的地區(qū)表現(xiàn)明顯[9,11]。
瀾湄流域下游地區(qū)即為地勢(shì)落差小的沖積平原，且耕地、人工開鑿水渠、城市化等人類活動(dòng)影響大。以DEM作為單一要素提取的河網(wǎng)水系與實(shí)際分布會(huì)存在較大差異，遂引入Agree算法修正DEM。該算法原理是通過疊加矢量河網(wǎng)，降低實(shí)際河網(wǎng)所在柵格的高程，控制水流方向，增加提取河網(wǎng)柵格的“匯流”能力。
算法主要步驟包括：①通過衛(wèi)星遙感影像解譯、Google earth清繪河道中線或國(guó)家基礎(chǔ)地理信息中心提供的地理信息數(shù)據(jù)等途徑獲取所需的實(shí)際主要河網(wǎng)數(shù)字化水系的矢量文件；②疊加矢量河網(wǎng)和填洼后的DEM，以矢量線要素為中心做緩沖區(qū)，緩沖區(qū)鄰域半徑的設(shè)置應(yīng)不小于1/2個(gè)DEM分辨率；③運(yùn)用轉(zhuǎn)換工具將矢量緩沖區(qū)范圍轉(zhuǎn)為柵格文件，進(jìn)而運(yùn)算降低重疊部分DEM柵格高程，獲取修正后DEM，命名為DEM-A。
2.1.2 流域水系提取
基于DEM的水系和流域邊界提取技術(shù)路線見圖2。工作原理是基于地表徑流漫流模型，關(guān)鍵步驟包括：迭代計(jì)算洼地填平；D8單流向算法確定水流流向；沿水流流向計(jì)算每個(gè)柵格單元的上游匯流能力即匯流累積量；設(shè)置匯流閾值，匯流量不小于閾值的柵格是潛在河網(wǎng)；根據(jù)Strahler水系分級(jí)法對(duì)河流進(jìn)行分級(jí)；捕捉傾斜點(diǎn)和計(jì)算分水嶺，獲取自然全流域和自然子流域出口和流域邊界；水系河網(wǎng)和流域邊界矢量化[14]。
圖2 水系和流域邊界提取技術(shù)路線
2.2 確定最佳集水閾值方法
在流域提取過程中有一個(gè)重要的變量參數(shù)是集水面積閾值。閾值的變化控制著水系河網(wǎng)的疏密程度，目視判讀的方式存在一定的主觀性，計(jì)算結(jié)果難以達(dá)成統(tǒng)一。于是，本文采用河網(wǎng)密度法，運(yùn)用擬合指數(shù)函數(shù)曲線的割線所對(duì)應(yīng)的斜率與擬合曲線相切的數(shù)學(xué)方法來確定切點(diǎn)，該切點(diǎn)的物理意義即為河網(wǎng)密度變化趨于平穩(wěn)時(shí)的最佳集水面積閾值。
擬合集水面積閾值與河網(wǎng)密度關(guān)系曲線，構(gòu)建集水面積閾值與河網(wǎng)密度關(guān)系方程f(x)，定義域?yàn)閇x1，x2]，值域?yàn)閇f(x2)，f(x1)]，令f(x)的一階導(dǎo)數(shù)為f′(x)。在定義域內(nèi)，曲線割線的斜率k=[f(x1)-f(x2)]/(x1-x2)，以該斜率作為擬合曲線的切線，即k=f′(x0)，求解x0，切線與擬合曲線的切點(diǎn)(x0，y0)代表了河網(wǎng)密度隨集水面積閾值由劇烈變化變?yōu)槠骄徸兓霓D(zhuǎn)折點(diǎn)。即，x0為目標(biāo)閾值[16]。
2.3 河網(wǎng)精度評(píng)價(jià)
本文對(duì)河網(wǎng)精度評(píng)價(jià)分為兩個(gè)方面：一方面，檢驗(yàn)提取的水系和數(shù)字化水系主要河道河長(zhǎng)的相對(duì)誤差，有效反映偏離真值的實(shí)際大?。涣硪环矫?，引入“河網(wǎng)套合差”的概念，指的是疊加提取水系和數(shù)字化水系兩部分，兩水系由于位置偏移會(huì)形成細(xì)碎的多邊形，計(jì)算這些細(xì)碎多邊形面積占流域總面積的比值即為套合差，該比值越小代表兩條水系的位置偏差越小[14,29](見圖3)。
圖3 河網(wǎng)套合差示意
2.4 流域面雨量計(jì)算
流域面雨量是指某一流域或區(qū)域整個(gè)面上的平均降雨量。它能客觀地反映流域的降水情況，在水文模型和水情預(yù)報(bào)等工作中應(yīng)用廣泛。本文的降水資料數(shù)據(jù)采用的是中央氣象臺(tái)5 km×5 km網(wǎng)格的智能網(wǎng)格預(yù)報(bào)降水產(chǎn)品，以算術(shù)平均法估算流域面雨量預(yù)報(bào)結(jié)果適用可行，計(jì)算簡(jiǎn)單。具體計(jì)算原理及公式不再詳述[30]。
3 結(jié)果與分析
3.1 集水面積閾值與水系河網(wǎng)的變化關(guān)系
本文選擇1 000、10 000、15 000、20 000、25 000、30 000、40 000、50 000、60 000、100 000共10個(gè)河網(wǎng)柵格數(shù)為閾值樣本，在提取水系河網(wǎng)的過程中以此判斷匯流累積量柵格數(shù)量的可變參數(shù)影響河網(wǎng)疏密程度。不同集水面積閾值對(duì)水系提取的影響特征見表1。
表1 受集水面積閾值影響提取水系河網(wǎng)特征變化
由表1可看出，隨著閾值面積由8.1 km2逐漸增加至810 km2的變化過程中，河道總數(shù)從44 429條減少為433條，河源數(shù)從22 261個(gè)減少為217個(gè)，總河長(zhǎng)從198 854.57 km減至23 363.51 km，河道級(jí)別從八級(jí)河流減少至五級(jí)河流；在集水面積閾值擴(kuò)大了100倍的情況下，河道總數(shù)、河源數(shù)相應(yīng)減少為1/100，總河長(zhǎng)則減少成了近原來的1/9，流域面積雖受到一定影響有一些減少但變化不大。另外，流域河網(wǎng)分級(jí)結(jié)果也存在一定變化規(guī)律。在同一集水面積閾值下，河道數(shù)和河長(zhǎng)隨河網(wǎng)級(jí)別增加逐漸減少且減小幅度逐漸變緩，一級(jí)河道數(shù)占河道總數(shù)50%左右，二級(jí)河道數(shù)和河長(zhǎng)是一級(jí)河道數(shù)和河長(zhǎng)的1/2，一級(jí)和二級(jí)的河長(zhǎng)和河道數(shù)在總數(shù)的占比達(dá)到75%左右。
3.2 流域最佳閾值的確定與精度評(píng)價(jià)
依據(jù)表1中河網(wǎng)閾值柵格數(shù)對(duì)應(yīng)的河網(wǎng)密度繪制圖4關(guān)系曲線。由圖4可知，隨閾值增大，河網(wǎng)密度逐漸減小，閾值范圍為103～104，最大減小幅度呈斷崖式；閾值范圍為104～4×104，仍呈減小趨勢(shì)但下降幅度逐漸變緩；閾值范圍為4×104～105，變化逐漸趨于穩(wěn)定。
圖4 河網(wǎng)柵格數(shù)閾值—河網(wǎng)密度關(guān)系曲線
本文引入多種函數(shù)關(guān)系擬合河網(wǎng)柵格數(shù)閾值和河網(wǎng)密度，最終確定采用擬合效果最好的冪函數(shù)關(guān)系，得到擬合方程
y=6.165 7x-0.465(R2=0999 3)
(1)
式中，y為河網(wǎng)密度，km/km2；x為河網(wǎng)柵格數(shù)閾值；R2為相關(guān)系數(shù)，代表擬合程度，R2=0.999 3表明擬合程度高，兩者相關(guān)性強(qiáng)。
如圖5所示，連接閾值樣本首尾兩端點(diǎn)做割線，以割線斜率為擬合曲線切線斜率做切線，求得切點(diǎn)即為擬合曲線由劇烈到平緩變化的轉(zhuǎn)折點(diǎn)。對(duì)擬合方程(1)求導(dǎo)，得到擬合曲線一階導(dǎo)數(shù)
圖5 最佳閾值的確定
k=y′=-2.867x-1.465
(2)
進(jìn)一步求解方程(2)與斜率k的關(guān)系，獲得切點(diǎn)坐標(biāo)為(15 933，0.068)，此時(shí)河網(wǎng)柵格數(shù)閾值15 933為提取水系的最佳閾值，即最佳集水面積閾值約為129.1 km2，以該閾值提取水系結(jié)果如圖6所示。切點(diǎn)縱坐標(biāo)0.068其物理意義為，擬合曲線上河網(wǎng)柵格數(shù)閾值15 933對(duì)應(yīng)的河網(wǎng)密度是0.068 km/km2；利用Arcgis水文分析工具以15 933為閾值提取水系，計(jì)算得到河網(wǎng)密度為0.069 km/km2，相對(duì)誤差為1.5%，擬合效果好。
圖6 90 m分辨率下瀾湄流域最佳集水面積閾值提取水系
為驗(yàn)證提取水系精度，反映提取水系與矢量化數(shù)字水系偏移程度，對(duì)河流總長(zhǎng)的相對(duì)誤差和河網(wǎng)“套合差”兩個(gè)指標(biāo)進(jìn)行計(jì)算分析。提取水系河流總長(zhǎng)52 061.6 km，數(shù)字水系河流總長(zhǎng)48 395.5 km，相對(duì)誤差為7.5%?；贏rcGIS圖層加載1∶1 000 000矢量化數(shù)字水系圖作為標(biāo)準(zhǔn)做檢驗(yàn)，對(duì)比數(shù)字化水系和提取水系之間的位置偏移，統(tǒng)計(jì)由此偏移產(chǎn)生的兩水系之間的細(xì)碎多邊形面積(見圖7)。統(tǒng)計(jì)結(jié)果顯示，整個(gè)瀾湄流域內(nèi)的細(xì)碎多邊形面積約為20 000 km2，河網(wǎng)套合差為2.5%，小于3%，吻合程度較好。由圖7a可見，上游地區(qū)海拔較高、地勢(shì)變化明顯且無(wú)較大湖泊，兩水系疊加形成的多邊形相對(duì)較少，多邊形面積約為2 100 km2，對(duì)應(yīng)的流域面積為166 000 km2，河網(wǎng)套河差為1.3%；圖7b顯示為流域中下游地區(qū)，該范圍內(nèi)地勢(shì)逐漸變緩，水系豐富，在柬埔寨境內(nèi)有較大湖泊洞里薩湖，兩水系之間的多邊形較上游明顯增多且面積增大，計(jì)算得中下游的河網(wǎng)套河差為2.9%，上游和中下游的河網(wǎng)套河差結(jié)果均在有效范圍內(nèi)，但中下游地區(qū)河網(wǎng)偏移程度較上游大，提取水系精度略差于上游地區(qū)?？傮w而言，本文提取水系結(jié)果與實(shí)際水系較為吻合，能有效、客觀地反映瀾湄流域水系整體特征，可為后續(xù)研究提供數(shù)據(jù)支撐。
圖7 提取水系與數(shù)字化水系河網(wǎng)套河差
3.3 河網(wǎng)分級(jí)和子流域劃分
由90 m分辨率高程數(shù)據(jù)提取最優(yōu)集水面積閾值的瀾湄流域水系情況如圖6所示，河網(wǎng)分為1～6級(jí)，一級(jí)為河道最低級(jí)別的河源水系，共計(jì)1 342條，占總河道數(shù)的50%，一級(jí)河網(wǎng)發(fā)育系數(shù)約為2.9，河系不均勻系數(shù)約為1.3。
基于Arcgis自動(dòng)化提取流域出口，共劃分小流域3 259個(gè)，流域面積普遍在50～900 km2范圍。為滿足業(yè)務(wù)應(yīng)用需要，依據(jù)水系上下游匯流關(guān)系，建立河道水系與小流域以及各小流域間的拓?fù)潢P(guān)系；參考流域內(nèi)地形特征以及湖泊范圍，滿足水系連續(xù)性、分水線完整性、出水口準(zhǔn)確性等關(guān)鍵條件；對(duì)流域面積在5 km2以下的微小流域或中間匯流區(qū)，以匯流關(guān)系為主要依據(jù)合并到相鄰小流域內(nèi)[31-32]；另外，又以流域內(nèi)干流上的8個(gè)主要水文站點(diǎn)作為流域出口對(duì)其進(jìn)行流域劃分；進(jìn)而完成人工合并小流域的工作。將瀾湄流域的子流域劃分出3個(gè)等級(jí)，根據(jù)子流域面積由小到大依次定義為一級(jí)、二級(jí)、三級(jí)流域。一級(jí)流域共計(jì)9個(gè)子流域，子流域面積最小的是昌都—舊州區(qū)間，約為1.6萬(wàn)km2，最大的是穆達(dá)漢—上丁區(qū)間，約為24萬(wàn)km2。二級(jí)流域共計(jì)61個(gè)子流域，子流域面積小于5 000 km2的占比25%，在5 000～10 000 km2的占23%，在10 000～20 000 km2的占33%，大于20 000 km2的占19%。三級(jí)流域共計(jì)328個(gè)子流域，子流域面積普遍控制在500～5 000 km2，占比達(dá)到90%。其中，小于500 km2的為不適合與相鄰合并的中間匯流區(qū)，另存在一個(gè)大于5 000 km2的子流域，受洞里薩湖水域限制，包含湖域區(qū)域，子流域面積約6 700 km2。
3.4 不同子流域尺度流域面雨量預(yù)報(bào)分析
選取2021年7月24日～25日的一次降水過程預(yù)報(bào)為例，即為24日8時(shí)起報(bào)的24 h面雨量，分析不劃分子流域的瀾湄流域以及劃分出一級(jí)流域、二級(jí)流域的面雨量結(jié)果對(duì)此次過程在空間分布上的響應(yīng)。未劃分子流域的全流域經(jīng)算術(shù)平均計(jì)算面雨量為25 mm，全流域面積大僅以一個(gè)數(shù)值表達(dá)面雨量，在降水分布均勻且降雨量較小時(shí)可起到一定參考作用，但在降雨空間分布不均時(shí)無(wú)法體現(xiàn)降雨中心。圖8a中一級(jí)流域面雨量計(jì)算結(jié)果顯示，昌都—舊州和穆達(dá)漢—上丁面雨量量級(jí)為小雨，舊州—允景洪量級(jí)為中雨，允景洪—清盛、清盛—瑯勃拉邦和瑯勃拉邦—萬(wàn)象段量級(jí)為大雨，舊州—允景洪和萬(wàn)象—穆達(dá)漢量級(jí)達(dá)到暴雨。圖8b顯示，二級(jí)流域進(jìn)一步縮小子流域面積，面雨量結(jié)果與智能網(wǎng)格預(yù)報(bào)降雨結(jié)果在暴雨中心的空間體現(xiàn)上更為一致，降雨較大范圍集中于中下游地區(qū)，泰國(guó)中東部和老撾中南部交界的位置，有6個(gè)子流域面雨量超過60 mm，量級(jí)為大暴雨，最大在泰國(guó)境內(nèi)達(dá)到143 mm，與一級(jí)流域萬(wàn)象—穆達(dá)漢段暴雨級(jí)面雨量相呼應(yīng)。
圖8 面雨量空間分布
4 結(jié)論與討論
本文分析了瀾湄流域集水面積閾值與水系特征變化的數(shù)量關(guān)系，建立了集水面積閾值與河網(wǎng)密度之間的冪函數(shù)，通過曲線割線斜率法客服主觀性，確定了90 m分辨率高程下的最佳集水面積閾值為129.1km2。并以1∶1 000 000的矢量化數(shù)字水系圖做驗(yàn)證依據(jù)，總河長(zhǎng)相對(duì)誤差為7.5%，河網(wǎng)套河差為2.5%，均說明提取水系與實(shí)際情況吻合程度較高，提取水系有效，進(jìn)一步劃分了瀾湄流域1～6級(jí)的水系和3個(gè)等級(jí)的子流域?；谥悄芫W(wǎng)格預(yù)報(bào)降水，以一次降水過程展現(xiàn)面雨量在不同等級(jí)子流域上的空間分布，子流域等級(jí)越大、劃分越細(xì)致對(duì)暴雨中心的體現(xiàn)越明顯。本文的研究成果，可在水文氣象業(yè)務(wù)上為瀾湄流域的面雨量監(jiān)測(cè)、預(yù)報(bào)和水文模型的構(gòu)建提供參考，為數(shù)字孿生流域建設(shè)與流域防洪減災(zāi)奠定基礎(chǔ)。
雖然本文對(duì)DEM進(jìn)行了修正工作，但是在處理較大湖泊、入?？谌侵薜忍厥獾匦紊吓c真實(shí)情況方面仍存在差異；研究區(qū)域面積大，在高精度高程數(shù)據(jù)情況下運(yùn)行速度慢，選取的集水閾值數(shù)據(jù)量有限。下一步工作可以考慮按照地形、水系類型等因素先將整個(gè)流域劃分為幾個(gè)不同的小區(qū)域，以求在保證高程精度情況下提高運(yùn)算效率，進(jìn)而綜合模擬流域水系。
猜你喜歡
水科學(xué)進(jìn)展(2023年4期)2023-10-07 11:23:44
農(nóng)業(yè)科學(xué)研究(2022年2期)2022-08-01 05:25:30
中國(guó)造紙(2022年2期)2022-04-04 16:40:13
河北果樹(2020年1期)2020-02-09 12:31:40
華東師范大學(xué)學(xué)報(bào)(自然科學(xué)版)(2019年3期)2019-06-24 05:29:09
成都信息工程大學(xué)學(xué)報(bào)(2019年1期)2019-05-20 09:14:26
中國(guó)港灣建設(shè)(2017年11期)2017-12-19 12:27:05
水利科技與經(jīng)濟(jì)(2017年5期)2017-04-22 02:39:28
西藏科技(2016年5期)2016-09-26 12:16:40
水利科技與經(jīng)濟(jì)(2016年9期)2016-04-22 01:07:22
