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

基于可解釋性分析的黃河流域生態系統 NPP 時空異質性研究

2025-09-28 00:00:00劉子華王海軍曹永強王金珂
人民黃河 2025年6期

關鍵詞:生態系統;凈初級生產力;碳匯;可解釋性分析;黃河流域中圖分類號:K903;TV882.1 文獻標志碼:A doi:10.3969 / j.issn.1000-1379.2025.06.016引用格式:,等.基于可解釋性分析的黃河流域生態系統 NPP 時空異質性研究[J].人民黃河,2025,47(6):103-109.

Exploring the Spatial?Temporal Evolution Characteristics and Driving Factors of Carbon Sequestration in the Yellow River Basin Based on Explainability

LIU Zihua1, WANG Haijun2, CAO Yongqiang1, XU Fang3, WANG Fei4, WANG Jinke1(1.Institute for Ecological Civilization Development in the Beijing?Tianjin?Hebei Region, Tianjin Normal University,Tianjin 300387, China; 2.Monitoring Center for Soil and Water Conservation, Ministry of Water Resources, Beijing 100053, China;3.Tianjin Institute of Geological Survey, Tianjin 300191, China; 4.School of Economic Geography (Institute of EconomicGeography of Hunan Province), Hunan University of Finance and Economics, Changsha , China)

Abstract: In order to enrich the research on Net Primary Productivity (NPP) and provide insights for the ecological sustainability and carbon sink enhancement of the Yellow River Basin, this study investigated the temporal and spatial patterns of NPP from 2000 to 2020. The basin was divided into three ecological zones of the semi?arid region of the mid?temperate zone (Zone I), the semi?humid region of the warm temperate zone (ZoneⅡ) and the semi?arid region of the plateau zone (Zone ). Based on multi?source remote sensing data, the Carnegie?Ames?Stanford Approach (CASA) model was used to estimate annual and monthly NPP values for the entire basin and its three sub?regions. To further inter? pret the impact of environmental factors, the Shapley additive explanations (SHAP) method was introduced to analyze the contributions of tem? perature, precipitation, NDVI and solar radiation to NPP variation. The main findings are as follows: a) Spatially, NPP shows significant het? erogeneity, with higher values concentrated in Zone II and lower values widely distributed in ZonesⅠand Ⅲ. The mean NPP values rank as fol? lows: ZoneⅡ> Zone Ⅲ> Zone I. Areas along the Yellow River generally exhibit higher NPP than other areas at the same latitude. b) Temporal? ly, NPP in the basin shows an overall upward trend with fluctuations and exhibits strong intra?annual seasonality. During the study period, 47.17% of the area experiences an increase in NPP, while 3.33% shows a decrease. ZoneⅡexperiences the most significant increase, and the growth rate of NPP gradually slows over time. c) At the basin scale, NDVI has the greatest influence on NPP, and its interaction with solar ra? diation is particularly strong. However, the dominant driving factors are varied by sub?region: NDVI in Zone I, air temperature in Zone II and NDVI again in Zone Ⅲ. d) The strategies for enhancing NPP in each ecological area of the Yellow River Basin should be adapted to local condi? tions, have their own focuses, and fully attach importance to the interaction effects among various natural factors.

Key words: ecosystem; NPP; carbon sink; interpretability analysis; Yellow River Basin

0 引言

聯合國政府間氣候變化專門委員會(IPCC)于2014 年發布的第五次評估報告指出工業革命后化石燃料產生的溫室氣體使得過去 100a 全球平均氣溫上升了 0.74°C ,2023 年發布的第六次評估報告指出全球升溫幅度可達 1.5°C 。 對此,全球正在致力于避免溫室氣體排放造成的最壞后果。 我國作為最大的發展中國家,正處于快速城市化進程中,具有較大的能源消耗和碳排放需求,盡管如此,依然彰顯大國擔當、提出“雙碳”目標(力爭2030 年前二氧化碳排放達到峰值、努力爭取2060 年前實現碳中和),并采取政策引導、能源結構調整、產業優化、碳市場建設、技術創新、生態保護等一系列措施[1],積極應對氣候變化。

