999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

干旱荒漠區人工綠洲土壤鹽堿化風險綜合評估與演變分析

2022-01-21 00:45:14劉子金徐存東朱興林周冬蒙田俊姣谷豐佑李智睿趙志宏
中國環境科學 2022年1期
關鍵詞:研究

劉子金,徐存東*,朱興林,周冬蒙,田俊姣,谷豐佑,黃 嵩,李智睿,趙志宏,王 鑫

干旱荒漠區人工綠洲土壤鹽堿化風險綜合評估與演變分析

劉子金1,2,3,徐存東1,2,3*,朱興林4,周冬蒙4,田俊姣1,2,谷豐佑1,2,黃 嵩1,2,李智睿1,2,趙志宏1,2,王 鑫1,2

(1.華北水利水電大學水利學院,河南 鄭州 450046;2.浙江省農村水利水電資源配置與調控關鍵技術重點實驗室,浙江 杭州 310018;3.河南省水工結構安全工程技術研究中心,河南 鄭州 450046;4.中國科學院西北生態環境資源研究院,甘肅 蘭州 730000)

為明晰干旱荒漠區人工綠洲水鹽時空分異特征與鹽堿化風險空間格局演變過程,以甘肅省景電灌區為研究區,以2002,2010及2018年為研究代表年,基于多級模糊理論從地質氣候驅動、水土環境驅動和自然-人類驅動3個驅動過程構建土壤鹽堿化風險評估體系,集成云發生器原理、黃金分割率法、組合賦權法以及排隊理論構建土壤鹽堿化空間風險綜合評價模型.將長序列監測數據、多時相空間數據以及經濟社會數據依據各驅動要素權重以ArcGIS10.2為技術平臺進行多源融合,對各空間風險狀態進行了空間可視化表達與流向追蹤.結果表明:1)研究區2002, 2010及2018年土地鹽堿化整體風險分別為“臨界狀態”、“臨界狀態—輕度風險”、“輕度風險”;2)地下水埋深及地下水礦化度是驅動土壤鹽堿化進程加劇的主導因素;3)灌區鹽堿化空間風險演變模式劇烈程度排序為持續變化型>前期變化型>后期變化型>持續穩定型>反復變化型.研究區鹽堿化整體風險狀態呈現出由“無風險—臨界狀態”以及“臨界狀態—輕度風險”過渡的趨勢;4)2002~2018年間研究區鹽堿化空間風險整體呈升級態勢,表征出明顯的地區差異性且總體呈現由西北向東南以弧射狀增高的空間格局.

干旱荒漠區;人工綠洲;鹽堿化風險;空間格局;演變分析

我國西北干旱荒漠區土地資源豐富、光熱條件充足,但水資源空間配置不均[1].通過提水灌溉的方式,極大緩解了區域內水土資源不協調的問題.但干旱荒漠區生態環境本地極端脆弱,對于外調水資源驅動背景下的水土資源開發響應十分劇烈[2].灌溉農業對于土地資源的掠奪性開發致使區域內水鹽平衡演變模式及原生驅動過程受到了潛在破壞,進一步推進了區域內土壤鹽堿化及次生鹽堿化的演變進程,使得鹽堿化空間風險逐漸升級[3-4].因此揭示干旱荒漠區水鹽驅動要素時空分異特征及鹽堿化空間風險演變模式,對于豐富和推進鹽堿化監測、防治的針對性與精準性研究具有重要的科學意義.

干旱荒漠區人工綠洲的土壤水鹽空間分異存在較為復雜的生態水文過程與水力聯系[5].準確揭示其驅動過程的動力學機制,探明其時空演變規律與參與介質、驅動過程之間的耦合關系是水土環境領域當前研究的熱點與難點[6].近年來對土地鹽堿化驅動機理、介質聯系、空間風險以及演變過程的認識已經取得一些進展.如針對灌溉背景下土地鹽堿化的演變趨勢、發生概率以及監測分區進行的分析與評價[7-9];運用SaltMod模型、多元地質統計原理及模糊隸屬函數就土地鹽堿化的演變態勢預測、空間變異過程及總鹽度風險評價進行分析[10-12];基于灰色系統理論分別就河西地區綠洲、蘇北海涂典型圍墾區以及于田綠洲的土壤鹽堿化風險進行的定量評價[13-16].但仍存在以下不足之處,其一就研究對象及研究目標而言,目前的研究主要針對于田間尺度下灌溉條件、方式、土質結構對于典型單元土壤鹽堿化的驅動機理及驅動過程研究.對于區域尺度背景下土壤鹽堿化整體空間風險狀態的評估研究相對較少;其二就研究方法而言,仍存在相對單一或偏主觀等問題,在一定程度上影響了鹽堿化風險狀態評價結果的客觀性和準確性;其三,對于鹽堿化不同風險狀態的空間變化以及風險地類特征演變流向考慮的仍不夠充分.云理論在刻畫模糊評價系統的隨機性及不確定性問題上更具優勢,使得評價結果更貼近實際情況[17-20].

鑒于此,以甘肅省景電灌區為例,基于多級模糊理論、云發生器原理、黃金分割率法、組合賦權法以及排隊理論構建了一套土壤鹽堿化空間風險評價的多級模糊評價模型,同時引入土地鹽堿化空間風險轉移矩陣與鹽堿化風險地類特征變化模式圖譜,通過多時相數據融合的方式,將2002~2018年間研究區土壤鹽堿化空間風險狀態進行了定量評價與可視化表達,對于不同風險地類特征變化模式與演變流向進行了靶向追蹤分析,揭示了灌區土壤鹽堿化風險在多時空尺度上的發展態勢與演變規律.以期為干旱荒漠區人工綠洲水土資源可持續發展提供理論依據與科學指導.

1 研究區及數據源

1.1 研究區概況

如圖1所示,景電灌區地處甘肅省中部(37°26′~ 38°41′N,103°20′~104°04′E),地理位置處于甘、寧、蒙3省的交界地帶,是為緩解我國西北干旱荒漠區土地沙化及生態環境惡化的生態建設工程.灌區總控制面積約1496km2,灌區內土地類型以耕地、鹽堿地、草地、旱地、沙地及戈壁為主.氣候條件為典型的溫帶大陸性氣候.灌區分為一期與二期,其中一期灌區于1974年基本建成,二期灌區于1994年基本建成,該地區干旱少雨,多年平均降雨量為185.6mm.多年平均蒸發量為2433.8mm.土壤類型以荒漠灰鈣土為主,植被結構單一,主要表現為旱生草本混合群落與超旱生小灌木特征.年日照時數長達2714h,無霜期約190d[5,21].

