<ol id="ebnk9"></ol>
    1. 工業開采誘發地震的應力降研究現狀和主要科學議題

      發布時間:2025-07-14 09:55:31   來源:心得體會    點擊:   
      字號:

      姜叢 尹欣欣 蔣長勝 翟鴻宇 王騫 張延保 來貴娟 尹鳳玲

      1)中國地震局地球物理研究所,北京 100081 2)甘肅省地震局,蘭州 730000

      在地殼淺部采用水力壓裂等方式進行的頁巖氣和干熱巖等非常規能源開采、廢水回注以及抽取方式進行的礦山開采,均可造成地殼應力狀態的擾動并產生誘發地震,給這些工業活動的順利實施以及地表社區與重大基礎設施的地震安全帶來較大挑戰。應力降Δσ是描述地震震源特征的重要參量,被定義為地震破裂前后斷層上的平均應力差(Hanks,1979)。由于應力降可被用來考察高頻地震動特征(Boore,1983)、區分所在地區的構造背景(Allmann et al,2009;
      Boyd et al,2017)、根據其自相似性給出關于地震危險性的認識(Abercrombie,1995;
      Oth et al,2010),因此其在天然地震的地球動力學研究和減輕地震災害風險研究中廣受青睞。在工業開采誘發地震活動中,應力降還可用來了解水力壓裂等工業活動與誘發地震震源性質之間因果關系(Yu et al,2020),以及作為重要參量來揭示斷層活化現象(Chen et al,2020),這使其成為目前誘發地震研究中的重點內容之一。尤其是開采活動中布設大量的密集地震觀測網絡,為更深層次利用應力降提供了新機遇(Shearer et al,2019)。

      (1)

      Δσ=cM0/S/L

      (2)

      從而實現了直接由地震波測量應力降的過程。由上述定義可見:①受到震源模型假設的影響,Δσ的具體表達形式依賴于震源模型,例如半徑為R的圓形斷層上的應力降為Δσ=7M0/16R3,長為L寬為w的矩形斷層上的應力降為Δσ=2M0/πw2L;
      ②由于斷層面積S≈L2,應力降Δσ與1/L3相關,其測定結果的不確定性很大(Allmann et al,2007);
      ③由于應力降的概念基于靜態位錯理論,僅考慮了斷層的相對移動,沒有考慮斷層的破裂,因此通常所說的靜態應力降要比實際的應力降小。

      鑒于上述關于應力降Δσ的定義和物理含義,以及其在新型工業活動和誘發地震災害風險管控中的重要價值,有必要系統全面地總結目前關于誘發地震應力降研究的新進展、技術發展趨勢和已獲得的科學認識。而在人為工業活動中的這些進展和認識,對天然地震的研究也同樣有重要啟發意義。本文嘗試對近20年的相關研究進行總結,以期為新型工業活動伴隨的誘發地震的減災研究以及天然地震的研究提供參考。

      包括應力降在內的震源參數的測定,需要分離或同時測定幾何擴散、介質衰減、場地響應等波形信號的耦合信息。新型工業活動誘發地震多為震級較小的微震事件,波形信號的低信噪比給測定工作帶來較大難度。目前已發展的測定方法有如下三類:

      (1)譜擬合方法是計算誘發地震的應力降的最為直接的方法。一般首先利用多錐形窗方法(Percival et al,1993)分別對三分量的地動速度譜進行平滑處理,并利用Boatwright點源模型擬合地動速度譜(Boatwright,1980),即

      (3)

      (2)單事件的經驗格林函數方法(EGF)是誘發地震應力降計算的重要方法,其優點是可一定程度提高拐角頻率fc和標量地震矩M0的測量精度(Kwiatek et al,2015)。EGF方法的原理是,對擬求解震源參數的地震,找到位置接近的小震級事件(震級差ΔM>1.0)作為經驗格林函數,利用譜閾的反卷積方法,消除待求解事件與EGF事件相同的場地響應、傳播路徑衰減項,從而獲得震源譜。其中擬求解事件i與EGF事件j的地動速度譜的比值為

      (4)

      對式(4)一般可采用最小二乘法擬合理論譜與實際的觀測比值譜的L2范數最小值來獲得參數結果(Holmgren et al,2019;
      Wang et al,2020)。在實際應用中,能否找到嚴格吻合EGF定義的小地震是制約經驗格林函數方法的重要因素(Huang et al,2019),此外對EGF方法計算結果分辨率的限度也有較多討論(Abercrombie,2015;
      Shearer et al,2019)。

      μ 為重力參數,ω =[00 0]探測器在軌道坐標系下相對于慣性坐標系的角速度,θ0為探測器的真近點角。設rc=[x y z]T探測器在軌道坐標系與小行星的相對位置,a=[axayaz]T是探測器各軸的控制加速度,是d=[dxdydz]T外部環境的擾動。

      (3)利用譜疊加和廣義反演方法(Castro et al,1990;
      Shearer et al,2006;
      Oth et al,2011)獲得更具普遍意義的經驗格林函數(EGF)是計算應力降的重要方法。例如,Warren等(2002)以及Shearer等(2006)發展了一種P波譜疊加方法,可以將震源譜與幾何擴散、傳播路徑的衰減、場地響應分離出來。該方法通過將同一臺站大量的地震記錄進行疊加,迭代后消除幾何擴散、路徑衰減等差異項,再利用震級分組的震源項與理論震源譜的擬合獲得修正后的經驗格林函數,從而更為完整地去除幾何擴散、路徑衰減和場地響應,獲得每個地震的震源譜,并利用Brune圓盤震源模型計算應力降等震源參數。Trugman等(2017)在該方法的基礎上,發展了非參數重采樣方法來估計震源參數的不確定性。為提高譜疊加和廣義反演方法中的震源譜參數反演的可靠性,Picozzi等(2017)還發展了一種基于數據驅動策略的無參數計算方法,其核心是設計了震源譜lg10(S(f))對各個震源參數的無量綱的偏導數作為頻域相依的敏感度函數,作為反演的權重函數。

      目前在應力降的數值上已取得一定共識。工業開采誘發地震的應力降通常較震源深度更深的天然地震的應力降值偏小(Fehler et al,1991;
      Abercrombie et al,1993;
      Hua et al,2013;
      Hough,2014)。已有的對全球天然地震的研究表明,應力降至少在三個數量級上變化(例如,Shearer et al,2006;
      Allmann et al,2007、2009),數值通常在1~100MPa之間(Abercrombie,1995),且符合正態分布(Allmann et al,2009)。而對工業開采誘發地震,應力降數值則通常在0.01~10MPa之間。本文羅列了一些工業開采區誘發地震應力降的測定結果,以便進行比較研究。

      在北美有大量的誘發地震應力降測定案例。例如,Wu等(2019)對2015年美國俄克拉荷馬州加斯里(Guthrie)MW4.0誘發地震應力降的測定結果為3.4MPa,這一數值接近Wu等(2018)測定的俄克拉荷馬州其他誘發地震應力降的中值(約2MPa)。Kwiatek等(2015)采用Madariaga震源模型觀測到美國加利福尼亞州蓋爾瑟斯地熱田Prati-9井附近的一組誘發地震,震級分布在MW2.0~3.8之間,這些地震的平均應力降約為7MPa,但明顯低于Viegas等(2011)采用經驗格林函數方法對同一水力壓裂注入井測定的應力降中值(28MPa)。Jeong等(2020)對美國得克薩斯州沃思堡盆地的90次誘發地震(ML2.0~3.9)的應力降測定結果顯示,平均應力降約為4.46MPa。在加拿大不列顛哥倫比亞省東北部的蒙特尼區塊,通過P波估算的應力降中值為0.08MPa,對數標準差為0.8(以lg為單位),同時利用S波估算的應力降中值為0.03MPa,對數標準差為0.7(Yu et al,2020)。在加拿大西部沉積盆地(WCSB),Holmgren等(2019)利用S波測定了87次震級為M2.3~4.2的誘發地震的應力降,發現中小震級事件的平均應力降為 7.5±0.5MPa,但測定結果在0.2~98MPa之間變化,且有一個極端事件的應力降測定值為370MPa。

      在歐洲和亞洲,部分頁巖氣和干熱巖開采誘發地震的應力降得到了測定和研究。其中,Kettlety等(2020)測定了英國蘭開夏郡普雷斯頓新路PNR-1z井2018年發生的震級為ML1.5的誘發地震,平均應力降值為1.0MPa。對于瑞士巴塞爾地熱田附近2006年發生的MW3.2±0.1誘發地震,以及圣加侖地區2013年發生的MW3.4±0.1誘發地震事件,Edwards等(2015)測得的平均應力降分別為3.5MPa和3.0MPa。在芬蘭赫爾辛基,Kwiatek等(2019)對OTN-3井附近震級為MW2.0的誘發地震測得的平均應力降為1.6MPa。在亞洲的韓國浦項地區,Song等(2019)對PH02井附近的誘發地震計算獲得的平均應力降約為2MPa,與Chai等(2020)得到的1.92MPa結果一致。

      表1給出了這些誘發地震應力降測定結果。由上述案例可見,對于多數震級在1.5~4.0的誘發地震,平均應力降分布在0.03~7MPa之間,且多數主要分布在1.0~7MPa之間。目前認為誘發地震應力降的數值相比天然地震更小的現象,主要是由誘發地震發生在地殼更淺部造成的。

      表1 全球部分地區工業開采誘發地震的平均應力降測定結果

      誘發地震的發生機制與天然地震明顯不同,至少包括注入流體影響巖石孔隙壓力、注入的較冷流體與高溫巖石相互作用引起熱彈性應變和裂縫滑動(Rawal et al,2014;
      尹欣欣等,2021)等較為獨特的機制,因此影響應力降數值的因素可能也有所不同。下面分別從主觀因素造成測定結果的不可靠性以及客觀因素造成測定結果的復雜性做出簡要總結。

      震源模型和常數參數的設置,是應力降計算的重要假設,可直接影響應力降的計算結果。例如,Yu等(2020)使用圓盤震源模型(Brune,1970)和取值S波的k=0.37,相比于使用Madariaga模型(Madariaga,1976)以及取值k=0.21,其計算得到的應力降較后者小約5.5MPa。Yamada等(2007)發現使用圓盤震源模型(Brune,1970)要比使用Madariaga模型(Madariaga,1976)獲得的應力降小,并分析認為使用Madariaga模型更為合理。破裂速度的取值也是應力降計算中的重要影響因素,較慢的破裂速度將導致震源半徑的估計值較小、應力降的估計值較大。而破裂速度的取值通常帶有明顯的主觀性,例如在Madariaga模型中常取剪切波速 3.9km/s的90%,這在采用中值計算平均應力降的研究中影響更為明顯(Allmann et al,2009)。地震波傳播介質的品質因子Q值的計算和取值,也會顯著影響應力降測定的可靠性(Ko et al,2012)。Q值對巖石的力學性質非常敏感,其空間變化可以作為水力壓裂影響局部應力狀態的指標(Yu et al,2020),但當對衰減校正不充分時常??傻玫捷^低的應力降測定結果(Tomic et al,2009)。由于上述震源模型和參數設置的人為主觀因素影響,Holmgren等(2019)、Yu等(2020)和Wang等(2020)建議,應力降的絕對數值只有在使用相同的震源模型和參數設置的條件下才可以進行比較。

      一些客觀因素也顯著影響了應力降的分布差異,包括構造環境的差異、存在的水量、破裂速度的變化、斷層面正應力的變化等(Tomic et al,2009)。這些因素導致應力降數值時空分布出現差異,另一方面也使得根據應力降數值反推地下施工區或儲層應力狀態變得更為復雜。這些影響因素至少包括:①水熱型地熱田和增強型地熱田(干熱巖)的高溫特征,對地殼的摩擦強度等力學行為有重要作用,例如Oth(2013)對日本地區的應力降的研究表明,低應力降的數值在空間上與高熱流區有關;
      ②誘發地震在深度分布上的差別以及流體介入地震破裂過程的程度也能明顯影響應力降的數值,這是由于震源深度較淺和流體注入時,可能會降低有效應力,從而導致剪切應力更容易克服莫爾-庫侖破壞準則進而產生誘發地震(Goertz-Allmann et al,2013);
      ③誘發地震的應力降在空間上具有方向性差異,且對中等大小的地震(Lui et al,2019)和微震(Folesky et al,2016)均可觀測到。而當地震記錄的臺站覆蓋較為稀疏時,應力降測量方式不當則可能帶來結果的較大不確定性,例如僅從斷層的前向或后向測量,可能會將拐角頻率fc的測定結果帶來2倍誤差,并導致應力降結果達到8倍偏差(Holmgren et al,2019)。

      4.1 誘發地震應力降的空間相關性

      探討誘發地震應力降的空間相關性問題具有重要的科學價值。這是由于斷層帶和裂隙上可能存在應力的不均勻性(Chen et al,2020;
      姜叢等,2022),應力降的空間變化可以揭示斷層強度的非均質性,以及反映斷層帶的水文地質性質,這有助于理解誘發地震的發生機制。

      應力降是否具有“近低、遠高”分布現象近年來在誘發地震研究中成為熱點問題。Yu等(2020)在對加拿大不列顛哥倫比亞省東北部的蒙特尼區塊水力壓裂井周邊3km范圍內的484次地震(M-1.0~3.0)的研究中發現,近井端的應力降(0.1~1MPa)要低于遠井端的應力降(1~10MPa),且約小一個數量級。類似的現象也在加拿大阿爾伯塔省克魯克湖(Crooked Lake)附近的注水誘發地震(Clerc et al,2016),以及瑞士巴塞爾和薩爾瓦多柏林地區的增強型地熱系統(EGS)開采區誘發地震研究(Goertz-Allmann et al,2011;
      Kwiatek et al,2014)中被證實。上述例子中應力降數值與井的距離相關的現象,在統計學意義上也得到了證明(Yu et al,2020)。Yu等(2020)認為,應力降的這種空間分布特點可能是由于井筒附近水力壓裂受損巖石的破裂速度較慢,裂隙密度高或孔隙壓力較高,阻止了井筒附近地殼巖石存儲和釋放較大應力,而遠端地震事件的應力擾動受水力壓裂的影響則相對有限或幾乎可以忽略。

      應力降是否具有空間分布的恒定性也是誘發地震研究的關注點之一。在對加拿大不列顛哥倫比亞省東北部的蒙特尼區塊水力壓裂誘發地震的研究中,Yu等(2020)還發現無論是發生在距離井筒近處還是遠處的地震,應力降的數值隨著時間變化而大致不變。Kwiatek等(2015)和Picozzi等(2017)分別采用不同的應力降計算方法,在對美國Geysers地熱田Prati-9井誘發地震的分析中,也表示并未觀察到靜態應力降隨時間的明顯變化。這種空間分布的恒定性,可能表明具有相似流變特性的巖石中的動態破裂通常具有相對恒定的應力降(Yu et al,2020)。還有一些研究發現,地震序列的應力降具有隨時間變化的特點,主要表現在主震后應力降隨時間增加,在余震序列的后期觀察到的應力降總體較高(Sumy et al,2017;
      Chen et al,2020;
      Goertz-Allmann et al,2011),Goertz-Allmann等(2011)認為這是孔隙壓力過高導致的。但這些序列的應力降時序變化是否在同一位置、余震序列活動的短期調整變化是否能影響空間分布的恒定性,還缺乏更深入的研究。

      4.2 誘發地震應力降與震源機制類型的相依性

      即使在天然構造地震研究中,應力降與震源機制類型的依賴性也是重要的研究課題。根據安德森斷層理論,逆斷層作用的剪應力最高,正斷層作用的剪應力最低,走滑型斷層的剪應力介于前兩者之間,盡管應力降和絕對應力之間的關系尚不清楚,但傳統觀點認為應力降和絕對應力在與震源機制類型的關系上相同(McGarr,1984;
      McGarr et al,2002)。但也正因為應力降和絕對應力之間的關系不明確,一些新的觀測現象在逐漸揭示其復雜性。例如Allmann等(2009)在對地震應力降的全球變化的研究中觀察到了應力降中值對震源機制類型有依賴性,其中走滑事件的應力降最高,中值約為10MPa,而正斷層和逆斷層事件的應力降較低,中值為2~3MPa,并未顯示逆沖型地震具有更高的應力降。

      對于誘發地震,Kwiatek等(2015)對美國蓋爾瑟斯地熱田Prati-9井水力壓裂誘發地震的研究表明,盡管正斷層型地震的應力降略低于走滑型地震的應力降,但由于走滑型地震在儲層深部分布的數量更多,因此這種應力降與震源機制類型的相關性難于驗證。研究這一問題的難點在于影響應力降數值大小的因素,例如在一些應力降的計算中假設斷層破裂速度是恒定的,而Allmann等(2009)采用可變的斷層破裂速度來計算應力降時,可以解釋走滑型地震具有較高的應力降可能是由更快的斷層破裂速度造成的。

      4.3 誘發地震應力降的非自相似性

      應力降是否具有自相似性在物理上具有特殊意義,其實質是大地震和小地震是否具有相似的動態破裂過程。在實際應用中,同一地區小地震的應力降空間分布是否與中等以上地震存在相關性、能否用平均的應力降代表這一地區的總體應力降水平,直接關系到是否可用數量更多的小震級的誘發地震來預測未來破壞性誘發地震的應力降,并減少概率性地震危險性評估和地震預警中的地面運動不確定性(Hardebeck,2020)。

      對應力降的自相似性/非自相似性研究,最重要的研究方式是考察應力降是否與震級或標量地震矩存在相關性。Holmgren等(2019)在對加拿大西部沉積盆地87次M2.3~4.2級誘發地震應力降的測定中,發現應力降隨震級增大呈上升趨勢。Novakovic等(2018)和Chen等(2020)在對美國俄克拉荷馬州以及Picozzi等(2017)對加利福尼亞州蓋爾瑟斯地熱田進行誘發地震研究中,也均發現較大震級的誘發地震具有相對更高的應力降。Chai等(2020)對韓國浦項誘發地震序列的應力降研究中也獲得了同樣結論。其他案例還包括Mayeda等(2005)和Venkataraman等(2006)等的研究。概括起來,目前支持工業開采誘發地震應力降具有自相似性的案例總體上較少,僅有Yamada等(2007)對南非金礦開采誘發的20次微震的靜態應力降的研究等少數案例,此外,盡管Picozzi等(2017)在美國加州蓋爾瑟斯地熱田誘發地震的研究中觀察到了大地震群體具有自相似趨勢,但對單個地震序列的研究往往發現強烈的非自相似行為。

      4.4 誘發地震應力降的深度依賴性

      在天然地震的應力降研究中,許多案例顯示出應力降與震源深度存在依賴性(Hardebeck et al,2009;
      Allmann et al,2007;
      Baltay et al,2017),并認為這可能與斷層的摩擦強度沿深度的變化有關。工業開采誘發地震多發生在地下3~10km深度的地殼淺部,表1給出的大量案例也表明誘發地震的應力降總體相比天然地震偏小,因此目前普遍持有的觀點是造成誘發地震應力降比天然地震偏小的原因也在于震源深度的影響(Satoh,2006;
      Hough,2014;
      Yenier et al,2015;
      Zhang et al,2016;
      Atkinson et al,2017)。

      但對無論是誘發地震還是天然地震是否存在應力降大小的深度依賴性目前還有較大爭議,例如Wang等(2020)認為在計算應力降時假設剪切波速恒定可能導致人為的深度依賴性,特別是發生在地殼淺部5km以上的地震,Allmann等(2007)也持有類似觀點。Holmgren等(2019)進一步解釋認為,由于應力降與速度的立方呈正相關,也認為由于剪切波速隨深度增加,應力降的估計對假定的震源深度很敏感。在具體的應力降測定研究中,Allmann等(2009)對全球范圍的天然地震研究表明,盡管存在較大的分散性,但平均應力降估計值隨深度的變化很小。進一步從不同震源機制類型的角度分析,Allmann等(2009)還發現正斷層事件和走滑事件的平均應力降均沒有深度依賴性,而逆斷層作用下平均應力降隨深度的增加與Bilek等(1998、1999)在俯沖帶地震中發現的結果一致,可以用深度相關的剛度變化來解釋。左可楨等(2021)使用譜比法分析四川長寧地區微震震源參數之間的相互關系以及應力降的時空分布特征,也得到應力降總體與震源深度關系并不明顯的結果。

      盡管正反兩種觀點的爭議較大,表1給出的大量案例至少給出了誘發地震的應力降數值分布范圍更窄、總體比天然地震偏小的初步認識。由于應力降測定結果的主觀可靠性、客觀影響因素較為復雜,這也使得對這一問題的爭論和研究在未來仍具有開放性。

      鑒于應力降在研究工業開采誘發地震的震源特征、高頻地震動特征、區分構造背景和地震危險性分析等諸多領域的重要價值,本文系統地總結了近20年來相關方向的研究進展。分別從計算所用理論方法、獲得的平均應力降數值、影響應力降數值的主客觀因素、圍繞應力降的重要科學議題4個方面做了總結。獲得的主要認識如下:

      (1)在工業開采誘發地震的理論計算方法上,針對誘發地震微震較多、信噪比較低、實時快速計算需求明顯等特點,目前主要發展了三類方法,包括直接采用譜擬合方法求取震源參數、利用小震級事件作為經驗格林函數(EGF)提高拐角頻率fc和標量地震矩M0測量精度以及利用譜疊加和廣義反演方法對多臺站多事件同時計算震源參數。不同方法尚存在明顯的不足之處,這也限制了對應力降數值本身和分布規律探索的過高期待。

      (2)對于工業開采誘發地震的應力降平均值,本文收集了相關案例,結果顯示多數震級在1.5~4.0的誘發地震平均應力降分布在0.03~7MPa之間,且多數主要分布在1.0~7MPa之間,應力降數值相比天然地震分布的范圍更窄、總體偏小。對這一現象的主要解釋集中在誘發地震往往發生在地殼更淺部,以及流體更多地參與誘發地震發生過程等因素上。

      (3)影響工業開采誘發地震測定數值的因素,在測定方法的可靠性上,主要受到采用相對簡單的圓盤模型等震源模型假定,以及對波速、斷層破裂速度、地殼介質品質因子Q值的具體取值差異等因素影響;
      此外,構造環境的差異、存在的水量、破裂速度的變化、斷層面正應力的變化等客觀因素也影響了應力降的數值本身。

      (4)圍繞工業開采誘發地震應力降的科學問題,總結了應力降是否具有“近低、遠高”分布特征以及是否存在空間分布恒定性等空間分布依賴性問題,以及誘發地震的應力降是否與震源機制類型存在相依性、是否存在相似性或非自相似性、是否存在深度依賴性等4個問題。其中對應力降存在距離依賴性、存在非自相似性有一定程度的共識,但鑒于應力降測定可靠性和多種客觀因素的影響,這些科學問題未來均有較大的開放性、需要持續探索。

      工業開采區通常布設更為密集、覆蓋程度更好的現代化地震臺網,這為誘發地震獲得更為可靠的包括應力降在內的震源參數提供了便利。但通過本文的總結可見,受觀測技術和測定理論方法的局限,誘發地震的應力降測定的可靠性問題并未獲得實質性突破。對誘發地震的應力降很難精確測量似乎成了為數不多的共識性問題(Shearer et al,2019),Stein等(2009)甚至認為當前的靜態應力降研究僅反映了觀測譜或震源譜的一般特征、無法直接反映震源的物理特性,甚至會誤導對真實的斷層活動的認識。由于應力降是由拐角頻率通過假設的震源模型和參數間接計算得到的,對此Holmgren等(2019)和Shearer等(2019)建議應該轉而關注拐角頻率fc,這是相比應力降更基本、有用且更接近實際數據的震源參數。

      致謝:中國地震局地球物理研究所博士研究生張琰在文獻檢索上予以幫助,郭祥云高級工程師等專家提供了諸多建議,在此一并表示感謝。

      猜你喜歡震級震源斷層多種震級及其巧妙之處*地震科學進展(2022年10期)2022-10-14如何跨越假分數的思維斷層數學教學通訊·小學版(2022年4期)2022-05-29基于累積絕對位移值的震級估算方法自然災害學報(2022年2期)2022-05-10嘛甸油田喇北西塊一區斷層修正研究西部探礦工程(2022年2期)2022-02-14X油田斷裂系統演化及低序級斷層刻畫研究非常規油氣(2021年5期)2021-11-13地震后各國發布的震級可能不一樣?奧秘(創新大賽)(2021年3期)2021-05-15Pusher端震源管理系統在超高效混疊采集模式下的應用*石油管材與儀器(2020年2期)2020-05-11一種改進的近斷層脈沖型地震動模擬方法震災防御技術(2019年3期)2019-06-02新震級標度ML和MS(BB)在西藏測震臺網的試用科技視界(2018年30期)2018-01-311988年瀾滄—耿馬地震前震源區應力狀態分析地震研究(2014年3期)2014-02-27
      国产另类无码专区|日本教师强伦姧在线观|看纯日姘一级毛片|91久久夜色精品国产按摩|337p日本欧洲亚洲大胆精

      <ol id="ebnk9"></ol>