?

氣候變化下水文模擬不確定性若干問題討論

2023-02-04 11:32趙自陽王紅瑞楊亞鋒李曉軍
水資源保護 2023年1期
關鍵詞:不確定性水文氣候變化

張 力,趙自陽,王紅瑞,楊亞鋒,2,李曉軍

(1.北京師范大學水科學研究院, 北京 100875; 2.華北理工大學理學院, 河北 唐山 063210;3.濟南市水利工程服務中心, 山東 濟南 250099)

氣候變化和水文不確定性是21世紀最大的挑戰之一,同時也是水文領域23個未決科學問題中的研究熱點[1-2]。聯合國政府間氣候變化專門委員會(IPCC)在第六次評估報告中指出,到21世紀末降水模式和溫度將發生顯著變化,進一步加劇水循環速度,與之相關的水資源供應和需求方面的不確定性不斷增加,為決策者統籌水資源管理系統帶來了更為嚴峻的挑戰[3-4]。然而,水資源系統內在的灰色、模糊以及隨機屬性為其組成和影響因素帶來廣泛不確定性[5],這種不確定性及其相互作用可能導致規劃工作的額外復雜性,并影響后續決策過程。水資源系統的多周期、多層次和多目標特征也可能進一步放大這些不確定性,引發用水者之間的潛在沖突,加劇缺水危機[6-7]。在氣候變化的背景下,預計到2025年將有18億人生活在缺水國家或地區[8],致使正確考慮不確定性已成為水資源系統規劃的基礎[9-10]。此外,聯合國教科文組織在2012年發布的第四版《世界水資源發展報告》中強調,歷史經驗將不足以估計用水量與未來需求變化之間的關系,水資源調控需充分考慮環境影響中的風險和不確定性[11]。

倘若與水文模型相關的預測不確定性結果沒有通過氣候變化影響分析得到正式處理和傳播,那么這些研究中得出的水文變化結論可能被夸大,從而誤導水資源管理決策者[12]?,F階段在全球和區域范圍內已發表許多有關水文系統不確定性主題的研究[13-14],學者們也開發了多種數學規劃方法用于處理水資源管理問題中的不確定性[15]。如金菊良等[16-17]分別采用投影尋蹤權重優化、聯系數和耦合協調度等方法,降低了水資源承載力評估中的不確定性風險;張金萍等[18]基于小波變換和集對分析研究了鄭州市降雨-徑流的不確定性關系,為該地區徑流預測提供了理論支撐;王京晶等[19]應用GLUE(generalized likelihood uncertainty estimation)算法討論了基流分割前后流量數據對城市雨洪模擬的不確定性影響。然而,不確定性分析在許多模型構建過程中仍然不是標準的做法,實際應用過程中向決策者們展示沒有充分考慮不確定性限制結果的方式仍屢見不鮮[20-21]。鑒于已有研究的不足,本研究基于不確定性本質,從氣候變化視角全面綜述氣候變化下水文系統模擬中的不確定性,進而展望該領域的未來研究方向,以期為深入探究區域水循環機理、指導防災減災、加快推進韌性城市建設等提供參考。

1 不確定性的來源與分類

廣義上,不確定性可以被視為信息的一種屬性[22]。Knight[23]將不確定性定義為由于決策者認知不足而不能將概率分配給結果,在水文科學背景下,這個定義為水文領域研究者提供了標準。Funtowicz等[24]將水資源系統中的不確定性描述為不準確、不可靠和近乎無知這3類信息不足的情況;Larry[25]將其定義為事件可能發生在可控范圍之外,是水文水資源系統的基本屬性。為了涵蓋水文不確定性的所有維度,IPCC第五次評估報告中將不確定性定義為:由于先驗信息不足或對已知事物存在分歧而導致的一種不完整的認知狀態。

