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

X射線熒光CT成像中熒光產額、退激時間、散射、偏振等關鍵物理問題計算與分析*

2021-11-01 06:10:26張芝振李亮
物理學報 2021年19期
關鍵詞:方向

張芝振 李亮?

1) (清華大學工程物理系,北京 100084)

2) (清華大學,粒子技術與輻射成像教育部重點實驗室,北京 100084)

X 射線熒光CT(X-ray fluorescence computed tomography,XFCT)是一種使用X 射線熒光(X-ray fluorescence,XRF)實現功能性成像的新技術,在生物醫學成像中表現出較大潛力.但是,X 射線穿過生物體的同時還會產生大量康普頓散射光子,對XRF 信號的采集形成很強的背景噪聲;因此,如何有效消除康普頓散射噪聲對于提高XFCT 成像質量至關重要.本文研究總結了XFCT 成像過程中涉及的物理過程,包括:熒光的產額、退激發時間、熒光發射角分布、熒光偏振態、康普頓散射角分布與散射光偏振態,并通過研究熒光與散射光物理性質的差異尋找去除康普頓散射噪聲的方法.經過物理過程推導和分析計算,發現:1) 高原子序數元素的K 層熒光退激發時間極短,在現有探測器的時間分辨率條件下,無法分辨散射光與熒光;2) K 層發射熒光的角分布各向同性,康普頓散射角分布在與入射光偏振方向附近取得最小值,而且入射光線偏振度越高,散射光的微分截面越小,偏振光源將有利于減少康普頓散射噪聲;3) K 層熒光線偏振度為零,而康普頓散射光子在一些散射方向上具有一定線偏振度,因此偏振態的差異可能用于降低康普頓散射噪聲.

1 引言

X 射線熒光CT(X-ray fluorescence computed tomography,XFCT)是一種針對高原子序數元素的功能成像技術.XFCT 通過探測X 射線與物體內高原子序數元素發生光電效應后產生的XRF 得到物質內該元素的濃度分布圖像,所以兼具透射型CT 結構成像和發射型CT 功能成像的優點.Boisseau 和Grodzins[1]利用美國國家同步輻射光源首次實現XFCT 實驗.同步輻射光源具有準單能和高線偏振度的特性,是XFCT 成像的理想光源,但是因為同步輻射占用場地和成本高昂的問題,所以基于同步輻射的XFCT 系統很難應用于實際醫學應用中.X 光機相比于同步輻射具有巨大的靈活性和經濟性優勢,近年來有很多基于X 光機的XFCT 系統設計的研究[2?9].因為X 光機光源采用電子軔致輻射的方式產生X 射線,所以其能譜是分布范圍較大的連續譜并且不具備線偏振特性.相比于同步輻射光源,X 光機光源會產生更大的散射背景,嚴重影響了XFCT 的成像質量.

去除散射背景的方法主要有兩種.第一種方法采用略高于和略低于熒光物質的K-edge 的單能入射光分別采取兩組數據,兩組數據之差可近似為熒光信號[10].第一種方法對光源單能性要求很高,一般采用同步輻射源.第二種方法利用散射光子能量連續分布的特點,通過三次多項式擬合的方式估計散射,進而在數據中減去散射的貢獻[2],該方法對探測器的能量分辨率要求高,一般使用單像素X 射線譜儀采集熒光信號,導致掃描時間很長.

逆康普頓散射源相比于同步輻射源不僅具有占地面積小和價格適中的優勢,而且具有同步輻射源準單能和線偏振的特性,適合應用于XFCT 成像系統設計中.逆康普頓散射源與同步輻射源降低XFCT 散射背景的原理都是應用了線偏振X 光康普頓散射的角分布特性.Chi 等[11]蒙卡模擬研究了逆康普頓光源線偏振特性對XFCT 成像質量的影響,發現逆康普頓光源能抑制XFCT 成像的康普頓散射背景約1.6 倍.