減少碳排放、增加碳匯是實現碳中和的兩條路徑[2]。 黨的二十大報告提出要提升生態系統碳匯能力。 2023 年4 月,自然資源部、國家發展改革委、財政部、國家林草局聯合印發了《生態系統碳匯能力鞏固提升實施方案》。 凈初級生產力(植被在單位時間、單位面積上通過光合作用累積的有機物總量,NPP)是衡量生態系統碳匯能力的一項重要指標[3-4],當前對生態系統碳匯的主流核算方法(模型)有統計模型、光能利用模型、生態過程模型等[5-6],碳匯量化方法是碳匯研究的核心課題之一。 遙感技術具有非侵入性且獲取圖像范圍大、頻率高等特點[7-8],成為快速獲取地面樣地土壤、土地利用、植被狀況等數據的有效途徑。 在NPP 時空分布特征研究中,大都采用數理統計分析方法[9]。 在 NPP 變化的影響因子探究中,大都選取生態系統 NPP 涉及的自然因子進行相關性分析,或在全球氣候變化的大背景下探究極端氣象事件對 NPP 的影響[10-13]。 雖然多數學者普遍認為氣候與水土環境對NPP 的變化有極大影響,但是對影響因子的作用效果進行分析時忽視了各影響因子間的相互關系、相互作用,將各因子視為獨立個體,這顯然不合理且容易造成偏差。 愈來愈多的學者已認識到這一問題,但還缺少針對性的研究方法及研究結論。

基于博弈論視角的 SHAP 可解釋性方法(模型),可以消除研究中依賴影響因子間線性關系造成的偏差[14-15],對預測結果提供透明和可理解的解釋,在高維度機器學習算法中具有不可替代的優勢。 本文基于SHAP 可解釋性分析方法,探討黃河流域 NPP 的時空異質性及各區域 NPP 變化的主導因子,以期豐富對NPP 的研究并為黃河流域生態系統良性發展、碳匯能力提升等提供參考。

1 研究區概況與研究方法

1.1 研究區概況

黃河流域總面積 79.5 萬 km2 ,橫跨我國東中西部,是我國重要的生態屏障,也是經濟社會發展的重要區域,在國家發展大局和社會主義現代化建設全局中具有舉足輕重的戰略地位[16]。 依據氣候條件、生態功能等[17],本研究將黃河流域分為 3 個生態區域(見圖1),其中:Ⅰ區為中溫帶半干旱地區,以草原、沙地和丘陵為主,有沙漠化現象,地勢較平坦但面臨水土流失和生態退化問題;Ⅱ區為暖溫帶半濕潤地區,以平原和低山丘陵為主,氣候溫暖、濕潤,植被類型豐富,森林和草原混合,農業發達;Ⅲ區為高原半干旱地區,多為高山、盆地,植被為較原始的戈壁草原和干旱草原,受人類干擾較為強烈。

圖 1 黃河流域生態系統分區

Fig.1 Ecosystem Division of Yellow River Basin

1.2 研究方法

1)NPP 計算方法。 計算 NPP 的常用方法(模型)有光合有效輻射法、植物生長模型、生態模型等。 本研究采用改進的 CASA 模型[18]計算 NPP ,該模型是對經典 CASA 模型( NPP=APAR×ε ,式中 APAR 為光合有效輻射吸收量、 ε 為光能利用率)的優化和擴展,根據植被類型和環境條件動態調整光能利用率、依據更精確的氣象數據(如地表溫度、降水量、太陽輻射強度)優化模型參數、更詳細地模擬碳循環過程,通過高分辨率數據和多源數據融合,提高 NPP 估算精度。