地球上具有廣闊的海洋面積和熱容,全球平均表面溫度的細微變化都需要大量的熱能來驅動。與1900—2000年的平均地表溫度相比,1880年以來地表溫度每10年上升0.08℃,1981年以來變暖速度進一步加快,每10年上升0.18℃,2020年是有史以來地表溫度第二高的年份,陸地面積創下歷史新高(圖1)。IPCC第四次評估報告中指出這些額外的熱量導致了區域性和季節性的極端溫度,加劇了強降雨,改變了區域水文循環規律。在氣候變化對水資源影響的量化中,研究者大多基于降尺度技術對某一特定情景下的全球氣候數據進行預處理,并應用其驅動水文模型來闡釋水文過程的響應機理。而由于自然氣候系統和流域水文循環過程具有不同時空尺度的非線性反饋,致使不確定性廣泛存在于氣候變化情景、評價模型以及評估過程中[26]。盡管氣候變化下水文系統中的不確定性本質上難以量化,但是了解不確定性的來源和影響對于理性決策以及提高對水資源風險和可靠性的認知仍然十分重要[27]。

圖1 1880—2020年地表溫度變化Fig.1 Surface temperature change from 1880 to 2020

1.1 氣候變化情景的不確定性

水文系統預測變化的不確定性源于氣候系統的內部可變性。氣候情景作為表征人類和環境相互作用系統中不確定性的手段,在結合水文模型和水資源評價模型的基礎上,已經成為評估氣候變化對水文過程影響的前沿方法[28]?,F有氣候情景可分為4類:增量情景、慣性情景、類比情景、GCMs(general circulation models)情景。增量情景是假定到未來某一特定時期氣候要素的變化量,特點是操作簡單,但不具備氣候、水文預測功能;慣性情景是假定未來氣候變化延續已有變化趨勢,特點是忽略了氣候的不同周期及其組合;類比情景又稱為古氣候比擬法,基于分析千萬年時間尺度下的氣候數據進行類比,特點是難以獲得可信資料,且人類活動對自然影響越來越大,難以準確衡量;GCMs情景是利用全球氣候模式模擬未來的氣候變化情況,能反映地球系統各圈層復雜關系和溫室氣體排放響應??紤]到前3個氣候情景的不足,將GCMs情景作為本文討論的核心,造成其不確定的關鍵來源可從氣候模式、情景設定以及降尺度方法方面考慮[29]。

1.1.1氣候模式不完善性

GCMs的結構、物質和能量交換過程和代表性濃度路徑是影響氣候模式模擬結果的主要不確定性因素[30]。GCMs結構由于情景的參考運行和擾動運行之間的差異多變,模型參數化處理難度增大,致使預測結果存在顯著差異性[31];海洋、大氣與大陸三者間的物質和能量交換機制為評估過程帶來不確定性,一些簡單模式為了簡化這些尺度上的物理過程,對海冰和大氣運動進行粗糙處理造成誤差;區域土地利用變化和氣溶膠也可能復雜化不同尺度上的歸因,造成一定誤差[32-33]。因此,完善區域氣候模式和全球海陸氣耦合模型已成為現階段精準預測全球或者區域變化的重中之重。在IPCC第六次評估報告中,盡管應用了更為全面的海洋-大氣耦合環流模型(coupled ocean-atmosphere general circulation model,CGCM),提高了其模擬精度和參數優化方案,但與實際氣候狀況相比其模擬水平仍需優化改進。

1.1.2情景設定不確定性

情景不確定性指由未來經濟社會和人類活動的不可預測性導致未來溫室氣體(greenhouse gas,GHG)排放預測的不確定[34]?,F階段應用較為廣泛的排放預測情景為1992年發布的IS92情景和后續提出的經濟社會共享情景(shared socioeconomic pathways,SSPs),其中SSPs情景應用于IPCC第五、第六次評估報告。使用預測情景進行影響研究代表了一種自上而下的方法來解決問題,其對于真實反映氣候變化下水文系統的脆弱性、敏感性來說至關重要。因此,目前影響氣候變化評估可靠性的最大問題,即精準可靠地研發未來幾十年甚至上百年的氣候情景。水文模擬不確定性的另一關鍵是選擇氣候情景,單一氣候變化情景下降尺度輸出只能得出許多現實可能性中的單一軌跡。然而,這種單一的軌跡本身并不能完全代表未來復雜水文過程下的水文情景,更不足以用于準確預估氣候變化帶來的水文影響。因此,在影響評估中常采用系列情景而非單一的最佳猜測或者平均水平情景。Simonovic等[35-36]指出,由于使用的未來排放情景存在差異,氣候對防洪影響的研究尚存在不確定性,該研究利用GCMs的輸出來評估防洪系統的有效性并得出結論,不同的氣候情景提供了不同的水文參數估計值。