圖1 研究區地理位置

1.2 多級模糊評價指標體系建立

引入多級模糊理論實現對評價目標多層次驅動過程的客觀揭示[22].結合已有針對水鹽分異驅動機理、過程以及土地鹽堿化時空演化格局的相關研究成果[2,8,16],參照研究區實際情況與歷史資料[5,20].將區域尺度土地鹽堿化空間風險狀態用“高度風險1”、“中度風險2”、“輕度風險3”、“臨界狀態4”、“無風險5”來描述.參與驅動的過程可以概述為地質氣候驅動、水土環境驅動和自然-人類驅動.將驅動因子列定為:地形、氣候、土壤、水資源、自然干擾及人類干擾8個因子.依據驅動狀態可進一步將驅動因子進行分解并表征為:海拔1、坡度2、年蒸發量1、年降水量2、土壤含鹽量1、土壤電導率2、地下水埋深1、地下水礦化度2、地表灌水量3、植被覆蓋度1、地表溫度2、人口密度1以及土地利用類型2.多級模糊評價指標體系見圖2.

圖2 多級模糊評價指標體系

1.3 數據來源

結合灌區建設運行過程中本底數據以及為滿足研究條件的充分性與代表性,選定研究代表年份分別為:2002年一期灌區續建配套全面完成,總灌溉面積為3.85萬hm2,年提水量為3.22億m3;2010年二期灌區續建配套建設2a,總灌溉面積為4.92萬hm2,年提水量為3.86億m3;2018年代表灌區運行現況,總灌溉面積6.05萬hm2,年提水量為4.60億m3.研究數據主要包含長序列監測數據、無人機航拍數據、遙感反演數據以及經濟社會數據.各數據獲取方式如下:(1)長序列監測對象包括年蒸發量、年降水量、土壤含鹽量、地下水埋深、地下水礦化度以及地表灌水量等鹽堿化驅動狀態因子,其主要由布設的空間監測設備以及現場樣本提取分析進行數據收集;(2)無人機航拍掃描數據來源于2017年4月25~26日在灌區采用Trimble UX5型固定機翼無人機(飛行高度500m,控制精度為15cm,采用“蛇型移動式”航線)搭載的SONY ILCE-5100相機(2400萬像素)配備影像傳感器對研究區東部的封閉型水文地質單元進行連續航拍所得;(3)遙感反演數據及氣象數據包括:由地理空間數據云(http://www. gscloud.cn)獲取的Landsat8及Landsat5系列衛星搭載數據;由Weather Underground(https://www. wunderground.com)所提供的歷史氣象數據;由通量監測設備所獲取的水熱通量數據;(4)經濟社會數據主要來源于景電灌區內各縣、鎮、村的統計年鑒、經濟社會歷史調查報告以及由中國科學院資源環境科學與數據中心(https://www.resdc.cn)提供的相關資料.各項評價指標獲取來源見表1,各指標監測點位空間分布見圖3.

表1 評價因子數據來源

圖3 長序列監測點位空間分布

2 研究方法

2.1 數據處理

為實現數據的空間面狀化,將源數據進行以下預處理:(1)長序列監測數據:通過Origin9.1正態QQ圖分析對數據先進行樣本源檢驗,結果表明,上述指標數據分布滿足正態性分布.基于ArcGIS- Spatial Analys-插值對年蒸發量1、年降水量2、土壤含鹽量1、土壤電導率2、地下水埋深1、地下水礦化度2、地表灌水量3進行空間插值.引入交叉驗證法[23]及誤差矩陣[24]進行插值精度驗證,最終優選插值方法及相應精度見表2.研究區年降水量及年蒸發量各研究節點無顯著變異性,故本文只列出2018年的相關結果.

表2 各評價指標因子所選插值方法及精度

注:MAE、MRE、RMSE分別表示平均絕對誤差、平均相對誤差以及均方根誤差.

遙感反演及航拍數據:對Landsat光學數據分別進行輻射定標、大氣校正[25]等預處理,基于ArcGIS10.2提取了研究區數字高程信息,地表溫度、植被覆蓋度分別采用單窗算法[26]、像元二分模型[27]進行解譯處理.土地利用在解譯過程中,為提高解譯精度,將航拍影像數據基于UASMaster進行影像集成,并利用POS數據完成數據定向及點云提取,最后通過配準、建模獲取整體航拍掃描數據.通過目視解譯及現場調查,確定感興趣區提取遙感影像解譯標志.并以此作為土地利用解譯的前置條件.通過反復對比建立訓練樣本,最終基于ENVI5.3采用人機交互型解譯方式進行解譯.經濟社會數據:將統計年鑒、經濟社會歷史調查報告以及中國科學院資源環境科學與數據中心提供的柵格數據進行比對驗證,確保滿足精度要求,對研究區不同時間節點的人口密度及分布進行空間表征.

2.2 云模型

2.2.1 云模型概念及云發生器原理 云模型能夠較好刻畫模糊系統隨機性與波動性,可以通過不確定性語言將研究系統的隨機性與模糊性進行有機結合,能夠更為準確揭示土地鹽堿化受多因素耦合影響發生演化這一復雜過程,彌補了以往研究的不足.云模型基于(期望)、(熵)、(超熵)3個云數字特征將多變量耦合關聯的模糊性進行表征,其轉換關聯系統CG包含正向與逆向兩個過程,可實現定量數值與定性概念的相互轉換.3個云數字特征分別代表云滴重心、云滴離散度以及云滴厚度[20].其概念特征及云發生器原理見圖4.

圖4 正態云模型數字特征

2.2.2 評價標準云 水鹽時空分異過程所表征出的具體狀態具有較強的模糊性,黃金分割率云生成法可以在一定程度上消除指標之間的模糊性,具有較好的理論價值與應用價值[28].以土地鹽堿化風險評語集={1,2,3,4,5}={高度風險,中度風險,輕度風險,臨界狀態,無風險}為基礎,引入黃金分割率,定義相鄰等級之間min=0.618max、min= 0.618max,設定“輕度風險”狀態云模型參數為3=0.5,3=0.005,由此確立各評價等級的數字分布特征分別為1(0,0.103,0.1031)、2(0.309,0.064, 0.0081)、3(0.50,0.031,0.005)、4(0.691,0.064, 0.0081)、5(1,0.103,0.0131),評語集系統通過MATLAB仿真得圖5.

圖5 土地鹽堿化風險狀態評估標準云

2.2.3 組合賦權法 為充分滿足指標賦權的客觀性與精度要求,避免與實際狀態不符以及個人經驗性對賦權結果的精度擾動.本文引入基于多目標加權的組合賦權法進行指標權重確定,該方法能夠避免傳統賦權方法的主、客觀性過強以及偏好性,能夠更加精準的評判各影響要素在復合系統中實際重要性程度,具有特有優勢[29].計算過程如式1所示.

式中: wi為第i個指標的權重; n為指標個數; i為指標重要性排序等級.

2.2.4 隸屬度云模型與綜合評價云模型 基于正、逆向云發生器原理,界定評價閾值為[0~1],聘請專家學者、灌區綜合管理人員以及灌區內運行工作人員依據評語層對狀態層評價指標打分,對評判數據基底進行有效性篩選、剔除及調整,確保后期計算所得的云數字特征滿足超熵<1/3及3原則[30],計算流程如圖6所示.

基于各指標權重及與隸屬度云模型依據式2計算綜合評價云.

式中:為指標個數;ω為各指標權重.

2.3 土地鹽堿化空間風險轉移矩陣

為揭示區域尺度土地鹽堿化風險在不同時期的空間分異過程與轉換關系,基于土地利用轉移矩陣[31],構建土地鹽堿化風險空間轉移矩陣對研究區水鹽分異態勢在數量上的分布特征以及不同風險等級地類流向進行靶向監督,其表達過程為:

式中:為研究區總面積,為土地鹽堿化風險等級類別,為研究期初土地鹽堿化風險等級地類特征,為研究期末土地鹽堿化風險等級地類特征.

2.4 土地鹽堿化風險地類特征變化模式圖譜

土地鹽堿化風險變化模式圖譜主要用于表征各研究節點對應的不同土地鹽堿化風險地域特征的空間-屬性-過程一體化數據,是探究時空動態變化下水鹽空間分異演變模式以及演變格局的基本信息單元.參考文獻[32-33],結合景電灌區水土資源的實際狀態,分別用1、2、3、4、5對高度風險區、中度風險區、輕度風險區、臨界狀態區以及無風險區進行編碼.將研究區土地鹽堿化風險演化模式定義為5種圖譜類型,分別為:1)持續變化型(1-2-3型),即研究區各階段鹽堿化風險地類特征均持續發生演化;2)反復變化型(1-2-1型),即不同階段鹽堿化風險地類特征在前期與后期演化模式相反,屬于正向與逆向交替演化型;3)前期變化型(1-2-2型),鹽堿化風險地類特征僅在前期演化,后期保持靜止穩定態;4)后期變化型(1-1-2型),鹽堿化風險地類特征僅在后期演化,前期保持靜止穩定態;5)持續穩定型(1-1-1型),各研究時間節點對應的鹽堿化風險地類特征均未發生演化.在圖譜特征定義基礎上,對已編碼的圖斑單元進行融合,將融合后的鹽堿化風險地類變化模式與圖譜模式比對分析從而得出區域尺度不同土地鹽堿化風險狀態的地類變化模式圖譜.計算公式為:

式中:分別表示研究期初、期中與期末土地鹽堿化風險空間地類圖譜柵格屬性;表示研究期內不同風險地類變化圖譜柵格屬性.

3 結果與分析

3.1 權重的求解

咨詢本領域專家3名,依據指標重要度確定其排隊等級.依據式1分別計算“驅動狀態層”、“驅動因子層”以及“驅動過程層”的權重值.結果見表3~5.

表3 “驅動狀態層”分層指標權重

表4 “驅動因子層”指標權重

表5 “驅動過程層”指標權重

3.2 隸屬度云模型計算與土地鹽堿化風險狀態綜合評估

聘請本領域專家學者3人,研究區管理人員3人,日常運行及維護工作人員4人,界定評語區間為[0,1],就驅動狀態要素對研究區各時間節點對應的土壤鹽堿化風險進行客觀評定.將所得評定分值基于云發生器原理進行數字表征,依據隸屬度云模型計算流程進行數據優選與正逆向反饋調整,確保結果客觀性及有效性,并依次計算得出各研究年“狀態層”“因子層”與“過程層”的云模型特征參數,最終依據式2計算得出綜合評語云數字特征,各研究節點過程層與綜合評語云數字特征見表6.2002、2010及2018年綜合評價云數字特征分別為1(0.7067, 0.0392,0.0107)、2(0.6162,0.0406,0.0111)、3(0.4768, 0.0491,0.0128).各時期云數字特征中,熵()與超熵()分別僅為0.0392, 0.0406, 0.0491, 0.0107, 0.0111, 0.0128.小于1/3,表明云滴離散度低,評價結果霧化度較低,更偏向實際狀態,可信性高.

表6 過程層與綜合評價云數字特征參數

基于MatLab將2002、2010及2018年土地鹽堿化風險云數字特征進行可視化表達,結果見圖7.由圖7可知,研究區2002, 2010及2018年土地鹽堿化整體風險分別為“4-臨界狀態”、“4-臨界狀態~3-輕度風險”、“3-輕度風險”.各時期鹽堿化風險整體表征出加劇態勢,表明灌區內土地鹽堿化風險等級不斷升級.究其原因,灌區內土地鹽堿化風險持續升級主要受控于地表水資源調控背景下大水漫灌式等傳統型、粗放式田間灌溉模式以及干旱荒漠型氣候條件的影響,水-熱-鹽輸移過程被不斷驅動,灌區內土地資源平衡演變模式受到了劇烈擾動,導致地表鹽分呈積聚型響應態勢.

圖7 隸屬度云模型計算過程

3.3 研究區土地鹽堿化空間風險分析

3.3.1 2002~2018年土地鹽堿化風險空間變化分析 參照上文中指標權重確定方式,通過專家咨詢,確定各指標因子的總體排隊等級,從驅動狀態層整體層面確定各指標因子權重,各驅動狀態指標因子權重見表7.

表7 “驅動狀態層”指標權重