2)Sen’s 斜率估計和 Mann-Kendall 檢驗。 Sen’ s斜率估計是一種非參數統計方法,用于分析時間序列的變化趨勢(斜率),具有對異常值不敏感、穩健性強等特點,廣泛應用于環境科學、水文、氣象等研究領域,公式為 (式中:median[ ]為中位數函數, i,j 為時間序列樣本序號, xj、xi 為時間序列樣本值),斜率 β 大于 0 表示時間序列呈上升趨勢、小于 0 表示時間序列呈下降趨勢。 Mann-Kendall 檢驗是一種非參數檢驗方法[19],用于檢驗時間序列是否存在顯著上升或下降趨勢,廣泛應用于環境、水文、氣象等領域。 本研究將 Sen’s 斜率估計和 Mann-Kendall檢驗結合起來(稱為 Sen’s-MK 檢驗),對黃河流域 3個生態區域 NPP 變化趨勢進行檢驗和分析。

3)SHAP 可解釋性分析。 NPP 受多因子影響,之前的研究通常將多源數據合并到一個復合向量(時間序列)中進行分析,且習慣性地將不同的變量(特征)視為獨立的,忽視了不同變量間的內在聯系,因而難以正確捕捉時間序列的變化趨勢。 本研究采用 SHAP 可解釋性分析方法對黃河流域不同區域 NPP 變化趨勢進行研究,該方法是基于博弈論和機器學習的可解釋性方法,用于解釋機器學習模型的預測結果,通過確定每個因子(變量)對模型預測的貢獻來幫助理解模型的預測結果,即將模型的預測值解釋為每個輸入變量(特征)的歸因值之和(這也是可加特征歸因方法的定義),用公式表示為

式中: g 為解釋模型預測值, φ0 為常數, ∴M 分別為變量(特征)編號、變量(特征)數量, Zj' 為表征變量 j 是否存在的變量(存在取 1,不存在取 0), φj 為變量 j 的歸因值(稱為 Shapley 值或 SHAP 值,其值越大表明對NPP 變化的貢獻越大)。

1.3 數據來源與處理

NPP 的計算通常依賴于氣象數據及土地利用數據等,本文研究時段為 2000—2020 年,采用的數據來源見表1(均為成熟的遙感圖像數據產品)。 對多源遙感圖像數據進行重采樣(合成或插值調整)至分辨率為 1km 。

表 1 數據來源

Tab.1 Data Sources

2 結果與分析

2.1 不同時間尺度 NPP 的時空演進特征

采用 CASA 模型計算年、月兩種時間尺度的 NPP值,采用 Mann - Kendall 檢驗方法對 2000—2007 年、2007—2014 年、 2014—2020 年 及 整 個 研 究 時 段(2000—2020 年)年尺度 NPP 變化趨勢進行顯著性檢驗,結合 Sen’s 斜率估計評估不同階段年尺度 NPP 變化趨勢和強度;提取各自然因子月尺度數據并計算逐月 NPP,進一步分析 NPP 受自然因子季節性影響呈現出的周期性變化特征。

2.1.1 年尺度 NPP 演變情況

2000 年、2007 年、2014 年、2020 年 4 個代表年份NPP 計算結果見圖 2,NPP 多年均值見圖 3。 由圖 2、圖3 可以看出,就整個流域而言,NPP 在研究期呈現逐漸上升趨勢,且表現出較為顯著的空間異質性,黃河干流沿線 NPP 較同緯度其他區域的高,流域整體高值出現在流域東南部;從 3 個生態區域比較來看,Ⅱ區NPP 明顯高于Ⅰ區和Ⅲ區、Ⅰ區最低;從各生態分區內部看,3 個生態區域內部均有明顯的空間異質性。

圖 2 代表年份 NPP 計算結果

Fig.2 NPP Calculation Results of Representative Years

?

圖 3 NPP 多年均值

Fig.3 NPP Multi?Year Average

