?

基于哨兵1號的臺風風場反演方法研究

2022-06-23 05:31于鵬鐘小菁耿旭樸
關鍵詞:臺風

于鵬 鐘小菁 耿旭樸

關鍵詞:哨兵1號;臺風;合成孔徑雷達;交叉極化

0引言

我國位于亞洲東部、太平洋西岸,東部和南部大陸海岸線超1.8萬km,特殊的地理位置決定了我國是世界上少數幾個受臺風影響最嚴重的國家之一.1949—2019年期間,共有491個臺風登陸我國,即平均每年約有7個臺風在我國登陸[1],給人民的生命安全和生產生活帶來了嚴重的危害,由其引起的極高風速、風暴潮以及風暴巨浪是我國海岸帶最主要的致災因素[2].臺風的風速是決定臺風潛在破壞性的重要指標,及時并準確地掌握臺風風場信息是深入研究風暴潮洪水、風暴浪和海岸侵蝕等臺風災害的基礎,對災害預警、風險評估和災后重建等具有重要意義.

相比于傳統的地面觀測方法和系統,衛星遙感具備在更大空間尺度上對海表面進行近實時觀測的優勢,這使得基于衛星遙感的手段非常適用于對像臺風這樣大尺度現象的監測.光學和紅外傳感器是探測臺風頂部云層性質的重要儀器,但由于云層的持續存在和阻擋,它們很難被用來監測臺風系統在海表面風速的大小.星載合成孔徑雷達(SyntheticApertureRadar,SAR)是一種非常獨特的工具,它具備全天時和全天候的對地觀測能力.且與其他微波傳感器和數值天氣預報估算結果(分辨率通常約10km以上)相比,SAR具有更高的空間分辨率(風速產品可達到百米量級),可以更好地用于估計臺風的風速和研究臺風的內部結構,在預報臺風路徑方面也具備極大的潛力.

海表面的風速可以通過同極化(Vertical-Vertical或Horizontal-HorizontalPolarization,VV或HH極化)SAR圖像以及經驗地球物理模型函數(GeophysicalModelFunctions,GMFs)[3-6]來獲取,但在高風速條件下,它們會受到信號飽和的影響,使反演結果小于實際風速[7-10].此外,如果沒有預先的外部風向輸入,基于同極化雷達數據的GMFs無法提供唯一的風速結果[11].近年來發射的星載C波段衛星傳感器(Radarsat-2、Sentinel-1以及Gaofen-3)提供了獲取來自地表交叉極化(Vertical-Horizontal或Horizontal-VerticalPolarization,VH或HV極化)雷達信號的新手段.與同極化信號相比,交叉極化雷達信號包含更多的海表面非布拉格散射的信息,這些信息與高風速條件下存在的波浪破碎現象有關[12-13].基于海表后向散射的物理模型[14]和GMFs模型[15]等的研究表明,C波段的交叉極化信號至少在觀測60m/s的海表風速時也不會出現信號飽和的情況.因此,利用交叉極化SAR信號對臺風風速進行反演的結果更為可靠,而且該反演方法并不像同極化信號一樣需要預輸入風向信息,具備在惡劣天氣條件下對海表風速進行獨立監測的能力.近些年來的研究提出了多種基于交叉極化SAR數據的高風速反演模型[16-22],但缺少對各類模型之間的比較研究.

本研究使用哨兵1號(Sentinel-1,S-1)SAR遙感觀測數據來研究臺風的風場信息.在交叉極化雷達數據的基礎上,以2020年的臺風“海高斯”及“莫拉菲”為例,首先應用噪音向量和去噪方法對超寬幅模式影像內部的加性和乘性噪音進行去除.然后利用近年來開發的幾種交叉極化GMFs對兩個臺風的風場進行估算,并與歐洲中期天氣預報中心(EuropeanCentreforMedium-rangeWeatherForecasts,ECMWF)模型風場信息和基于同極化數據的哨兵1號風場產品做比較.分析和研究適合哨兵1號的臺風風場反演方法.

1數據及預處理

1.1 哨兵1號SAR遙感數據