此外,熒光與散射光的物理性質的差別可能用于降低散射背景,比如是否有可能利用熒光與散射光飛行時間、角分布和偏振態的差別減少散射背景或尋找散射背景最小的探測方位.本文分析了XFCT中熒光與散射光相關物理過程的性質:熒光與散射光的角分布,熒光退激發時間與熒光產額,熒光與散射光的偏振態.這些物理過程決定了XFCT 熒光成像的極限,通過這些物理過程性質能夠指導XFCT 的系統設計,減少XFCT 的康普頓散射背景和提高XRF 信號強度.

2 XFCT 成像方法和物理原理

在XFCT 所涉及的能量范圍內,X 射線與物質相互作用方式有光電效應、相干散射和康普頓散射.在X 射線能量為幾十keV 時,相干散射截面與康普頓散射截面相當,所以在處理XFCT 散射問題時相干散射和康普頓散射一般都要考慮,其中相干散射主要沿著小角方向散射,對光子的傳播會有影響,但通過后面的討論可以看出,在大角度方向康普頓散射是主要的.

XFCT 成像系統幾何設計有很多種方式[2?12],其中基于小孔成像、扇束X 射線源和XRF 探測器陣列的成像系統設計能大幅提高掃描速度[5],其XFCT 實驗布置示意圖如圖1 所示.小孔準直器作用是使XRF 探測器具有空間分辨能力.示意圖中被掃描物體是一個包含四個小圓柱的大圓柱模型,四個小圓柱里含有不同濃度的某種高原子序數元素的溶液.圖1 中描繪了某一條光路以及該路徑上可能產生的XRF 光子和散射光子,入射光與溶液中的重金屬元素原子發生光電效應,重金屬原子發射熒光(特征X 射線),熒光經過小孔后,被XRF探測器陣列探測,同時散射光子也會進入XRF 探測器,形成散射背景.XRF 探測器可通過設置能窗閾值得到入射到每個探測器單元的相應能窗的光子計數.X 射線探測器陣列用于傳統的透射型CT成像,其得到的衰減系數分布圖可用于XFCT 成像的衰減校正.

圖1 基于小孔成像方式的XFCT 實驗設置Fig.1.Experimental setup based on pinhole imaging.

為描述入射和散射X 射線的偏振態,引入Stokes矢量[13]:

其中I表示光強度;P1和P2表示在與傳播方向垂直的平面內成 4 5?夾角的兩對正交方向下分別測量得到的平面偏振度;P3表示圓偏振度.光束的線偏振度和圓偏振度定義分別為[14]

其中PL表示線偏振度;PC表示圓偏振度;Imax和Imin分別表示經過理想線偏振片后光強的最大與最小值;IR和IL分別表示光束中右旋光和左旋光的強度.定義歸一化Stokes 矢量為

Stokes 矢量與正交方向選擇有關,選定ξ1?1 對應的線偏振方向為x軸,則ξ1+1 對應的線偏振方向為y軸,傳播方向為z軸正向,并規定:迎著光子傳播方向(z軸正向)看時,繞z軸順時針旋轉角度為正,繞z軸逆時針旋轉角度為負,則在旋轉ψ角度后的新正交方向下的Stokes 矢量 [ 1,ζ]T與原正交方向下的Stokes 矢量 [ 1,ξ]T由變換矩陣M聯系,即 [ 1,ζ]TM[1,ξ]T.M的非零元為

一束光經過某個相互作用后,其Stokes 矢量會改變,通過一個 4×4 變換矩陣T聯系相互作用前后的Stokes 矢量:

3 K 層熒光物理性質

3.1 熒光產額

某個殼層或子殼層的熒光產額定義為該殼層產生一個空穴后,通過發射熒光退激發的概率.K 殼層產生空穴后,不同躍遷方式產生的熒光的能量和強度不同.一般地,殼層越高的電子通過退激發填充K 層空穴,產生的熒光的能量越高,某種熒光的強度正比于該熒光對應躍遷的分支比.熒光產額的半經驗公式[15]:

式中,ωK表示K 層熒光產額;Z表示原子序數;Ci表示擬合參數.(6)式中參數Ci可以由實驗得到的不同元素熒光產額數據擬合得到.Bambynek 等[15]和Hubbell 等[16]結合理論計算與實驗給出擬合參數和元素熒光產額數據表.