1.1.3降尺度技術不確定性

早期情景不確定性的來源僅通過GCMs的結構差異或排放情景設定來解釋[37],氣候模式降尺度方法的選擇作為另一不確定性的重要來源,很少受到關注[38-39]。一般來說GCMs的空間尺度可達幾百千米,時間上也是以年或月作為基本單位,相比較之下水文模型的時空尺度更為精細化。因此從時空分辨率的角度,GCMs的輸出數據與水文模型需求不一致,導致其不能用于模型驅動來源,無法科學支撐未來水資源評價,亟須降尺度處理以滿足氣候模式和水文模型之間接口匹配性問題。目前已經開發了基于動態降尺度的區域氣候模型(regional climate models,RCMs)和統計降尺度(statistical downscaling,SD)方法來滿足這一要求[40]。RCMs的開發基于動態公式,使用GCMs的初始和與時間相關的橫向邊界條件,以犧牲有限區域建模為代價實現更高的空間分辨率,主要缺陷是計算成本較大,僅適用于部分研究區。盡管有所改進,RCMs的輸出尺度對于某些實際應用來說仍然過于粗糙,例如小流域水文響應機制分析和田間農業影響研究,需要當地和特定地點的精細化氣候情景。此外,根據不確定性爆炸理論,RCMs建模中的每一部分都會對后續級聯產生影響[41]。SD方法通過將代表大尺度一些變量的狀態與其他預測變量的狀態聯系起來,具有靈活性高、計算方法相對容易實現的優點,因此應用較為廣泛[42]。然而,無論使用哪一種降尺度方法,計算結果中始終存在誤差,如何對數據進行降尺度處理也是氣候變化情景應用的不確定性根源之一。

1.2 評價模型的不確定性

在氣候變化對水文過程影響的研究中,常用手段是借助率定好的水文模型,結合插值、區域氣候因子空間降尺度轉換等方法將氣溫、降水作為輸入數據,獲取水文模型的產匯流結果,進而分析氣候變化對徑流的貢獻[43-44]。然而,評價過程中與水文模型相關的不確定性往往被忽略[45]。事實上,流域水文過程不僅受到降水、氣溫、日照、相對濕度等氣候因子的影響,同時也受到地形地貌、植被覆蓋等下墊面因子的影響。而流域水文模型通過模型結構、參數和輸入數據的數學描述方式對這一復雜系統形成機理和變化規律進行概化,其計算過程和結果必然會衍生出一系列的不確定性。

1.2.1模型結構不確定性

模型結構不確定性反映了仿真模型或設計程序無法精確表示系統真實的物理機制或過程(圖2),主要原因有兩方面[46]:一是缺乏對變化環境下水文物理循環機制的足夠認知,導致無法構建反映真實水循環機理的流域或區域水文模型。以下墊面變化顯著的黃土高原地區為例,賀振等[47]認為該地區下墊面條件較差,無法有效滯留地表徑流,因此得出黃土高原地區產流機制以超滲地面徑流為主導這一共識。然而胡彩虹等[48]基于水文學原理和數理統計等方法,對黃河中游典型流域場次洪水的降雨徑流特征進行了分析,初步表明蓄滿產流以及混合產流模式可能存在于黃土高原地區,那么在該區域水文模型產流模塊構建時,僅考慮一種產流機制將無法保證模擬結果的可靠性。二是由于水文過程極其復雜,包括截留、填洼、下滲、蒸散發等多個過程[49],同時在氣候因子和下墊面因子的共同影響下,加劇了精細刻畫水文過程形成機理和認知變化規律的難度[50]。盡管部分學者認為,模型結構不確定性影響的重要性遠不如氣候變化[51-53],然而這一結論可以部分歸因于所使用的小樣本量。通過量化對比結構和參數的不確定性,H?jberg等[54]表明,結構不確定性對模型性能的影響占主導地位,并且無法通過優化參數不確定性補償;Rojas等[55]也得出相似結論,結構不確定性可能占預測結果不確定性的30%以上。

圖2 模型結構不確定性示意圖Fig.2 Schematic diagram of model structure uncertainty

1.2.2模型參數不確定性

