分裂窗算法完整版本_第1頁
分裂窗算法完整版本_第2頁
分裂窗算法完整版本_第3頁
分裂窗算法完整版本_第4頁
分裂窗算法完整版本_第5頁
已閱讀5頁,還剩47頁未讀 繼續(xù)免費(fèi)閱讀

下載本文檔

版權(quán)說明:本文檔由用戶提供并上傳,收益歸屬內(nèi)容提供方,若內(nèi)容存在侵權(quán),請進(jìn)行舉報(bào)或認(rèn)領(lǐng)

文檔簡介

第5講:分裂窗算法基于NOAA-AVHRR數(shù)據(jù)的分裂窗算法一、分裂窗算法:Introduction二、分裂窗算法的推導(dǎo)三、大氣透過率的確定四、地表比輻射率估計(jì)(略)五、SensitivityanalysisBi(Ti)=

i(

)[

iBi(Ts)+(1-

i)Ii

]+Ii

WhereBi(Ti) observedradiance

Bi(Ts) groundradiance

Ii

downwellatmosphericradiance

Ii

upwellatmosphericradiance

i(

)

atmospherictransmittance

i groundemissivity

i channelThermalradiancetransferequationNOAA-14:ch1–visible0.58-0.68

m

Channel2–nearinfrared(NIR)0.725-1.10

mChannel3–middleinfrared(MIR)3.55-3.93

m

Channel4–thermalinfrared10.5-11.3

m

Channel5–thermalinfrared11.5-12.5

m

NOAA-AVHRRNOAA--NationalOceanicandAtmosphericAdministration

AVHRR--advancedveryhighresolutionradiometerAVHRRisabroad-band,four(NOAA6,7,8)orfive(NOAA7,9,11,12,14-16)channelscanner,sensinginthevisible,near-infrared,andthermalinfraredportionsoftheelectromagneticspectrum.ThissensoriscarriedontheNationalOceanicandAtmosphericAdministration's(NOAA's)PolarOrbitingEnvironmentalSatellites(POES),beginningwithTIROS-Nin1978.Ithas2399kmwideswath.ThesatelliteorbitstheEarth14timeseachdayfrom833kmaboveitssurface.Theaverageinstantaneousfield-of-view(IFOV)is1.4milliradians,yieldingaLAC/HRPTGroundFieldofView(GFOV)ofapproximately1.1kmatthesatellitenadirfromthenominalorbitaltitudeof833km.NOAA-1412/30/94-PresentNOAA-1505/13/98-Present

NOAA-11一、分裂窗算法:Introduction分裂窗方法是到目前為止應(yīng)用最廣泛的地表面溫度反演方法,尤其是用來分析NOAA-AVHRR數(shù)據(jù)。從1981年NOAA-7衛(wèi)星發(fā)射以后,我們就能通過8-14

m范圍內(nèi)的兩個(gè)狹窄的熱紅外波段來對向行輻射進(jìn)行常規(guī)觀測。分裂窗方法的理念可追逆到AndingandKauth(1970)的研究,它是以11-12

m附近的兩個(gè)相鄰的熱紅外波段的不同大氣透過特征為基礎(chǔ)進(jìn)行地表面溫度的反演。由于它的理論是以同一大氣窗口分裂為兩個(gè)相鄰的熱紅外波段來觀測,所以稱為分裂窗方法。這種方法充分利用了AVHRR通道4和5的兩個(gè)相鄰波段的同步觀測的優(yōu)勢來進(jìn)行地表面溫度的反演。由于NOAA系列衛(wèi)星上都安裝有這兩個(gè)通道的遙感器,并且至少有兩個(gè)NOAA衛(wèi)星同時(shí)在軌道上運(yùn)行,所以每天至少能獲得4次覆蓋全球的遙感觀測數(shù)據(jù),從而使NOAA-AVHRR數(shù)據(jù)獲得非常廣泛的應(yīng)用。D.AndingandR.Kauth,1970.Estimationofseasurfacetemperaturefromspace.RemoteSens.Environ.,1,217-220.提出大氣影響C.Prabhakara,G.Dalu,andV.G.Kunde,1974.Estimationofseasurfacetemperaturefromremotesensinginthe11-13