本研究所使用的SAR遙感影像由歐洲航天局(EuropeanSpaceAgency,ESA)哥白尼計劃的哨兵1號任務提供,該數據是全世界首個向用戶免費開放的SAR遙感影像數據,可在赤道地區提供重訪周期為6d的海表面風場數據.哨兵1號任務由兩顆衛星(Sentinel-1A和Sentinel-1B)組成,配備有中心頻率為5.405GHz的C波段SAR傳感器,可在4種特有的成像模式下工作,即:條帶(stripmap,SM)、波(wave,WV)、干涉寬幅(InterferometricWide,IW)和超寬幅(ExtraWide,EW)模式.

本文所使用影像的具體信息如表1所示,臺風“海高斯”于北京時間(下同)2020年8月18日18:25左右被影像1和2捕捉到,兩景影像的成像模式為IW模式,空間分辨率為20m×22m.臺風“莫拉菲”于2020年10月27日18:39左右被影像3和4捕捉到,兩景影像的成像模式為EW模式,空間分辨率為93m×87m.

1.2 臺風及最大風速數據

2007號臺風“海高斯”(Higos)是2020年登陸我國的最強臺風之一,也是該年登陸我國廣東的最強臺風.該臺風于2020年8月16日在菲律賓呂宋島以東的西北太平洋洋面生成,8月17日以熱帶低壓的強度通過了呂宋海峽后到達中國南海,8月18日十幾個小時內完成熱帶低壓、熱帶風暴、強熱帶風暴和臺風的四級連跳(圖1(a)),并最終于2020年8月19日6時左右在廣東省珠海市金灣區沿海登陸.“海高斯”所攜帶的極具破壞力的強風以及極端的降水給我國東南沿海地區帶來了嚴重的洪澇災害,導致了大面積的農田、房屋、公路及水利工程等受損.

2018號臺風“莫拉菲”(Molave)是2020年通過我國南海的最強臺風之一.該臺風于2020年10月21日在西太平洋的法斯島南部海域生成,10月25日完成熱帶風暴、強熱帶風暴及臺風的三級連跳(圖1(b)),此后總計在菲律賓登陸了5次,于10月26日進入中國南海,并最終于10月28日中午在越南廣義省附近沿海登陸.該臺風總計致使71人死亡,46人失蹤,并造成超過6.6億美元的財產損失,鑒于該臺風所導致的嚴重自然災害,臺風委員會已將其從臺風命名表中進行除名,以后將不再使用“莫拉菲”這一臺風名稱.