由表7可知,研究區區域尺度水鹽分異進程驅動影響因素排序為地下水埋深>地下水礦化度>土壤含鹽量>地表灌水量>土壤電導率>地表溫度>年蒸發量>年降水量>植被覆蓋度>土地利用類型>人口密度>海拔>坡度.可見地下水埋深、地下水礦化度、土壤含鹽量以及地表灌水量是驅動灌區土壤水鹽運移的控導性要素.主要由于地下水作為干旱區人工綠洲水鹽運移的主要載體,也是地表鹽分積聚的主要來源;土壤含鹽量作為灌區土壤鹽堿化的終點因素,表征土壤鹽堿化進程中的屬性及具體狀態;地表灌水量主要通過在微觀尺度對土壤、地下水中的鹽分及礦化物質進行輸移,以此實現對土壤鹽分的再分布作用.除上述主要控導因素外,水文氣候要素、植被覆蓋、土地利用、人口密度、地形要素等主要通過潛在影響水鹽輸移過程從而對水鹽分異態勢起到微觀驅動作用,這種作用的表征較為緩慢且呈隱性,故對土地鹽堿化演化過程影響相對較弱.

根據研究區土壤鹽堿化風險因子插值、反演及解譯結果,將各鹽堿化驅動因子進行時空疊加,考慮到各驅動要素的量綱不一致,故在進行多因素耦合疊加分析時,先基于ArcGIS10.2對各驅動要素的空間數據圖件進行重分類處理,重分類過程基于“正向累加原則”,即各指標最終累加數值結果與土壤鹽堿化風險呈正相關性.結合所構建的多級模糊評價系統評語集和層,將分類標準劃定為5類,即“高度風險1”、“中度風險2”、“輕度風險3”、“臨界狀態4”及“無風險5”.結合各驅動要素權重,通過柵格疊合計算,分別得出2002、2010及2018年的土壤鹽堿化空間風險格局見圖8.基于ArcGIS10.2對灌區3個時期的土地鹽堿化空間風險面積進行統計.

圖8 研究區土壤鹽堿化風險空間分布示意

由圖8及表8可知,2002~2018年灌區土地鹽堿化空間風險總體表征出無風險區、臨界狀態區面積減少,輕度風險、中度風險以及高度風險區域面積增加的演變態勢.究其原因,結合研究區地表灌水量及土地利用方式的空間演變過程可知,隨著灌水量的不斷增高,灌區內土地利用方式發生了很大變遷,耕地面積呈持續增長狀態,沙地、未耕地以及草地面積均明顯下降,表明隨著時間推進,研究區水土資源開發度不斷升高,這也進一步加劇了灌區內人類活動對區域內水土資源的掠奪性開發過程,催動了土地鹽漬化的地表生態演變過程.從土壤鹽堿化空間風險面積轉移過程來看,2002、2010和2018年灌區內土壤鹽堿化無風險區域面積占比分別為37.04%、21.81%和19.39%,臨界狀態由2002年的41.44%下降至31.18%,面積減少了16057.2hm2,無風險及臨界狀態區面積呈縮減態勢,主要歸因于隨著灌區內水利設施以及生產方式的不斷轉變,增大了土地資源的利用度,使得灌區運行過程中農業生產活動對原生水土資源的破壞度不斷提升.輕度鹽堿化風險空間占比從2002年的9.62%提升至2018年的21.64%,中度鹽堿化風險空間占比從2002年的10.14%提升至2018年的16.14%,重度鹽堿化風險空間占比從2002年的1.76%提升至2018年的11.65%,灌區內鹽堿化整體呈現出持續演化態勢.結合研究區土壤鹽堿化風險空間分布可知,灌區內鹽堿化風險呈現出明顯的“地區差異性”,并呈現出以弧射狀由西北向東南增加的態勢,低風險區則主要集中分布在海子灘鎮、大靖鎮、裴家營鎮、紅水鎮等區域.中高風險區主要集中于研究區東部的蘆陽鎮、草窩灘鎮以及五佛鄉等區域.結合研究區的水文地質條件可知,低風險區所對應水文地質結構為開敞型水文地質單元,其水文地質單元由灌溉入滲帶(入滲主導、鹽隨水移)、溶質活躍遷移帶(水量多變,運移遲滯)以及水鹽耗散帶(鹽隨水移、深層耗散)3部分組成,地下水鹽運移的流向與地勢變化相吻合,這為高礦化度地下水的排泄提供了良好的地質基礎,水鹽易耗散、難積聚.而高風險區對應的封閉型水文地質單元,該地質結構單元從外圍向盆地中心逐漸分別形成了潛水交替相對流暢的灌溉入滲帶、溶質活躍遷移帶,并逐漸轉化為水鹽交替遲緩的匯水聚鹽帶(蒸發主導,水散鹽聚),受該地質層結構所影響,高風險區地下水儲水層上通下阻,加上多年開放式漫灌等不合理灌溉方式,地下水水文及礦化度被不斷抬升,隨著水—熱—鹽驅動轉換過程的不斷進行,導致該地質區土壤鹽堿化風險不斷升級.在這種水鹽驅動模式下,灌區水鹽分異態勢呈激變式,最終導致研究區土壤鹽漬化及次生鹽漬化空間風險不斷升高.結合上述分析內容,研究區在后續應進一步對灌溉背景下水資源對土壤鹽堿化的驅動過程進行調控,可通過大力推廣滴灌、噴灌等節水灌溉方式,一方面可以有效減少灌溉用水量,節約水資源;另一方面灌溉水量的減少對于地下水埋深等水鹽驅動要素也可以進行反向調控.此外,在高風險區可通過開挖排堿溝等方式,加快土壤中集聚鹽分的耗散過程.

3.3.2 2002~2018 年土地鹽堿化風險變化圖譜分析 為研究景電灌區土地鹽堿化風險在2002~ 2010~2018年的流向,依據式4,基于ArcGIS10.2對2002~2018年的土地鹽堿化風險變化圖譜進行計算分析,結果見表9所示.

由表9可知,景電灌區在2002~2018年間土地鹽堿化空間風險演變類型中,持續變化型所占比重最大,像元數為872978個,演變面積為73924.64hm2;其次為前期變化型,像元數為422096個,演變面積為35743.57hm2;后期變化型像元數為377623個,演變面積為31977.49hm2;持續穩定型像元數為110165個,演變面積為9328.95hm2;反復變化型像元數為65284個,演變面積為5528.35hm2;結合各風險變化模式及其演變面積可知,灌區鹽堿化空間風險演變模式劇烈程度排序為持續變化型>前期變化型>后期變化型>持續穩定型>反復變化型.