mwindowregion.J.Geographys.Res.,79,5039-5044.,分析兩個(gè)波段的大氣影響差異對大氣校正的作用L.M.McMillin,1975.Estimationofseasurfacetemperaturesfromtwoinfraredwindowmeasurementswithdifferentabsorption.J.Geographys.Res.,80,5113-5117.提出利用兩波段大氣影響差異建立溫度反演算法P.Y.DeschampsandPhulpin,1980.Atmosphericcorrectionofinfraredmeasurementsofseasurfacetemperatureusingchannelsat3.7,11and12m.BoundaryLayerMeteorol.,18,131-143.進(jìn)一步完善海面溫度的分裂窗算法E.P.McClain,W.G.Pichel,C.C.Walton,Z.Ahmad,andJ.Sutton,1983.Multichannelimprovementstosatellite-derivedglobalseatemperatures.Adv.SpaceRes.,2,43-47.建立分裂窗算法在全球海面溫度監(jiān)測方法HistoryofSplitwindowmethods(a)Planck函數(shù)的線性化。這一假定要求用兩個(gè)非常鄰近的波段來進(jìn)行觀測。這就是說,一個(gè)波段的觀測值可能表達(dá)為另一個(gè)波段的線性函數(shù)。(b)所有相關(guān)溫度(亮度溫度、氣溫和表面溫度)的相同量級。在大多數(shù)大氣條件下,這一假定是正確的,從而讓我們能消除掉大氣的上行輻射貢獻(xiàn)。(c)地表比輻射率在兩個(gè)波段之間的差異性不大,存在線性關(guān)系。(d)大氣水分的吸收作用較弱,并且在兩個(gè)波段之間不相同。這一假定需求大氣的總水分含量不要太大。大氣水分含量的吸收作用與這兩個(gè)波段間的亮度溫度差異之間存在線性關(guān)系,只有在這些條件下,才能通過這種線性關(guān)系來進(jìn)行近似估計(jì)。在潮濕的熱帶大氣條件并有大視角的情況下,這個(gè)假定可以打破。

分裂窗算法的基本假設(shè):(a)水在熱紅外波段范圍內(nèi)的比輻射率已為人們所熟知,并且該比輻射率接近于1(~0.99)。此外,在8-14m范圍內(nèi),水的比輻射率隨光譜和視角的變化很小,因而可以忽略不計(jì)。(b)海表面的溫度與近水面的氣溫相差常常很小。(c)海表面溫度在空間上和時(shí)間上的變化都非常緩慢,從而可以從低空間分辨率的觀測值中獲得非常高精度的溫度反演結(jié)果。

SplitwindowmethodforLSToverseasurface因素影響著陸地表面溫度反演的精度

(a)陸地的比輻射率在空間上和時(shí)間上都可能有很大的變化。在8-14

m窗口范圍內(nèi),其絕對值可從0.91變化到0.98。地表面類型、土壤水分含量、以及植被物候階段和真實(shí)植被覆蓋結(jié)構(gòu)的不同都可能有顯著的影響。(b)陸地的比輻射率可能有顯然的光譜依賴性。盡管陸地的比輻射率在AVHRR通道4和5之間的差異通過小于0.01,但這些微小差異對LST反演也有顯著影響(Becker1987,Colletal.1994a,Colletal.1994b)。(c)有效比輻射率(effectiveemissivity),即遙感器從某個(gè)視角實(shí)際觀測物體時(shí)所獲得的比輻射率,取決于視角和地表面的各向異性(Gossmann1987,Choudhury1989,WanandDozier1989,CassellesandSobrino1989,LabedandStoll1991)。

(d)陸面溫度可能大大高于海面溫度。因此,針對海表面建立進(jìn)來的Planck方程的線性化,可能不適用于陸面溫度的變化范圍。此外,AVHRR遙感器在320K(47C)達(dá)到飽和,在觀測一些很熱環(huán)境的LST時(shí),這可能會引起一些問題。e)在一個(gè)AVHRR像元范圍內(nèi),LST可能會有很大的變化。這就會引起所觀測到的溫度值的確切含義有難以理解和實(shí)際意義問題,因?yàn)檫@種溫度觀測是在至少1.2km2的面積范圍內(nèi)進(jìn)行取值。類似地,“表面”在植被地區(qū)可能很難確切地定義。(f)LST有強(qiáng)烈的晝夜起伏變化。這就會使觀測時(shí)間成為一個(gè)重要問題。同時(shí)會引起這種溫度的實(shí)際應(yīng)用價(jià)值(g)在陸地上,(近地)氣溫可能與表面溫度有很大的差異。(h)由于地形起伏和高程變化,大氣的路徑長度可能有很大差別。根據(jù)某個(gè)高度來建立的大氣影響與亮度溫度之間的線性回歸可能不適合于其它高程的地面。Ts=T4+A(T4-T5)+BTs=A0+A1T4+A2T5