K 層熒光產額ωK是原子序數Z的增函數,對Z >60 的元素,ωK>90%,所以高原子序數原子的K 層空穴態通過發射熒光退激的概率最大.因此XFCT 應該選擇高原子序數元素作為熒光元素,在高Z元素范圍內熒光產額基本飽和,提高Z對增加熒光產額的收益很低,但提高Z可增大熒光能量,能量高的光子更容易穿透物體被探測器采集到,從而降低統計噪聲.

K 殼層空穴態不同躍遷方式產生不同能量的熒光,熒光按末態空穴所在層分成兩組,即Kα代表來自L 層電子的躍遷產生的熒光,Kβ表示來自M 層及更高層電子躍遷產生的熒光,同一組內的熒光能量相近.Ertu?ral 等[17]給出了原子序數16≤Z≤92 范圍內59 種元素的Kβ/Kα熒光強度比測量值,Scofield[18]給出了10≤Z≤98 原子序數范圍內Kβ/Kα熒光強度比理論計算值,從實驗和理論計算結果可以看出,熒光強度比Kβ/Kα隨原子序數增大而增大,對于Z>60 的元素,Kβ/Kα熒光強度比在0.3 左右.通過熒光強度比和K 層熒光產額ωK可以得到K 層不同能量熒光的產額.

入射X 射線也會與重金屬原子的L,M 和N等更高電子殼層的電子發生光電效應,且發射的熒光能量隨層數增高而減小,其中熒光能量最高的L 層熒光的能量一般為十幾keV,屬于低能區的光子,幾乎無法穿透模體,所以L,M 和N 等更高電子殼層的熒光在XFCT 熒光成像中可以忽略.

3.2 熒光退激發時間

如果熒光退激發半衰期很長,且長于探測器的時間分辨率,那么可以采用極短時間的X 光脈沖作為光源,在散射光子到達探測器后,延遲一段時間采集熒光光子,這樣就可以在時間上分開散射和熒光從而去除康普頓散射背景;如果探測器的時間分辨率大于熒光退激發半衰期,那么熒光和散射光子到達探測器的平均時間差小于探測器時間分辨率,從而無法分辨熒光和散射光子.為了證明這種方法是否可行,需要計算熒光退激發時間.原子發生光電效應后處于激發態,激發態是不穩定的,假設t=0 時刻原子處于激發態,則在以后的t時刻體系仍然停留在這個態的概率是

式中,τ是激發態平均存在時間,亦稱激發態壽命,其與激發態能級寬度Γ關系為

K層熒光分寬度ΓR與K層熒光產額的關系是[15]

對于Z>40 的元素,K層空穴態能級寬度半經驗公式為[15]

通過(8)式和(9)式得到K層熒光壽命τK為

可見,K 層熒光壽命隨原子序數增加快速下降,對已知的XFCT 應用的最低原子序數的元素釓(Z64),K 層熒光壽命為0.03 ps.當前對伽馬光子探測器時間分辨率要求最高的PET 設備的符合時間分辨率在百皮秒量級[19],遠大于中高原子序數元素的K 層熒光退激發時間,因此,無法通過時間差別來分辨熒光與散射光子進而減少散射背景.

3.3 熒光發射角分布與線偏振度

原子內殼層熒光發射角分布和偏振態在理論和實驗上已經有很多研究[20?26]:若空穴態的總角動量量子數J1/2,則熒光發射角分布各向同性且線偏振度為零[20];若J>1/2,則熒光發射一般各向異性并且部分偏振[22].因為K 層空穴態J1/2[22],所以K 層空穴態發射熒光的角分布各向同性且線偏振度為零.L,M 等不同子層具有不盡相同的總角動量量子數,而且某些子層的空穴態會通過Coster-Kronig 躍遷將空穴轉移到其他子層,因此L,M 等殼層的空穴態發射熒光更加復雜多樣,剔除Coster-Kronig 效應的影響后,不同子層空穴態發射熒光角分布與線偏振度也滿足上述角動量規則[24].