表8 研究區土地鹽堿化風險空間變化

在持續變化型土地鹽堿化演變模式中,無風險—臨界狀態—輕度風險的面積占比最高,為43937.59hm2,59.43%.這一驅動模式表征研究區土壤鹽堿化的主要演變態勢,即表明研究區內鹽漬化風險正處于一個不斷加劇的演變狀態,反映人類活動與土地資源空間適應性不斷降低,同時也對反復變化型模式占比最小進行了反向驗證;前期變化型與后期變化型演變模式對應的最大演變像元分別是臨界狀態—輕度風險—輕度風險,無風險—無風險—臨界狀態,兩者變化面積分別為10128.19, 9039.73hm2,占各對應變化模式比重分別為28.34%與28.27%.結合3個研究時期所對應的鹽堿化風險空間分布狀態及各水鹽分異驅動因子的空間分布狀態,這兩種變化模式主要集中分布在研究區的中部及西北部區域,該區域受地質單元及人類活動影響,水鹽驅動過程相對滯緩,此外,受地下水礦化度、地下水埋深的動態擾動較東部地區較弱,故該地區鹽堿化空間風險對水鹽驅動要素的耦合響應程度相對較弱,致使這兩種模式對應區域鹽堿化進程相對緩慢,并在局部地區表征出驅動停滯的狀態.

表9 研究區2002~2018年土地鹽堿化空間風險演變圖譜分析

表10 灌區2002~2010年土地鹽堿化空間風險轉移面積矩陣(hm2)

3.3.3 2002~2018年土地鹽堿化空間風險演變流向分析 由上述研究內容可知,景電灌區在2002~2018年間土地鹽堿化不同空間風險面積發生了較大演變,基于ArcGIS10.2的空間分析工具對研究區內各類鹽堿化風險之間的轉移過程用鹽堿化風險轉移矩陣的方式進行定量計算.鹽堿化空間風險轉移矩陣中,行代表灌區期初的鹽堿化風險狀態,列表示期末的鹽堿化風險狀態,灌區內2002~2010年及2010~2018年的鹽堿化空間風險轉移面積矩陣如表10,表11所示.

表11 灌區2010~2018年土地鹽堿化空間風險轉移面積矩陣(hm2)

由表10可知,在2002~2010年景電灌區土壤鹽堿化空間風險轉移矩陣中,轉移面積最大的是無風險狀態向臨界狀態的轉移,演變轉移面積為25477.48hm2,占無風險轉出面積的43.9%,占臨界狀態轉入面積的37.5%.其次是臨界狀態向輕度風險的轉移,轉移面積為17289.2hm2,占臨界狀態轉出面積的26.7%,占輕度風險轉入面積的73.1%.究其原因,主要是由于二期灌區在2002~2010年間,隨著灌溉提水及配水設備的不斷完善,灌區內提水量及灌溉用水量不斷增加,這一過程使得灌區內水土資源演化態勢由原生態向次生態演繹,土地資源對水資源的敏感性不斷增強所致.同理由表11可知,在2010~2018年土壤鹽堿化空間風險轉移矩陣中,轉移面積最大的是臨界狀態向輕度風險的轉移,演變轉移面積為24695.8hm2,占臨界狀態轉出面積的36.4%,占輕度風險轉入面積的72.9%.其次是輕度風險向臨界狀態的逆向演變過程,面積為9395.59hm2,占輕度風險轉出面積的39.7%,占臨界狀態轉入面積的19.3%.究其原因,上述正向環境惡化驅動過程的演變以及局部地區出現逆向健康化恢復過程的原因主要有以下兩點:一是灌區內水鹽運移受“灌水—蒸發—耗散”過程以及地下水礦化度、地下水埋深等鹽堿化驅動要素的循環式驅動所影響,鹽堿化風險整體格局仍處于增長態勢.同時,隨著土地鹽堿化過程的不斷演進,灌區內環境狀態對灌溉背景下水土資源的適配性有改善的趨勢,但受正向驅動過程劇烈性所制約,這種逆向健康化演變過程演繹較為緩慢,且演繹速度遠小于正向惡化驅動速度.

4 討論

本文從“地質氣候”、“水土環境”以及“自然—人類”3個主要驅動過程介入,基于多級模糊理論、云發生器原理、黃金分割率法、組合賦權法以及排隊理論構建了一套適用于干旱荒漠區人工綠洲土壤鹽堿化空間風險評價的多級模糊評價模型,并以云數字特征的形式對研究區3個時期的土壤鹽堿化空間整體風險進行了定量評價.在此基礎上,基于長序列監測數據、空間數據以及經濟社會數據的空間信息提取處理與多源信息融合,對景電灌區2002, 2010年及2018年的土壤鹽堿化空間風險格局進行了可視化表達,同時建立了土地鹽堿化空間風險轉移矩陣與風險地類特征變化模式對2002~2018年間的土壤鹽堿化風險特征以及演變流向進行了分析.可為本研究區及相似研究區水土環境狀態評估與演變流向靶向追蹤提供一種新思路.

雖然在一定程度上克服了傳統土地鹽堿化風險評價的單一性、局限性、模糊性以及可視化效果不突出等弊端.但仍存在一定的不足之處,在構建多層模糊評價系統時,對水鹽運移的驅動力,即熱力條件驅動的系統性考慮不夠充分,僅引入了地表溫度以及蒸發量兩個主導要素,對于地表輻射等潛在熱力要素的微觀驅動作用考慮不足.其二在進行多源數據融合時,由于研究區面積相對較大,對于鹽堿化空間風險的尺度效應考慮不足,使得多源數據融合結果與實際土壤鹽堿化空間風險狀態存在一定差距.針對上述不足,在后續研究中,可進一步豐富與完善鹽堿化驅動要素體系,并將柵格單元做更進一步的劃分,以期提高鹽堿化空間風險的評價精度.

5 結論

5.1 研究區2002、2010及2018年綜合評價云數字特征分別為1(0.7067,0.0392,0.0107)、2(0.6162, 0.0406,0.0111)、3(0.4768,0.0491,0.0128).各時期土地鹽堿化整體風險分別為“4-臨界狀態”、“4-臨界狀態~3-輕度風險”、“3-輕度風險”.

