?

基于GA.GRNN地應力場反演的超長深埋隧洞圍巖穩定性研究

2024-02-18 18:50郭沖王翰盛鄧偉杰趙順利杜衛長
人民珠江 2024年1期
關鍵詞:遺傳算法神經網絡

郭沖 王翰盛 鄧偉杰 趙順利 杜衛長

摘要:基于超長距離引調水工程地質條件及實測地應力資料,建立了工程引水線路初始地應力場三維反演計算模型。通過遺傳算法優化廣義回歸神經網絡方法(GA.GRNN),較為準確地反演了工程引水線路的初始地應力場?;诔跏嫉貞?,對長距離引水線路地應力進行分區,并在各分區選擇典型斷面進行隧洞圍巖穩定分析。研究結果表明:工程區應力場以多隆溝斷裂帶為界,西側為水平構造應力區,極高地應力區和高地應力區多分布于地質構造帶附近,需注意應力釋放導致的巖爆發生;多隆溝斷裂帶東側為自重應力區,中等地應力區多分布于此,圍巖穩定性相對較好。

關鍵詞:初始地應力;神經網絡;遺傳算法;圍巖穩定

中圖分類號:TV672+.1? 文獻標識碼:A? 文章編號:1001.9235(2024)01.0103.11

Research on Stability of Utralong and Deep.Buried Tunnel Surrounding Rocks Based on GA.GRNN Geostress Field Inversion

GUO Chong1,2,WANG Hansheng3,DENG Weijie1,2,ZHAO Shunli1,2,DU Weichang1,2

(1.Jianghe Engineering Inspection and Testing Co.,Ltd.,Zhengzhou 450003,China;

2.Yellow River Engineering Consulting Co.,Ltd.,Zhengzhou 450003,China;

3.Hangzhou Huajin Lianzhe Intellectual Property Agency Co.,Ltd.,Hangzhou 310051,China)

Abstract: According to the geological condition and measured geostress data of the utralong.distance water diversion project,a 3D inversion calculation model of the initial geostress field of the water diversion route was established.The initial geostress field of the water diversion route was exactly simulated through the method of generalized regression neural network (GRNN) optimized by genetic algorithms (GA.GRNN).Based on the initial geostress field,the geostress of the water diversion route was partitioned,and the stability of the tunnel surrounding rocks in typical sections was analyzed.The research results show that the geostress field is bounded by Duolong gully fault belt,and the west of the fault belt is a horizontal tectonic stress area.Extremely high and high geostress areas are distributed in geological structural belts where attention should be paid to rock bursts caused by stress release.The east of the Duolong gully fault belt is the self.weight stress area,and the moderate geostress area is distributed here.The stability of the surrounding rocks here is relatively good.

Keywords:initial geostress;neural networks;genetic algorithms;stability of surrounding rocks

地應力是賦存于巖體內的天然應力,對工程的設計、施工具有重要意義[1-3]。工程中地應力原位測量是提供巖體初始地應力最直接、有效的手段[4]。但由于地應力場成因復雜多變,測點相對分散,地應力實測結果往往具有一定的離散性,難以體現整個工程區地應力場的宏觀規律。特別是對于超長距離引調水工程,由于場地、經費等客觀條件限制,不能進行大規模、全范圍的原位測量。因此,需要在實測地應力資料的基礎上,結合實際工程地質條件,通過有效的分析方法進行反演計算,以獲得范圍更大、適用范圍更廣的初始地應力場[5-7]。

目前,相關科研工作者針對初始地應力場的反演分析開展了一定的研究工作,以多元回歸法為主,如張建國等[8]、裴啟濤等[9]采用多元回歸法,對工程地應力場反演進行了深入研究。多元回歸法地應力狀態表達式較為單一,能較好地適用于小范圍的地應力場反演,如水電站地下廠房、穿山公路隧洞、邊坡、煤層巷道開挖的地應力反演。而對于超長距離的引調水工程,引水隧洞穿越多地層和應力區,應力狀態錯綜復雜,單一的應力表達式無法完全表示工程區全段的應力狀態,適用性有限。隨著計算機技術的發展,灰色理論、神經網絡法、遺傳算法等智能方法被應用于巖體初始地應力場的反演計算。如戚藍等[10]將系統工程灰色理論引入到地應力場分析中;裴書鋒等[11]、馬玉巖等[12]、陳正林等[13]將人工神經網絡應用于地應力場的反演計算;張金[14]、楊志強等[15]、謝學斌等[16]基于遺傳算法開展了巖體初始應力場反演研究。但傳統的神經網絡和遺傳算法在反演計算中存在精度和效率上的不足,因此有必要建立一種適用于長距離引調水工程的考慮多重地質因素、計算效率高的反演計算方法。