對黃河流域年尺度 NPP 的 Sens-MK 檢驗結果(見圖4)表明:在研究期內黃河流域 NPP 整體呈現較為明顯的上升趨勢,NPP 上升的面積占比為 47.17% 、下降的面積占比僅 3.33% ;NPP 上升最為顯著的是Ⅱ區(其 13.81% 的面積為顯著上升、 55.93% 的面積為極顯著上升)、其次為Ⅰ區(上升面積占比為 28.61% )、再次為Ⅲ區(上升的面積占比為 19.86% ),值得說明的是Ⅰ區、Ⅲ區的上升面積占比雖然不高但是下降的面積占比也較小。 研究時段初期(2000—2007 年)3個生態區域 NPP 均明顯上升,但Ⅲ區在研究時段中后期(2007—2014 年、2014—2020 年) NPP 上升的面積占比減小甚至在2014—2020 年有 10.79% 的面積出現了下降。 出現上述現象的原因可能是研究時段初期因NPP 較低而上升較為容易,研究時段后期因 NPP 相對較高甚至達到了峰值而進一步上升的難度較大。

圖 4 Sen’s-MK 檢驗結果Fig.4 Sen’s?MK Test Results

綜上所述,2000—2020 年黃河流域NPP 呈現上升趨勢但上升速度逐漸減緩,顯著上升區域主要集中在Ⅱ區;黃河流域 NPP 在空間上分異明顯,最高值出現在Ⅱ區南部,Ⅰ區 NPP 整體較低,Ⅲ區 NPP 介于Ⅱ區與Ⅰ區之間,各生態區域黃河干流沿線 NPP 均相對較高。

2.1.2 月尺度自然因子及 NPP 演變情況

鑒于影響 NPP 的自然因子具有明顯的季節性,根據月尺度自然因子數據計算了逐月NPP。 圖5 為3 個生態區域地表溫度、太陽輻射強度、降水量、NDVI 及NPP 各月份多年平均值年內變化情況。

?

圖 5 各因子及 NPP 各月份多年平均值年內變化情況Fig.5 Variation of Each Factor and Npp in EachMonth of Multi?Year Mean Annual

由圖 5 可以看出,3 個生態區域月尺度 4 個因子及 NPP 年內季節性變化均較為明顯,反映出受季節變化強烈影響的共性,各生態區域 NPP 及其影響因子差異顯著。 從地表溫度來看,年內變化幅度大,夏季高、冬春低,溫差Ⅰ區最大、Ⅱ區和Ⅲ區相對較小,各生態區域地表溫度高低順序為Ⅱ區 gt; Ⅲ區>Ⅰ區;從太陽輻射強度來看,具有較大的季節性波動,以冬春換季時最強、其他時段較為平穩,空間上Ⅰ區 gt; Ⅲ區>Ⅱ區,即從西北向東南遞減,這可能與地勢(海拔)、植被狀況、大氣濕度等的影響有關;從降水量來看,呈現顯著的季節集中性,主要集中于 6—9 月,空間上由西北向東南增加,降水量Ⅰ區最小、Ⅱ區與Ⅲ區大致相當,Ⅱ區降水的季節性變化最為顯著;從 NDVI 來看,呈現強烈的季節性波動尤其Ⅲ區波動最為劇烈,夏季大、冬春小,空間上呈現由西北向東南遞增的分異規律,Ⅱ區>Ⅲ區>Ⅰ區,原因是Ⅱ區和Ⅲ區植被覆蓋度較高、Ⅰ區植被覆蓋度較低;從 NPP 來看,與地表溫度、降水量、NDVI 一致的季節波動性,即“春生夏旺秋衰冬歇”,空間上Ⅱ區 gt; Ⅲ區 gt;I 區,3 個生態區域間表現出明顯的梯度特征,這與各季節及各生態區域的氣候、植被條件的差異性相吻合,整體來看Ⅱ區是黃河流域自然環境條件相對較好、NPP 相對較高的區域。

2.2 黃河流域 NPP 空間分異的驅動因子