5.2 各驅動要素對土壤鹽堿化演進的驅動影響排序為地下水埋深>地下水礦化度>土壤含鹽量>地表灌水量>土壤電導率>地表溫度>年蒸發量>年降水量>植被覆蓋度>土地利用類型>人口密度>海拔>坡度,表明地下水埋深及地下水礦化度是研究區內土壤鹽堿化進程加劇的主導因素.

5.3 2002~2018年間,研究區土地鹽堿化空間風險變化趨勢呈現無風險區域和臨界狀態區域急劇減少,輕度風險面積顯著增加,中度風險與高度風險面積增長速度相對緩慢等特點.灌區鹽堿化空間風險演變模式劇烈程度排序為持續變化型>前期變化型>后期變化型>持續穩定型>反復變化型.土地鹽堿化空間風險主體演變流向為“無風險—臨界狀態”以及“臨界狀態—輕度風險”.

5.4 研究區2002~2018年土壤鹽堿化空間風險格局整體呈增高趨勢,且呈現出明顯的空間差異性,整體空間風險表征出從西北部向東南部弧射狀升高的演化態勢.高風險區主要集中分布在研究區東南部的封閉型水文地質單元,而西北部的開敞型水文地質單元以及中部地區鹽堿化風險等級相對較低.

[1] 鄧銘江.中國西北“水三線”空間格局與水資源配置方略[J]. 地理學報, 2018,73(7):1189-1203.

Deng M J. "ThreeWater Lines" strategy: Its spatial patterns and effectson water resources allocation in northwest China [J]. Acta Geographica Sinica, 2018,73(7):1189-1203.

[2] 郭 勇,尹鑫衛,李 彥,等.農田-防護林-荒漠復合系統土壤水鹽運移規律及耦合模型建立[J]. 農業工程學報, 2019,35(17):87-101.

Guo Y, Yin X W, Li Y, et al. Soil water and salt dynamics and its coupling model at cropland-treebelt-desert compound system [J]. Transactions of the Chinese Society of Agricultural Engineering, 2019,35(17):87-101

[3] 楊勁松.中國鹽漬土研究的發展歷程與展望[J]. 土壤學報, 2008, 45(5):837-845.

Yang J S. Development and prospect of the research on salt-affected soils in China [J]. Acta Pedologica Sinica, 2008,45(5):837-845.

[4] 張 帥,董會忠,曾文霞.土地生態系統脆弱性時空演化特征及影響因素—以黃河三角洲高效生態經濟區為例 [J]. 中國環境科學, 2019,39(4):1696-1704.

Zhang S, Dong H Z, Zeng W X. The time-space evolution characteristics of the vulnerability of land ecosystems and influencing factors: A case study of theYellow River Delta Efficiency Eco- economic Zone [J]. China Environmental Science, 2019,39(4):1696- 1704.

[5] 徐存東.景電灌區水鹽運移對局域水土資源影響研究[D]. 蘭州:蘭州大學, 2010.

Xu C D. Research on the effect of local water and soil environment caused by water-salt transportation in Jing-Dian Irrigation district [D]. Lanzhou University, 2010.

[6] 王全九,鄧銘江,寧松瑞,等.農田水鹽調控現實與面臨問題[J]. 水科學進展, 2021,32(1):139-147.

Wang Q J, Deng M J, Ning S R, et al. Reality and problems of controlling soil water and salt in farmland [J]. Advances In Water Science, 2021,32(1):139-147.

[7] 胡宏昌,田富強,張 治,等.干旱區膜下滴灌農田土壤鹽分非生育期淋洗和多年動態[J]. 水利學報, 2015,46(9):1037-1046.

Hu H C, Tian F Q, Zhang Z, et al. Soil salt leaching in non-growth period and salinitydynamics under mulched drip irrigation in arid area [J]. Journal of Hydraulic Engineering, 2015,46(9):1037-1046.

[8] 毛海濤,黃慶豪,吳恒濱.干旱區農田不同類型土壤鹽堿化發生規律[J]. 農業工程學報, 2016,32(S1):112-117.

Mao H T, Huang Q H, Wu H B. Soil salinization features in arid areas farmland [J]. Transactions of the Chinese Society of Agricultural Engineering, 2016,32(S1):112-117.

[9] 楊厚翔,雷國平,徐 秋,等.基于危險度與風險格局的土地鹽堿化監測區優先級評價[J]. 農業工程學報, 2019,35(7):238-246.

Yang H X, Lei G P, Xu Q, et al. Priority evaluation of land salinization monitoring area based on danger degree and risk pattern [J]. Transactions of the Chinese Society of Agricultural Engineering, 2019, 35(7):238-246.

[10] Chang X, Gao Z, Wang S, et al. Modelling long-term soil salinity dynamics using SaltMod in Hetao Irrigation District, China [J]. Computers and Electronics in Agriculture, 2019,156:447-458.

[11] Castrignanò A, Buttafuoco G, Puddu R. Multi-scale assessment of the risk of soil salinization in an area of south-eastern Sardinia (Italy) [J]. Precision Agriculture, 2008,9(1/2):17-31.

[12] Buchanan S, Huang J, Triantafilis J. Salinity risk assessment using fuzzy multiple criteria evaluation [J]. Soil Use and Management, 2017,33(2):216-232.

[13] 李自珍,李維德,石洪華,等.生態風險灰色評價模型及其在綠洲鹽漬化農田生態系統中的應用 [J]. 中國沙漠, 2002,22(6):617-622.

Li Z Z, Li W D, Shi H H, et al. Gray model for ecological risk assessment and its application in salinization oasis agroecosystem [J]. Journal of Desert Research, 2002,22(6):617-622.

[14] 姚榮江,楊勁松,陳小兵,等.蘇北海涂典型圍墾區土壤鹽漬化風險評估研究[J]. 中國生態農業學報, 2010,18(5):1000-1006.

Yao R J, Yang J S, Chen X B, et al. Evaluating soil salinization risk in typical coastal reclaimed regions in North Jiangsu Province [J]. Chinese Journal of Eco-Agriculture, 2010,18(5):1000-1006.

[15] 李冬順,楊勁松,姚榮江.生態風險分析用于蘇北灘涂土壤鹽漬化風險評估研究 [J]. 土壤學報, 2010,47(5):857-864.

Li D S, Yang J S, Yao R J, et al. Application of ecological risk analysis to soil salinization risk assessment of coastal tidal flat in north Jiangsu province [J]. Acta Pedologica Sinica, 2010,47(5):857-864.