針對隧洞圍巖穩定性研究,科研工作者建立了Mohr.Coulomb、Drucker.Prager、概率強度、非線性統一等圍巖強度準則[17-20]和容許位移、容許位移速率、變形速率比等變形準則[21-23],以及圍巖松動圈理論、圍巖局部失穩監控指標體系、圍巖開挖位移預警值等圍巖破裂區判據等[24-26]。針對高水平構造應力區深埋超長隧洞的圍巖穩定性研究相對較少,當前研究的重點大多集中在自重應力的影響上。同時地應力分區標準和研究斷面的選取較為粗略,并不完全適用于超長引調水工程。因此,有必要對工程沿線地應力場進行細化分區研究,并基于地應力分區結果和地層巖性對隧洞典型斷面進行精準圍巖穩定分析。

基于此,針對長距離引水工程地質條件及實測地應力資料,建立引水線路初始地應力場三維反演計算模型,通過遺傳算法優化神經網絡方法(GA.GRNN),反演計算工程引水線路的初始地應力場?;诔跏嫉貞?,對長距離引水線路地應力進行細化分區,并在各分區選擇典型斷面進行隧洞圍巖穩定分析,為工程的設計和施工提供重要依據。

1 高水平構造應力區地應力實測

1.1 工程概況

青海省某超長距離引調水工程由引水工程和供水工程組成,工程區位于青藏高原東北部,海拔高度多在2 500~3 000 m以上。日月山、野牛山、拉脊山自西北向東南斜貫工程區。區域內地形地貌主要有構造剝蝕、侵蝕高山-中高山、構造侵蝕中山丘陵和侵蝕堆積河谷三大地貌單元。工程位于青藏高原東北部構造帶,地處祁連、青海南山-拉脊山兩大山系的交接處,區域大地構造上屬祁連中間隆起帶、南祁連褶皺帶及秦嶺褶皺帶內。受多條構造帶的影響,工程區域內地應力場分布較為復雜。

1.2 水壓致裂法地應力測試

在工程引水線路開展了SZK03、SZK13和SZK30的水壓致裂法地應力測試工作,地應力實測結果見表1。

水壓致裂法地應力實測結果表明,工程測試區域內SH范圍為9.85~35.48 MPa,Sh范圍為8.01~22.45 MPa,最大水平主應力方向為N32°~79°E,應力場狀態以NE~NEE向擠壓為主。

SZK13和SZK30工程區域地應力分布規律為SH>Sv≈Sh,區域內應力以水平構造應力為主,區域構造作用比較強烈;SZK03工程區域地應力分布規律為Sv>SH>Sh,區域內應力以垂直應力為主,為自重應力區,表明工程引水線路地應力場的復雜性。

2 基于GA.GRNN的地應力反演分析模型

2.1 初始地應力場反演計算原理

地應力反演原理[9]是將地應力回歸計算值k作為因變量,把數值計算求得的自重應力場和構造應力場相應于實測點的應力計算值σik作為自變量,回歸方程的形式為:

n——工況數。

對于m個觀測點,最小二乘法殘差平方和為:

式中 σ*jk——k觀測點j應力分量的觀測值;σijk——i工況下k觀測點j應力分量的數值計算值;j=1~6,對應初始應力的6個分量。

解此方程,得n個待定回歸系數L=(L1,L2,…,Ln)T,則計算域內任一點P的回歸初始應力σjp,可由該點各工況計算值迭加而得:

2.2 遺傳算法優化神經網絡(GA.GRNN)