將特定的參數集作為概念簡化結果,對于每個水文模型來說都是必要的[56],但不是所有模型參數都能通過觀測試驗直接獲取,參數不確定性是指在集成和概念化過程中無法估計這些有效參數的結果。在構建水文模型模擬時由于缺乏相應的觀測數據而不得不借助參數估計的方法基于已有數據間接確定,給出滿足模擬要求的模型參數取值[57]。參數不確定性也可能是水文系統空間異質性和觀測誤差的結果,例如某些參數(如水力傳導率)可以在點尺度上進行測量,但它們的值在流域尺度上差異很大,盡管在操控模型時采用一套參數的模擬效果能獲得令人滿意的精度要求[58],但是對整個流域做均一化處理,在物理意義方面勢必會有所欠缺。有模型構建者提出基于現場觀測試驗和先驗信息能夠直接率定出合適的模型參數[59],然而,縱使在高密度試驗條件下得到某一特定時空尺度下的參數結果,往往也難以達到令人滿意的模擬效果。另外,即便模型能夠精細刻畫水文循環過程,由于校準數據中的錯誤,也可能產生參數不確定性。

1.2.3模型輸入不確定性

水文系統中模型的輸入數據也是高度不確定的(圖3),這種不確定性的影響分為3類:①測量誤差。例如,水文模型中影響模擬結果最重要的輸入資料是降水數據[60],然而已有長時間序列降雨徑流模擬研究中通常存在一個問題:由于20世紀60年代及以后的一段時間內,對降雨的觀測技術還很落后,水文站一般是利用自制雨量筒結合刻度尺對降雨進行觀測后編入當地水文年鑒,致使原始數據與實際降雨過程有較大出入。因此,基于這類數據得出的模擬結果可靠性不高。②觀測站點空間布局無序,數據樣本缺乏代表性。準確計算水文氣象數據的空間分布值是復雜的,為了提高這些數據的質量,最好的辦法是增加監測站點的密度,然而在高程、坡度等不規則的區域往往是難以實現的[61-62],同樣構成輸入不確定性問題。③數據處理和轉錄誤差[63]。水文建模通常涉及使用觀測流量數據評估模型精度,然而實際河道中流量數據難以直接連續測量,而連續監測水位高度的方法較容易實現,所以通常只觀測水位,流量則由水位-流量關系曲線轉換得來[64]。這種轉換不僅傳播隨機誤差,還受到水位-流量關系曲線插值和外推等方法的影響[65]。已有研究表明,外推誤差在傳統水位-流量關系曲線計算方法的不確定性中占據主導,導致流量數據的不確定范圍為5%~25%[66-67]。

圖3 模型輸入不確定性示意圖Fig.3 Schematic diagram of model input uncertainty

1.3 評價過程的不確定性

1.3.1人類活動影響不確定性

人類活動為氣候變化提供了更多不穩定因素[68-69]。通常,氣候變化會導致大氣溫度上升和降水模式的改變,而人類活動對水循環過程有直接(如取水)或間接(如土地利用變化和水利工程建設)的影響[70]。隨著人口膨脹和經濟社會發展,人類活動在水文過程中的重要性很可能超過氣候變化[71],因此,區分它們對水循環的驅動效應對于為水資源管理決策提供科學支撐來說非常重要。人類活動主要從3個方面影響區域水循環機理:①水土保持措施改變下墊面狀況從而使原有產匯流范圍和路徑發生改變,坡面產流過程的糙率系數也隨之改變,縱使氣候因子不發生顯著突變,未來的水文情勢也會或多或少受其影響[72-73];②跨流域調水、開采地下水等引水工程直接影響區域用水結構;③修建水利工程導致河網縱向連通性降低,影響河道匯流過程[74]。在實際過程中這三者共同作用且相互影響,評估其作用時應適當將其分割開來單獨考慮,而不是將所有人類活動歸為一個影響因子。

1.3.2概率分布函數不確定性

評估極端氣候條件下流量的演變情勢通常需要進行水文頻率計算,因此需要考慮選取合適概率分布函數的不確定性。其不確定性主要來源于概率分布函數線型選擇、分布參數的確定以及水文資料樣本選取等。目前水文頻率計算常用的概率分布函數有皮爾遜Ⅲ型分布、廣義極值分布、三參數對數正態分布等[75](表1)。由于洪水產生的主導因素不同,各流域選取的分布線型也存在差異。參數估計方法包括矩法、適線法以及權函數法等,方法的選擇會對結果產生顯著影響[76]。此外,水文序列的長度、系列代表性及樣本抽樣方法對水文頻率計算結果的影響也不容忽視。杜鴻等[77]應用4種抽樣方法和2種極值統計模型,結合K-S(Kolmogorov-Smirnov)法分析淮河流域極端流量時空演變規律,表明兩種統計模型均能較好地擬合極端流量,百分位閾值法的模擬精度優于其他超閾限峰值法。