4 XFCT 中的散射

4.1 散射光角分布

實驗室參考系下單個自由電子與線偏振光的康普頓散射截面Klein-Nishina 公式為[25]

式中,r0是經典電子半徑;k和k0分別是出射和入射光子能量,以電子靜質量能為單位;e'和e分別是出射光偏振方向和入射光偏振方向; d?是立體角元.假設在出射方向n放置偏振器,該偏振器透光軸沿e'方向,一個偏振態為e的入射光子經過康普頓散射到出射方向n,則該光子通過偏振器的概率正比于n方向上的 dσKN/d?[27].散射方向和偏振方向如圖2 所示,入射光沿Z軸正向入射,用矢量k0(|k0|k0)表示;散射光沿方向出射,用矢量k(|k|k)表示.e是描述入射光偏振方向的單位矢量,θ∠SOZ是散射角,φ∠AOX是散射的方位角.y和y'表示垂直于散射平面SOX的線偏振方向,x和x'表示平行于平面SOX的線偏振方向,位于平面x′Sy′內的e'表示散射光的偏振方向,e'與y'夾角為β,顯然,x',y'和e'垂直kf.從圖2描述的幾何關系得到 (e·e′)2cos2β(1?sin2θcos2φ).

圖2 偏振康普頓散射示意圖Fig.2.Schematic diagram of polarized Compton scattering.

對出射光相互垂直的兩個偏振方向x'和y'的分截面求和得到入射光為線偏振光且不測量散射光偏振的實驗中的截面[25]為