廣義回歸神經網絡(GRNN)[11]是徑向基網絡的一種變形形式,可自動建立輸入和輸出的高維非線性映射,在解決復雜的非線性問題時,不需要假定具體的函數形式,從而大大簡化了求解的難度。廣義回歸神經網絡預測結果的好壞依賴于基函數光滑因子的選取,由于地應力變化規律不確定,無法用函數定向地確定光滑因子的最佳值,且試算法計算量較大,無法保證精度要求。遺傳算法[14]是一種全局最優化方法,適用于多極值點的優化問題,將選擇、交叉、變異等概念引入到算法中,通過構成一組初始可行解群體并對其進行操作,使其逐漸移向最優解。為此,引入遺傳算法,對光滑因子可能的取值進行全局搜索,自動搜索最優光滑因子參數值。通過把遺傳算法和廣義回歸神經網絡結合到一起,構建出GA.GRNN方法,計算流程見圖1。

采用遺傳算法搜索最優光滑因子,建立GA.GRNN模型,其具體步驟如下:①對樣本數據進行歸一化預處理;②確定遺傳算法參數;③遺傳算法初始化,生成初始種群;④GRNN讀入學習樣本進行網絡訓練,按給定的適應度函數進行適應度評價;⑤遺傳算法按照每個個體的適應度大小進行選擇、交叉和變異操作,得到新的種群;⑥判斷是否達到最大進化代數,若已達到,則停止計算,返回適應度最高的個體;否則轉至步驟④,直到達到最大進化代數;⑦輸出適應度最高的個體對應的實值數,即為最優的光滑因子;⑧用最優光滑因子建立GRNN模型,對測試樣本進行預測,得到預測結果;⑨對預測數據進行反歸一化,對GRNN網絡性能進行評估,并存儲數據。

2.3 計算范圍及計算模型

通過分析該工程區水文地質、工程地質條件和實測地應力測點的分布情況,建立反演分析計算模型。將引水隧洞線路方向設為X軸,垂直隧洞線路方向設為Y軸,豎直向上為Z軸。計算范圍為:X軸方向取30 000 m,Y軸方向取2 500 m,Z軸方向為海拔2 300 m至地表。

計算模型考慮引水隧洞沿線穿越的8條斷層,包含4條較大的斷裂:多隆溝斷裂(F1)、青海南山南緣斷裂(F2)、青海南山北緣斷裂(F3)、倒淌河—循化斷裂(F4)。工程區主要地層巖性分別為:印支期中三疊世花崗閃長巖(γδT2)、印支期中三疊世花崗閃長巖(ηγT2)、新近系貴德群(Ng)灰褐桔紅色砂礫巖、早中三疊系隆務河組(T1l)青灰色砂巖板巖互層夾灰巖、加里東期(γS.O3)灰白色花崗巖、中三疊系香阿洞組(T2x)灰-深灰色砂巖板巖夾灰巖、第四系上更新統洪積(Qpl3)角礫、砂、黃土狀亞砂土及黃土。數值模型采用四面體單元,模型網格單元數為426 182,節點數為80 543,區域三維模型和實體模型見圖2。

2.4 模型參數選取及邊界條件

反演計算過程中采用線彈性模型,根據工程勘察試驗結果,同時參考工程地質勘察報告中第四系松散層孔隙潛水、基巖裂隙水及相關鉆孔水位特征等信息,在選取計算參數時,有針對性地將部分參數參考試驗結果中的飽和狀態進行參數取值。計算采用的巖石基本參數見表2。

工程引水線路較長,跨越多個地質構造帶,地應力場分布規律較為復雜,回歸巖體初始應力場的7種基本因素見圖3。

通過查閱均勻設計表,確定包括自重影響系數在內的邊界條件一共7個參數,每個參數分為28個水平。根據遺傳算法求得廣義回歸神經網絡的最優光滑因子為0.553 8,分別代入實測點SZK30、SZK03和SZK13歸一化處理后的數據,得到的模型邊界條件見表3。

3 初始地應力反演結果

3.1 反演結果比較分析

將邊界條件代入數值分析模型進行計算,獲得模型各地應力分量,與實測結果對比見表4和圖4。

由表4和圖4可以看出,大部分測點地應力分量的反演值與實測值在數值上較為接近,應力分量平均絕對誤差為5.94 MPa,誤差相對較小。另外,相較于水平應力分量,垂直應力分量的誤差相對較大,且回歸值均高于實測值,這是由于垂直應力受地形地貌、歷史應力、構造作用等綜合因素的影響,其值一般為巖體自重應力和地質構造垂直應力疊加的結果,而水壓致裂法僅能測量水平主應力,表中所列垂直應力分量實測值為根據上覆巖體厚度進行的估算值,其值均小于實際垂直應力??偟膩碚f,反演結果與現場試驗測試結果和歷史地質構造資料中的地應力場一般規律相符,可較好地模擬工程引水線路的初始地應力場。