表1 水文頻率計算常用的概率分布函數Table 1 Commonly used probability distribution functions for hydrological frequency calculation

1.3.3氣象水文耦合技術不確定性

現有氣象-陸面水文全耦合研究技術還有長足發展空間[78],傳統的氣候水文耦合研究大多只考慮了大氣和陸面過程單向反饋過程,采用不同類別的模型對大氣或陸面過程單獨進行研究,操作過程較為單一且應用簡單,應用范圍相對較廣[79-80]。但流域水文要素變化與區域氣候存在互饋機制,缺乏二者間雙向耦合的理論方法與關鍵技術研究,對水文和氣象兩個基礎學科的交叉融合造成阻隔,對研究成果的完整性和精確性也產生了一定誤差??紤]到大氣與水文之間存在互饋作用,一些學者通過將氣候模式與水文模型協同編譯運行并且不破壞二者間的獨立性,進一步完善了陸面水文與氣候全耦合技術。如Senatore等[81]以意大利南部的一個流域為研究區,通過集成氣象研究預報模型(weather research and forecasting,WRF)和一套水文物理參數構建陸面水文-氣候全耦合模型,并將模擬結果與原始非耦合模型橫向對比,表明全耦合模型在觀測降雨以及水循環模擬方面性能更好;Wagner等[82]采用氣象與水文全耦合的方式,將WRF和HMS水文模型擴展為中尺度模型系統,在鄱陽湖流域進行應用和檢驗,結果表明全耦合模型系統能夠實現地表和地下的側向徑流以及地下水、非飽和帶、地表和大氣之間的雙向相互作用。

2 氣候變化對水文過程影響不確定性的研究進展

2.1 敏感性分析法

敏感性分析用于量化模型參數調整或輸入數據誤差對運行結果的影響[83],可分為局部和全局敏感性分析[84]。局部敏感性分析的核心思想是控制變量,即每次只修改單個輸入參數,通過計算模擬結果的變化進而確定所選參數對分析結果不確定性的貢獻?;诖?,水文敏感性被定義為模型對已知氣候變化因子的響應,即降水和溫度的變化導致作物蒸散、徑流和地下水補給的變化。這種方法簡單易行,例如Bao等[85]通過假定氣候影響下的降雨變化(以2%為步長,變化范圍為-20%~20%)和溫度變化(以0.2℃為步長,變化范圍為-2~2℃),對淮河流域4個子匯水區的徑流、蒸發和土壤濕度進行了敏感性分析,結果表明水文變量對降雨更為敏感,其次是溫度,但其缺陷在于忽略了模型參數之間的內在相關性。全局敏感性使用不同GCMs的集合或者在參數略微改變的情況下多次運行來確定,優勢在于可同時分析模擬結果對多個參數或輸入因子的響應[86]。然而,GCMs設計的結構差異,如大氣和海洋的垂直和水平分辨率以及各種過程的參數化,為計算引入了更多的不確定性。

2.2 概率分析法

概率是研究者最為熟悉和應用最為廣泛的量化氣候變化下水文不確定性的手段。目前應用較廣的概率分析方法有:基于貝葉斯理論的蒙特卡羅(Monte Carlo)隨機模擬[87]、廣義似然不確定性估計法(GLUE)[88]和方差分析(analysis of variance,ANOVA)法[89]。蒙特卡羅法不僅可以生成多種情景,同時可以從概率分析的角度來衡量不同情景的權重。Guo等[90]曾使用該方法生成具有不同參數集的徑流序列,并通過非參數方法估計徑流的概率密度函數,然后針對不同的GCM輸出和假設氣候情景定量估計模型參數和徑流不確定性。隨后研究者又結合隨機過程中的馬爾可夫過程衍生出馬爾可夫鏈-蒙特卡洛方法(Markov chain Monte Carlo,MCMC),極大地提高了蒙特卡洛方法的采樣效率,成為水文氣候科學領域不可或缺的工具[91]。GLUE方法在似然函數的選取與閾值確定方面有著更高的主觀性和靈活性,是目前最常用于不確定性估計的經驗頻率方法。ANOVA法將總的集合不確定性進一步分解為來自不同來源(情景、模型結構等)的貢獻率以及它們之間的相互作用[92],比其他概率分析方法所需要的假設更少,并且可以很容易地用于量化與模型-情景耦合相關的不確定性。Bosshard等[93]基于ANOVA法定量計算了氣候模型、統計處理方法和水文模型對萊茵河近景和遠景下水文模擬不確定性的貢獻率,結果表明,總體不確定性不是由3個因素的個體不確定性累加而來,而三者之間的交互作用占總體不確定性的5%~40%。