SplitwindowalgorithmsWhereTs Landsurfacetemperature

T4andT5 Brightnesstemperatures

AandB

theparameters A0,

A1,and

A2

theparametersThegeneralforms17LSTalgorithmsinliterature1.Kerretal.[1992)2.OttléandVidal-Madjar[1992]3.Price[1984]4.BeckerandLi[1990]5.PrataandPlatt[1991]6.Vidal[1991]7.Uliveriretal.[1996]8.Colletal.[1994]9.Sobrinoetal.[1991]10.Prata[1993]11.Fran?aandCracknell[1994]12.Qinetal.[2001]13.SobrinoandRaissouni[2000]14.Fran?oisandOttlé[1996]/Q15.Fran?oisandOttlé[1996]/W16.BeckerandLi[1995]17.Sobrinoetal.[1994]T4

brightnesstemperatureatchannel4 knownT5

brightnesstemperatureatchannel5 known

i(

) atmospherictransmittance estimated

i groundemissivity estimatedIi

,

Ii

=f(

i(

),Bi(Ta))Ta

effectivemeanatmospherictemperature二、Qinetal.(2001)分裂窗算法的推導(dǎo)B4(T4)=

4(

)[

4B4(Ts)+(1-

4)I4

]+I4

B5(T5)=

5(

)[

5B5(Ts)+(1-

5)I5

]+I5

Bi(Ti)=

i

i(

)Bi(Ts)+[1-

i(

)](1-

i)

i(

)Bi(Ta

) +[1-

i(

)]Bi(Ta)Ii

=[1-

i(

)]Bi(Ta)Ii

=(1-

i(

))Bi(Ta

)關(guān)于大氣平均作用溫度的分析

根據(jù)這一替代,NOAA-AVHRR遙感器所接收到的熱輻射強(qiáng)度可以近似地表示為Bi(Ti)=

i

i(

)Bi(Ts)+[1-

i(

)][1+(1-

i)

i(

)]Bi(Ta)把公式(12)應(yīng)用到AVHRR的通道4和5數(shù)據(jù)上,得到B4(T4)=

4

4(

)B4(Ts)+[1-

4(

)][1+(1-

4)

4(

)]B4(Ta)B5(T5)=

5

5(

)B5(Ts)+[1-

5(

)][1+(1-

5)

5(

)]B5(Ta)通道4和5的輻射傳導(dǎo)方程TaylorexpansiontoPlanckfunctionBi(Tj)=Bi(T)+(Tj-T)

Bi(T)/

T=(Li+Tj-T)

Bi(T)/

T,式中i表示通道4或5。Tj表示溫度;若j=4或5,則表示通道4或5的亮度溫度;若j=s,則表示我們所要求解的地表溫度;若j=a,則表示大氣溫度。參數(shù)Li是中間變量,由下式給出Li=Bi(T)/[

Bi(T)/

T],

因此,對于通道5,在亮度溫度為T5的情況下,Planck輻射函數(shù)的Taylor展開式可表示為B5(T5)=B5(T4)+(T5-T4)

B5(T4)/

T=(L5+T5-T4)

B5(T4)/

T類似地,對于通道4和5所對應(yīng)的各個(gè)溫度,我們有B4(T4)=(L4+T4-T4)

B4(T4)/

T=L4

B4(T4)/

TB4(Ts)=(L4+Ts-T4)

B4(T4)/

TB4(Ta)=(L4+Ta-T4)

B4(T4)/

TB5(Ts)=(L5+Ts-T4)

B5(T4)/

TB5(Ta)=(L5+Ta-T4)

B5(T4)/

TB4(T4)=

4

4(

)B4(Ts)+[1-

4(

)][1+(1-

4)

4(

)]B4(Ta)B5(T5)=

5