各生態區域 NPP 因自然環境條件不同而具有明顯的空間分異性,同時 NPP 與各自然因子存在不同程度的聯系(即各因子對 NPP 的變化具有不同程度的影響)。 圖6、圖 7 展示了對 3 個生態區域 NPP 變化的SHAP 可解釋性分析部分結果(圖 7 中各因子與自身的交互 SHAP 值稱為該因子的主效應,不同因子交互的 SHAP 值稱為副效應)。

圖 6 各因子 SHAP 值分布范圍

Fig.6 SHAP Value Distribution Range of Each Factor

Ⅰ區各因子SHAP 值的平均值(絕對值的平均值,下同)大小排序為 NDVI(6.605) >降水量 (2.358)gt; 地表溫度(0.858) >太陽輻射強度(0.698),其中 NDVI 的SHAP 值平均值最大、分布范圍最廣、主效應較強(8.240),且與其他因子有較強的交互效應(與降水量的交互效應最強,二者交互的SHAP 值為0.864),說明NDVI 對于該區 NPP 變化的影響最大(即 NDVI 對于該區 NPP 的變化有關鍵作用)。

圖 7 因子間的交互效應

Fig.7 Interaction Effect Between Factors

Ⅱ區各因子 SHAP 值的平均值大小排序為地表溫度 (15.43)gt;NDVI(3.25)gt; 降水量(2.35) gt; 太陽輻射強度(0.89),與Ⅰ區相比,各因子對 NPP 變化的貢獻均明顯提升、各因子間的交互效應相較Ⅰ區顯著提升,其中主效應最大的是地表溫度(20.20),其與 NDVI 的交互效應最強(SHAP 值為 1.55),說明該區域自然因子的交互更加緊密且對 NPP 影響顯著。

Ⅲ區各因子 SHAP 值的平均值大小排序為 NDVI(9.57) gt; 地表溫度(5.02) gt; 降水量(2.44) gt; 太陽輻射強度(1.03),NDVI 的主交互性最強( 達到 13.80),各因子間交互效應最大的是地表溫度與 NDVI 的交互效應(SHAP 值為2.03),說明該區各自然因子間交互也較為密切且對 NPP 變化的影響顯著。

從流域整體來看,NDVI 對 NPP 變化的影響最大且其與太陽輻射強度的交互效應較大,地表溫度對Ⅱ區 NPP 變化的作用遠大于其他因子(原因是該區植被較其他區域豐富,地表溫度對植被生長、光合作用與呼吸作用的影響集中反映在 NPP 這一量化指標上)。 黃河流域幅員遼闊、氣候類型多樣、植被類型豐富等,主要由自然環境條件決定的 NPP 值存在顯著的空間差異性,從上述 SHAP 可解釋性分析結果可以看出各因子在不同生態分區的作用有顯著差異,因此為了提升黃河流域 NPP,應因地制宜,根據各區域 NPP 的主要驅動因子制定相應的策略,如Ⅱ區應充分考慮地表溫度與其他因子的交互效應。

2.3 討論

上述分析表明,黃河流域 NPP 的時空分異性均較為顯著,不同生態分區 NPP 變化的驅動因子及其驅動作用也有顯著差異。 SHAP 可解釋性分析結果驗證了前人關于自然因子對 NPP 變化作用效果的研究結論,如氣候因子對 NPP 變化的作用巨大、夏季是 NPP 的高峰時期[20]。 SHAP 可解釋性分析方法是研究視角和方法的拓展與豐富,其基于各自然因子相互聯系的系統觀點得出的研究結論,與將各自然因子作為獨立的個體得出的研究結論“ 降水是決定 NPP 的關鍵因子”“溫度升高提升區域 NPP” [21-22]有一定區別。

鑒于高程、土壤類型、土地利用類型、人口等因子在月時間尺度上變化甚微,在后續的研究中,可引入這些因子對面狀 NPP 的年尺度變化進行 SHAP 可解釋性分析,以探究包括自然因子和社會因子的 NPP 變化歸因,進而診斷區域碳盈虧的關鍵因子。 此外,可開展其他因子如城鄉景觀差異等對 NPP 日尺度變化的影響等,以豐富對 NPP 的研究。