用于驗證模型結果中臺風最大風速的數據來源于國家氣象中心研制的全球再分析數據產品(ChinaMeteorologicalAdministration’sglobalatmosphericRe-Analysis,CRA)[23](http://data.cma.cn),該數據集基于集合-變分混合同化技術,融合了探空、地面等中國特有的常規觀測資料和風云衛星資料,可提供逐小時的臺風中心位置、強度和移動方向等信息,能夠為臺風強度及其他氣象預報提供有效檢驗.與時間分辨率為6h的臺風最佳路徑(BestTrack)數據相比,采用CRA最大風速信息能夠有效降低因時間插值導致的誤差.圖1(a)顯示了國家氣象中心提供的“海高斯”從2020年8月16—20日的最大風速信息,圖1(b)為“莫拉菲”從2020年10月25—29日的最大風速信息.

1.3 ECMWF模型數據

大氣模型的風速預報結果是研究臺風強度和參數的重要數據來源之一.本文將利用ECMWF提供的模式風場輸出信息與本研究中其他方法獲取的臺風最大風速和結構進行比較.ECMWF預報模型提供了空間分辨率為0.125°×0.125°、時間分辨率為3h的海表面10m的風場數據.

1.4 SAR數據預處理

S-1衛星SAR數據可由歐洲航天局提供的SNAP(SentinelApplicationPlatform)軟件進行預處理,處理過程包括輻射定標、幾何校正、陸地掩膜和均值濾波等.由于“海高斯”登陸地點在珠江口附近,該海域的航運較為繁忙,進行風速反演時,結果容易受到船只等信號的干擾.考慮到在交叉極化SAR影像的多視處理中,隨著視數的增大,尤其是當視數達到6以上時,船只與背景海面雷達信號的對比度會顯著降低[24],因此“海高斯”的SAR遙感圖像的像素間距首先被重采樣為10m×10m,后被統一采用多視及平均處理為600m×600m以降低船只等噪聲對風速反演的影響.而對于“莫拉菲”的SAR影像,由于采用的是EW成像模式,扇狀效應及熱噪音的存在嚴重影響了雙極化SAR數據的成像結果[25],尤其限制了交叉極化數據在第1子條帶上風場反演的應用[26].因此,針對“莫拉菲”的兩景交叉極化影像,本研究應用了S-1中的噪音向量以及Sun等[27]提出的噪聲去除方法對EW模式的交叉極化SAR影像進行處理,該算法除對已有的加性噪音去除方法進行了改進外,還提出了一種乘性噪音的去除方法以降低EW模式影像子條帶邊界的噪音.圖2展示了應用該算法前后的SAR影像,通過兩圖的對比可以發現,該算法能夠在各子條帶之間的邊界區域獲得較好的處理結果,在第1子條帶中明顯存在的殘差噪音也得到了有效去除,經過處理后的整幅SAR影像不再出現明顯的不連續特征.

圖3(a)顯示了2020年8月18日IW模式的哨兵1號VH極化方式捕捉到的經預處理和陸地掩膜后的“海高斯”的SAR影像,表1所示的兩幅時間間隔約為1min的SAR圖像被合成在一起以提供“海高斯”更完整的結構.圖3(b)展示的則是2020年10月27日EW模式的哨兵1號VH極化方式捕捉到的經處理后的臺風“莫拉菲”的SAR影像,同樣由兩幅影像合成以更完整展示“莫拉菲”的臺風結構,由圖上的經緯度范圍及表1信息可知,EW模式影像的幅寬更大,能夠展示比“海高斯”SAR影像更大的范圍.

2基于交叉極化SAR數據的風速反演方法

與同極化雷達信號相比,交叉極化信號對風向不敏感,反演風場不會因為風向變化而出現多個風速解的情況[28].目前,針對C波段的交叉極化SAR遙感數據,研究人員已經建立了多個交叉極化海洋模型以將海表面的粗糙度與較高風速聯系起來,這些模型并不需要外部風向的輸入[16-17,19,21-22],其中主要有基于Radarsat-2衛星數據提出的C波段交叉極化海洋模型(C-bandCrossPolarizationOceanmodel,C-2PO)[16-17,19]和C波段交叉極化耦合參數海洋模型(C-bandCross-PolarizationCoupled-ParametersOceanmodel,C-3PO)[20],以及基于我國發射的Gaofen-3衛星提出的全極化條帶交叉極化模型(Quad-polarizationStripmapCross-polarizationmodel,QPS-CP)[21-22].

2.1 C-2PO和QPS-CP模型

本文采用的C-2PO模型有4種,分別以C2011[16]、C2012[17]、C2014Z[18]及C2014V[19]進行表示,其公式可以表示為

上式中:sVH表示VH極化的雷達后向散射系數(單位為dB);U10代表海表面以上10m標準高度處的風速(單位為m/s);b1和b2是用來擬合風速估算值的經驗系數.參數b1在C2011、C2012、C2014Z及C2014V模型中分別選取為0.592、0.58、0.332和0.218;參數b2分別為35.6、35.652、30.143和29.07.

QPS-CP基于Gaofen-3衛星的全極化條帶模型SAR產品提出,其表達形式與C-2PO的公式相同,都是關于風速和交叉極化SAR信號的線性公式,但參數大小不同.本文采用了C2019[21]和C2021[22]兩種模型,參數b1分別為0.6683和0.4273;參數b2分別為37.3732和34.3875.

2.2 C-3PO模型

雷達的入射角是風速反演中的另外一個影響因素,但在上述基于交叉極化的模型中沒有體現.因此,本研究應用了C-3PO以考慮雷達入射角對估算高風速信息的影響.C-3PO模型的形式為

上式中:q表示雷達入射角(單位為°).

本研究基于S-1的交叉極化SAR數據對以上7種交叉極化模型結果、ECMWF模型結果及S-1海洋二級產品提供的風場進行比較,分析出最適合進行臺風風場反演的方法.

3基于SAR數據估算的臺風風場結果

3.1 模型估算的最大風速結果

直接利用同極化數據的GMFs進行風速反演已有諸多應用,但與該系列模型相比,基于同極化數據和變分方法的模型能夠將更多的誤差因素考慮在內,可以在中低風速條件下獲得更高精度的風場信息[29].如S-1提供的風速產品就是利用Portabella等[30]提出的基于二維變分的方法,將C-bandModelFunctionIfremer2(CMOD-IFR2)同極化模型以及ECMWF模型的風矢量輸出結果進行調整以獲取最佳風場結果.因此,盡管CMOD-IFR2模型能在小于20m/s的風速條件下得到很好的結果[31-32],但歐空局仍然采用了基于二維變分方法的風場數據作為S-1的風場產品.表2列出了通過ECMWF氣象預報模型、S-1二級海洋產品以及采用不同交叉極化模型獲得的臺風“海高斯”及“莫拉菲”的海表面最大風速結果.其中,S1代表哨兵1號提供的由ECMWF和CMOD-IFR2兩者同化得來的風速產品,從其值(25.5m/s和42.5m/s)遠大于ECMWF風速大?。?8.3m/s和35.8m/s)可知,其最大風速會受到CMOD模型的影響,并且與最佳路徑數據相比,其具有更小的平均誤差,更適合進行臺風的風速反演.需要注意的是,國家氣象中心CRA產品提供的最大風速為2min最大平均風速值,而SAR影像反演獲得的最大風速為影像獲取時刻的最大風速值.