5(

)B5(Ts)+[1-

5(

)][1+(1-

5)

5(

)]B5(Ta)Taylor展開式的物理含義參數(shù)Li與溫度Ti之間也呈很明顯的接近于線性的關(guān)系

L4=-62.239281+0.430589T4,

R2=0.9991SEE=0.1893L5=-66.540666+0.465845T5, R2=0.9993SEE=0.1860Li=Ti/ni為了推導(dǎo)分裂窗算法,

我們可以用如下近似式來表示參數(shù)Li,Li=ai+biTi.

一般地講,該式中的系數(shù)ai和bi可以用公式(17)中的回歸系數(shù)來表示:對于通道4,a4=-62.23928,b4=0.43059對于通道5,a5=-66.54067,b5=0.46585分裂窗算法的推導(dǎo)

為了簡化起見,我們定義:Ci=

i

i(

),

(19)Di=[1-

i(

)][1+(1-

i)

i(

)]

(20)因此,輻射傳導(dǎo)方程可重寫成如下形式:B4(T4)=C4B4(Ts)+D4B4(Ta)

(21a)B5(T5)=C5B5(Ts)+D5B5(Ta)

(21b)把Planck函數(shù)的開展式代入,

我們得L4

B4(T4)/

T=C4(L4+Ts-T4)

B4(T4)/

T

+D4(L4+Ta-T4)

B4(T4)/

T(L5+T5-T4)

B5(T4)/

T=C5(L5+Ts-T4)

B5(T4)/

T

+D5(L5+Ta-T4)

B5(T4)/

T消去方程中的

B4(T4)/

T

B5(T4)/

T

得L4=C4(L4+Ts-T4)+D4(L4+Ta-T4)L5+T5-T4=C5(L5+Ts-T4)+D5(L5+Ta-T4)分裂窗算法的推導(dǎo)

從上述聯(lián)立方程中消去Ta,得D5L4-D4(L5+T5-T4)=D5C4(L4+Ts-T4)-D4C5(L5+Ts-T4) +D5D4(L4-T4)-D4D5(L5-T4)對方程求解

Ts,我們得到一個(gè)新的分裂窗算法,其公式與分裂窗算法的一般形式相同: Ts=T4+A(T4-T5)+B式中的參數(shù)分別定義如下: A=D4/E0,

B=E1L4-E2L5 E1=D5(1-C4-D4)/E0,

E2=D4(1-C5-D5)/E0 E0=D5C4-D4C5DerivationofsplitwindowalgorithmDerivationofsplitwindowalgorithm由于參數(shù)Li

是一個(gè)溫度函數(shù),我們還可進(jìn)一步推導(dǎo),以便得到一個(gè)更簡化的演算地表溫度的算法。代入?yún)?shù)Li

,我們得到B=E1(a4+b4T4)-E2(a5+b5T5)重新組織上式,得到B=E1b4T4-E2b5T5+E1a4-E2a5分裂窗算法公式中并合并同類項(xiàng),我們得到:Ts=A0+A1T4-A2T5式中系數(shù)

A0,A1,和

A2

分別定義為A0=E1a4-E2a5 A1=1+A+E1b4

A2=A+E2b5.Qinetal.[2001]:

Ts=A0+A1T4+A2T5

where A0=[66.54067D4(1-C5-D5) -62.23928D5(1-C4-D4)]/(D5C4-D4C5) A1=1+[D4+0.43059D5(1-C4-D4)]/(D5C4-D4C5)

A2=-D4[1+0.46585(1-C5-D5)]/(D5C4-D4C5)

Ci=

i

i(

)

Di=(1-

i(

))(1+(1-

i)

i(

))