3 結論

1)從空間維度來看,黃河流域 NPP 分異性顯著,高值出現在Ⅱ區、低值廣泛分布在Ⅰ區和Ⅲ區,各生態區域NPP 均值大小排序為Ⅱ區 gt; Ⅲ區>Ⅰ區,黃河干流沿線 NPP 較同緯度其他區域的高。

2)從時間維度來看,黃河流域 NPP 在研究期波動上升、在年內具有明顯的季節性變化,研究期 NPP 上升的面積占比為 47.17% 、下降的面積占比為 3.33% ,3個生態區域中 NPP 上升最顯著的是Ⅱ區,研究期黃河流域 NPP 上升速度逐漸減緩。

3)從流域整體來看,NDVI 對 NPP 變化的影響最大且其與太陽輻射強度的交互效應較大,但各生態區域 NPP 變化的主要驅動因子有所不同,Ⅰ區為 NDVI、Ⅱ區為地表溫度、Ⅲ區為 NDVI。

4)黃河流域各生態區域提升 NPP 的策略應因地制宜、各有側重點,并充分重視各自然因子間的交互效應。

參考文獻:

[1] 中國應對氣候變化的政策與行動 2023 年度報告(摘編)[J].環境保護,2023,51(21):50-58.

[2] 左其亭,邱曦,鐘濤.“雙碳”目標下我國水利發展新征程[J].中國水利,2021(22):29-33.

[3] 陳舒婷,郭兵,楊飛,等.2000—2015 年青藏高原植被 NPP時空變化格局及其對氣候變化的響應[J].自然資源學報,2020,35(10):2511-2527.

[4] 劉洋洋,章釗穎,同琳靜,等.中國草地凈初級生產力時空格局及其影響因素[J].生態學雜志,2020,39(2):349-363.

[5] 張振東,常軍.2001—2018 年黃河流域植被 NPP 的時空分異及生態經濟協調性分析[J].華中農業大學學報,2021,40(2):166-177.

[6] 焦偉,陳亞寧,李稚.西北干旱區植被凈初級生產力的遙感估算及時空差異原因[ J]. 生態學雜志,2017,36( 1):181-189.

[7] TURNER W,RONDININI C,PETTORELLI N,et al.Free andOpen?Access Satellite Data are Key to Biodiversity Conservation[J].Biological Conservation,2015,182:173-176.

[8] BROWN S.Measuring Carbon in Forests:Current Status andFuture Challenges [ J]. Environmental Pollution, 2002, 116(3):363-372.

[9] 王菲,曹永強,周姝含,等.黃河流域生態功能區植被碳匯估算及其氣候影響要素[J].生態學報,2023,43(6):2501-2514.

[10] 王娟,何慧娟,董金芳,等.黃河流域植被凈初級生產力時空特征及自然驅動因子[ J].中國沙漠,2021,41(6):213-222.

[11] 施亞林,曹艷萍,苗書玲.黃河流域草地凈初級生產力時空動態及其驅動機制[J].生態學報,2023,43(2):731-743.

[12] 田智慧,張丹丹,赫曉慧,等.2000—2015 年黃河流域植被凈初級生產力時空變化特征及其驅動因子[J].水土保持研究,2019,26(2):255-262.

[13] 夏冰,馬鵬宇,徐聰,等.近 20 年黃河流域植被凈初級生產力時空分布及其對極端天氣變化的時空響應[J].水土保持研究,2023,30(2):256-266,284.

[14] HUANG H B,PAN Y,WANG C P,et al.Nonlinear FloodResponses to Tide Level and Land Cover Changes in SmallWatersheds[J].Land,2023,12(9):1743.

[15] FRANZKE C L E. Nonlinear Climate Change [ J]. NatureClimate Change,2014,4(6):423-424.

[16] 中共中央國務院印發《黃河流域生態保護和高質量發展規劃綱要》[J].中華人民共和國國務院公報,2021(30):15-35.