[16] 依力亞斯江?努爾麥麥提,師慶東,阿不都拉?阿不力孜,等.灰色評估模型定量評價于田綠洲土壤鹽漬化風險[J]. 農業工程學報, 2019,35(8):176-184.

Ilyas Nurmemet, Shi Q D, Abdulla Abliz, et al. Quantitative evaluation of soil salinization risk in Keriya Oasis based on grey evaluation model [J]. Transactions of the Chinese Society of Agricultural Engineering, 2019,35(8):176-184

[17] 劉華新,苑一鳴,周 沛,等.基于融合理論的風電機組狀態評價正態云模型[J]. 太陽能學報, 2018,39(10):2891-2900.

Liu H X, Yuan Y M, Zhou P, et al. Normal cloid model for condition evaluation of wind turbines based on fusion theory [J]. Acta Energiae Solaris Sinice, 2018,39(10):2891-2900.

[18] 冶運濤,梁犁麗,曹 引,等.基于可變集和云模型的河湖水系連通方案優選決策方法[J]. 農業機械學報, 2018,49(12):211-225,313.

Ye Y T, Liang L L, Cao Y, et al. Optimal decision method for interconnected river system network schemes based on variable sets and cloud model [J]. Transactions of the Chinese Society for Agricultural Machinery, 2018,49(12):211-225,313.

[19] Wang F, Zhong D H, Yan Y L, et al. Rockfill dam compaction quality evaluation based on cloud-fuzzy model [J]. Journal of Zhejiang University-SCIENCE A, 2018,19(4):289-303.

[20] 徐存東,程 慧,王 燕,等.灌區土壤鹽漬化程度云理論改進多級模糊評價模型[J]. 農業工程學報, 2017,33(24):88-95.

Xu C D, Cheng H, Wang Y, et al. Improved multi-level fuzzy evaluation model based on cloud theory for evaluation of soil salinization degree [J]. Transactions of the Chinese Society of Agricultural Engineering, 2017,33(24):88-95.

[21] 王榮榮.景電灌區土壤鹽漬化特征及水鹽運移規律研究[D]. 鄭州:華北水利水電大學, 2017.

Wang R R. Characteristics of soil salinization and law of water and salt transport in Jingdian Irrigation Area [D]. Zhengzhou: North China University of Water Resources and Electric Power, 2017.

[22] 汪培莊,李洪興.模糊系統理論與模糊計算機[M]. 北京:科學出版社, 1996.

Wang P Z, Li H X. Fuzzy system theory and fuzzy computer [M]. Beijing: Science Press, 1996.

[23] 李 佳,段 平,呂海洋,等.基于改進的逐點交叉驗證的RBF形態參數優化方法及其空間插值實驗[J]. 地理與地理信息科學, 2016, 32(3):39-42,48.

Li J, Duan P, Lv H Y, et al. YRBF shape parameter optimization approach basd on iloocve and its spatial interpolation experiments [J]. Geography and Geo-Information Science, 2016,32(3):39-42,48.

[24] 馬宏宏,余 濤,楊忠芳,等.典型區土壤重金屬空間插值方法與污染評價[J]. 環境科學, 2018,39(10):4684-4693.

Ma H H, Yu T, Yang Z F, et al. Spatial interpolation methods and pollution assessment of heavy metals of soil in typical areas [J]. Environmental Science, 2018,39(10):4684-4693.

[25] 彭繼達,張春桂.基于高分一號遙感影像的植被覆蓋遙感監測——以廈門市為例[J]. 國土資源遙感, 2019,31(4):137-142.

Peng J D, Zhang C G. Remote sensing monitoring of vegetation coverage by GF-1 Satellite: A case study of Xiamen City [J]. Remote Sensing for Land & Resource, 2019,31(4):137-142.

[26] 胡德勇,喬 琨,王興玲,等.單窗算法結合Landsat8熱紅外數據反演地表溫度[J]. 遙感學報, 2015,19(6):964-976.

Hu D Y, Qiao K, Wang X L, et al. Land surface temperature retrieval from Landsat 8thermal infrared data using mono-window algorithm [J]. National Remote Sensing Bulletin, 2015,19(6):964-976.

[27] 江 淼,張顯峰,孫 權,等.不同分辨率影像反演植被覆蓋度的參數確定與尺度效應分析[J]. 武漢大學學報(信息科學版), 2011,36(3): 311-315.

Jiang M, Zhang X F, et al. Vegetation coverage retrieval scale effect analysis using multi sensor data [J]. Geomatics and Information Science of Wuhan University, 2011,36(3):311-315.

[28] 傅鶴林,黃 震,黃宏偉,等.基于云理論的隧道結構健康診斷方法[J]. 工程科學學報,2017,39(5):794-801.

Fu H L, Huang Z, Huang H W, et al. Health diagnosis method of shield tunnel structure based on theory [J]. Chinese Journal of Engineering, 2017,39(5):794-801.

[29] 魏巍賢,馮 佳.多目標權系數的組合賦值方法研究[J]. 系統工程與電子技術, 1998,(2):14-16.

Wei W X, Feng J. Study on multiobjective weights Combination Assigning Method [J]. Systems Engineering and Electronics, 1998,(2): 14-16.

[30] Wang J Q, Peng L, Zhang H Y, et al. Method of multi-criteria group decision-making based on cloud aggregation operators with linguistic information [J]. Information ences, 2014,274:177-191.

[31] 張曉娟,周啟剛,王兆林,等.基于MCE-CA-Markov的三峽庫區土地利用演變模擬及預測[J]. 農業工程學報, 2017,33(19):268-277.

Zhang X J, Zhou Q G, Wang Z L, et al. Simulation and prediction of land use change in Three Gorges Reservoir Area based on MCE- CA-Markov [J]. Transactions of the Chinese Society of Agricultural Engineering, 2017,33(19):268-277.

[32] 龔文峰,袁 力,范文義.基于CA-Markov的哈爾濱市土地利用變化及預測[J]. 農業工程學報, 2012,28(14):216-222.

Gong W F, Yuan L, Fan W Y. Dynamic change and prediction of land use in Harbin city based on CA-Markov model [J]. Transactions of the Chinese Society of Agricultural Engineering, 2012,28(14):216-222.

[33] 王友生,余新曉,賀康寧,等.基于CA-Markov模型的藉河流域土地利用變化動態模擬[J]. 農業工程學報, 2011,27(12):330-336,442.