Qin,Z.,G.Dall'Olmo,A.Karnieli,andP.Berliner.2001,DerivationofsplitwindowalgorithmanditssensitivityanalysisforretrievinglandsurfacetemperaturefromNOAA-advancedveryhighresolutionradiometerdata.J.Geophys.Res.,106(D19):22655-22670Splitwindowalgorithms三、大氣透過率的確定大氣透射率受許多大氣因素影響。氣壓、氣溫、氣溶膠含量、大氣水分含量、O3、CO2、CO、NH4等對熱輻射傳導(dǎo)均有不同程度的作用,從而使地表的熱輻射在大氣中的傳導(dǎo)產(chǎn)生衰減。因此,準(zhǔn)確的大氣透射率的求算比較復(fù)雜,需要較詳細(xì)的大氣剖面數(shù)據(jù)。一般來說,準(zhǔn)確地求算大氣透射率需要進(jìn)行大氣模擬。目前較普遍使用的大氣模擬程序有LOWTRAN、MODTRAN和6S等。然而,這種大氣模擬需要很詳細(xì)的大氣剖面數(shù)據(jù)。正如上面指出,在大多數(shù)情況下,實(shí)時(shí)大氣剖面數(shù)據(jù)并不具備,從而使大氣模擬法難以實(shí)施,雖然也可以使用這些程序提供的標(biāo)準(zhǔn)大氣來替代,但模擬過程較為復(fù)雜,很難在實(shí)際研究中普遍使用。在大氣各影響因素中,大氣水分含量的變化較快。研究表明,大氣透射率的變化主要取決于大氣水分含量的動態(tài)變化,其它因素因其動態(tài)變化不大而對大氣透射率的變化沒有顯著影響,并且大氣水分對熱輻射的吸收作用較大。因此,它就成為大氣透射率估計(jì)的主要考慮因素。根據(jù)這一特征,我們運(yùn)用大氣模擬程度LOWTRAN7來模擬大氣水分含量的變化與大氣透射率的變化之間的關(guān)系,然后建立相關(guān)方程,以便用來進(jìn)行大氣透射率的近似估計(jì)。在這一模擬中,我們考慮了大氣水分在0.4-6.4g/cm2區(qū)間范圍內(nèi)變動,這一區(qū)間基本上代表了天空晴朗條件下的大氣水分含量變化幅度。對于沙漠地區(qū)干燥的氣候,大氣水分含量一般較低,

只有0.5-1.5g/cm2左右,

而在較為濕潤的地區(qū),大氣水分含量一般較高,可達(dá)2.0-3.5g/cm2。大氣模擬還需要假定一個(gè)地面附近的氣溫所對應(yīng)的大氣剖面溫度分布。為此,我們考慮兩種情形:夏季和冬季。對于夏季,

我們假定地面附近的氣溫為35

C,

而冬季則為18

C。由于AVHRR的圖像寬度達(dá)2400km,遙感器視角是影響大氣透過率的重要因素。因此,在大氣模擬中我們考慮了10

的天頂視角。其它大氣狀態(tài)如CO2等使用中緯度夏季平均大氣的剖面數(shù)據(jù)。夏季大氣剖面模擬結(jié)果冬季大氣剖面模擬結(jié)果wgcm-2 EstimationEquations

R2 SEE 0.4-1.6

4(10)=0.979160

0.062918w

0.99425

0.002266

5(10)=0.968144

0.098942w

0.99716 0.0025011.6-3.0

4(10)=1.035378

0.097514w

0.99746 0.002602

5(10)=1.026468

0.135133w

0.99879 0.0024863.0-5.0

4(10)=1.098068

0.118847w

0.99999 0.000825

5(10)=1.034865

0.139598w

0.99947 0.0022465.0-6.4

4(10)=1.060569

0.111773w

0.99954 0.001269

5(10)=0.866450

0.105842w

0.99747 0.002821大氣透過率的估計(jì)方程:Summerprofile大氣透過率的估計(jì)方程:Winterprofile0.4-1.6

4(10)=0.983311-0.072444w

0.99469 0.002215

5(10)=0.981868-0.121979w

0.99679 0.0028961.6-3.0

4(10)=1.058059-0.121354w

0.99817 0.002749

5(10)=1.048364-0.163678w

0.99948 0.0019733.0-4.4

4(10)=1.111348-0.140080w

0.99998 0.000262

5(10)=1.033785-0.160330w

0.99934 0.0021844.4-5.4

4(10)=1.071156-0.131166w

0.99969 0.000966

5(10)=0.879292-0.125053w

0.99821 0.0022145.4-6.4

4(10)=0.964106-0.111430w

0.99907 0.001424

5(10)=0.681518-0.088431w

0.99649 0.002195wgcm-2EstimationequationRSEEDependenceoftransmittanceonzenithangleofviewing

4(

)=(-2.3994

10-3)+(2.2976

10-5)

2 R2=0.998487SEE=0.00043

5(

)=(-3.2766

10-3)+(3.1454

10-5)