對比交叉極化模型的結果可以發現,C-3PO模型提供了最小的平均誤差(3.6m/s),與“海高斯”的最大風速相比,其誤差只有0.8m/s,在估算最大風速較高的莫拉菲(45m/s)時,其誤差(6.3m/s)較大.與預期結果相似,ECMWF預報模型結果較差,尤其是在估算“海高斯”的風速時僅能提供最大風速為18.3m/s的8級風力,其遠遠低于CRA提供的接近30m/s的11級風力結果.

3.2 基于交叉極化模型的“海高斯”臺風風場結果

圖4(a)展示了2020年8月18日18時25分ECMWF模型預報的“海高斯”臺風的風場結果,圖4(b)為S-1二級海洋產品提供的基于ECMWF和CMOD-IFR2模型的同化結果.圖4(c)—(i)分別展示了基于交叉極化模型C2011、C2012、C2014Z、C2014V、C-3PO、C2019及C2021的風場結果.本研究中,基于SAR數據的風場反演結果(圖4(b)—(i))的空間分辨率為600m,遠高于ECMWF預報模型(圖4(a))的結果(約10km),因此,高分辨率的SAR影像能夠提供其他監測和模型手段無法獲得的更精準的臺風內部結構.交叉極化模型(圖4(c)—(i))與S-1海洋產品(圖4(b))相比,可以提供更為完整的臺風結構,尤其是在臺風眼的西側區域,海洋產品只提供了較低的風速結果.與其他模型結果相比,C-3PO模型(圖4(g))除可以獲得更準確的最大風速外,還能夠很好地還原臺風內部的高風速信息.如表2,交叉極化模型中除了C-3PO外,C2021模型的平均誤差(4.1m/s)最小,但其對臺風眼區域的風速估算誤差達到10m/s以上,無法準確獲取該區域的風速信息.

3.3 基于交叉極化模型的“莫拉菲”臺風風場結果

圖5(a)—(i)分別展示了臺風“莫拉菲”的ECMWF預報模型、S-1二級海洋產品以及不同交叉極化GMFs反演的臺風風場結果.盡管與臺風“海高斯”的結果(圖4(a))相比,ECMWF的模型風場結果更為對稱(圖5(a)),但其臺風眼區域位于北緯14°左右,與SAR影像呈現出的臺風眼區域(圖5(g)中間藍色區域,北緯13.5°左右)存在很大差距.數據同化產品能夠顯示出更準確的臺風眼位置(圖5(b)),但與交叉極化模型的結果相比(圖5(c)—(i)),其在臺風眼附近無法呈現出強臺風一般所具有的較為對稱的風圈結構.在所有交叉極化GMFs估算的最大風速結果中,C2014V模型(圖5(f))的平均誤差相對較小,但與圖4(f)類似,其對臺風的最大風速具有明顯的高估現象.

3.4 基于交叉極化模型的中低風速反演結果驗證