Wang Y S, Yu X X, He K N, et al. Dynamic simulation of land use change in Jihe watershed based on CA-Markov model [J]. Transactions of the CSAE, 2011,27(12):330-336,442.

Comprehensive assessment and evolution analysis of soil salinization in artificial oasis in arid desert area.

LIU Zi-jin1,2,3, XU Cun-dong1,2,3*, ZHU Xing-lin4, ZHOU Dong-meng4, TIAN Jun-jiao1,2, GU Feng-you1,2, HUANG Song1,2, LI Zhi-rui1,2, ZHAO Zhi-hong1,2, WANG Xin1,2

(1.School of Water Conservancy, North China University of Water Resources and Electric Power, Zhengzhou 450046, China;2.Key Laboratory for Technology in Rural Water Management of Zhejiang Province,Hangzhou 310018, China;3.Henan Provincial Hydraulic Structure Safety Engineering Research Center, Zhengzhou 450046, China;4.Northwest Institute of Eco-Environment and Resources, Chinese Academy of Sciences, Lanzhou 730000, China)., 2022,42(1):367~379

To clarify the spatiotemporal variation characteristics of water and salt and the spatial evolution of salinization in artificial oasis in arid desert area, Jingtaichuan electric power irrigation area in Gansu Province was taken as the research area, 2002, 2010 and 2018 were selected as representative years. The soil salinization risk assessment system was constructed based on multi-level fuzzy theory from three driving processes: geological-climatic driven, soil-water environmental driven and natural-human driven. The cloud generator principle, golden ratio method, combined assignment method and queuing theory were integrated to construct a comprehensive spatial risk evaluation model of soil salinization. The long-sequence monitoring data, multi-temporal spatial data, and economic and social data were fused based on the element weights using ArcGIS10.2. The spatial risk state of salinization in each period was visually expressed and flow traced. The results of the study showed that: 1) The overall risks of land salinization in 2002, 2010 and 2018 were "critical", "critical-mild risk" and " mild risk"; 2) The dominant factors in the soil salinization were the depth of groundwater and the mineralization of groundwater; 3) In the irrigation area, the intensity of spatial evolution pattern of salinization risk: continuous change type>early change type>late change type>continuous stable type>repeated change type. The overall salinization risk in the study area showed a transition trend from "risk-free to critical state" and "from critical state to mild risk"; 4) From 2002 to 2018, the salinization in the study area was aggravating with obvious regional differences and generally presented an arcing increase from northwest to southeast.

arid desert areas;artificial oasis;salinization risk;spatial pattern;evolution analysis

X171.4

A

1000-6923(2022)01-0367-13

劉子金(1996-),男,甘肅白銀人,華北水利水電大學碩士研究生,主要從事水利工程環境效應評估與灌區水鹽調控方面的研究.發表論文8篇.

2021-05-18

國家自然科學基金資助項目(51579102);河南省高校科技創新團隊支持計劃(19IRTSTHN030);中原科技創新領軍人才支持計劃(204200510048);浙江省重點研發計劃項目(2021C03019);浙江省聯合基金項目(LZJWD22E090001)

* 責任作者, 教授, xcundong@126.com

猜你喜歡
研究
FMS與YBT相關性的實證研究
2020年國內翻譯研究述評
遼代千人邑研究述論
視錯覺在平面設計中的應用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
關于遼朝“一國兩制”研究的回顧與思考
EMA伺服控制系統研究
基于聲、光、磁、觸摸多功能控制的研究
電子制作(2018年11期)2018-08-04 03:26:04
新版C-NCAP側面碰撞假人損傷研究
關于反傾銷會計研究的思考
焊接膜層脫落的攻關研究
電子制作(2017年23期)2017-02-02 07:17:19
主站蜘蛛池模板: 亚洲成人精品久久| 欧美一级高清视频在线播放| 国产精品三区四区| 中字无码av在线电影| 波多野结衣AV无码久久一区| 中文字幕亚洲专区第19页| 久久免费精品琪琪| 免费一级α片在线观看| 国产成人在线无码免费视频| 亚洲中文精品久久久久久不卡| 亚洲国产成人精品一二区| 国产视频a| 亚洲伊人天堂| 国产成人无码播放| 乱人伦视频中文字幕在线| AV熟女乱| 国产激爽大片高清在线观看| 国产在线日本| 国产精品一区在线麻豆| 国产在线精品99一区不卡| 日韩中文欧美| 国产不卡网| 亚洲无码日韩一区| 欧美精品综合视频一区二区| 亚洲国产一区在线观看| 日韩成人在线一区二区| 久久9966精品国产免费| 亚洲女人在线| 国产特级毛片aaaaaa| 精品99在线观看| 黄色免费在线网址| 免费无码AV片在线观看中文| 欧美a在线| 不卡色老大久久综合网| 无码视频国产精品一区二区| 亚洲天堂网在线播放| 欧美一区二区三区不卡免费| 国产精品流白浆在线观看| 国产麻豆精品久久一二三| 国产成人亚洲无吗淙合青草| 欧美一区二区福利视频| 国产精品久久久久鬼色| 奇米影视狠狠精品7777| 久久综合干| www.亚洲天堂| 久久久久88色偷偷| 国产精品va免费视频| 激情综合网激情综合| 久久动漫精品| 任我操在线视频| 免费99精品国产自在现线| 国产乱人伦精品一区二区| 成人午夜网址| 一边摸一边做爽的视频17国产| 国产精品自在自线免费观看| 女同国产精品一区二区| 国产粉嫩粉嫩的18在线播放91| 国产理论最新国产精品视频| 91精品视频在线播放| 国产成人av一区二区三区| 久久久久久久久亚洲精品| 国产精品手机在线观看你懂的| 免费观看成人久久网免费观看| 丝袜无码一区二区三区| 毛片网站观看| 激情在线网| 日本a∨在线观看| 免费看a毛片| 国产精品无码一区二区桃花视频| 2021亚洲精品不卡a| 日韩人妻无码制服丝袜视频| 美女一级免费毛片| 国产乱人免费视频| 99热这里只有免费国产精品| 亚洲色图狠狠干| 18禁色诱爆乳网站| 国产91小视频| 日韩欧美91| 久久精品最新免费国产成人| 国产91成人| 国产综合另类小说色区色噜噜| 制服丝袜一区二区三区在线|