2 R2=0.998729SEE=0.00054大氣透過率的遙感視角校正

=

0*|Di-D0|Where:

0=2arctg(hps/H)=0.0757°

Di

locationofthepixel D0locationofthenadirpixelhps=halfpixelsize=0.55kmH=satellitealtitude=833km

4(

)=

4(10)-

4(

)

5(

)=

5(10)-

5(

)

大氣透過率的估計(jì)四、地表比輻射率估計(jì)(略)參考:LandsatTM中的地表比輻射率估計(jì)ForconvenienceofexpressionwecomputetheprobableLSTestimationerror

Tsasfollows:

Ts=

Ts(x+

x)

Ts(x)

wherexisthevariablethatsensitivityanalysisorientsto(transmittanceoremissivity),

xispossibleerrorofthevariablex,andTs(x+

x)andTs(x)aretheLSTsimulatedbyouralgorithmforx+

xandx,respectively.

五、SensitivityanalysisFortheanalysisweuseassumptionofaveragegroundemissivity0.97andfollowthemethodcreatedby

SobrinoandCaselles[1991]todeterminetheemissivityofthetwochannelsas

4=0.967and

5=0.971.ForbrightnesstemperatureweassumethatT4isgreaterthanT5andthattheirdifferenceisT4-T5=0.7.FormostnaturalsurfacesoftheEarth,groundemissivityisbetween0.95and0.98[Humesetal.,1994].AsshowninFigure3,theT4-T5=0.7assumptionisrationalformostcases.ConsideredthepossibleLSTchangeoftheEarth,0

-70

CisselectedastherangeofbrightnesstemperatureT4change。ConditionsforsensitivityanalysisHistogramofT4-T5,computedfromtheIsrael-Sinai(Egypt)peninsulaimageacquiredonJuly18,1998.

Probablelandsurfacetemperature(LST)estimationerrorwithtransmittanceduetopossibletransmittanceerrorinchannel4Probablelandsurfacetemperature(LST)estimationerrorwithtransmittanceduetopossibletransmittanceerrorinbothchannels4and5.

ProbableLSTestimationerrorwithbrightnesstemperaturedifferenceduetopossibletransmittanceerrors(a)

4=0.01,(b)

5=0.01,(c)

4=

5=0.01,and(d)

4=

5=0.05.

AverageLSTestimationerrorduetotransmittanceerrorThisaverageLSTerroriscomputedunder

T4-T5

=1

Cwithtransmittancerange0.70-0.90andemissivityrange0.95-0.98.ProbableLSTestimationerrorduetosimultaneousemissivityerrorsinbothchannels4and5,illustratingthechangeof

Tsagainstemissivity.ProbableLSTestimationerrorduetosimultaneousemissivityerrorsinbothchannels4and5,illustratingthech

溫馨提示

  • 1. 本站所有資源如無特殊說明,都需要本地電腦安裝OFFICE2007和PDF閱讀器。圖紙軟件為CAD,CAXA,PROE,UG,SolidWorks等.壓縮文件請下載最新的WinRAR軟件解壓。
  • 2. 本站的文檔不包含任何第三方提供的附件圖紙等,如果需要附件,請聯(lián)系上傳者。文件的所有權(quán)益歸上傳用戶所有。
  • 3. 本站RAR壓縮包中若帶圖紙,網(wǎng)頁內(nèi)容里面會有圖紙預(yù)覽,若沒有圖紙預(yù)覽就沒有圖紙。
  • 4. 未經(jīng)權(quán)益所有人同意不得將文件中的內(nèi)容挪作商業(yè)或盈利用途。
  • 5. 人人文庫網(wǎng)僅提供信息存儲空間,僅對用戶上傳內(nèi)容的表現(xiàn)方式做保護(hù)處理,對用戶上傳分享的文檔內(nèi)容本身不做任何修改或編輯,并不能對任何下載內(nèi)容負(fù)責(zé)。
  • 6. 下載文件中如有侵權(quán)或不適當(dāng)內(nèi)容,請與我們聯(lián)系,我們立即糾正。
  • 7. 本站不保證下載資源的準(zhǔn)確性、安全性和完整性, 同時(shí)也不承擔(dān)用戶因使用這些下載資源對自己和他人造成任何形式的傷害或損失。

最新文檔

評論

0/150

提交評論