通過表1、圖4和圖5的對比可以發現,與其他模型結果相比,C-3PO模型無論是在臺風眼內外區域的風速估算上,還是在臺風結構的對稱性上,都能呈現出較好的結果.但C-3PO及C2014V等交叉極化模型在臺風外圍低風速區域的風速估算上存在問題,如在圖4(f)和圖4(g)中存在大范圍的藍色低風速區域,與海表面實際的氣象情況不符.考慮到交叉極化模型對低風速估算存在一定的問題,而S-1的海洋產品在小于20m/s的中低風速條件下能夠獲得較為理想的結果,本文利用該產品對幾種交叉極化模型的中低風速結果進行了評估,如表3所示.盡管與其他模型結果相比,C-3PO反演的中低風速結果與風場產品的相關系數(0.41)較高,但也僅為中低程度相關,并且存在較大的均方根誤差.C2011、C2012及C2019能夠獲得小于3m/s的均方根誤差,但與風場產品僅為弱相關或低度相關.這可能與低風速區域對應的雷達后向散射系數過于接近S-1數據的最大噪聲等效散射系數(MaximumNoiseEquivalentSigmaZero,NESZ,大約–23dB)值有關.總體而言,對于哨兵1號SAR數據,交叉極化模型在中低風速條件下的表現較差.

3.5 臺風風場的融合結果

基于同極化SAR數據的S-1海洋產品能夠提供精度較高的中低風速結果,并且已經應用于常規的海表面風場監測,但在臺風的高風速反演上則存在一定的問題.相比之下,C-3PO等交叉極化模型能夠獲得較為合理的臺風風場的對稱結構以及較為可靠的最大風速結果,但又無法獲取準確的中低風速信息.因此,無法采用單一波段(同極化或交叉極化)的哨兵1號SAR數據對整體的臺風風場進行準確估算.綜合以上因素,本研究對S-1二級海洋產品和交叉極化模型結果進行了融合,其中背景的中低風速由S-1海洋產品提供,當C-3PO反演的風速結果大于20m/s且大于海洋產品提供的風速時,融合數據的風速由C-3PO模型提供.

圖6展示了臺風“海高斯”及“莫拉菲”的融合風場結果.結果表明,該融合結果綜合了同極化和交叉極化SAR數據的優勢,在很好地再現臺風內部的較強風圈的同時,還能在臺風外圍得到較為合理的風速結果.對于臺風“海高斯”,融合后的臺風外圍的中低風速結果(圖6(a))比C-3PO模型的結果(圖4(g))更為可信,圖4(g)中臺風東側的藍色低風速值的區域結果得到了較大改善.而對于臺風“莫拉菲”,S-1海洋產品對C-3PO模型的反演結果進行了有效補充,其在東側區域的高風速估算上有了一定改善(圖6(b)和圖5(g)),獲得了42.5m/s的最大海表面風速,與CRA提供的45m/s的最大風速信息更為接近.總之,該融合方法無論在空間分辨率還是在極高風速的估算方面,都能夠提供比ECMWF模型(圖4(a)和5(a))更好的臺風風場結果.

4討論

SAR衛星通常采用掃描模式(ScanSAR),即通過天線在距離向的周期性轉向以獲得更大范圍的觀測空間,但在波束邊緣會獲得比波束中心更大的回波強度,導致在每個子條帶之間會出現較大的噪音.在S-1的IW和EW模式下SAR數據進一步采用了漸進式地形掃描(TerrainObservationsbyProgressiveScansSAR,TOPSAR)技術,即引入天線波束在方位向的轉變以提高SAR圖像質量,但噪音對于海表面這類后向散射信號較弱的地物的影響仍然十分顯著,尤其是對于擁有5個子條帶的EW模式交叉極化海表面影像(圖2(b)).圖7(a)展示了未經噪音去除而應用C-3PO交叉極化模型的臺風“莫拉菲”的風場結果,可以看出在各個子條帶之間形成了不合理的高風速條帶,該現象在第1和第2子條帶之間較為明顯,這種現象限制了哨兵1號EW模式交叉極化數據在臺風風場反演等領域的應用.本文以臺風“莫拉菲”的EW模式VH極化SAR影像為例,應用了噪聲去除方法[27]及歐空局提供的噪音向量以降低影像中的加性和乘性噪音,對比應用該方法前、后的風場結果(圖7),發現各子條帶之間的風場反演結果得到了極大改善,尤其是在第1和第2子條帶之間的高風速結果得到了有效優化,第1條帶內的東西向的條形結果也得到了有效去除.此外,背景噪音也會對反演結果造成影響,如目前采用的交叉極化模型多是基于Radarsat-2和Gaofen-3數據提出的,而Radarsat-2數據(如雙極化掃描模式產品,約–29dB)擁有優于Gaofen-3數據(如寬幅掃描模式產品,優于–22dB)和S-1數據(如IW模式產品,約–23dB)的背景噪音水平.由于海表同極化后向散射信號要強于交叉極化信號,背景噪音水平對同極化數據影響較小,但對交叉極化數據影響較大,尤其是在海表風速較小時,這也導致在將基于Radarsat-2提出的C2014Z、C2014V和C-3PO(圖4(e)—(g))模型直接應用于S-1交叉極化數據時會出現不合理的低風速結果,這也是本文反演小于20m/s的臺風風速時未采用交叉極化模型的原因.因此,未來可以應用本文采用的針對哨兵1號SAR數據提出的類似噪音去除方法,開發適合S-1交叉極化數據的反演模型,改善基于S-1衛星數據的風速反演效果.