3.2 工程引水線路初始地應力場分布規律

基于GA.GRNN方法反演結果,對工程引水線路區域地應力進行反演合成計算,截取引水隧洞軸線縱剖面的最大主應力、中間主應力和最小主應力,見圖5。引水隧洞沿線區域初始地應力場最大主應力范圍為15~30 MPa,中間主應力范圍為10~20 MPa,最小主應力范圍為10~20 MPa。

引水隧洞軸線縱剖面的地應力矢量分布見圖6。由6a可以看出,構造應力方向與隧洞軸線方向近平行,引水隧洞沿線最大主應力矢量以多隆溝斷裂為界分為2個區域:第一區域為多隆溝斷裂以西的閃長巖一區,以水平構造應力為主,構造應力方向與隧洞軸線方向近平行,為構造應力發育區;第二區域為多隆溝斷裂以東的砂巖、板巖互層區和青海南山北緣斷裂帶(F2)東部閃長巖二區,以自重應力為主,為自重應力區。反演所用測點SZK30位于閃長巖一區,測點SZK13位于閃長巖二區的F3斷裂區,二者均以水平構造應力為主,測點SZK03位于閃長巖二區,以自重應力為主,反演計算地應力方向與實測地應力方向基本吻合。

4 高水平構造應力區隧洞圍巖穩定分析

4.1 地應力分級

初始地應力的判別與劃分,采用DL/T 5419—2009《水電水利工程地下建筑物工程地質勘查技術規程》中考慮地應力大小與強度應力比的綜合分級方案,強度應力比計算見式(4)。

式中 Rb——巖石飽和單軸抗壓強度,MPa;σm——最大主應力,MPa。

巖體初始地應力分級方案見表5。

基于地應力反演所得初始地應力場,提取引水隧洞沿線地應力數值,統計各測點巖體的強度應力比,并根據地應力分級方案對引水隧洞沿線地應力進行分區,結果見圖7。

根據圖7的統計結果,將隧洞沿線地應力按地層進行分類,結果見表6。

由表6可知,引水線路極高地應力區主要集中在斷層附近,典型斷面如7.8~8.4 km的閃長巖段,局部初始地應力達59.6 MPa。高地應力區分布范圍較廣,典型斷面如2.2~7.6 km和20.6~25.2 km的閃長巖段、8.8~11.6 km的花崗巖段、12.2~20.0 km和25.4~29.0 km的砂巖板巖互層段,初始地應力水平主要集中于20.0~50.0 MPa。中等地應力區主要集中分布于1.6~2.0 km和21.8~23.4 km的閃長巖段,其初始地應力水平相對較低。

4.2 隧洞圍巖穩定計算模型

隧洞開挖對應力場的擾動主要集中在一定范圍內,通常不超過3倍洞徑。利用反演計算的初始地應力場評價隧洞圍巖穩定時,在數值建模過程中,應建立足夠大的邊界,按3倍洞徑考慮,提取典型斷面處的初始應力場應力值,作為圍巖穩定分析的應力邊界條件。

選取3種地應力分區下隧洞的典型斷面(8.4 km;20.8、11.0、15.0 km;22.0 km)進行隧洞圍巖穩定分析。計算模型尺寸為100 m×50 m×50 m,共112 000個單元,119 700個節點。計算采用莫爾-庫倫力學模型和應力邊界條件,加載時初始化巖體內部的賦存應力,并施加邊界條件。根據工程勘察試驗結果和初始應力場反演結果,圍巖基本參數和初始地應力見表7。

4.3 隧洞典型斷面圍巖穩定分析

根據地應力分區結果,在引水線路中選取典型斷面進行圍巖穩定分析,以此評價施工期圍巖的安全性,并為保障引水隧洞的安全與快速施工提供數據支撐。

4.3.1 極高地應力區(8.4 km)

本區域所選的典型斷面為測點8.4 km處閃長巖斷面。隧洞開挖后圍巖的位移場、最大主應力見圖8。隧洞開挖后,圍巖最大變形為0.22 m,圍巖最大主應力高達65.0 MPa,未出現拉應力,受偏應力的作用,最大主應力出現偏轉,普遍集中于拱頂偏左和拱底偏右位置。鑒于圍巖應力較高,建議隧洞開挖時注意應力釋放導致的巖爆發生,同時在巖體應力釋放后及時支護,防止圍巖塑性區向深部擴展。