2.3 氣候模型和情景模擬分析方法

并非所有的不確定性都能以概率的形式描述。情景分析是目前氣候影響研究中最依賴的技術,其通過組合能夠模擬可能發生的不同情景,例如采用不同的氣候模式、排放情景、降尺度方法和水文模型均可視為不同的情景。田燁[94]采用4種未來排放情景、3種氣候模式、3種水文模型模擬未來情景下的最大徑流,發現水文模型參數對未來極端徑流不確定性的貢獻作用最大,GCMs的貢獻作用僅次于水文模型參數,水文模型結構的貢獻率最低;Jung等[95]考慮了GCMs結構、未來GHG排放情景、土地利用變化情景、自然變化和水文模型參數這5種不確定性,采用兩種溫室氣體排放情景下的兩種未來土地利用條件來量化城市洪水頻率變化的不確定性,結果表明流域洪水頻率變化對氣候變化的敏感性高于土地利用變化。然而,在情景分析中會不可避免地遇見兩個問題:一是目前情景設定的發展略滯后于氣候模型,研究者們應該致力于解決這一發展不平衡問題,以提供最新的影響評估情景;二是大多數氣候情景構建方法是基于模型的氣候變化估計與觀測氣候數據相結合,而由于觀測數據集難以涵蓋年代際尺度氣候變化的全部范圍,以及由于使用不同的方法將模式和觀測氣候數據結合起來,進一步將不確定性引入到氣候情景中。在氣候情景中,這些與使用觀測氣候數據有關的不確定性通常被忽略。

3 研究展望

不確定性仍然是可持續水資源管理的挑戰,本研究從氣候變化情景、水文模型和評價過程3方面討論氣候變化下水文模擬中的不確定性,概述了目前量化不確定性的主流方法,未來對于氣候變化下水文模擬的不確定性應著重開展以下研究:

a.極端水文事件研究。極端氣候事件伴隨著全球變暖頻頻發生,如2021年7月河南遭受強暴雨,鄭州及周邊地區多個站點的日降雨量突破歷史極值,造成嚴重的暴雨洪澇災害。鑒于目前氣候模式對極端氣候事件的模擬效果不盡如人意,而復雜網絡可對此類非線性演變特征進行深入剖析,可考慮應用其增強對極端氣候事件預估的可靠性。

b.輸入數據冗余性處理。水文系統模擬中需要大量輸入數據驅動模型,是否可以從不確定性人工智能角度中的貝葉斯網絡出發,通過局部信息采集、特征識別、目標修正等技術達到整體認知,科學處理各種輸入數據時間窗問題和冗余性,從而為無資料地區徑流預測提供科學支撐是亟待解決的問題。

c.水文過程異方差性。變化環境下流域水文序列存在波動集群效應,具有較強的異方差性和波動性。徑流模擬過程中眾多動態不確定性因素的累積傳導效應會造成預測結果的不可靠,然而其在水文領域中并沒有得到應有的關注。因此,基于異方差性波動模型揭示變化環境下非平穩異方差性水文序列的發生規律是解決當下水科學問題的現實需求。

猜你喜歡
不確定性水文氣候變化
法律的兩種不確定性
《應對氣候變化報告(2022)》發布
繼往開來 守正創新——河北省水文工程地質勘查院
水文
水文水資源管理
英鎊或繼續面臨不確定性風險
氣候變化與環保法官
氣候變化:法官的作用
水文
具有不可測動態不確定性非線性系統的控制
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合