如表2所示,交叉極化模型獲得的“海高斯”的最大風速與CRA的較為接近,但在最大風速達到45m/s的“莫拉菲”的風速估算上存在一定誤差,除模型導致的誤差外,這可能與臺風期間伴隨的強降水導致的C波段SAR回波信號的減弱有關[33-35].Alpers等[34]基于15年以上的實驗及觀測結果發現,當風速小于10m/s時,降雨會導致C波段海表雷達后向散射信號增強;當風速大于10m/s時,降雨會降低雷達信號強度.因此,對于臺風的不同區域,其反演精度可能也會存在差異.如在臺風眼區,風速和降雨量較小,通過同極化數據能夠獲得較好的風速反演結果;在云墻區,風速和降雨量最大,受到兩者的雙重影響,誤差可能最大;而對于螺旋雨帶區,其風速和降雨強度會介于臺風眼和云墻區之間,其反演精度也會介于兩者之間.今后,需要搜集更多的臺風期間的實測數據,增加對臺風不同區域風速估算誤差的分析研究,并在交叉極化GMFs中引入降雨的貢獻以提高對臺風風場的反演精度.

5結論與展望

本文以哨兵1號SAR數據的交叉極化影像為數據源,應用噪音向量和噪音去除方法對EW模式影像進行了預處理,影像中的加性和乘性噪音得到了有效去除,各子條帶之間的扇狀效應也得到了很大改善.利用2011—2021年10年來開發的C-2PO,C-3PO及QPS-CP等7種交叉極化模型對2020年的“海高斯”和“莫拉菲”兩個臺風的風場進行估算,結果表明,C-3PO模型在高風速反演方面能夠獲得較小的相對誤差,且與ECMWF模型以及S-1二級海洋風場產品相比,能夠在臺風眼附近獲得較為合理的對稱風場結構.與風場產品的比較發現,交叉極化模型在中低風速的估算上存在一定誤差,通過與已有風場產品的融合,可以兼顧同極化和交叉極化SAR數據的優勢,很好地再現臺風外部中低風速以及內部的海表面高風速結構.

交叉極化的SAR數據在高風速研究方面具有極高的價值,但目前的交叉極化模型在反演風速時仍存在一定的問題,如在本研究中需要與同極化產品進行融合以獲得較為可靠的中低風速信息,未來需要對基于S-1交叉極化數據的GMFs進行優化.另外,在軌SAR衛星數量較少導致很難實現對同一個臺風的多次重復觀測,這極大地限制了SAR數據在臺風監測與預報方面的應用.目前,國內已批準及正在規劃一系列的SAR衛星任務,尤其是低成本的小衛星星座計劃,如已于2020年12月發射首顆輕小型SAR衛星的廈門大學“海絲”衛星星座計劃[36],有望為臺風等海洋觀測和研究提供更短的重訪周期以及更高分辨率的SAR觀測數據,這也將加強和提高SAR衛星遙感在應對諸如臺風“海高斯”和“莫拉菲”這類極端災害事件的能力.

猜你喜歡
臺風
臺風過韓
英國皇家空軍臺風戰斗機
臺風,你到底是什么?
瓶子里的臺風怪
EF2000臺風戰斗機
臺風收割機
臺風寶寶
呼呼!臺風來了
臺風天外出小心
臺風地震連襲日本
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合