• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    不確定性條件下地下水污染監(jiān)測(cè)井網(wǎng)優(yōu)化設(shè)計(jì)——基于XGBoost替代模型

    2022-06-02 02:20:02董廣齊盧文喜潘紫東
    中國(guó)環(huán)境科學(xué) 2022年5期
    關(guān)鍵詞:污染監(jiān)測(cè)井網(wǎng)污染源

    董廣齊,盧文喜*,范 越,潘紫東

    不確定性條件下地下水污染監(jiān)測(cè)井網(wǎng)優(yōu)化設(shè)計(jì)——基于XGBoost替代模型

    董廣齊1,盧文喜1*,范 越2,潘紫東1

    (1.吉林大學(xué)新能源與環(huán)境學(xué)院,吉林 長(zhǎng)春 130012;2.長(zhǎng)江科學(xué)院,湖北 武漢 430010)

    在地下水污染監(jiān)測(cè)井網(wǎng)優(yōu)化設(shè)計(jì)中,應(yīng)用模擬優(yōu)化方法時(shí)客觀存在的模型參數(shù)不確定性往往會(huì)影響到設(shè)計(jì)監(jiān)測(cè)井網(wǎng)的可靠性.針對(duì)該問(wèn)題,重點(diǎn)考慮了滲透系數(shù)和污染源釋放強(qiáng)度的不確定性,應(yīng)用模擬優(yōu)化方法和蒙特卡羅方法求解上述不確定性參數(shù)影響下的最優(yōu)監(jiān)測(cè)井布設(shè)方案.為緩解蒙特卡羅方法多次調(diào)用模擬模型所產(chǎn)生的巨大計(jì)算負(fù)荷,本研究建立了XGBoost替代模型,代替模擬模型與優(yōu)化模型進(jìn)行耦合.為提高監(jiān)測(cè)井網(wǎng)對(duì)實(shí)際污染羽的監(jiān)測(cè)精度,污染監(jiān)測(cè)井網(wǎng)優(yōu)化模型以監(jiān)測(cè)空間矩誤差極小化為優(yōu)化目標(biāo).此外,本次研究還考慮了監(jiān)測(cè)井網(wǎng)設(shè)計(jì)中污染源釋放強(qiáng)度的動(dòng)態(tài)變化過(guò)程.最后,以撫順市某煤矸石堆放場(chǎng)地為基礎(chǔ)建立假想例子,驗(yàn)證所提方法的有效性.結(jié)果表明:1.XGBoost能夠有效近似模擬模型的輸入輸出關(guān)系,顯著降低了計(jì)算負(fù)荷.2.空間矩能夠有效評(píng)估監(jiān)測(cè)井網(wǎng)插值污染羽和實(shí)際污染羽的逼近程度,優(yōu)化設(shè)計(jì)后的監(jiān)測(cè)井網(wǎng)能夠較為準(zhǔn)確地捕捉到實(shí)際污染羽的狀態(tài).3.模擬優(yōu)化方法結(jié)合蒙特卡羅方法能有效求解不確定性條件下最優(yōu)監(jiān)測(cè)井網(wǎng)的設(shè)計(jì)問(wèn)題.本文為地下水污染監(jiān)測(cè)井網(wǎng)設(shè)計(jì)提供了一種穩(wěn)定可靠的方法.

    地下水污染;監(jiān)測(cè)井布設(shè);模擬優(yōu)化方法;不確定性;XGBoost替代模型

    地下水污染具有存在的隱蔽性、發(fā)現(xiàn)的滯后性、修復(fù)困難等特點(diǎn),在污染治理和修復(fù)的過(guò)程中建立相應(yīng)的地下水污染監(jiān)測(cè)網(wǎng)絡(luò),可以及時(shí)地獲取地下水污染資料,從而為污染修復(fù)工作提供數(shù)據(jù)支持和檢驗(yàn)保障.然而,由于地下水監(jiān)測(cè)井網(wǎng)的建設(shè)成本高昂,監(jiān)測(cè)工作的展開(kāi)需在有限的資金條件下獲得更多的地下水有效信息,因此對(duì)地下水污染監(jiān)測(cè)井網(wǎng)進(jìn)行優(yōu)化設(shè)計(jì)就顯得十分必要.其基本目的是在一定的人力,物力,資金限制下,獲取最具有代表性的地下水水質(zhì)時(shí)空分布信息,找到最經(jīng)濟(jì)有效的監(jiān)測(cè)井布設(shè)方案.

    模擬優(yōu)化方法將地下水?dāng)?shù)值模擬與運(yùn)籌學(xué)中的優(yōu)化方法相結(jié)合,使得其既能夠以地下水系統(tǒng)固有的物理規(guī)律為基礎(chǔ)來(lái)預(yù)測(cè)地下水位、水質(zhì)的分布趨勢(shì),又能夠在滿足各種環(huán)境、經(jīng)濟(jì)、技術(shù)等要求的前提下獲得最優(yōu)的井位布設(shè)方案[1].國(guó)內(nèi)外基于模擬優(yōu)化方法對(duì)監(jiān)測(cè)井網(wǎng)優(yōu)化設(shè)計(jì)的研究有很多.Meyer等首次考慮了不確定性條件下監(jiān)測(cè)井網(wǎng)的設(shè)計(jì),運(yùn)用模擬-優(yōu)化方法以最大化檢出概率為目標(biāo),確定了監(jiān)測(cè)井網(wǎng)的位置[2].Reed以提高監(jiān)測(cè)井網(wǎng)監(jiān)測(cè)精度為目的對(duì)優(yōu)化模型的目標(biāo)函數(shù)進(jìn)行改進(jìn),以全局質(zhì)量評(píng)估誤差最小,監(jiān)測(cè)成本最小為目標(biāo)函數(shù),結(jié)合反距離權(quán)重法和普通克里格法進(jìn)行插值,建立了地下水污染的長(zhǎng)期監(jiān)測(cè)網(wǎng)絡(luò)[3].Wu基于Reed的研究,在目標(biāo)函數(shù)中增加了污染羽一階矩陣和二階矩陣評(píng)估,提高監(jiān)測(cè)井網(wǎng)插值對(duì)實(shí)際污染羽的擬合精度,并建立了更為有效的監(jiān)測(cè)井網(wǎng)絡(luò)[4].駱乾坤整合前人在目標(biāo)函數(shù)上的研究成果,建立了包括監(jiān)測(cè)費(fèi)用評(píng)估、質(zhì)量評(píng)估和空間矩評(píng)估在內(nèi)的多目標(biāo)優(yōu)化模型,并用改進(jìn)小生境Pareto遺傳算法(NPGA)進(jìn)行求解,充分揭示了監(jiān)測(cè)費(fèi)用和監(jiān)測(cè)精度之間的權(quán)衡關(guān)系[5].范越以最大覆蓋高污染區(qū)域?yàn)楸O(jiān)測(cè)方案優(yōu)化目標(biāo)并應(yīng)用蒙特卡羅方法,在污染監(jiān)測(cè)井優(yōu)化設(shè)計(jì)中考慮了滲透系數(shù)不確定性帶來(lái)的影響[6].羅建男應(yīng)用模擬優(yōu)化法并結(jié)合0-1整數(shù)規(guī)劃模型,對(duì)某垃圾填埋場(chǎng)地下水質(zhì)監(jiān)測(cè)井進(jìn)行優(yōu)化設(shè)計(jì),取得了很好的結(jié)果[7].

    盡管前人在監(jiān)測(cè)井網(wǎng)的優(yōu)化設(shè)計(jì)上已取得了很大的進(jìn)展,在優(yōu)化模型的形式,優(yōu)化目標(biāo)的選取,優(yōu)化模型的求解等方面均有較為廣泛的研究,但長(zhǎng)期以來(lái)研究者對(duì)于不確定性的研究往往只局限于污染源信息或者模型參數(shù)的單一角度,缺少對(duì)于二者的綜合分析,針對(duì)該問(wèn)題本文同時(shí)對(duì)污染源釋放強(qiáng)度和模擬模型滲透系數(shù)的不確定性進(jìn)行了研究,運(yùn)用蒙特卡洛方法對(duì)兩者共同影響下的地下水污染監(jiān)測(cè)井網(wǎng)進(jìn)行優(yōu)化設(shè)計(jì).此外,過(guò)去的文獻(xiàn)一般僅針對(duì)釋放強(qiáng)度恒定的污染源,而對(duì)釋放強(qiáng)度隨時(shí)間變化的污染源少有涉及,因此本文對(duì)釋放強(qiáng)度非恒定的污染源情形進(jìn)行了研究.最后,在應(yīng)用蒙特卡羅方法時(shí),需要多次調(diào)用模擬模型,勢(shì)必會(huì)產(chǎn)生極大的計(jì)算負(fù)荷,而替代模型能夠以較小的計(jì)算量近似擬合模擬模型的輸入輸出關(guān)系,在蒙特卡洛分析中直接調(diào)用替代模型而不是模擬模型可以大幅減小計(jì)算負(fù)荷,減少計(jì)算時(shí)間.傳統(tǒng)的替代模型構(gòu)建方法(如Kriging、支持向量機(jī)、BP神經(jīng)網(wǎng)絡(luò)等[8-9])屬于淺層機(jī)器學(xué)習(xí),在簡(jiǎn)單問(wèn)題上擬合效果較好,但對(duì)于復(fù)雜條件下模擬模型存在的非線性映射關(guān)系擬合精度還有待提高.學(xué)習(xí)能力更強(qiáng)的深度神經(jīng)網(wǎng)絡(luò)方法(如卷積神經(jīng)網(wǎng)絡(luò)[10]、循環(huán)神經(jīng)網(wǎng)絡(luò)、深度置信神經(jīng)網(wǎng)絡(luò)等)的非線性表達(dá)能力更強(qiáng),但缺點(diǎn)是依賴數(shù)據(jù)驅(qū)動(dòng),需要建立較大的樣本集,同時(shí)模型設(shè)計(jì)復(fù)雜,不夠輕量化,容易出現(xiàn)過(guò)擬合現(xiàn)象.XGBoost(極限梯度提升樹(shù))是基于決策樹(shù)的集成機(jī)器學(xué)習(xí)算法,通過(guò)組合多個(gè)簡(jiǎn)單的弱分類器形成了一個(gè)學(xué)習(xí)能力更強(qiáng)的強(qiáng)分類器并且構(gòu)造形式簡(jiǎn)單,其次引入了葉節(jié)點(diǎn)懲罰函數(shù)有效避免了模型的過(guò)擬合,最后它對(duì)數(shù)據(jù)量的要求較低,在中小型數(shù)據(jù)集中也能有較好的表現(xiàn).目前為止尚未見(jiàn)到應(yīng)用XGBoost方法建立地下水?dāng)?shù)值模擬模型的替代模型的報(bào)道,亟需進(jìn)行實(shí)例驗(yàn)證.因此,本研究采用XGBoost方法擬合多維不確定變量下的模擬模型,研究其在多維輸入多維輸出情況下的擬合效果.

    綜上,本文基于模擬優(yōu)化方法建立了以最大化污染監(jiān)測(cè)精度為優(yōu)化目標(biāo)的地下水污染監(jiān)測(cè)井網(wǎng)優(yōu)化模型,采用蒙特卡羅方法研究滲透系數(shù)和污染源釋放強(qiáng)度不確定性對(duì)監(jiān)測(cè)井網(wǎng)優(yōu)化設(shè)計(jì)的影響,并同時(shí)考慮了污染源釋放強(qiáng)度的動(dòng)態(tài)變化過(guò)程.為減小求解過(guò)程中的計(jì)算負(fù)荷,應(yīng)用XGBoost方法構(gòu)建模擬模型的替代模型.最后以撫順市某煤矸石堆放場(chǎng)為基礎(chǔ)建立假想例子驗(yàn)證了方法的有效性.研究結(jié)果可為地下水污染監(jiān)測(cè)井網(wǎng)的布設(shè)提供方法依據(jù).

    1 模擬優(yōu)化方法

    模擬優(yōu)化方法的本質(zhì)是將模擬模型和優(yōu)化模型進(jìn)行耦合.模擬部分即建立地下水污染質(zhì)運(yùn)移模擬模型,反映地下水系統(tǒng)本身固有的物理規(guī)律;優(yōu)化部分即根據(jù)實(shí)際的優(yōu)化目標(biāo)構(gòu)建合適的目標(biāo)函數(shù),進(jìn)而建立優(yōu)化模型.本次研究采用嵌入法對(duì)兩者進(jìn)行耦合,將模擬模型(或其替代模型)作為約束條件嵌入優(yōu)化模型.而建立模擬模型是應(yīng)用模擬優(yōu)化方法的基礎(chǔ).

    1.1 地下水污染質(zhì)運(yùn)移模擬模型

    1.1.1 地下水水流模型的建立 根據(jù)研究區(qū)水文地質(zhì)概念模型,建立研究區(qū)二維潛水模型,如下所示:

    1.1.2 地下水污染質(zhì)運(yùn)移模型的建立

    在地下水水流數(shù)值模型基礎(chǔ)之上,建立二維地下水污染質(zhì)運(yùn)移數(shù)學(xué)模擬模型:

    1.2 模擬模型的XGBoost替代模型

    替代模型能夠以較小的計(jì)算量近似逼近模擬模型的輸入輸出關(guān)系,通過(guò)對(duì)一定數(shù)量的已知樣品特征的輸入-輸出關(guān)系進(jìn)行擬合(訓(xùn)練樣本),來(lái)構(gòu)建模擬模型的擬合函數(shù),用于預(yù)測(cè)未知樣品的特征輸出響應(yīng)[11],在應(yīng)用蒙特卡洛方法處理不確定性問(wèn)題時(shí),替代模型能夠大幅度降低計(jì)算負(fù)荷提高計(jì)算效率.

    XGBoost是基于決策樹(shù)的集成機(jī)器學(xué)習(xí)算法,在非結(jié)構(gòu)數(shù)據(jù)(圖像、文本等)的預(yù)測(cè)問(wèn)題中,人工神經(jīng)網(wǎng)絡(luò)數(shù)據(jù)的表現(xiàn)要優(yōu)于其他算法或框架,但在處理中小型結(jié)構(gòu)數(shù)據(jù)或表格數(shù)據(jù)中,基于決策樹(shù)的算法表現(xiàn)更為優(yōu)異[12].XGBoost近年來(lái)在機(jī)器學(xué)習(xí)領(lǐng)域表現(xiàn)出眾,是機(jī)器學(xué)習(xí)賽事Kaggle中大多數(shù)冠軍方案的主體算法,已經(jīng)應(yīng)用在多個(gè)工業(yè)前沿中[12-14].

    XGBoost[15]以CART(分類回歸樹(shù))為基礎(chǔ). CART通過(guò)建立迭代過(guò)程不斷地將特征空間通過(guò)二叉樹(shù)劃分為若干單元,每一個(gè)特征單元(葉節(jié)點(diǎn))對(duì)應(yīng)于一個(gè)輸出(一個(gè)簡(jiǎn)單的二層CART切分過(guò)程如圖1所示).通過(guò)不斷優(yōu)化特征空間的切分方式(切分特征和切分點(diǎn))和各特征單元的輸出值,使回歸樹(shù)輸出和已知輸出誤差最小,進(jìn)而對(duì)數(shù)據(jù)集進(jìn)行擬合.CART有一個(gè)較為顯著的缺點(diǎn),對(duì)于一個(gè)訓(xùn)練集特征空間易劃分的過(guò)于詳細(xì),導(dǎo)致過(guò)擬合的現(xiàn)象,泛化能力較差.

    圖1 二層CART切分示意

    XGBOOST算法彌補(bǔ)了CART的不足,采用集成思想,線性組合多個(gè)CART,將上一棵樹(shù)的擬合誤差作為下一棵樹(shù)的擬合目標(biāo),通過(guò)不斷加入新的回歸樹(shù)使預(yù)測(cè)值不斷逼近已知值.

    該過(guò)程如下所示:

    ●初始化(模型中沒(méi)有樹(shù)):

    ●在模型中插入第1棵樹(shù):

    ●在模型中插入第2棵樹(shù):

    ●在模型中插入第棵樹(shù):

    在明確了XGBoost組合多個(gè)CART的過(guò)程后,還需要進(jìn)一步確定模型中各個(gè)CART的構(gòu)建方法.XGBoost采用以下目標(biāo)函數(shù)對(duì)決策樹(shù)的切分特征、切分節(jié)點(diǎn)和各葉節(jié)點(diǎn)的輸出值進(jìn)行優(yōu)化.

    本文以5個(gè)不確定性變量作為輸入(各應(yīng)力期污染源釋放強(qiáng)度1、2、3,兩分區(qū)滲透系數(shù)1、2),以各潛在監(jiān)測(cè)井的污染質(zhì)濃度作為輸出,建立XGBoost替代模型,算法基于Python語(yǔ)言實(shí)現(xiàn).

    1.3 優(yōu)化模型

    地下水污染監(jiān)測(cè)井網(wǎng)的優(yōu)化設(shè)計(jì)問(wèn)題可以近似概化為運(yùn)籌學(xué)中經(jīng)典的離散選址問(wèn)題.在研究區(qū)中根據(jù)實(shí)際情況預(yù)先設(shè)置可以建設(shè)監(jiān)測(cè)井的潛在監(jiān)測(cè)井位置,而后從中挑選出信息量較大的監(jiān)測(cè)井,這可以應(yīng)用0-1整數(shù)規(guī)劃模型進(jìn)行表示,潛在井被選中取1,否則取0.本次研究的優(yōu)化目標(biāo)是最大化污染監(jiān)測(cè)信息量.如果設(shè)計(jì)監(jiān)測(cè)井網(wǎng)的監(jiān)測(cè)數(shù)據(jù)經(jīng)過(guò)克里格插值方法得到的插值污染羽和實(shí)際污染羽越相近,那么監(jiān)測(cè)井網(wǎng)的選址越好,地下水污染監(jiān)測(cè)信息量越大,否則監(jiān)測(cè)數(shù)據(jù)包含無(wú)效信息,應(yīng)進(jìn)行監(jiān)測(cè)井優(yōu)化,刪除冗余井,增加有效的控制井.但事實(shí)上,實(shí)際污染羽的形態(tài)很難得到,通常可以根據(jù)已有地質(zhì)和水文地質(zhì)信息建立數(shù)值模擬模型,求得模擬模型污染羽并用其代替實(shí)際污染羽.因此,優(yōu)化問(wèn)題轉(zhuǎn)化為怎樣甄選出信息量大的潛在監(jiān)測(cè)井,組合成合適的監(jiān)測(cè)井網(wǎng),使監(jiān)測(cè)井網(wǎng)插值污染羽和模擬模型污染羽更為相近.

    在進(jìn)行實(shí)際優(yōu)化求解過(guò)程中,需要將兩種污染羽的逼近程度進(jìn)行數(shù)學(xué)化表達(dá),建立評(píng)估監(jiān)測(cè)有效信息量的數(shù)學(xué)表達(dá)式.污染羽狀態(tài)以濃度為基礎(chǔ)可以處理成一個(gè)灰度圖像,兩種污染羽逼近程度的評(píng)估可以轉(zhuǎn)化為圖像的相似性評(píng)估.空間矩可以較為準(zhǔn)確地描述圖像的空間特征,它包括零階矩、一階矩和二階矩.零階矩可以代表污染羽的質(zhì)量,一階矩代表污染羽的重心位置,二階矩代表污染羽的形狀特征.只有同時(shí)知道這3種不同的物理特征,才能準(zhǔn)確地描述污染羽流的分布[5].

    污染羽的零階矩(全局質(zhì)量)、一階矩(重心位置、二階矩(污染羽形狀特征)公式如下所示:

    對(duì)于模型中滲透系數(shù)和釋放強(qiáng)度的不確定性問(wèn)題,應(yīng)用蒙特卡羅方法,通過(guò)多次運(yùn)轉(zhuǎn)替代模型將模型參數(shù)的不確定性轉(zhuǎn)化為了實(shí)際污染羽形態(tài)的不確定性.而后進(jìn)一步計(jì)算不同污染羽形態(tài)下的空間矩誤差,優(yōu)化模型的優(yōu)化目標(biāo)為各污染羽形態(tài)下空間矩誤差的總和最小,由此加強(qiáng)了監(jiān)測(cè)井網(wǎng)優(yōu)化結(jié)果的魯棒性,保證了其在不確定參數(shù)的整個(gè)取值范圍內(nèi)都能表現(xiàn)較好.在以上基礎(chǔ)上,建立污染監(jiān)測(cè)井網(wǎng)優(yōu)化模型,數(shù)學(xué)表達(dá)式如下:

    2 案例應(yīng)用

    2.1 問(wèn)題概述

    研究區(qū)域?yàn)閾犴樖心趁喉肥瘽B漏區(qū),東西長(zhǎng)約9700m,南北約6100m,總面積約為52.35km2.由于長(zhǎng)期進(jìn)行煤生產(chǎn)活動(dòng),大量堆放的煤矸石隨著降雨入滲向地下水中泄漏污染物,在當(dāng)?shù)匾鹆藝?yán)重的地下水污染,依據(jù)《地下水質(zhì)量標(biāo)準(zhǔn)》(GB/T 14848- 2017)[16],研究區(qū)地下水水質(zhì)已達(dá)V類.本文將硫酸根離子作為典型污染物,模擬研究區(qū)概況如圖2所示,本次研究的主要含水層為一層8m厚的潛水含水層,研究區(qū)根據(jù)水文地質(zhì)參數(shù)的非均質(zhì)性可劃分為兩個(gè)區(qū)域,概化為非均質(zhì)各項(xiàng)同性含水層,地下水流可概化為二維非穩(wěn)定流.煤矸石滲漏區(qū)概化為面狀污染質(zhì)源匯項(xiàng).研究區(qū)北部G1為渾河邊界,概化為給定水頭邊界;東北部邊界G2和南部邊界G4為隔水巖層,概化為隔水邊界;將東部邊界G3概化為側(cè)向徑流補(bǔ)給邊界,西南部邊界G5概化為側(cè)向徑流排泄邊界.在溶質(zhì)運(yùn)移模型中,根據(jù)實(shí)際情況以硫酸根作為模擬污染質(zhì),研究區(qū)北部邊界G1和西南部邊界G5概化已知對(duì)流-彌散通量邊界;東北部邊界G2和南部邊界G4概化為零通量的水動(dòng)力彌散邊界;東部邊界G3概化為已知濃度邊界,地下水總體由東南向西北方向流動(dòng).

    圖2 研究區(qū)概況

    在數(shù)值模型建立后,采用2017年10月的水位、水質(zhì)數(shù)據(jù)進(jìn)行模型校正,采用2018年10月的水位、水質(zhì)數(shù)據(jù)進(jìn)行模型檢驗(yàn).在模型的校正和檢驗(yàn)過(guò)程中,水頭的觀測(cè)值和模擬值誤差小于0.5m,污染質(zhì)濃度觀測(cè)值和模擬值的相對(duì)誤差小于10%,模型校正和檢驗(yàn)的準(zhǔn)確度均滿足要求,兩分區(qū)滲透系數(shù)校正結(jié)果分別為82.7m/d和97.4m/d.此外,根據(jù)污染物的釋放時(shí)間和觀測(cè)濃度,運(yùn)用溶質(zhì)運(yùn)移模型進(jìn)行反演識(shí)別,獲得煤矸石向地下水泄漏硫酸根的釋放強(qiáng)度大約為7000mg/d.研究區(qū)各水文地質(zhì)參數(shù)校正值見(jiàn)表1.

    在以上實(shí)際場(chǎng)地的基礎(chǔ)上建立假想例子.相關(guān)部門對(duì)煤矸石堆污染物進(jìn)行階梯性修復(fù),10內(nèi)其污染源釋放強(qiáng)度逐漸降低,預(yù)計(jì)釋放強(qiáng)度為1~ 47000mg/d,5~74000mg/d,8~101000mg/d,本次研究為對(duì)10后的修復(fù)效果檢驗(yàn)提供數(shù)據(jù)支持,設(shè)計(jì)監(jiān)測(cè)井網(wǎng)擬合10后地下水污染羽狀態(tài).

    表1 模型校正后研究區(qū)水文地質(zhì)參數(shù)表

    煤矸石堆污染源的釋放強(qiáng)度受多種因素影響,如降水強(qiáng)度、降水持續(xù)時(shí)間、人類活動(dòng)等,滲透系數(shù)在空間上同樣具有一定的變異性,因此污染物釋放強(qiáng)度和滲透系數(shù)均具有一定的不確定性,這兩者都會(huì)對(duì)地下水污染井網(wǎng)設(shè)計(jì)產(chǎn)生較大影響.

    因此,本文同時(shí)考慮污染源3個(gè)應(yīng)力期釋放強(qiáng)度和2個(gè)分區(qū)滲透系數(shù)的不確定性,采用了蒙特卡洛方法對(duì)10后的地下水污染監(jiān)測(cè)井網(wǎng)進(jìn)行設(shè)計(jì),根據(jù)實(shí)際場(chǎng)地條件確定50口潛在監(jiān)測(cè)井,其在煤矸石滲漏區(qū)的地下水下游區(qū)域分布較密,在上游區(qū)域分布較稀疏,從中選取合適的監(jiān)測(cè)井構(gòu)建地下水污染監(jiān)測(cè)井網(wǎng),最大化地下水污染監(jiān)測(cè)精度.

    不確定性變量大致變化范圍依據(jù)專家經(jīng)驗(yàn)如表2所示.

    表2 不確定參數(shù)分布范圍

    在本案例中,研究區(qū)污染物濃度場(chǎng)使用GMS中MODFLOW和MT3DMS模塊進(jìn)行計(jì)算.

    2.2 XGBoost替代模型的建立

    本文在研究區(qū)內(nèi)設(shè)置50口潛在監(jiān)測(cè)井,從中選取合適的監(jiān)測(cè)井構(gòu)建地下水污染監(jiān)測(cè)井網(wǎng).因此替代模型的輸入向量為污染源3個(gè)應(yīng)力期的釋放強(qiáng)度以及兩個(gè)分區(qū)的滲透系數(shù)(1、2、3、1、2),輸出向量為50口潛在監(jiān)測(cè)井處的硫酸根濃度.

    本文采用R2(決定系數(shù))和MRE(平均相對(duì)誤差)兩個(gè)指標(biāo)衡量替代模型的擬合效果,其中MRE的數(shù)學(xué)表達(dá)式如下:

    式中:為樣本組數(shù);y為模擬模型輸出值;y為替代模型輸出值.

    表3 5折交叉驗(yàn)證中5次訓(xùn)練的評(píng)價(jià)結(jié)果

    為了更好的評(píng)估模型的擬合能力,減少隨機(jī)誤差帶來(lái)的影響,提高樣本信息利用率,采用折交叉驗(yàn)證方法對(duì)模型進(jìn)行交叉驗(yàn)證,該方法將樣本集平均分割為個(gè)互斥子集,每次選擇(-1)個(gè)子集作為訓(xùn)練集,1個(gè)子集作為測(cè)試集;重復(fù)多輪計(jì)算,保證每個(gè)子集都被且僅被作為一次測(cè)試集,最終獲得個(gè)測(cè)試結(jié)果,并輸出其平均值作為擬合精度[18].在本次研究中,取折數(shù)為5,5次訓(xùn)練后所得評(píng)價(jià)指標(biāo)及平均值見(jiàn)表3.其中Fold-4表現(xiàn)最好,該次訓(xùn)練替代模型和模擬模型的擬合結(jié)果如圖3所示.5次訓(xùn)練2和MRE的平均值分別為0.995和0.87%,說(shuō)明XGBoost替代模型的擬合精度較高.

    圖3 替代模型與模擬模型的擬合結(jié)果

    2.3 優(yōu)化模型的建立和求解

    在對(duì)XGBoost模型進(jìn)行交叉驗(yàn)證后,選取訓(xùn)練表現(xiàn)最好的一次作為模擬模型的替代模型.而后運(yùn)用拉丁超立方抽樣(LHS)在不確定變量分布范圍內(nèi)抽取1000組輸入樣本,將其帶入到替代模型而不是模擬模型中進(jìn)行求解,得到各潛在監(jiān)測(cè)井的硫酸根濃度,從而獲得不確定條件下的1000組污染羽.

    根據(jù)上文中優(yōu)化模型構(gòu)建方法,建立最大化污染監(jiān)測(cè)精度的監(jiān)測(cè)井網(wǎng)優(yōu)化模型,優(yōu)化目標(biāo)為最小化1000組污染羽的空間矩誤差總和,優(yōu)化模型具體的數(shù)學(xué)形式如下:

    式中:第一個(gè)約束條件是模擬模型的替代模型,該約束條件使整個(gè)優(yōu)化模型滿足地下水系統(tǒng)本身固有的物理規(guī)律;第二個(gè)約束條件為監(jiān)測(cè)井的限制個(gè)數(shù);第三個(gè)約束條件表示決策變量是0-1整數(shù)變量.

    圖4 不同監(jiān)測(cè)井總數(shù)下的最優(yōu)布設(shè)方案

    該優(yōu)化模型運(yùn)用MATLAB平臺(tái)的遺傳算法工具箱進(jìn)行求解,獲得了不同監(jiān)測(cè)井總數(shù)下的最優(yōu)布設(shè)方案,如圖4所示.此外,運(yùn)行GMS求解地下水溶質(zhì)運(yùn)移模型平均單次耗時(shí)約為2min,如直接使用模擬模型進(jìn)行不確定性分析,需要計(jì)算模擬模型1000次,總計(jì)耗時(shí)33h.而運(yùn)用替代模型代替模擬模型則只需運(yùn)行模擬模型100次獲取輸入-輸出樣本集,總耗時(shí)3.3h,替代模型的輸出約需2~3s,可忽略不計(jì).因此,利用替代模型進(jìn)行1000次蒙特卡洛實(shí)驗(yàn),可節(jié)省90%的計(jì)算時(shí)間.

    2.4 優(yōu)化結(jié)果的檢驗(yàn)

    本模型的目的是在滲透系數(shù)和污染源釋放強(qiáng)度不確定的條件下最大化地下水污染監(jiān)測(cè)精度,最后應(yīng)當(dāng)對(duì)優(yōu)化結(jié)果進(jìn)行檢驗(yàn).

    選取各輸入變量的標(biāo)準(zhǔn)值作為輸入對(duì)優(yōu)化結(jié)果進(jìn)行檢驗(yàn),即污染源釋放強(qiáng)度為7000(1~4)、4000(5~7)、1000(8~10)mg/d,兩個(gè)分區(qū)的滲透系數(shù)分別為82.7m/d(區(qū)1)和97.4m/d(區(qū)2).其余水文地質(zhì)條件不變,運(yùn)用GMS進(jìn)行求解,獲得10后的模擬污染羽.而后將模擬污染羽與監(jiān)測(cè)井布設(shè)方案得到的插值污染羽進(jìn)行相似性評(píng)估,評(píng)估指標(biāo)采用零階矩、一階矩、二階矩,兩污染羽空間矩的相對(duì)誤差如表4所示,模擬污染羽與監(jiān)測(cè)井網(wǎng)插值污染羽的對(duì)比如圖5所示.

    可以看出,隨著監(jiān)測(cè)井?dāng)?shù)量的不斷增加,插值污染羽與模擬污染羽的誤差評(píng)估指標(biāo)都在不斷降低,其中二階矩敏感性最強(qiáng),下降最為明顯.當(dāng)監(jiān)測(cè)井?dāng)?shù)量為18口時(shí),全局質(zhì)量誤差(m)為2.27%,質(zhì)心相對(duì)誤差(1)為1.07%,二階矩相對(duì)誤差(E2)為10.91%,擬合效果較好,結(jié)合污染羽對(duì)比圖(圖5)也可以看出,監(jiān)測(cè)井網(wǎng)設(shè)計(jì)方案的插值污染羽能夠相對(duì)精確地刻畫實(shí)際污染羽的濃度、位置、形狀.在污染質(zhì)濃度大的區(qū)域,監(jiān)測(cè)井布設(shè)密集,反之稀疏,且隨著監(jiān)測(cè)井?dāng)?shù)目的增加,刻畫精度越來(lái)越高.本次研究為污染監(jiān)測(cè)井網(wǎng)的優(yōu)化設(shè)計(jì)提供了一種穩(wěn)定可靠的方法.

    表4 各個(gè)監(jiān)測(cè)方案模擬污染羽與插值污染羽的空間矩相對(duì)誤差

    圖5 不同監(jiān)測(cè)方案下模擬污染羽和插值污染羽的對(duì)比

    3 結(jié)論

    3.1 運(yùn)用XGBoost方法建立模擬模型的替代模型精度較高,能有效近似模擬模型的輸入輸出關(guān)系,顯著降低了計(jì)算負(fù)荷.因此,在運(yùn)用蒙特卡羅方法分析不確定性時(shí),使用XGBoost替代模型代替模擬模型與優(yōu)化模型進(jìn)行耦合是可行的.

    3.2 空間矩圖像評(píng)估方法能有效評(píng)估監(jiān)測(cè)井網(wǎng)插值污染羽和實(shí)際污染羽的逼近程度,在以空間矩誤差極小化為目標(biāo)函數(shù)進(jìn)行監(jiān)測(cè)井網(wǎng)優(yōu)化設(shè)計(jì)時(shí),最優(yōu)監(jiān)測(cè)方案能夠較為準(zhǔn)確地捕捉到實(shí)際污染羽的形態(tài).

    3.3 將0-1整數(shù)規(guī)劃和模擬模型相結(jié)合的模擬優(yōu)化方法,能夠從一系列潛在監(jiān)測(cè)井中選擇合適位置構(gòu)建監(jiān)測(cè)井網(wǎng).將蒙特卡羅方法和模擬優(yōu)化方法相結(jié)合,可以用來(lái)求解不確定性條件下污染監(jiān)測(cè)井網(wǎng)的最優(yōu)布設(shè)方案.

    [1] 盧文喜.地下水系統(tǒng)的模擬預(yù)測(cè)和優(yōu)化管理[M]. 北京:科學(xué)出版社, 1999:104-105.

    Lu W X. Groundwater system simulation prediction and optimization management [M]. Beijing: Science Press, 1999:104-105.

    [2] Meyer P D, Brill Jr E D. A method for locating wells in a groundwater monitoring network under conditions of uncertainty [J]. Water Resources Research, 1988,24(8):1277-1282.

    [3] Reed P, Barbara M, Albert J V. Cost-effective long-term groundwater monitoring design using a genetic algorithm and global mass interpolation [J]. Water Resources Research, 2000,36(12):3731-3741.

    [4] Wu J F, Zheng C M. Cost-effective sampling network design for contaminant plume monitoring under general hydrogeological conditions [J]. Journal of Contaminant Hydrology, 2005,77:41-65.

    [5] 駱乾坤,吳劍鋒,林 錦,等.地下水污染監(jiān)測(cè)網(wǎng)多目標(biāo)優(yōu)化設(shè)計(jì)模型及進(jìn)化求解[J]. 水文地質(zhì)工程地質(zhì), 2013,40(5):97-103.

    Luo Q K, Wu J F, Lin J, et al. An evolutionary-based multi-objective optimization model for groundwater monitoring network design [J].Hydrogeology & Engineering Geology, 2013,40(5):97-103.

    [6] 范 越,盧文喜,歐陽(yáng)琦,等.基于Kriging替代模型的地下水污染監(jiān)測(cè)井網(wǎng)優(yōu)化設(shè)計(jì)[J]. 中國(guó)環(huán)境科學(xué), 2017,37(10):3800-3806.

    Fan Y, Lu W X, Ouyang Q, et al. Optimum design of groundwater pollution monitoring well network based on Kriging surrogate model [J]. China Environmental Science, 2017,37(10):3800-3806.

    [7] 羅建男,李多強(qiáng),范 越,等.某垃圾填埋場(chǎng)地下水質(zhì)監(jiān)測(cè)井網(wǎng)優(yōu)化設(shè)計(jì)——基于模擬優(yōu)化法[J]. 中國(guó)環(huán)境科學(xué), 2019,39(1):196-202.

    Luo J N, Li D Q, Fan Y, et al. Optimization of groundwater quality monitoring network at a landfill——based on simulation optimization method [J]. China Environmental Science, 2019,39(1):196-202.

    [8] 邢貞相,曲睿卓,趙 瑩,等.不同替代模型在地下水污染源釋放歷史反演中適用性研究[J]. 東北農(nóng)業(yè)大學(xué)學(xué)報(bào), 2018,49(12):59-68.

    Xing Z X, Qu R Z, Zhao Y, et al. Applicability study of different surrogate models in retrieving release history of groundwater contaminant sources [J]. Journal of Northeast Agricultural University, 2018,49(12):59-68.

    [9] 潘紫東,盧文喜,范 越,等.基于模擬-優(yōu)化方法的地下水污染源溯源辨識(shí)[J]. 中國(guó)環(huán)境科學(xué),2020,40(4):1698-1705.

    Pan Z D, Lu W X, Fan Y, et al. Inverse identification of groundwater pollution source based on simulation-optimization approach [J]. China Environmental Science, 2020,40(4):1698-1705.

    [10] 莫紹星.基于深度學(xué)習(xí)的地下水模擬高維不確定性分析和反演[D]. 南京:南京大學(xué), 2019.

    Mo SX. Towards Efficient high-dimensional uncertainty quantification and inverse analysis in groundwater modeling using deep learning [D]. Nanjing: Nanjing University, 2019.

    [11] 侯澤宇,盧文喜,王 宇.基于替代模型的地下水DNAPLs污染源反演識(shí)別[J]. 中國(guó)環(huán)境科學(xué), 2019,39(1):188-195.

    Hou ZY, Lu WX, Wang Y. Surrogate-based source identification of DNAPLs-contaminated groundwater [J]. China Environmental Science, 2019,39(1):188-195.

    [12] 何 龍.深入理解XGBoost:高效機(jī)器學(xué)習(xí)算法與進(jìn)階[M]. 北京:機(jī)械工業(yè)出版社, 2020:6-8.

    He L. Deep understanding of XGBoost: efficient machine learning algorithms and advanced problem [M]. Beijing: China Machine Press, 2020:6-8.

    [13] 陳振宇,劉金波,李 晨,等.基于LSTM與XGBoost組合模型的超短期電力負(fù)荷預(yù)測(cè)[J]. 電網(wǎng)技術(shù), 2020,44(2):614-620.

    Chen Z Y, Liu J B, Li C, et al. Ultra Short-term Power Load Forecasting Based on Combined LSTM-XGBoost Model [J]. Power System Technology, 2020,44(2):614-620.

    [14] 趙洪山,閆西慧,王桂蘭,等.應(yīng)用深度自編碼網(wǎng)絡(luò)和XGBoost的風(fēng)電機(jī)組發(fā)電機(jī)故障診斷[J]. 電力系統(tǒng)自動(dòng)化, 2019,43(1):81-86.

    Zhao H S, Yan X H, Wang G L, et al. Fault diagnosis of wind turbine generator based on deep autoencoder network and XGBoost [J].Automation of Electric Power Systems, 2019,43(1):81-86.

    [15] Chen T, Guestrin C. XGBoost: a scalable tree boosting system [C]// Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining. San Francisco, California, USA: Association for Computing Machinery, 2016:785–94.

    [16] GB/T 14848-2017 地下水質(zhì)量標(biāo)準(zhǔn) [S].

    GB/T 14848-2017 Standard for groundwater quality [S].

    [17] Jin Y. A comprehensive survey of fitness approximation in evolutionary computation [J]. Soft Computing, 2005,9(1):3-12.

    [18] 張鈞博,何 川,嚴(yán) 健,等.基于交叉驗(yàn)證的XGBoost算法在巖爆烈度分級(jí)預(yù)測(cè)中的適用性探討[J]. 隧道建設(shè), 2020,40(S1):247-253.

    Zhang J B, He C, Yan J, et al. Discussion on the applicability of XGBoost algorithm based on cross validation in prediction of rockburst intensity classification [J]. Tunnel Construction, 2020,40 (S1):247-253.

    [19] Fan Y, Lu W X. Optimal design of groundwater pollution monitoring network based on the SVR surrogate model under uncertainty [J]. Environmental Science and Pollution Research, 2020,27:24090–24102.

    [20] 黃滿紅,顧國(guó)維.地下水環(huán)境監(jiān)測(cè)系統(tǒng)的設(shè)計(jì)[J]. 環(huán)境監(jiān)測(cè)管理與技術(shù), 2003,15(1):13-15.

    Huang M H, Gu G W. The design about environmental monitoring system of underground water [J].The Administration and Technique of Environmental Monitoring, 2003,15(1):13-15.

    [21] Philip D M, Albert J V, Eheart W. Monitoring network design to provide initial detection of groundwater contamination [J]. Water Resources Research, 1994,30(9):2647-2659.

    Optimal design of groundwater pollution monitoring network under uncertainty.

    DONG Guang-qi1, LU Wen-xi1*, FAN Yue2, PAN Zi-Dong1

    (1.College of New Energy and Environment, Jilin University, Changchun 130012, China;2.Yangtze River Scientific Research Institute, Wuhan 430010, China)., 2022,42(5):2144~2152

    When applying the simulation-optimization method, objective parameter uncertainty will usually affect the reliability of the design result of groundwater pollution monitoring network. For this problem, the study simultaneously considered the uncertainty of hydraulic conductivity and emission intensity of pollution source, applying Monte Carlo method to design the optimal monitoring wells scheme under the influence of model uncertainty. But Monte Carlo method need to invoke the simulation model many times which will cause a huge amount of calculation. To reduce the calculation load, the study proposed to use Extreme Gradient Boosting (XGBoost) method to construct the surrogate model replacing the simulation model to couple the optimization model in the optimal design of GPMN. To sufficiently improve the monitoring precision of GPMN, the optimization model applied error of spatial moment as objective function. Besides, the dynamic change of emission intensity of pollution source was also considered. Finally, we proposed a hypothetical example based on a coal gangue pile in Fushun City to verify the validity of the method. The results are demonstrated: 1.the XGBoost surrogate method can fit the input-output relationship of the simulation model to a high degree with less computation. 2.spatial moment can effectively assess the approximation degree between interpolation pollution plume of GPMN and actual pollution plume, through which the optimized monitoring network can accurately depict actual pollution plume 3.the simulation-optimization method combines Monte Carlo method can solve the problem of the design of GPMN under uncertainty. In conclusion, the paper provides a stable and reliable method for the design of GPMN.

    groundwater pollution;monitoring wells network design;simulation-optimization method;uncertainty;XGBoost surrogate model

    X523,X830

    A

    1000-6923(2022)05-2144-09

    董廣齊(1998-),男,河南信陽(yáng)人,吉林大學(xué)碩士研究生,主要從事地下水污染監(jiān)測(cè)井網(wǎng)優(yōu)化設(shè)計(jì)研究.

    2022-10-27

    國(guó)家自然科學(xué)基金資助項(xiàng)目(41972252);國(guó)家科技部重點(diǎn)研發(fā)計(jì)劃項(xiàng)目(2018YFC18000405)

    * 責(zé)任作者, 教授, lwx999@163.com

    猜你喜歡
    污染監(jiān)測(cè)井網(wǎng)污染源
    大氣環(huán)境污染監(jiān)測(cè)與環(huán)境保護(hù)對(duì)策探究
    持續(xù)推進(jìn)固定污染源排污許可管理全覆蓋
    全球污染監(jiān)測(cè)站搜尋隱秘殺手
    超低滲透油藏水平井注采井網(wǎng)設(shè)計(jì)優(yōu)化研究
    青海省人民政府辦公廳轉(zhuǎn)發(fā)省環(huán)境保護(hù)廳省氣象局關(guān)于青海省大氣污染監(jiān)測(cè)預(yù)報(bào)預(yù)警工作方案的通知
    基于污染源解析的空氣污染治理對(duì)策研究
    十二五”期間佳木斯市污染源排放狀況分析
    看不見(jiàn)的污染源——臭氧
    各向異性油藏菱形反九點(diǎn)井網(wǎng)合理井排距研究
    現(xiàn)代工業(yè)技術(shù)在農(nóng)田重金屬污染監(jiān)測(cè)中的應(yīng)用研究
    国产av一区在线观看免费| 色吧在线观看| 3wmmmm亚洲av在线观看| 亚洲片人在线观看| 久久久久久久久久黄片| 99riav亚洲国产免费| 禁无遮挡网站| 操出白浆在线播放| 日日夜夜操网爽| 国产精品亚洲一级av第二区| 一区二区三区高清视频在线| 国产真实伦视频高清在线观看 | 国产 一区 欧美 日韩| 成人永久免费在线观看视频| 欧美日本亚洲视频在线播放| 午夜福利免费观看在线| 国产美女午夜福利| 小蜜桃在线观看免费完整版高清| 麻豆久久精品国产亚洲av| 色综合亚洲欧美另类图片| 午夜亚洲福利在线播放| www国产在线视频色| 一区二区三区国产精品乱码| 国产高清视频在线播放一区| 99在线人妻在线中文字幕| 99久久无色码亚洲精品果冻| 内地一区二区视频在线| 久久伊人香网站| 熟女少妇亚洲综合色aaa.| 少妇高潮的动态图| 国产探花在线观看一区二区| 亚洲人成网站在线播| 久久性视频一级片| 在线观看日韩欧美| 国内精品久久久久久久电影| 狂野欧美激情性xxxx| 内地一区二区视频在线| 日本一本二区三区精品| 少妇的逼水好多| 国内精品一区二区在线观看| 亚洲中文字幕日韩| 大型黄色视频在线免费观看| 久久精品人妻少妇| 欧美性猛交黑人性爽| 国产淫片久久久久久久久 | 国产成人欧美在线观看| 欧美日韩瑟瑟在线播放| 国产精品av视频在线免费观看| 国产精品爽爽va在线观看网站| tocl精华| 看免费av毛片| 中文字幕av成人在线电影| 国产欧美日韩精品亚洲av| 高清毛片免费观看视频网站| 九色国产91popny在线| 小说图片视频综合网站| 美女黄网站色视频| 欧美乱妇无乱码| 国产又黄又爽又无遮挡在线| 亚洲一区二区三区色噜噜| 亚洲,欧美精品.| 99riav亚洲国产免费| av黄色大香蕉| 精品一区二区三区人妻视频| 国产一区二区激情短视频| 麻豆国产av国片精品| 国产高清三级在线| 两个人的视频大全免费| 淫秽高清视频在线观看| www国产在线视频色| 搞女人的毛片| 热99在线观看视频| 成人无遮挡网站| 国内久久婷婷六月综合欲色啪| 成人国产综合亚洲| 日韩欧美一区二区三区在线观看| 一区二区三区免费毛片| 国产 一区 欧美 日韩| 99视频精品全部免费 在线| 久久婷婷人人爽人人干人人爱| 啪啪无遮挡十八禁网站| 观看免费一级毛片| 精品电影一区二区在线| 极品教师在线免费播放| 国产精品久久电影中文字幕| 国产探花极品一区二区| 美女大奶头视频| 两个人视频免费观看高清| 色综合欧美亚洲国产小说| 熟妇人妻久久中文字幕3abv| 国产伦精品一区二区三区视频9 | 久久精品国产亚洲av涩爱 | 国产成人系列免费观看| 精品人妻1区二区| 久久婷婷人人爽人人干人人爱| 搡女人真爽免费视频火全软件 | 亚洲国产欧美网| 黄色日韩在线| 国产乱人视频| 国产高清激情床上av| 少妇的逼好多水| 麻豆成人午夜福利视频| 国产亚洲精品av在线| 亚洲真实伦在线观看| 俄罗斯特黄特色一大片| 亚洲avbb在线观看| 婷婷精品国产亚洲av| 九九久久精品国产亚洲av麻豆| 中国美女看黄片| 中文字幕人妻熟人妻熟丝袜美 | 久99久视频精品免费| 国产伦精品一区二区三区四那| 在线观看av片永久免费下载| 国产伦人伦偷精品视频| 国产精品久久久久久久久免 | 免费av不卡在线播放| 精品国产美女av久久久久小说| 18禁裸乳无遮挡免费网站照片| 中文在线观看免费www的网站| 欧美区成人在线视频| 美女高潮喷水抽搐中文字幕| 成人性生交大片免费视频hd| 国产探花极品一区二区| 97超视频在线观看视频| 午夜精品久久久久久毛片777| 日日摸夜夜添夜夜添小说| 嫩草影视91久久| 天天添夜夜摸| 99热这里只有精品一区| 一级黄色大片毛片| 久久久久久大精品| 久久久国产成人精品二区| 免费看光身美女| 国产精品av视频在线免费观看| 99riav亚洲国产免费| 精品一区二区三区av网在线观看| ponron亚洲| 国产一区二区三区在线臀色熟女| 亚洲精品乱码久久久v下载方式 | e午夜精品久久久久久久| 麻豆一二三区av精品| 国产精品1区2区在线观看.| 在线a可以看的网站| 少妇人妻一区二区三区视频| 女生性感内裤真人,穿戴方法视频| 好男人在线观看高清免费视频| 精品日产1卡2卡| 成人永久免费在线观看视频| 久久久久久九九精品二区国产| 亚洲国产欧洲综合997久久,| 18禁美女被吸乳视频| 国产老妇女一区| 美女免费视频网站| 欧美午夜高清在线| 国产精品女同一区二区软件 | 老熟妇乱子伦视频在线观看| 日本黄大片高清| 亚洲自拍偷在线| 天堂动漫精品| 欧美zozozo另类| 国产精品嫩草影院av在线观看 | 日韩欧美精品免费久久 | 两个人看的免费小视频| 在线观看66精品国产| 十八禁网站免费在线| 1024手机看黄色片| 欧美日韩瑟瑟在线播放| 19禁男女啪啪无遮挡网站| 夜夜夜夜夜久久久久| 亚洲最大成人中文| 草草在线视频免费看| 久久久色成人| 日韩国内少妇激情av| www.色视频.com| tocl精华| 国产午夜精品论理片| 亚洲国产精品999在线| 国产成年人精品一区二区| 夜夜爽天天搞| tocl精华| 久久草成人影院| xxx96com| 搞女人的毛片| 国产精品久久久人人做人人爽| 亚洲五月婷婷丁香| 欧美绝顶高潮抽搐喷水| 亚洲av美国av| 欧美一级毛片孕妇| 91久久精品电影网| 一区二区三区免费毛片| 国产av麻豆久久久久久久| 国产真实伦视频高清在线观看 | 18禁在线播放成人免费| 国产成人av激情在线播放| 午夜福利在线在线| 日本一二三区视频观看| 国产精品野战在线观看| 国产午夜精品论理片| 日本黄大片高清| 露出奶头的视频| 国产精品久久久久久久久免 | 久久久精品大字幕| 99国产极品粉嫩在线观看| 两个人视频免费观看高清| 国产在线精品亚洲第一网站| 丁香欧美五月| 男人和女人高潮做爰伦理| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 久久久久久人人人人人| 母亲3免费完整高清在线观看| 精品国内亚洲2022精品成人| 成人三级黄色视频| 国产精华一区二区三区| 色老头精品视频在线观看| 88av欧美| 成年女人毛片免费观看观看9| 在线观看免费视频日本深夜| 国产伦人伦偷精品视频| 午夜视频国产福利| 日韩欧美三级三区| 久久久国产成人精品二区| 一区二区三区激情视频| 丝袜美腿在线中文| 久久久久免费精品人妻一区二区| 免费高清视频大片| 精品久久久久久久久久久久久| 久久欧美精品欧美久久欧美| 午夜免费激情av| 91九色精品人成在线观看| 他把我摸到了高潮在线观看| 亚洲一区高清亚洲精品| 日本免费一区二区三区高清不卡| www.www免费av| 99久久99久久久精品蜜桃| 男女做爰动态图高潮gif福利片| 一级作爱视频免费观看| 在线观看免费视频日本深夜| 岛国在线免费视频观看| 91九色精品人成在线观看| 国产私拍福利视频在线观看| 欧美成人免费av一区二区三区| 九九在线视频观看精品| 亚洲成人中文字幕在线播放| 日本 av在线| 毛片女人毛片| 一个人看的www免费观看视频| 亚洲精品国产精品久久久不卡| 内地一区二区视频在线| 久久国产乱子伦精品免费另类| 亚洲成av人片免费观看| 变态另类成人亚洲欧美熟女| 亚洲av第一区精品v没综合| 日韩大尺度精品在线看网址| 国产精品久久久人人做人人爽| 免费在线观看成人毛片| 老熟妇仑乱视频hdxx| 老鸭窝网址在线观看| 午夜日韩欧美国产| 国产av不卡久久| 人妻久久中文字幕网| 日本黄大片高清| 国产精品日韩av在线免费观看| 操出白浆在线播放| 在线播放无遮挡| 97碰自拍视频| 色综合站精品国产| 国产麻豆成人av免费视频| av片东京热男人的天堂| 夜夜躁狠狠躁天天躁| av在线蜜桃| 丁香欧美五月| 九九热线精品视视频播放| 久久久国产成人精品二区| 国产毛片a区久久久久| 国产真人三级小视频在线观看| 老汉色av国产亚洲站长工具| 可以在线观看的亚洲视频| 99热精品在线国产| 色哟哟哟哟哟哟| 99精品欧美一区二区三区四区| 亚洲av成人av| 99在线视频只有这里精品首页| 国产又黄又爽又无遮挡在线| 免费在线观看影片大全网站| 国产精品一区二区三区四区久久| 在线观看免费午夜福利视频| 亚洲中文日韩欧美视频| 午夜福利高清视频| 国产91精品成人一区二区三区| 日韩亚洲欧美综合| 小蜜桃在线观看免费完整版高清| 午夜福利18| 搡老妇女老女人老熟妇| 黑人欧美特级aaaaaa片| 国产精品一及| 欧美成人免费av一区二区三区| 97碰自拍视频| 好男人在线观看高清免费视频| 黄色日韩在线| 啪啪无遮挡十八禁网站| 欧美激情久久久久久爽电影| 超碰av人人做人人爽久久 | 亚洲成人久久性| 国模一区二区三区四区视频| 内地一区二区视频在线| netflix在线观看网站| 一级a爱片免费观看的视频| 男女视频在线观看网站免费| 国产成人av激情在线播放| 两个人的视频大全免费| 长腿黑丝高跟| 在线免费观看不下载黄p国产 | 久久天躁狠狠躁夜夜2o2o| 天堂影院成人在线观看| 精品一区二区三区人妻视频| 噜噜噜噜噜久久久久久91| 免费看日本二区| 在线a可以看的网站| 99热只有精品国产| 99热6这里只有精品| 又黄又粗又硬又大视频| 日本 av在线| tocl精华| 搞女人的毛片| 美女黄网站色视频| 亚洲成人久久爱视频| 两性午夜刺激爽爽歪歪视频在线观看| 国产一区二区在线观看日韩 | 女人被狂操c到高潮| 国产午夜精品论理片| av国产免费在线观看| 久久精品91无色码中文字幕| 国产av一区在线观看免费| 国语自产精品视频在线第100页| 国产91精品成人一区二区三区| 精品一区二区三区视频在线 | 久久久精品欧美日韩精品| 亚洲久久久久久中文字幕| 亚洲专区中文字幕在线| 国产精品女同一区二区软件 | 在线观看一区二区三区| 一个人观看的视频www高清免费观看| 尤物成人国产欧美一区二区三区| 国产精品久久久久久久电影 | 精品久久久久久久久久免费视频| 国产三级中文精品| 日韩av在线大香蕉| 国产又黄又爽又无遮挡在线| 亚洲专区中文字幕在线| 亚洲欧美日韩卡通动漫| 国产成人av激情在线播放| 99国产极品粉嫩在线观看| 天天躁日日操中文字幕| 国产欧美日韩精品一区二区| 日韩欧美精品v在线| 黄片大片在线免费观看| 欧美日本视频| 久久精品国产亚洲av香蕉五月| 国产真实乱freesex| 又黄又爽又免费观看的视频| a级毛片a级免费在线| 国产在线精品亚洲第一网站| 亚洲国产欧洲综合997久久,| 99热精品在线国产| 久久精品综合一区二区三区| 99久久99久久久精品蜜桃| 九九热线精品视视频播放| 真人做人爱边吃奶动态| 色综合站精品国产| 国产亚洲精品av在线| 国产精品爽爽va在线观看网站| 麻豆一二三区av精品| 国产亚洲欧美在线一区二区| 欧美一级a爱片免费观看看| 高清日韩中文字幕在线| 最近最新中文字幕大全电影3| 国产高清激情床上av| 欧美中文综合在线视频| 亚洲激情在线av| 亚洲第一电影网av| 国产精品免费一区二区三区在线| 中文字幕人妻熟人妻熟丝袜美 | 真人做人爱边吃奶动态| 久久久精品欧美日韩精品| 亚洲成av人片在线播放无| 国产精品爽爽va在线观看网站| 一级毛片女人18水好多| 国产真人三级小视频在线观看| 欧美日韩乱码在线| 天天添夜夜摸| 禁无遮挡网站| 国产免费男女视频| 中文字幕av成人在线电影| 综合色av麻豆| www国产在线视频色| 久久欧美精品欧美久久欧美| 午夜精品久久久久久毛片777| 国产亚洲精品久久久久久毛片| 欧美黑人巨大hd| 一区二区三区国产精品乱码| 亚洲激情在线av| 老司机午夜十八禁免费视频| 白带黄色成豆腐渣| 久久精品综合一区二区三区| 老司机午夜福利在线观看视频| 亚洲 欧美 日韩 在线 免费| 日日干狠狠操夜夜爽| 看片在线看免费视频| 欧美bdsm另类| 夜夜夜夜夜久久久久| 国产亚洲精品久久久久久毛片| 少妇裸体淫交视频免费看高清| 国产高清三级在线| 国产高清视频在线播放一区| 女人十人毛片免费观看3o分钟| 成熟少妇高潮喷水视频| av在线蜜桃| 亚洲狠狠婷婷综合久久图片| 国产 一区 欧美 日韩| 97人妻精品一区二区三区麻豆| 禁无遮挡网站| 一区二区三区国产精品乱码| 午夜福利欧美成人| 51国产日韩欧美| 久久精品91蜜桃| 女人十人毛片免费观看3o分钟| 色综合婷婷激情| avwww免费| 国产成人av教育| 亚洲av二区三区四区| 极品教师在线免费播放| 国产成人a区在线观看| 五月玫瑰六月丁香| 在线国产一区二区在线| 51午夜福利影视在线观看| 免费无遮挡裸体视频| 最后的刺客免费高清国语| 成年女人看的毛片在线观看| 国产成人aa在线观看| 国产av在哪里看| 久久草成人影院| 国产高清激情床上av| 观看免费一级毛片| 国产成人福利小说| 国产午夜精品论理片| 69av精品久久久久久| 女同久久另类99精品国产91| 精品熟女少妇八av免费久了| 村上凉子中文字幕在线| svipshipincom国产片| 亚洲,欧美精品.| 可以在线观看毛片的网站| 嫩草影院精品99| 国产精品乱码一区二三区的特点| 国模一区二区三区四区视频| 一进一出抽搐gif免费好疼| 国产美女午夜福利| 精华霜和精华液先用哪个| 蜜桃久久精品国产亚洲av| 日韩欧美三级三区| 午夜免费男女啪啪视频观看 | 国内精品一区二区在线观看| а√天堂www在线а√下载| 国语自产精品视频在线第100页| 男插女下体视频免费在线播放| 老熟妇仑乱视频hdxx| 亚洲精华国产精华精| 精品国产超薄肉色丝袜足j| 久久久色成人| 国产野战对白在线观看| 丰满人妻熟妇乱又伦精品不卡| 天堂动漫精品| 国产熟女xx| 一个人观看的视频www高清免费观看| 全区人妻精品视频| 国产成人av教育| 日本五十路高清| 性欧美人与动物交配| 久久精品国产亚洲av香蕉五月| 淫妇啪啪啪对白视频| 国产精品日韩av在线免费观看| 韩国av一区二区三区四区| 亚洲国产精品成人综合色| 国产老妇女一区| 在线观看免费午夜福利视频| 一级a爱片免费观看的视频| 中文字幕人妻熟人妻熟丝袜美 | 亚洲国产欧美网| 岛国在线免费视频观看| 国产精品98久久久久久宅男小说| 久久午夜亚洲精品久久| 久久人人精品亚洲av| 99久国产av精品| 久久久久久久久大av| 免费一级毛片在线播放高清视频| 午夜福利免费观看在线| 岛国在线观看网站| 好男人在线观看高清免费视频| 人人妻人人澡欧美一区二区| 欧美中文综合在线视频| 十八禁网站免费在线| 午夜免费男女啪啪视频观看 | 国产精品乱码一区二三区的特点| 免费观看的影片在线观看| 神马国产精品三级电影在线观看| 欧洲精品卡2卡3卡4卡5卡区| 国产黄a三级三级三级人| 国产免费男女视频| 午夜福利视频1000在线观看| av天堂在线播放| 成人18禁在线播放| 久久精品综合一区二区三区| 国产亚洲精品综合一区在线观看| 欧美成人a在线观看| 成人av一区二区三区在线看| 久久草成人影院| 日韩欧美三级三区| 亚洲国产欧美网| 亚洲国产精品成人综合色| 两个人的视频大全免费| 特大巨黑吊av在线直播| 国产蜜桃级精品一区二区三区| 脱女人内裤的视频| 最近视频中文字幕2019在线8| 欧美最新免费一区二区三区 | 日本撒尿小便嘘嘘汇集6| 黄色女人牲交| 国产激情偷乱视频一区二区| 免费人成视频x8x8入口观看| 黄色日韩在线| 午夜激情福利司机影院| 国产精品久久久人人做人人爽| 看片在线看免费视频| 亚洲aⅴ乱码一区二区在线播放| 日本黄大片高清| 成年女人毛片免费观看观看9| 亚洲av美国av| 老司机深夜福利视频在线观看| 亚洲精品久久国产高清桃花| 午夜免费成人在线视频| 老司机午夜福利在线观看视频| 欧美av亚洲av综合av国产av| 老司机午夜福利在线观看视频| 91av网一区二区| 日韩欧美在线二视频| 中文字幕人妻丝袜一区二区| 久久精品国产综合久久久| 中国美女看黄片| 亚洲中文字幕一区二区三区有码在线看| 欧美成狂野欧美在线观看| 午夜久久久久精精品| 欧美乱码精品一区二区三区| 琪琪午夜伦伦电影理论片6080| 亚洲在线观看片| 五月玫瑰六月丁香| 国产色婷婷99| 亚洲熟妇熟女久久| 精品久久久久久久人妻蜜臀av| 少妇的逼好多水| 麻豆一二三区av精品| a级毛片a级免费在线| 国产亚洲精品av在线| 性色av乱码一区二区三区2| 成人国产综合亚洲| 亚洲精品影视一区二区三区av| 变态另类成人亚洲欧美熟女| 国产成人影院久久av| 色哟哟哟哟哟哟| 黄色片一级片一级黄色片| av天堂在线播放| 蜜桃亚洲精品一区二区三区| 色噜噜av男人的天堂激情| 国产av一区在线观看免费| 深夜精品福利| a级一级毛片免费在线观看| av视频在线观看入口| 亚洲专区中文字幕在线| 国产高清videossex| 免费看美女性在线毛片视频| 国内精品美女久久久久久| 欧美极品一区二区三区四区| 欧美性猛交黑人性爽| 熟女电影av网| 国产精品久久久久久久久免 | 亚洲激情在线av| 久久精品亚洲精品国产色婷小说| 精品欧美国产一区二区三| 亚洲国产精品sss在线观看| 国产精品电影一区二区三区| 亚洲精华国产精华精| 欧美日本亚洲视频在线播放| 怎么达到女性高潮| 老鸭窝网址在线观看| 亚洲真实伦在线观看| 日韩欧美在线二视频| 国产成年人精品一区二区| 亚洲久久久久久中文字幕| 在线播放无遮挡| 日本 av在线| 高潮久久久久久久久久久不卡| 九色成人免费人妻av| 日韩欧美在线乱码| а√天堂www在线а√下载| 亚洲av二区三区四区| 青草久久国产| 伊人久久大香线蕉亚洲五| 免费看日本二区| 亚洲激情在线av| 国产伦精品一区二区三区四那| 极品教师在线免费播放| 国内精品一区二区在线观看| 午夜免费激情av| 亚洲 国产 在线| 成熟少妇高潮喷水视频| 精品乱码久久久久久99久播|