一般地,如果入射光是線偏振度為PL(0

式中,pp 表示部分偏振(partially polarized),相應的圖2 中X軸方向取為入射光透過理想偏振片后的透射強度最大方向.注意到公式中k/k01/[1+k0(1?cosθ)] 與 方位角φ無關,所以公式對φ積分后得到截面沿散射角θ的分布與線偏振度PL無關,進而總截面與PL無關.

下面討論(14)式的最小值,做變換xk/k0,則

式中,x∈[1/(1+2k0),1].記(15)式右邊方括號內的項為函數g(x,φ),則

?g/?x在x1/(1+2k0) 處的值為

令f(k0,a)0,得到k0關于a的唯一解,記為k0K(a).易知,當f(k0,a)≤0 時,方程h(x,a)0 在定義域內存在唯一解,記為xX(a).討論函數(15)的最小值得到如下結論:若k0

由于K(1+PL)≥K(1)≈1.74,熒光CT 所用X 射線能量一般小于0.2 個電子靜能(100 keV),遠小于K(1),所以最終得到截面最小值點對應的θmin:

不同入射光能量下的截面最小點θmin隨偏振度P的變化如圖3 所示,可見熒光CT 所用的X光能量范圍內,微分截面最小值對應散射角在90°附近.偏振度越高且入射能量越低,則θmin越接近90°.由于部分偏振光可以統計描述為完全偏振光(PL1 )和一個完全非偏振光(PL0)的疊加,所以最小微分截面方向接近入射光完全偏振成分的電矢量方向.

圖3 不同入射能量時,θmin 隨線偏 振度變化Fig.3.θmin varying with polarization for different incident energy.

圖4 描述了不同入射能量下,微分截面最小值隨偏振度的變化.可見在熒光CT 成像的能量范圍內,隨入射光偏振度的增大,最小微分截面值減小.這為熒光CT 去散射提供了思路,即入射光偏振度越高,在θθmin,φ0 or π 處的康普頓散射本底越低.

圖4 不同入射能量時最小微分截面隨偏振度的變化Fig.4.Minimum differential cross section varying with polarization for different incident energy.

前面討論的情形是單個自由電子與入射光的康普頓散射截面,這種情形很理想,方便分析.實際上物質中的電子都不是自由的,此時,一個原子的總散射截面與其各個電子在自由狀態下散射截面的總和一般不相等,引入非相干散射函數S(q,Z),則一個原子的非相干散射微分截面為[29]

式中,Z是原子序數;xsin(θ/2)/λ0,

λ0是入射光子波長(單位為?).非相干散射函數的理論計算方法和數值表可在相關文獻中查尋[30?33].非相干散射函數滿足 0 ≤S(x,Z)/Z≤1,對于低Z原子的散射,只有在x較小(軟X 射線或小角度散射,x<1 )時,S(x,Z)/Z才顯著小于1.例如100 keV 入射光子,波長約為0.124 ?,入射到氧原子(Z8 ),散射角滿足θ<4?才使S(x,Z)/Z <0.5,實際實驗中,入射光束總有發散角和寬度,在小散射角處測量熒光信號存在很強的入射光背景,所以實驗設置應該在較大角度處測量熒光信號,因此在這個測量角度范圍內單個自由電子微分截面最小值的分析依然是有效的,即考慮非相干散射函數后,入射光線偏振優勢方向是依然是近似的散射強度最小點.

單個自由電子的Klein-Nishina 散射微分截面公式在kk0→0 時得到單個自由電子的Thomson散射截面公式:

若入射光線偏振度為PL,則相應的散射截面公式為

由(22)式易知,最小微分截面方向為θπ/2 且φ0,π,所以最小微分截面方向為入射光完全偏振成分的電矢量方向.從經典電動力學角度出發也可以得到相同的散射截面公式,即電子在入射光電場驅動下振動輻射電磁波.

一個原子的相干散射截面為[25]

其中F(x,Z) 是原子形狀因子,可通過文獻[30?32]查詢原子形狀因子和理論計算方法.原子形狀因子隨x的變化趨勢與非相干散射函數相反,x越大F(x,Z) 越接近零,相反地,x趨于零時F(x,Z) 趨于Z.

根據(20)式和(23)式,一個原子的康普頓散射與相干散射微分截面的總和為

(20)式、(23)式和(24)式忽略光的偏振對原子形狀因子和散射因子的影響[26],所以根據(14)式和(22)式的最小值的討論,一個原子的總微分截面的最小值點一定滿足φ0,π.

一個分子的散射截面角分布為[34]

其中m(x) 表示分子干涉函數;s(x) 和f(x) 表示采用獨立原子模型(independent atomic modelling,IAM)計算得到的非相干散射函數和分子形狀因子[33],即

其中Fi(x,Zi) 表示分子中第i種元素原子的形狀因子;Si(x,Zi) 表示第i種元素的非相干散射函數;ni表示分子中第i種元素的原子個數;m(x) 表示分子干涉函數[34].

計算得到水分子在方位角φ0,π 處的微分散射截面隨散射角θ的分布如圖5 所示.從圖5(a)可以看出,相干散射微分截面隨散射角增大快速減小,在大角度處,微分散射截面由主要由非相干散射貢獻.由圖5(b)中總微分截面角分布隨線偏振度的變化可知,線偏振度越高,垂直方向散射的微分截面越小.PL0.5 時不同入射光能量下的總微分截面如圖5(c)所示,可見在XFCT 的光源能量范圍內,不同能量光子的水分子散射截面最小值均在入射光線偏振優勢方向附近,且該方向上微分截面相近.

圖5 φ=0,π 處微分截面隨散射角變化 (a) PL=1 時相干、非相干和總微分截面;(b) 不同線偏振度時的總微分截面;(c) PL=0.5 時不同入射能量下的總微分截面Fig.5.Differential cross section varying with θ at φ=0,π :(a) Incoherent,coherent and total differential cross section at PL=1 ;(b) total differential cross section varying with θ for different PL ;(c) total differential cross section varying with θ for different incident energy at PL=0.5.

4.2 散射光的偏振

聯系康普頓散射前后光束Stokes 矢量的變換矩陣T為[25]

其中k0k0n0,kkn,S是電子初始自旋方向,n0是入射光動量方向,n是散射光動量方向.在各向同性材料中,電子初始自旋方向的平均值為0,所以矩陣T中含S的矩陣元為零.在變換矩陣T成立的坐標系中要求垂直散射平面方向完全線偏振光的歸一化Stokes 矢量為即y和y'軸垂直于散射平面,由圖2 易知,坐標系xOy可 由XOY繞OZ軸旋轉ψ?φ得 到,則 由(4a)式、(4b)式和(4c)式可知,在XOY坐標系內Stokes 矢量為的入射光轉換到坐標系xOy后歸一化Stokes 矢量變為

在各向同性材料某點發生康普頓散射后,散射光Stokes 矢量變為 [If,ζ′]T,即:

由(29)式亦可得到散射光偏振的主方向與散射平面法向夾角β,符號規定為:迎著光子動量方向,在表示散射光偏振的局部坐標系x′Sy′內,y'軸以最小角度旋轉到e',若其旋轉方向為順時針,則β為正,否則為負.規定β的符號后,在局部坐標x′Sy′旋轉β角得到的新坐標系中,散射光Stokes矢量的第二個分量應達到正的最大值,據此條件可求得

其中函數 a tan2(y,x) 為雙參數反正切函數.

在φ0 or π 散射平面上散射光的歸一化Stokes矢量、線偏振度和偏振主方向角為β分別為

特別地,在θπ/2,φ0,π 散射方向上散射光的歸一化Stokes 矢量和線偏振度為

圖6 展示了在不同入射光能量下φ0,π 散射平面內隨散射角θ的變化.圖7 展示了不同能量下θπ/2,φ0,π 散射方向上散射光線偏振度隨入射光線偏振度的變化.若入射光為完全線偏振光(1 ),則θπ/2,φ0,π 散射方向上散射光的線偏振度為零,這與(12)式是相符的.從圖6 可以看出,在XFCT 所涉及的能量范圍內,相同散射角處不同能量下散射光線偏振度差距較小,而在θπ/2 附近,部分偏振的入射光,其散射光也具有較高偏振度.從圖7 可知,只有當入射光線偏振度接近完全線偏振時,θπ/2 的散射光的偏振度才顯著降低.

圖6 φ=0 or π 處 隨散 射角變化Fig.6. varying with θ at φ=0, π.

圖7 不同能量下 θ =π/2,φ=0, π 處散射光線偏振度隨入射光線偏振度變化Fig.7.The linear polarization of scattering photons varying with at θ =π/2,φ=0, π for different incident energy.

從(30)式、(31)式和(32a)式可知:完全非偏振入射光的散射光的偏振主方向始終垂直于散射平面,部分偏振入射光的散射光偏振;部分偏振入射光在φ0,π 平面內的散射光偏振度是散射角的函數,且存在兩個零點,左零點左側和右零點的右側偏振優勢方向平行于散射平面,在兩個零點之間,偏振優勢方向垂直于散射平面,在兩個零點處散射光完全非偏振;在φ0,π 散射平面內,完全線偏振入射光的散射光在θπ/2 時的偏振優勢方向均平行于散射平面,在θπ/2 時散射光完全非偏振.

根據以上關于散射光的線偏振度及熒光線偏振度的分析,并假設存在理想的X 光偏振片,因為K 層熒光的線偏振度為零,K 層熒光經過理想偏振片后減少50%,而散射光具有一定線偏振度,比如當入射光完全非偏振時,在單個自由電子這個理想的情況下,垂直于入射光方向的散射光的線偏振度可達到90%,經過理想偏振片后散射光可減少95%,因此,在經過偏振片后,雖然熒光強度減半,但是熒光強度與散射強度比值即信噪比在理想情況下可提高十倍,所以在原理上,散射光的線偏振度和K 層熒光線偏振度的差異可用于提高熒光CT 成像的信噪比.

5 結論

本文分析了XFCT 的物理過程及其性質.計算熒光退激發時間和熒光產額發現,中高Z元素熒光產額基本飽和,而熒光退激發時間遠遠小于現有探測器時間分辨率,因此,熒光退激發時間無法用于分辨散射光與熒光,提高原子序數可以增加熒光的能量,熒光能量越高,其穿透性越強,有助于提高探測到的熒光光子數,減小統計噪聲,提高信噪比.偏振光源有助于減少康普頓散射背景,入射光偏振度不影響康普頓散射總截面,但對截面的角分布有影響,部分偏振入射光在偏振方向附近具有最小的散射截面,在此方向上探測熒光光子有助于降低康普頓散射背景,并且對于完全非偏振或部分偏振的入射光,在大部分方向上的散射光仍然具有一定的線偏振度,而K 層熒光線偏振度為零,因此熒光和散射光線偏振度差異可以用于進一步降低熒光探測中的康普頓散射背景,由于實際成像過程的復雜性,通過偏振性質來減少散射背景需要進一步模擬和實驗驗證.

猜你喜歡
方向
2023年組稿方向
計算機應用(2023年1期)2023-02-03 03:09:28
方向
青年運動的方向(節選)
2022年組稿方向
計算機應用(2022年2期)2022-03-01 12:33:42
2022年組稿方向
計算機應用(2022年1期)2022-02-26 06:57:42
2021年組稿方向
計算機應用(2021年4期)2021-04-20 14:06:36
如何確定位置與方向
2021年組稿方向
計算機應用(2021年3期)2021-03-18 13:44:48
2021年組稿方向
計算機應用(2021年1期)2021-01-21 03:22:38
大自然中的方向
主站蜘蛛池模板: 伊大人香蕉久久网欧美| 精品亚洲欧美中文字幕在线看| 五月婷婷欧美| 超碰精品无码一区二区| 久久情精品国产品免费| 国产欧美视频在线观看| 国产激爽大片在线播放| 成年女人a毛片免费视频| 国产激情无码一区二区APP| 亚洲视频四区| 无码在线激情片| 亚洲色大成网站www国产| 亚洲无码熟妇人妻AV在线| 青青草原偷拍视频| 黄色网址手机国内免费在线观看| 97se亚洲综合| 欧美午夜网站| 国产一区二区三区在线观看视频| 亚洲成人手机在线| 毛片在线播放网址| 亚洲中文在线看视频一区| 无码精品福利一区二区三区| 精品视频91| 国产亚洲高清视频| 92精品国产自产在线观看| 无码丝袜人妻| 五月六月伊人狠狠丁香网| 欧美日韩北条麻妃一区二区| 国产又粗又爽视频| 在线看片中文字幕| 五月综合色婷婷| 国产精品蜜臀| 国产凹凸视频在线观看 | 97久久超碰极品视觉盛宴| 99尹人香蕉国产免费天天拍| 色综合天天视频在线观看| 最新国产网站| 91精品免费高清在线| 欧美啪啪一区| 午夜三级在线| 在线观看国产精品第一区免费 | 91麻豆国产视频| 在线看片免费人成视久网下载| 亚洲久悠悠色悠在线播放| 亚洲色图在线观看| 9cao视频精品| 成人在线欧美| 免费人成视网站在线不卡| 久久香蕉国产线看观看精品蕉| 米奇精品一区二区三区| 久996视频精品免费观看| 青青网在线国产| 国产精品女在线观看| 综合天天色| 蜜臀av性久久久久蜜臀aⅴ麻豆| 26uuu国产精品视频| 91精品人妻一区二区| 国产乱人免费视频| 精品91自产拍在线| 尤物成AV人片在线观看| 美女扒开下面流白浆在线试听 | 美女免费精品高清毛片在线视| 中美日韩在线网免费毛片视频| 色一情一乱一伦一区二区三区小说| 99re这里只有国产中文精品国产精品| 欧美天天干| 一级一级特黄女人精品毛片| 91麻豆精品视频| 免费女人18毛片a级毛片视频| 少妇高潮惨叫久久久久久| 99尹人香蕉国产免费天天拍| 国产一区成人| 夜精品a一区二区三区| 播五月综合| 69综合网| 性色一区| 国产91九色在线播放| 免费在线看黄网址| 67194亚洲无码| 国产精品福利一区二区久久| 国产内射一区亚洲| 精品伊人久久久香线蕉 |