[17] TANG Han,WEN Tong,SHI Peng,et al.Analysis of Char?acteristics of Hydrological and Meteorological Drought Evolutionin Sou?thwest China[J].Water,2021,13(13):1846.

[18] 王菲.“雙碳”目標下黃河流域碳源/碳匯時空分布及碳盈虧分析[C]//張弛,彭勇.水科學與智慧水務:第十九屆中國水論壇論文集.:中國水利水電出版社,2022:342-348.

[19] 王保良,王鴻翔,張海濤,等.基于 IHA-RVA 法的黃河水文情勢演變分析[J].人民黃河,2024,46(5):33-39.

[20] CAO Y, WEN X, WANG Y, et al. The Analysis of NPPChanges Under Different Climatic Zones and UnderDifferent Land Use Types in Henan Province,2001 - 2020[J].Sustainability,2024,16(18):8096.

[21] ZHANG Meiling,LAL Rattan, ZHAO Youyi, et al. Spatialand Temporal Variability in the Net Primary Production ofGrassland in China and Its Relation to Climate Factors[J].Plant Ecology,2017,218(9):1117-1133.

[22] BHATTARAI P,ZHENG Z,BHATTA K P,et al. Climate?Driven Plant Response and Resilience on the TibetanPlateau in Space and Time:A Review[ J].Plants,2021,10(3):480.

【責任編輯 張智民】

主站蜘蛛池模板: 欧美性色综合网| 91综合色区亚洲熟妇p| 国产精品成人一区二区| 在线色综合| www亚洲天堂| 香蕉精品在线| 国产福利小视频高清在线观看| 日韩欧美中文字幕在线精品| 亚洲国产成人久久77| 亚洲一级毛片免费观看| 国产爽歪歪免费视频在线观看| 色呦呦手机在线精品| 亚洲αv毛片| 国产高潮流白浆视频| 亚洲男人天堂2018| 91探花国产综合在线精品| 国产制服丝袜无码视频| 91精品国产自产在线观看| 超碰91免费人妻| 亚洲最大综合网| 东京热av无码电影一区二区| aⅴ免费在线观看| 久久久久夜色精品波多野结衣| 五月天综合网亚洲综合天堂网| 中文字幕欧美日韩| 日韩欧美国产精品| 九九精品在线观看| a欧美在线| 美女一区二区在线观看| 91人妻在线视频| 高清免费毛片| 亚洲精品黄| 色偷偷av男人的天堂不卡| 九色国产在线| 亚洲美女一级毛片| 国产在线精彩视频二区| 狠狠ⅴ日韩v欧美v天堂| 久热中文字幕在线| 国产正在播放| 内射人妻无码色AV天堂| 亚洲婷婷在线视频| 亚洲天堂精品在线观看| 久久久久免费精品国产| 91久久精品国产| 日本午夜三级| jizz国产视频| 不卡无码网| 成人午夜视频免费看欧美| 日韩毛片免费视频| 特黄日韩免费一区二区三区| 丁香婷婷久久| 亚洲欧州色色免费AV| 日韩一区二区三免费高清| 九九这里只有精品视频| 精品视频91| 99久久免费精品特色大片| 国产成人一区免费观看| 99国产精品国产| 日韩一级二级三级| 亚洲色欲色欲www网| 99精品国产自在现线观看| 亚洲人成在线精品| 欧美激情福利| 91九色视频网| 欧美三级日韩三级| 国产网站免费观看| AV熟女乱| 91原创视频在线| 国内毛片视频| 波多野结衣无码AV在线| 欧美精品1区| 国产无吗一区二区三区在线欢| 这里只有精品国产| 国产一区在线观看无码| 国产女人在线视频| 伊人久久大香线蕉aⅴ色| 欧美特黄一免在线观看| 中国国产一级毛片| 久久免费看片| 亚洲免费三区| 精品自窥自偷在线看| 亚洲国产91人成在线|