亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        基于TRMM數(shù)據(jù)提取降雨侵蝕力的校正方法
        ——以薊運(yùn)河上游地域?yàn)槔?/h1>
        2015-06-15 18:31:44靳秋桐史明昌張舉濤王珊胡影
        中國(guó)水土保持科學(xué) 2015年4期
        關(guān)鍵詞:雨型降雨量分辨率

        靳秋桐,史明昌,?,張舉濤,王珊,胡影

        (1.北京林業(yè)大學(xué)水土保持學(xué)院,100083,北京;2.北京林業(yè)大學(xué)水土保持與荒漠化防治教育部重點(diǎn)實(shí)驗(yàn)室,100083,北京;3.北京地拓科技有限責(zé)任公司,100084,北京;4.高德軟件有限公司,102200,北京)

        基于TRMM數(shù)據(jù)提取降雨侵蝕力的校正方法
        ——以薊運(yùn)河上游地域?yàn)槔?/p>

        靳秋桐1,2,史明昌1,2,3?,張舉濤1,2,王珊4,胡影3

        (1.北京林業(yè)大學(xué)水土保持學(xué)院,100083,北京;2.北京林業(yè)大學(xué)水土保持與荒漠化防治教育部重點(diǎn)實(shí)驗(yàn)室,100083,北京;3.北京地拓科技有限責(zé)任公司,100084,北京;4.高德軟件有限公司,102200,北京)

        由于TRMM數(shù)據(jù)的時(shí)間分辨率較低,導(dǎo)致用TRMM數(shù)據(jù)計(jì)算某一區(qū)域指定時(shí)間段內(nèi)的降雨侵蝕力存在較大的系統(tǒng)誤差。為解決該問(wèn)題提出通過(guò)統(tǒng)計(jì)研究區(qū)氣象站實(shí)際降雨過(guò)程的數(shù)據(jù),或者通過(guò)模擬研究區(qū)典型降雨過(guò)程,確定不同降雨過(guò)程中不同時(shí)間分辨率下降雨侵蝕力的差異程度,并以此對(duì)每一場(chǎng)降雨的降雨侵蝕力進(jìn)行修正,進(jìn)而校正指定時(shí)間段的降雨侵蝕力。為驗(yàn)證該方法的有效性,以薊運(yùn)河上游地區(qū)為研究區(qū),在Matlab 2010b中對(duì)當(dāng)?shù)氐湫徒涤赀^(guò)程進(jìn)行模擬,計(jì)算當(dāng)?shù)夭煌瑫r(shí)間分辨率下降雨侵蝕力差異程度,并利用典型降雨年份2008年的TRMM 3B42數(shù)據(jù)(0.25°×0.25°空間分辨率及3 h時(shí)間分辨率)通過(guò)ENVI-IDL語(yǔ)言程序進(jìn)行校正計(jì)算。結(jié)果表明,經(jīng)過(guò)校正后,研究區(qū)年降雨侵蝕力分布約為900~2 900 MJ·mm/(hm2·h)。該結(jié)果與這一地區(qū)多年平均降雨侵蝕力的量級(jí)及分布一致,說(shuō)明該校正方法是一種切實(shí)可行的、能有效提高TRMM數(shù)據(jù)應(yīng)用精度的方法;但是受TRMM數(shù)據(jù)空間分辨率的限制,該方法只適用于研究大中尺度區(qū)域降雨侵蝕力。

        降雨侵蝕力;TRMM;系統(tǒng)誤差;校正方法;薊運(yùn)河流域

        近些年來(lái)全球氣候變化導(dǎo)致極端天氣頻繁發(fā)生,我國(guó)很多地區(qū)旱澇災(zāi)害事件劇增,這加重了我國(guó)土壤侵蝕嚴(yán)重程度[1]。降雨是引起土壤侵蝕的重要因子,由降雨引起土壤侵蝕的潛在能力,即為降雨侵蝕力[2]5-8[3]。準(zhǔn)確計(jì)算和評(píng)估降雨侵蝕力對(duì)各種土壤侵蝕模型的設(shè)計(jì)、模型精度的評(píng)價(jià)及準(zhǔn)確定量預(yù)報(bào)土壤侵蝕工作都有重要的意義[4-6]。

        降雨侵蝕力是用降雨的總動(dòng)能和降雨強(qiáng)度估算的多年平均指標(biāo)[2]50-51[6]。指定時(shí)間長(zhǎng)度內(nèi),每場(chǎng)降雨的總動(dòng)能E和對(duì)應(yīng)場(chǎng)次最大降雨強(qiáng)度I的乘積EI的累計(jì)值作為該段時(shí)間內(nèi)的降雨侵蝕力[3,7]。通常降雨資料從氣象站點(diǎn)獲取[8-9],但氣象站點(diǎn)數(shù)量有限且分布不均,這使獲取空間面域上的降雨侵蝕力分布情況變得極為困難[5,10];此外,及時(shí)準(zhǔn)確地評(píng)價(jià)和預(yù)報(bào)土壤侵蝕,需要及時(shí)準(zhǔn)確的降雨數(shù)據(jù)去估算降雨侵蝕力,這些也是在現(xiàn)階段技術(shù)條件下氣象站難以滿足的[9,11]。1997年11月美國(guó)和日本合作開(kāi)發(fā)的全球第一顆搭載降水雷達(dá)的熱帶降雨觀測(cè)衛(wèi)星(Tropical Rainfall Measuring Mission,TRMM)發(fā)射成功,使獲取準(zhǔn)全球空間面域上的降雨分布即時(shí)數(shù)據(jù)成為現(xiàn)實(shí),彌補(bǔ)了氣象站觀測(cè)數(shù)據(jù)的不足[11-13]。

        學(xué)者們對(duì)TRMM數(shù)據(jù)在測(cè)算降雨強(qiáng)度和降雨量時(shí)的精度問(wèn)題爭(zhēng)議較大。C.M.Mannaerts等[14]認(rèn)為T(mén)RMM數(shù)據(jù)的精度可用于大尺度空間降雨的研究;而T.Y.Pham[15]認(rèn)為T(mén)RMM數(shù)據(jù)由于空間分辨率太低,并不能得到較為準(zhǔn)確的降雨分布。研究區(qū)尺度越小,用TRMM數(shù)據(jù)計(jì)算結(jié)果的誤差越大。TRMM數(shù)據(jù)不宜用于小尺度區(qū)域的研究,用其研究中等尺度空間范圍內(nèi)降雨侵蝕力也很少,通常只用其研究大尺度空間范圍內(nèi)的降雨量和降雨侵蝕力的分布[9-11,13,16]?,F(xiàn)有TRMM數(shù)據(jù)計(jì)算降雨侵蝕力的方法分為2種,第1種是根據(jù)TRMM數(shù)據(jù)估算的降雨量結(jié)合研究區(qū)的經(jīng)驗(yàn)公式間接計(jì)算,第2種是根據(jù)TRMM數(shù)據(jù)利用RUSLE的公式直接計(jì)算。第1種方法因計(jì)算過(guò)程簡(jiǎn)單,被廣泛使用但是其可靠性差,誤差不能估計(jì)。如Fan Jianrong等[9]選用經(jīng)驗(yàn)公式利用TRMM計(jì)算的降雨量數(shù)據(jù)去估算降雨侵蝕力分布。第2種方法計(jì)算過(guò)程較為復(fù)雜,但是能夠準(zhǔn)確卻反映研究區(qū)降雨侵蝕力變化趨勢(shì);為解決TRMM數(shù)據(jù)精度帶來(lái)的誤差,Zhu Qiang等[11]利用TRMM數(shù)據(jù)3 h統(tǒng)計(jì)的平均降雨強(qiáng)度和最大降雨強(qiáng)度代替30 min統(tǒng)計(jì)的平均降雨強(qiáng)度和最大降雨強(qiáng)度。A.Vrieling等[10,13,16]也簡(jiǎn)單認(rèn)為2種時(shí)間分辨率下的最大降雨強(qiáng)度相等;但是這些處理方式與實(shí)際情況不符(詳見(jiàn)第2節(jié)),從而降低了研究的可信度。

        筆者基于以往利用TRMM數(shù)據(jù)計(jì)算降雨侵蝕力方法的研究,找到了一種有效提高TRMM數(shù)據(jù)在中等尺度上應(yīng)用精度的方法,這既擴(kuò)大TRMM數(shù)據(jù)的應(yīng)用范圍,且彌補(bǔ)部分地區(qū)地面氣象觀測(cè)的空白。此外,筆者在用TRMM數(shù)據(jù)計(jì)算降雨侵蝕力的方法基礎(chǔ)上分析了該方法的系統(tǒng)誤差,提出了一種可行的校正方法,并以薊運(yùn)河流域上游地區(qū)為研究區(qū),利用典型降雨年份2008年的TRMM 3B42數(shù)據(jù),對(duì)該校正方法進(jìn)行了驗(yàn)證。

        1 降雨侵蝕力計(jì)算方法

        降雨強(qiáng)度是一個(gè)實(shí)時(shí)變化的物理量,時(shí)間分辨率越高,監(jiān)測(cè)到的降雨強(qiáng)度變化越快,變化的幅度越大,統(tǒng)計(jì)到的降雨強(qiáng)度極值越大;反之,時(shí)間分辨率越低,變化幅度越小,統(tǒng)計(jì)到的極值越小。降雨侵蝕力是在平降雨強(qiáng)度和最大降雨強(qiáng)度的基礎(chǔ)上計(jì)算得到。與降雨強(qiáng)度一樣,降雨侵蝕力也受到時(shí)間分辨率的影響。為此RUSLE手冊(cè)建議,用每30 min統(tǒng)計(jì)的平均降雨強(qiáng)度和最大降雨強(qiáng)度計(jì)算降雨侵蝕力[3,7],在降雨侵蝕力的研究計(jì)算中,只針對(duì)具有侵蝕性的降雨進(jìn)行研究和計(jì)算(通常將降雨量>12 mm的次降雨認(rèn)定為侵蝕性降雨[17])。

        1.1 地面氣象站降雨資料計(jì)算降雨侵蝕力

        地面氣象站降雨監(jiān)測(cè)數(shù)據(jù)的時(shí)間分辨率通常為30 min。RUSLE手冊(cè)建議,利用這種數(shù)據(jù)計(jì)算降雨侵蝕力的公式[3,18]如下:

        式中:R為年平均降雨侵蝕力,MJ·mm/(hm2·h);n為計(jì)算R值得總年數(shù);j為參與計(jì)算的年份序數(shù);m為第j年內(nèi)的暴雨次數(shù);k為參與計(jì)算的次暴雨序數(shù);E為次降雨總動(dòng)能,MJ/hm2;I30為30 min時(shí)間分辨率資料統(tǒng)計(jì)的次降雨最大30 min降雨強(qiáng)度,mm/h;i30為30 min時(shí)間分辨率資料統(tǒng)計(jì)的次降雨平均降雨強(qiáng)度,mm/h;ΔV為對(duì)應(yīng)時(shí)間段內(nèi)的降雨量,mm。

        1.2 TRMM數(shù)據(jù)計(jì)算降雨侵蝕力

        TRMM3B42數(shù)據(jù)集由TRMM科學(xué)資料信息系統(tǒng)和TRMM辦公室聯(lián)合制作(在http:∥trmm.gsfc. nasa.gov/或 http:∥old-cdc.cma.gov.cn/免費(fèi)獲取),是基于紅外亮溫資料采用3B-42算法得到的準(zhǔn)全球的降雨估量數(shù)據(jù),其圖像的一個(gè)像素點(diǎn)代表0.25°×0.25°經(jīng)緯度范圍,像素值為每3 h平均降雨強(qiáng)度資料,單位mm/h。即TRMM數(shù)據(jù)得到降雨強(qiáng)度資料的時(shí)間分辨率為180 min。

        基于TRMM的降雨侵蝕力計(jì)算方法是根據(jù)RUSLE手冊(cè)建議采用的降雨侵蝕力計(jì)算公式推導(dǎo),采用TRMM 3B42中的3 h降雨量資料,篩除次降雨量小于12 mm的非侵蝕性降雨,采用K.G.Renard等[3]和Zhu Qiang等[11]的計(jì)算方法計(jì)算次降雨侵蝕力、月降雨侵蝕力和年降雨侵蝕力,計(jì)算公式如下:

        式中:e為次降雨動(dòng)能,MJ/(mm·hm2);Rs為次降雨侵蝕力,MJ·mm/(hm2·h);i180為180 min時(shí)間分辨率下的次降雨平均降雨強(qiáng)度,mm/h;I180為T(mén)RMM資料中的次降雨最大180 min降雨強(qiáng)度,mm/h;p為次降雨持續(xù)3 h段數(shù),il為該次降雨第l個(gè)3 h段的平均降雨強(qiáng)度。計(jì)算時(shí),將該次降雨中il的最大值作為次降雨最大180 min降雨強(qiáng)度,將il的平均值作為180 min時(shí)間分辨率下的次降雨平均降雨強(qiáng)度。

        對(duì)于月降雨侵蝕力,則是將一個(gè)月內(nèi)的次降雨侵蝕力疊加得到:

        式中:Rm為月降雨侵蝕力,MJ·mm/(hm2·h);q為該月次降雨的次數(shù)。

        對(duì)于年降雨侵蝕力,則是將一年內(nèi)的次降雨侵蝕力疊加而得到:

        式中:Ry為年降雨侵蝕力,MJ·mm/(hm2·h);r為該年次降雨的次數(shù)。

        1.3 TRMM數(shù)據(jù)計(jì)算降雨侵蝕力分布

        編寫(xiě)ENVI-IDL語(yǔ)言程序一次性讀入研究范圍內(nèi)的TRMM 3B42數(shù)據(jù),利用式(3)~(6)計(jì)算出每個(gè)網(wǎng)格(0.25°×0.25°)站點(diǎn)的降雨量以及降雨侵蝕力。這樣得到的降雨侵蝕力可以表述大尺度上的降雨分布,但結(jié)果粗糙。通常利用地統(tǒng)計(jì)學(xué)方法對(duì)格網(wǎng)上的降雨侵蝕力進(jìn)行克里格插值,得到更精細(xì)平滑的降雨侵蝕力分布。

        為了準(zhǔn)確構(gòu)建插值模型,利用ArcGIS中探索性空間數(shù)據(jù)分析(ESDA)工具對(duì)數(shù)據(jù)進(jìn)行檢測(cè)及篩選。其中:使用自然間斷點(diǎn)的分類符號(hào)方法對(duì)TRMM數(shù)據(jù)初步計(jì)算結(jié)果進(jìn)行繪制和檢查;利用趨勢(shì)面分析任務(wù),將經(jīng)緯度坐標(biāo)(x,y)與降雨侵蝕力的值(z)構(gòu)成的點(diǎn)對(duì)投影在南北和東西2個(gè)平面上,得到的投影點(diǎn)分別進(jìn)行趨勢(shì)線回歸擬合,得到南北和東西2個(gè)方向平面上的降雨侵蝕力變化趨勢(shì)線,即得到趨勢(shì)面分析圖。根據(jù)趨勢(shì)分析的結(jié)果決定克里金插值是否進(jìn)行趨勢(shì)剔除處理。

        協(xié)同克里格方法可以給出有限區(qū)域內(nèi)區(qū)域變量的最佳線性無(wú)偏估計(jì)量,在遙感圖像處理中也得到了越來(lái)越廣泛的應(yīng)用[19]。因降雨受地形影響較大[20],所以結(jié)合研究區(qū)域的數(shù)字高程模型(DEM)可以對(duì)格網(wǎng)上降雨侵蝕力進(jìn)行協(xié)同克里格插值。在對(duì)TRMM數(shù)據(jù)進(jìn)行克里格插值之前,需要將上一節(jié)中計(jì)算出的每個(gè)網(wǎng)格的降雨侵蝕力值作為該網(wǎng)格中心點(diǎn)位置的降雨侵蝕力,然后進(jìn)行插值[10-11,13,16]。

        2 降雨侵蝕力計(jì)算方法的校正

        2.1 TRMM數(shù)據(jù)獲取降雨侵蝕力時(shí)的誤差來(lái)源

        由遙感影像得到的TRMM數(shù)據(jù)源計(jì)算的降雨侵蝕力結(jié)果必然存在一定程度的誤差。該誤差分為2類:隨機(jī)性誤差和系統(tǒng)性誤差。隨機(jī)性誤差由TRMM衛(wèi)星遙感獲取數(shù)據(jù)時(shí)降雨強(qiáng)度大小及取樣時(shí)間的隨機(jī)性引起,也由不同尺度凸顯的地學(xué)要素(地形、地貌及地理位置等)影響的隨機(jī)性引起。隨機(jī)性誤差不可避免,尤其地學(xué)要素對(duì)降雨侵蝕力總體趨勢(shì)的影響會(huì)隨著尺度的縮小被逐步放大。系統(tǒng)性誤差由TRMM衛(wèi)星搭載的測(cè)雨雷達(dá)的精度水平、監(jiān)測(cè)時(shí)間間隔,及TRMM數(shù)據(jù)計(jì)算降雨侵蝕力的計(jì)算方法差異導(dǎo)致。由設(shè)備精度帶來(lái)的誤差隨著技術(shù)水平的改進(jìn)可以提高,由計(jì)算方法帶來(lái)的系統(tǒng)誤差可以估計(jì)并消除。

        利用不同時(shí)間分辨率氣象站降雨資料得到的降雨侵蝕力大小不一樣[8,18]。研究TRMM數(shù)據(jù)計(jì)算的降雨侵蝕力精度及誤差值,需要選定評(píng)判的標(biāo)準(zhǔn),以往研究中通常采用 30 min降雨資料作為基準(zhǔn)[3,7]。因?yàn)闅庀笳緮?shù)據(jù)只是針對(duì)特定空間點(diǎn)的降雨數(shù)據(jù),所以對(duì)于面域的降雨分布來(lái)說(shuō),氣象站定點(diǎn)數(shù)據(jù)作為T(mén)RMM數(shù)據(jù)計(jì)算結(jié)果的基準(zhǔn)并不合適;故不能對(duì)由隨機(jī)性造成的誤差進(jìn)行準(zhǔn)確的估計(jì)。筆者只分析由計(jì)算方法引起的系統(tǒng)誤差。

        2.2 不同分辨率數(shù)據(jù)計(jì)算降雨侵蝕力方法的差異

        受TRMM 3B42數(shù)據(jù)時(shí)間分辨率3 h的限制,不能獲取式(1)和式(2)中所需的30 min時(shí)間分辨率統(tǒng)計(jì)的次降雨最大降雨強(qiáng)度I30和平均降雨強(qiáng)度i30;因此,在利用TRMM數(shù)據(jù)計(jì)算降雨侵蝕力時(shí),采用180 min分辨率下統(tǒng)計(jì)的次降雨平均降雨強(qiáng)度i180和最大降雨強(qiáng)度I180來(lái)替代較難獲得的30 min時(shí)間分辨率統(tǒng)計(jì)的次降雨最大降雨強(qiáng)度I30和平均降雨強(qiáng)度i30,利用式(3)和式(4)來(lái)計(jì)算次降雨侵蝕力。

        不同的計(jì)算方法應(yīng)用于同一地區(qū)會(huì)存在不同(甚至相差數(shù)倍)的降雨侵蝕力值[6]。降雨量、歷時(shí)及降雨時(shí)程(雨量隨歷時(shí)分配過(guò)程)的不同對(duì)一場(chǎng)降雨的降雨強(qiáng)度有顯著的影響。對(duì)于一場(chǎng)歷時(shí)恰好為30 min的侵蝕性降雨來(lái)說(shuō),如果在降雨開(kāi)始時(shí)刻計(jì)時(shí)統(tǒng)計(jì),最大降雨強(qiáng)度I30是I180的6倍,而平均降雨強(qiáng)度i30也是i180的6倍。由此可見(jiàn),氣象站數(shù)據(jù)和TRMM數(shù)據(jù)計(jì)算降雨侵蝕力所得結(jié)果的差異,其實(shí)質(zhì)上是不同時(shí)間分辨率下(最大和平均)降雨強(qiáng)度的差異[18]。在這種情況下,確定30 min與180 min時(shí)間分辨率下降雨強(qiáng)度的關(guān)系是校正在TRMM數(shù)據(jù)提取降雨侵蝕力中系統(tǒng)誤差的前提。

        確定同一場(chǎng)降雨在不同時(shí)間分辨率下降雨強(qiáng)度之間關(guān)系的途徑有2種:1)通過(guò)研究區(qū)氣象站大量的實(shí)際降雨過(guò)程數(shù)據(jù),統(tǒng)計(jì)推算出研究區(qū)不同類型降雨的降雨歷時(shí)、降雨量和降雨強(qiáng)度之間的關(guān)系式;2)通過(guò)研究區(qū)域的暴雨強(qiáng)度公式并結(jié)合設(shè)計(jì)雨型,進(jìn)行研究區(qū)不同典型降雨過(guò)程的時(shí)程分配,以這些時(shí)程分配統(tǒng)計(jì)出不同時(shí)間分辨率下不同類型降雨的降雨歷時(shí)、降雨量和降雨強(qiáng)度之間的關(guān)系式[21][22]3-5[23]14-17。 顯然,在缺少精細(xì)降雨資料的情況下,第二種途徑更為可行。根據(jù)得到不同時(shí)間分辨率下降強(qiáng)度的關(guān)系,計(jì)算出對(duì)于同一場(chǎng)降雨不同時(shí)間分辨率下降雨侵蝕力之間的關(guān)系:

        式中:μ為降雨侵蝕力差異系數(shù);R30和R180分別為30 min和180 min時(shí)間分辨率降雨資料計(jì)算的次降雨侵蝕力,MJ/(mm·hm2)。

        2.3 TRMM提取降雨侵蝕力校正方法

        基于以上分析,對(duì)于1.2節(jié)中使用TRMM數(shù)據(jù)計(jì)算大中尺度區(qū)域的降雨侵蝕力的方法,提出了以下校正方法:

        1)建立不同歷時(shí)和不同降雨量的典型降雨過(guò)程降雨侵蝕力之間的關(guān)系表(μ值表);

        2)對(duì)于實(shí)際降雨,根據(jù)該場(chǎng)降雨的歷時(shí)和降雨量,在步驟1得到的差異值表中找出對(duì)應(yīng)的降雨侵蝕力差異值μ;

        3)利用該差異值μ乘以TRMM數(shù)據(jù)計(jì)算的初步結(jié)果,得到每一場(chǎng)次降雨降雨侵蝕力的修正值R′s,

        4)通過(guò)累加每日、每月、每年內(nèi)的對(duì)應(yīng)降雨場(chǎng)次的修正值,得到修正后日、月、年的降雨侵蝕力。

        得到修正后的降雨侵蝕力之后,根據(jù)1.3節(jié)所述利用協(xié)同克里格插值計(jì)算降雨侵蝕力分布。

        3 應(yīng)用案例

        3.1 研究區(qū)概況

        薊運(yùn)河流域上游地區(qū)位于河海流域北部,東臨灤河,西臨潮白河,地跨北京、天津及河北3省市, E116°50′~118°15′,N39°45′~40°45′之間(圖1)。薊運(yùn)河主要支流有泃河、州河,泃、州兩支于天津市寶坻區(qū)九王莊匯流后稱薊運(yùn)河。九王莊以上地區(qū)為薊運(yùn)河的主要匯水區(qū)域,即為薊運(yùn)河上游地區(qū)。該區(qū)域銜接燕山山脈與華北平原,其北部地區(qū)為起伏的中低山、丘陵地形,南部地區(qū)以平原區(qū)為主。該地區(qū)屬于溫帶半濕潤(rùn)大陸性季風(fēng)氣候,季風(fēng)氣候顯著。

        3.2 數(shù)據(jù)來(lái)源

        根據(jù)中國(guó)氣象數(shù)據(jù)共享服務(wù)網(wǎng)年降水量數(shù)據(jù),薊運(yùn)河上游地區(qū)所在的北京市、天津市和承德市2003—2013年的11年間的降水量相關(guān)性很強(qiáng),且各市2008年的降雨量相當(dāng)于其11年平均值(圖2);因此選取2008年為薊運(yùn)河上游地區(qū)的典型降雨年份,并以此計(jì)算薊運(yùn)河上游地區(qū)降雨侵蝕力分布特征和分析計(jì)算結(jié)果的精度及誤差。

        圖1 研究區(qū)域及區(qū)內(nèi)主要?dú)庀笳军c(diǎn)分布示意圖Fig.1 Location of the study area with distribution of meteorological stations

        圖2 天津市、北京市和承德市近10年間年降雨量分布Fig.2 Annual precipitation distribution of Tianjin,Beijing and Chengde in the recent decade

        研究區(qū)降雨數(shù)據(jù)采用2008年的TRMM 3B42數(shù)據(jù)(數(shù)據(jù)獲取見(jiàn)1.2節(jié)),在ArcGIS中裁剪范圍在E116°26′57″~118°8′10″,N39°16′35″~40°47′33″的數(shù)據(jù)進(jìn)行研究。

        研究區(qū)數(shù)字高程(DEM)數(shù)據(jù)來(lái)自日本METI與美國(guó)NASA于2011年聯(lián)合開(kāi)發(fā)的ASTER GDEM數(shù)據(jù)(可在http:∥gdem.ersdac.jspacesystems.or.jp/ index.jsp注冊(cè)下載)。選用DEM數(shù)據(jù)格式為Geo-TIFF,圖像的空間分辨率為30 m。

        3.3 薊運(yùn)河上游地區(qū)μ值計(jì)算

        因缺少研究區(qū)以往精細(xì)降雨資料,采用設(shè)計(jì)典型降雨過(guò)程去計(jì)算薊運(yùn)河上游地區(qū)μ值表。典型降雨過(guò)程的設(shè)計(jì)需要確定雨型和暴雨強(qiáng)度公式。

        國(guó)內(nèi)外常用的雨型有芝加哥雨型、模式雨型、Huff雨型、三角形雨型、SCS雨型等等。各種雨型之間存在差異,目前還沒(méi)有一種公認(rèn)的雨型作為設(shè)計(jì)的依據(jù)[23]46-57。本研究區(qū)內(nèi),缺乏唐山及承德雨型的相關(guān)研究,而天津市區(qū)用Haff雨型的研究比較多,但天津市區(qū)屬于海洋性氣候,降雨過(guò)程與內(nèi)陸有別,故不采用天津市的雨型公式。北京市降雨雨型多用芝加哥雨型或在其基礎(chǔ)上改進(jìn)的雨型,并結(jié)合北京市暴雨公式進(jìn)行相關(guān)設(shè)計(jì)研究[21][22]20-21,筆者選用北京市用過(guò)的芝加哥雨型。

        通過(guò)薊運(yùn)河上游地區(qū)暴雨強(qiáng)度公式結(jié)合常用的芝加哥雨型計(jì)算μ值表,但缺少薊運(yùn)河上游地區(qū)暴雨強(qiáng)度公式的研究和記載;而天津市相關(guān)暴雨強(qiáng)度公式是根據(jù)天津市區(qū)氣象站資料求得[24]。天津市區(qū)臨海,其降雨量受海洋影響較大,而薊運(yùn)河上游地區(qū)距其較遠(yuǎn),且在內(nèi)陸,地形地貌差距較大,選用天津市暴雨強(qiáng)度公式不合適;因此采用緯度相近、地理范圍相互重疊、地形地貌相似的北京市第二暴雨分區(qū)的暴雨強(qiáng)度公式作為替代。筆者所用北京市第二區(qū)設(shè)計(jì)暴雨強(qiáng)度公式[21][22]3為

        式中:i為平均降雨強(qiáng)度,mm/min;td為降雨歷時(shí), min;P為設(shè)計(jì)降雨重現(xiàn)期,a;在選定總降雨量下,迭代試算出P值。

        降雨歷程的分配使用芝加哥雨型,雨峰系數(shù)r取0.25,則峰前時(shí)間序列i(tb)和峰后時(shí)間序列i (tc)分別用如下公式[21,25]計(jì)算:

        為了簡(jiǎn)化運(yùn)算,在校正降雨侵蝕力時(shí),不考慮每場(chǎng)降雨的降雨量大小不同,只考慮侵蝕性降雨的降雨歷時(shí)變化。選取降雨量為25 mm(國(guó)家氣象局規(guī)定的中雨和大雨的分界線),歷時(shí)分別為0.5、1、2、3、4、5、6 h的7場(chǎng)降雨,利用式(8)~(11)在Matlab 2010b中合成降雨過(guò)程,模擬的步長(zhǎng)為1 min。

        對(duì)上述擬合所得不同降雨過(guò)程曲線進(jìn)行統(tǒng)計(jì)計(jì)算,以降雨開(kāi)始時(shí)刻為統(tǒng)計(jì)起點(diǎn),得到30 min與180 min分辨率下最大、平均降雨強(qiáng)度及對(duì)應(yīng)的降雨侵蝕力(表1)。

        表1 不同降雨過(guò)程特征物理量及μ值Tab.1 Characteristics of physical quantities including μ value in different rainfall processes

        平均降雨強(qiáng)度隨時(shí)間分辨率的差異隨降雨歷時(shí)而不同(表1),但以此計(jì)算降雨動(dòng)能E相差只有16.5%左右,這相比于最大降雨強(qiáng)度隨時(shí)間分辨率的差異較小,此結(jié)果與前人研究結(jié)果[26]一致。最大降雨強(qiáng)度隨時(shí)間分辨率的差異均在135%以上,即I30≈2.35I180,這表明最大降雨強(qiáng)度受時(shí)間分辨率的影響程度與降雨歷時(shí)的關(guān)系不緊密。同時(shí)證明Zhu Qiang等[11]和 A.Vrieling等[10,16]研究時(shí)簡(jiǎn)單認(rèn)為I30=I180是不對(duì)的。

        TRMM數(shù)據(jù)只能反映出每3 h內(nèi)是否發(fā)生降雨,降雨持續(xù)了幾個(gè)3 h段,并能得出實(shí)際降雨的歷時(shí)。根據(jù)李建等[27]研究,20世紀(jì)90年代以后北京市降雨以6 h以下降雨為主,且貢獻(xiàn)最大的為2~4 h降雨。并且筆者通過(guò)TRMM數(shù)據(jù)發(fā)現(xiàn)>6 h的降雨非常少,只占到6.7%左右,沒(méi)有連續(xù)>9 h的降雨。另外,根據(jù)表1,隨著降雨歷時(shí)的增加,μ值變化幅度較小。因此我們采用了2、4、6 h的降雨代表<3 h、介于3~6 h的降雨、>6 h的降雨。即校正時(shí),歷時(shí)≤3 h的降雨,其降雨侵蝕力乘以μ2h(2.85);降雨介于3~6 h的降雨,其降雨侵蝕力乘以 μ4h(2.75);歷時(shí)>6 h的降雨,其降雨侵蝕力乘以μ6h(2.62)。本文中所模擬的2、4、6 h的降雨時(shí)程過(guò)程如圖3所示。

        3.4 結(jié)果與分析

        表2 薊運(yùn)河上游地區(qū)7個(gè)縣市的主氣象站點(diǎn)所在網(wǎng)格的降雨侵蝕力校正前后對(duì)比Tab.2 Rainfall erosivity on grids with main meteorological sites on seven counties of Jiyun River upriver before and after the calibration

        3.4.1 降雨侵蝕力校正前后數(shù)據(jù)對(duì)比 表2說(shuō)明,對(duì)每場(chǎng)降雨的降雨侵蝕力經(jīng)過(guò)校正之后,不同網(wǎng)格處累計(jì)得到2008年降雨侵蝕力的值均提高了2倍以上,使研究區(qū)2008年的降雨侵蝕力的量級(jí)與多年平均降雨侵蝕力[28-29]一致;故這樣的校正顯著提高了TRMM數(shù)據(jù)使用時(shí)的精度。

        圖3 設(shè)計(jì)25 mm降雨量不同歷時(shí)降雨過(guò)程Fig.3 Rainfall process in different designed durations at precipitation of 25 mm

        3.4.2 降雨侵蝕力分布趨勢(shì) 圖4表明較大的侵蝕力出現(xiàn)在該區(qū)域南部地區(qū),中部地區(qū)降雨侵蝕力居中,北部地區(qū)較小,即該地區(qū)降雨侵蝕力自東南到西北逐漸減弱的趨勢(shì)。圖5中可以直觀地看出,東西、南北方向變化趨勢(shì)近似。自東向西,數(shù)據(jù)點(diǎn)分布較為分散,東西方向擬合趨勢(shì)線有微弱的先升再降的趨勢(shì);而南北方向擬合趨勢(shì)線數(shù)據(jù)點(diǎn)自南向北有微弱的減弱趨勢(shì)。以上2方向上的趨勢(shì)與薊運(yùn)河上游地區(qū)地形相符合,北部山勢(shì)影響山前平原降雨比較多。

        3.4.3 降雨侵蝕力插值結(jié)果及分析 協(xié)同克里格方法是根據(jù)半變異函數(shù)分析所提供的變量空間自相關(guān)程度的信息來(lái)進(jìn)行插值,插值準(zhǔn)確程度與半變異函數(shù)擬合的好壞直接相關(guān)[19]。表3示出,指數(shù)模型在平均值、標(biāo)準(zhǔn)平均值、標(biāo)準(zhǔn)均方根3個(gè)指標(biāo)上是誤差最低。在均方根指標(biāo)上,指數(shù)模型的結(jié)果與該指標(biāo)最低值差異不大。只有在平均標(biāo)準(zhǔn)誤差上,指數(shù)模型的結(jié)果高于其他模型。綜合其他2個(gè)模型的誤差結(jié)果來(lái)看,指數(shù)模型模擬效果最好,故本研究采用指數(shù)模型進(jìn)行后面的插值參數(shù)改進(jìn)。

        圖4 2008年降雨侵蝕力自然間斷點(diǎn)分類符號(hào)圖Fig.4 Classification symbols of natural discontinuities of rainfall erosivity in 2008

        塊金值與基臺(tái)值的比值C0/(C0+C)稱為塊金效應(yīng),它表示隨機(jī)部分引起的空間變異占系統(tǒng)總變異的比例[30]。如果該比值較低,表明結(jié)構(gòu)性因素引起的空間變異起主要作用;如果該值較高,說(shuō)明隨機(jī)部分引起的空間變異起主要作用。通常認(rèn)為:塊金效應(yīng)<25%,空間相關(guān)性強(qiáng);25%~75%之間,空間相關(guān)性中等;>75%,空間相關(guān)性弱。本文采用指數(shù)模型計(jì)算的塊金值為0,基臺(tái)值為72 241,塊金效應(yīng)<25%,表明TRMM像素值具有強(qiáng)烈的空間相關(guān)性,其結(jié)構(gòu)性因素引起的空間變異起主要作用。

        圖5 2008年降雨侵蝕力趨勢(shì)面分析Fig.5 Trend surface analysis of rainfall erosivity in 2008

        依據(jù)擬合的半變異函數(shù)模型協(xié)同DEM進(jìn)行克里格插值。研究區(qū)年降雨侵蝕力分布約為900~2 900 MJ·mm/(hm2·h)。圖6表明,該區(qū)域降雨侵蝕力東南到西北逐漸減弱,南部地區(qū)大于北部地區(qū),該變化趨勢(shì)和章文波等[28]、劉斌濤等[29]研究結(jié)果一致。

        表3 不同半變異模型的誤差評(píng)價(jià)結(jié)果Tab.3 Error analysis of different fitting models of semi-variance function

        圖6 薊運(yùn)河上游各縣域降雨侵蝕力分布圖Fig.6 Rainfall erosivity distribution on the upriver counties of Jiyun River

        4 結(jié)論與討論

        以典型降雨年份2008年薊運(yùn)河上游地區(qū)的降雨侵蝕力計(jì)算為研究案例,結(jié)果顯示校正后的降雨侵蝕力接近了多年平均降雨侵蝕力的量級(jí)及分布。這說(shuō)明通過(guò)確定不同降雨過(guò)程中不同時(shí)間分辨率下降雨侵蝕力差異的μ值表,并以此對(duì)每一場(chǎng)降雨進(jìn)行修正,進(jìn)而減少TRMM數(shù)據(jù)計(jì)算降雨侵蝕力是系統(tǒng)誤差的方法,是一種切實(shí)可行的、有效提高TRMM應(yīng)用結(jié)果精度的方法。

        受到TRMM數(shù)據(jù)時(shí)空分辨率的影響,該方法只能用于大中尺度區(qū)域降雨侵蝕力的研究,不能用于小尺度區(qū)域的研究。該方法并不限制如何選擇降雨雨型和下墊面條件的影響。降雨雨型不同選擇,會(huì)使影響校正結(jié)果的精度,但是這種影響遠(yuǎn)比未校正之前的誤差小。下墊面條件引起的誤差屬于隨機(jī)性誤差,在現(xiàn)有研究的基礎(chǔ)上,暫時(shí)無(wú)法估計(jì)下墊面條件的影響規(guī)律和引起誤差大小。

        確定不同降水過(guò)程對(duì)應(yīng)的μ值是應(yīng)用該方法的前提和基礎(chǔ)。這一步驟需要以實(shí)際降雨過(guò)程數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析為基礎(chǔ),或者以更多典型降水過(guò)程的模擬數(shù)據(jù)為基礎(chǔ)。由于各地降雨量不同,降雨過(guò)程隨機(jī)變化,應(yīng)用該方法是需要因地制宜,確定出適合當(dāng)?shù)氐摩讨当?。?dāng)然在誤差允許的條件下,也可以制訂各氣候區(qū)或降雨類型區(qū)統(tǒng)一的μ值表,方便以后隨時(shí)調(diào)取使用。從本文案例的校正過(guò)程來(lái)看,影響校正結(jié)果效果的因素有雨型確定的方法、雨型種類、不同降雨歷時(shí)的降雨所占的比例、不同降雨歷時(shí)及不同降雨量的降雨對(duì)應(yīng)的雨型種類等。追求更高精度的校正效果時(shí),須對(duì)這些問(wèn)題進(jìn)行深入的研究,這也是今后進(jìn)一步研究和改進(jìn)完善該校正方法途徑。

        [1] 王占禮.中國(guó)土壤侵蝕影響因素及其危害分析[J].農(nóng)業(yè)工程學(xué)報(bào),2000,16(4):32- 36

        [2] Wischmeier W H,Smith D D.Predicting rainfall erosion losses:a guide to conservation planning[J].Agriculture Handbook 537,1978:5- 8,50- 51

        [3] Renard K G,Foster G R,Weesies G A,et al.Predicting soil erosion by water:aguide to conservation planning with the revised universal soil loss equation(RUSLE) [J].Agriculture Handbook 703,1997:22- 38

        [4] 葉芝苗,劉寶元,章大波,等.北京市降雨侵蝕力及其空間分布[J].中國(guó)水土保持科學(xué),2003,1(1): 16- 20

        [5] 章文波,付金生.不同類型雨量資料估算降雨侵蝕力[J].資源科學(xué),2003,25(1):35- 41

        [6] 卜兆宏,董勤瑞,周伏建,等.降雨侵蝕力因子新算法的初步研究[J].土壤學(xué)報(bào),1992,29(4):408- 417

        [7] Yu Baofu,Rosewell C J.A robust estimator of the R-factor for the Universal Soil Loss Equation[J].Transactions of the ASAE,1996,39(2):559- 561

        [8] 章文波,謝云,劉寶元.用雨量和雨強(qiáng)計(jì)算次降雨侵蝕力[J].地理研究,2002,21(3):384- 390

        [9] Fan Jianrong,Chen Yang,Yan Dong,et al.Characteristics of rainfall erosivity based on tropical rainfall measuring mission data in Tibet,China[J].Journal of Mountain Science,2013,10(6):1008- 1017

        [10]Vrieling A,Hoedjes J C B,van der Velde M.Towards large-scale monitoring of soil erosion in Africa:Accounting for the dynamics of rainfall erosivity[J].Global and Planetary Change,2014,115:33- 43

        [11]Zhu Qiang,Chen Xiuwan,Fan Qixiang,et al.A new procedure to estimate the rainfall erosivity factor based on Tropical Rainfall Measuring Mission(TRMM)data[J]. Science China Technological Sciences,2011,54(9): 2437- 2445

        [12]Kummerow C,Barnes W,Kozu T,et al.The tropical rainfall measuring mission(TRMM)sensor package[J]. Journal of atmospheric and oceanic technology,1998,15 (3):809- 817

        [13]Vrieling A,Sterk G,de Jong S M.Satellite-based estimation of rainfall erosivity for Africa[J].Journal of hydrology,2010,395(3):235- 241

        [14]Mannaerts C M,Saavedra C P.Regional scale erosion modeling and monitoring using remotely sensed data: some spatial data scale issues[C]∥Proceedings of the international symposium on 25 years of assessment of erosion.Ghent,Belgium:Ghent University,2003:433- 440

        [15]Pham T Y.Soil erosion risk modeling within upland landscapes in Vietnam using remotely sensed data and the RUSLE model[D].Nova Scotia:Dalhousie University, 2008:79- 84

        [16]Vrieling A,de Jong S M,Sterk G,et al.Timing of erosion and satellite data:a multi-resolution approach to soil erosion risk mapping[J].International Journal of Applied Earth Observation and Geoinformation,2008,10, 267- 281

        [17]Xie Yun,Liu Baoyuan,Nearing M A.Practical thresholds for separating erosive and non-erosive storms[J]. Transactions of the ASAE,2002,45(6):1843- 1847

        [18]Yin Shuiqing,Xie Yun,Nearing M A,et al.Estimation of rainfall erosivity using 5-to 60-minute fixed-interval rainfall data from China[J].Catena,2007,70(3): 306- 312

        [19]馮益明,雷相東,陸元昌.應(yīng)用空間統(tǒng)計(jì)學(xué)理論解譯遙感影像信息“缺失”區(qū)[J].遙感學(xué)報(bào),2004,8(4): 317- 322

        [20]傅抱璞.地形和海拔高度對(duì)降水的影響[J].地理學(xué)報(bào),1992,59(4):302- 314

        [21]張大偉,趙冬泉,陳吉寧,等.芝加哥降雨過(guò)程線模型在排水系統(tǒng)模擬中的應(yīng)用[J].給水排水,2008,34 (增刊2):354- 357

        [22]北京市城市規(guī)劃設(shè)計(jì)研究院.城市雨水系統(tǒng)規(guī)劃設(shè)計(jì)暴雨徑流計(jì)算標(biāo)準(zhǔn):DB11/T 969—2013[S].北京:中國(guó)計(jì)劃出版社,2013:3- 5,20- 21

        [23]范澤華.天津市降雨趨勢(shì)分析及設(shè)計(jì)暴雨研究[D].天津:天津大學(xué),2012:14- 17,46- 57

        [24]黃津輝,向文艷,戶超,等.天津市設(shè)計(jì)暴雨方法比較及公式修正[J].天津大學(xué)學(xué)報(bào),2013,46(4): 354- 360

        [25]Keifer C J,Chu H H.Synthetic storm pattern for drainage design[J].Journal of the hydraulics division,1957,83 (4):1- 25

        [26]殷水清,謝云,王春剛.用小時(shí)降雨資料估算降雨侵蝕力的方法[J].地理研究,2007,26(3):541- 547

        [27]李建,宇如聰,王建捷.北京市夏季降水的日變化特征[J].科學(xué)通報(bào),2008,53(7):829- 832

        [28]章文波,謝云,劉寶元.中國(guó)降雨侵蝕力空間變化特征[J].山地學(xué)報(bào),2003,21(1):33- 40

        [29]劉斌濤,陶和平,宋春風(fēng),等.1960—2009年中國(guó)降雨侵蝕力的時(shí)空變化趨勢(shì)[J].地理研究,2013,32(2): 245- 256

        [30]Schabenberger O,Pierce F J.Contemporary statistical models for the plant and soil sciences[M].Boca Ration: CRC press,2002:563- 599

        (責(zé)任編輯:宋如華)

        Calibration of rainfall erosivity calculation based on TRMM data: A case study of the upriver basin of Jiyun River,North China

        Jin Qiutong1,2,Shi Mingchang1,2,3,Zhang Jutao1,2,Wang Shan4,Hu Ying3
        (1.School of Soil and Water Conservation,Beijing Forestry University,100083,Beijing,China;2.Key Lab.of Soil&Water Conservation and Desertification Combating,Beijing Forestry University,Ministry of Education,100083,Beijing,China;3.Beijing Datum Science and Technology Development Co.,Ltd.,100084,Beijing,China;4.AutoNavi Software Co.,Ltd.,102200,Beijing,China)

        The low temporal resolution of TRMM data could cause systematic errors in the rainfall erosivity calculation within the specified time period.For this reason,we propose a method of calibration in order to reduce the errors in this paper.Firstly,a relationship table of rainfall erosivity at different events and duration with typical rainfall process is established with historical rainfall data or numerical simulation;then,the corresponding difference is calculated from the above table based on the precipitation and duration of each rainfall in the study area;afterwards,the adjusted value of rainfall erosivity in each rainfall is calculated by multiplying the figure found in the pervious step by the original value of rainfall erosivity correspondingly;finally,the calibrated value of daily,monthly or yearly rainfall erosivity isobtained by adding daily, monthly oryearly adjusted value ofrainfallerosivity correspondingly.To verify the effectiveness of this calibration method,the typical rainfall process in the study area,the upriver regions of Jiyun River,was simulated within Matlab 2010b to obtain therelationship table,and then the rainfall erosivity in the study area was calculated using TRMM 3B42 data in 2008 with the spatial resolution of 0.25°×0.25°and the time resolution of 3 h.The calculation was conducted in the ENVI-IDL Language program and calibrated by the relationship table.The results showed that:the rainfall erosivity in the upriver of Jiyun River ranged from 900 to 2 900 MJ·mm/(hm2· h)with the highest values in central area,followed by southern and northern parts,which was consistent with the local average annual values and distribution of rainfall erosivity.The results confirm that this method can effectively improve the precision of TRMM data's application.Due to the limitation of spatial resolution of TRMM data,this calibration method is only applied in the rainfall erosivity calculation at large and middle scales.In order to improve the precision of this calibration method,the typical rainfall process,rain type design and distribution of rainfall duration should be detailed.

        rainfall erosivity;TRMM;systematic error;calibration method;Jiyun River

        S157.1

        A

        1672-3007(2015)04-0094-09

        2014- 11- 15

        2015- 05- 03

        項(xiàng)目名稱:林業(yè)公益性行業(yè)科研專項(xiàng)經(jīng)費(fèi)項(xiàng)目課題“基于生態(tài)安全的水土保持措施空間配置技術(shù)”(201404209)

        靳秋桐(1990—),男,碩士研究生。主要研究方向:3S技術(shù)集成與系統(tǒng)開(kāi)發(fā)。E-mail:jinqiutong@bjfu.edu.cn

        ?通信作者簡(jiǎn)介:史明昌(1969—),男,博士,教授。主要研究方向:地理信息科學(xué)。E-mail:shimc@bjfu.edu.cn

        猜你喜歡
        雨型降雨量分辨率
        降雨量與面積的關(guān)系
        概化的累計(jì)暴雨量百分?jǐn)?shù)法在太湖區(qū)域設(shè)計(jì)暴雨雨型研究的應(yīng)用
        江蘇水利(2020年9期)2020-10-09 02:53:52
        天津市設(shè)計(jì)暴雨雨型的演變
        EM算法的參數(shù)分辨率
        深圳市流域暴雨雨型及變化趨勢(shì)分析
        原生VS最大那些混淆視聽(tīng)的“分辨率”概念
        上海市設(shè)計(jì)雨型對(duì)雨水管網(wǎng)模擬的影響研究
        山西建筑(2017年21期)2017-09-03 10:29:20
        基于深度特征學(xué)習(xí)的圖像超分辨率重建
        一種改進(jìn)的基于邊緣加強(qiáng)超分辨率算法
        洞庭湖區(qū)降雨特性分析

        亚洲sm另类一区二区三区| 国产成人美女AV| 中文字幕精品亚洲二区| 国产内射一级一片内射高清视频1| 夜夜夜夜曰天天天天拍国产| 成 人 免费 黄 色 视频| 亚洲欧美在线观看一区二区| 精品日本免费观看一区二区三区| 一区二区三区人妻少妇| 欧美内射深喉中文字幕| 国产在线不卡AV观看| 一本久久伊人热热精品中文| 青青草大香蕉视频在线观看| 国产成人涩涩涩视频在线观看| 亚洲国产精品国自产电影| 91九色国产在线观看| 91成人自拍国语对白| 亚洲成av人片一区二区| 国产熟女亚洲精品麻豆| 国产白浆大屁股精品视频拍| 一边捏奶头一边高潮视频| 久久人人爽人人爽人人片av麻烦| 久久久久久久久国内精品影视| 大陆少妇一区二区三区| 亚洲欧美色一区二区三区| 色偷偷88888欧美精品久久久 | 国产熟女自拍av网站| 国产乱妇无乱码大黄aa片| 中文乱码人妻系列一区二区| 亚洲国产精品第一区二区三区 | 国产成人精品成人a在线观看 | 久久国产自偷自免费一区100| 狼人av在线免费观看| 国产一级内射视频在线观看| 欧美 变态 另类 人妖| 亚洲中文欧美日韩在线| 中文字幕日韩高清乱码| 特黄做受又硬又粗又大视频小说| 久久这里只精品国产99热| 黄色大片国产精品久久| 国内精品视频一区二区三区八戒|