4.3.2 高地應力區(20.8、11.0、15.0 km)

高地應力區閃長巖段(20.8 km)隧洞開挖后圍巖的位移場、最大主應力見圖9。洞段開挖后,圍巖最大變形為0.13 m,圍巖最大主應力為52.0 MPa左右,未出現拉應力。受偏應力的作用,最大主應力的位置出現近45°的偏轉,普遍集中于拱頂偏右和拱底偏左位置。

高地應力區花崗巖段(11.0 km)隧洞開挖后圍巖的位移場、最大主應力見圖10。洞段開挖后,圍巖最大變形為0.12 m,該洞段圍巖最大主應力為36.0 MPa左右,未出現拉應力。

高地應力區砂巖與板巖互層段(15.0 km)隧洞開挖后圍巖的位移場、最大主應力見圖11。洞段開挖后,圍巖最大變形為0.20 m,洞段圍巖最大主應力為47.0 MPa左右,未出現拉應力。受偏應力影響,拱頂和拱底的最大主應力小于拱腰。

高地應力區圍巖應力較高,建議隧洞開挖時注意應力釋放導致的巖爆發生,并及時支護,防止塑性區擴展和應力松弛引發的掉塊。

4.3.3 中等地應力區(22.0 km)

中等地應力區閃長巖段(22.0 km)隧洞開挖后圍巖的位移場、最大主應力見圖12。洞段開挖后,圍巖最大變形為0.11 m,圍巖最大主應力為46.0 MPa左右,未出現拉應力。最大主應力變化明顯,受水平擠壓作用,隧洞拱頂偏左和拱底偏右處應力逐漸遞減至20.0 MPa。鑒于圍巖變形較小,測點22.0 km處閃長巖區隧洞圍巖安全性較好,巖體基本能夠實現自穩。

5 結論

結合高水平構造應力區長距離引水工程地質條件及實測地應力資料,建立引水線路初始地應力場反演計算模型,通過遺傳算法優化神經網絡方法反演得出工程引水線路的初始地應力場。對引水線路地應力進行細化分區,并在各分區選擇典型隧洞斷面進行圍巖穩定分析。通過研究得到如下結論。

a)地應力反演結果表明,引水隧洞沿線區域初始地應力場最大主應力范圍為15~30 MPa,約為最小主應力的1.5倍。構造應力方向與隧洞軸線方向近平行,工程區應力場以多隆溝斷裂帶為界,西側以水平構造應力為主,東側以自重應力為主。地應力反演回歸值與實測值較為吻合,為引水隧洞的設計和施工提供了合理的三維初始地應力場。

b)地應力分區研究表明,極高地應力區多出現在地質構造帶附近;由于隧洞埋深大,地質構造作用強烈,隧洞大部分區域處于高地應力區;中等地應力區主要分布于工程東側的自重應力區,初始地應力水平相對較低。

c)根據地應力分區結果,開展了引水隧洞典型斷面的圍巖穩定分析,斷裂帶附近的極高地應力區和分布范圍最廣的高地應力區段隧洞圍巖應力較高,變形較大,需注意應力釋放導致的巖爆發生,并及時支護,防止塑性區擴展和應力松弛引發的掉塊。中等地應力區段圍巖變形較小,巖體基本能夠實現自穩,圍巖穩定性相對較好。

參考文獻:

[1]周亞萍,姜海波.深埋高地應力水工隧洞地應力特征及巖爆脆性破壞深度預測研究[J].隧道建設(中英文),2022,42(S1):321-330.

[2]王琛濤.深埋隧洞高地應力巖體破壞型式研究及防治措施[J].水利規劃與設計,2022(7):81-86,103,138.

[3]王超,王益騰,韓增強,等.垂直孔應力解除法地應力測試技術及工程應用[J].巖土力學,2022,43(5):1412-1421.

[4]王成虎.地應力主要測試和估算方法回顧與展望[J].地質論評,2014,60(5):971-991,996,992-995.

[5]周朝,尹健民,董志宏,等.考慮邊界荷載作用方向的特長隧道初始應力場分區反演方法[J].巖石力學與工程學報,2022,41(S1):2725-2734.

[6]王金安,李飛.復雜地應力場反演優化算法及研究新進展[J].中國礦業大學學報,2015,44(2):189-205.

[7]趙雨,白金朋.基于FLAC3D的多元線性回歸法地下廠房初始地應力場反演重構[J].水電能源科學,2022,40(3):149-152,185.

[8]張建國,張強勇,楊文東,等.大崗山水電站壩區初始地應力場反演分析[J].巖土力學,2009,30(10):3071-3078.

[9]裴啟濤,李海波,劉亞群.南水北調西線工程壩區初始地應力場反演分析[J].巖土力學,2012,33(S2):338-344.

[10]戚藍,崔溦,熊開智,等.灰色理論在地應力場分析中的應用[J].巖石力學與工程學報,2002,21(10):1547-1550.

[11]裴書鋒,趙金帥,于懷昌,等.考慮洞室巖體應力型破壞特征的局部地應力反演方法及應用[J].巖土力學,2020,41(12):4093-4104.

[12]馬玉巖,沈陽,侯東奇.基于人工神經網絡和地層剝蝕原理的地應力場反演研究[J].水電站設計,2021,37(3):1-5,12.

[13]陳正林,何國志,張劼超,等.小樣本數據下三維地應力反演分析[J].科學技術與工程,2022,22(22):9822-9829.

[14]張金.遺傳算法在地應力反演中的應用[D].青島:中國石油大學,2010.

[15]楊志強,高謙,翟淑花,等.復雜工程地質體地應力場智能反演[J].哈爾濱工業大學學報,2016,48(4):154-160.

[16]謝學斌,羅海霞,楊承祥,等.基于遺傳單純形算法與RBF網絡的地應力場反演方法[J].鐵道科學與工程學報,2015,12(1):72-78.

[17]甘玉葉,鄭小燕.Hoek.Brown準則中參數不同選取值對Mohr.Coulomb參數值的影響[J].科學技術與工程,2012,12(21):5379-5383.

[18]王彬,榮傳新,程樺,等.基于DP準則雙排管凍結壁力學特性理論分析[J].科學技術與工程,2016,16(25):44-50,62.

[19]喬建剛,彭瑞,李景文,等.基于修正莫爾庫倫的基坑開挖對隧道安全影響研究[J].中國安全生產科學技術,2022,18(2):177-183.

[20]周鳳璽,邵彥平,甘東彪.基于廣義非線性統一強度理論的隧道塌落拱分析[J].應用力學學報,2020,37(2):682-688,935-936.

[21]黃海.九綿高速公路可溶巖隧道施工期圍巖位移控制基準研究[D].成都:成都理工大學,2020.

[22]鐘放平,張偉,宮鳳強,等.山嶺公路隧道圍巖穩定性的位移判別分析方法[J].公路工程,2008,33(6):91-94.

[23]李世煇,宋軍.變形速率比值判據與貓山隧道工程驗證[J].中國工程科學,2002(6):85-91.

[24]蘇士龍,高海海,周康樂.基于統一強度理論的巷道圍巖松動圈計算方法[J].科學技術與工程,2020,20(27):11045-11050.

[25]李美晨,趙啟峰,程志恒,等.巷道圍巖失穩垮冒風險等級分類與差異性分級管控[J].煤炭技術,2021,40(11):106-109.

[26]沈才華,古文博,李鶴文,等.基于損傷擴容理論的圓形隧洞圍巖松動圈位移計算方法[J].隧道建設(中英文),2019,39(1):40-47.

猜你喜歡
遺傳算法神經網絡
神經網絡抑制無線通信干擾探究
遺傳算法對CMAC與PID并行勵磁控制的優化
基于自適應遺傳算法的CSAMT一維反演
一種基于遺傳算法的聚類分析方法在DNA序列比較中的應用
基于遺傳算法和LS-SVM的財務危機預測
協同進化在遺傳算法中的應用研究
基于神經網絡的拉矯機控制模型建立
基于改進的遺傳算法的模糊聚類算法
復數神經網絡在基于WiFi的室內LBS應用
基于支持向量機回歸和RBF神經網絡的PID整定
91香蕉高清国产线观看免费-97夜夜澡人人爽人人喊a-99久久久无码国产精品9-国产亚洲日韩欧美综合