張雙圣,劉漢湖,強 靜,劉喜坤,朱雪強
基于貝葉斯公式的地下水污染源及含水層參數(shù)同步反演
張雙圣1,3,劉漢湖1,強 靜2*,劉喜坤3,朱雪強1
(1.中國礦業(yè)大學(xué)環(huán)境與測繪學(xué)院,江蘇 徐州 221116;2.中國礦業(yè)大學(xué)數(shù)學(xué)學(xué)院,江蘇 徐州 221116;3.徐州市城區(qū)水資源管理處,江蘇 徐州 221018)
監(jiān)測方案優(yōu)化;貝葉斯公式;信息熵;Kriging替代模型;差分進(jìn)化自適應(yīng)Metropolis算法;參數(shù)后驗均值偏離率
地下水污染源及含水層參數(shù)反演是指通過建立地下水水流及溶質(zhì)運移模型,運用已有的污染物濃度監(jiān)測數(shù)據(jù)反求污染源位置、污染源釋放強度、釋放時間等污染源信息以及含水層參數(shù)等.其本質(zhì)是運用監(jiān)測數(shù)據(jù)對地下水水流及溶質(zhì)運移模型的源匯項及控制方程參數(shù)進(jìn)行反演.地下水污染源及含水層參數(shù)的反演識別方法主要包括解析法、模擬優(yōu)化方法和不確定性分析方法等.解析法形式簡單,而且計算效率較高,但該方法主要適用于水文地質(zhì)條件比較單一,而且必須滿足一定的假定條件的地下水系統(tǒng).由于反演問題均可轉(zhuǎn)化為優(yōu)化問題,因此模擬優(yōu)化方法成為地下水污染源識別及含水層參數(shù)反演的常用方法.求解模擬優(yōu)化問題的算法主要有模擬退火算法[1]、遺傳算法[2]、禁忌搜索算法[3]等,這些算法均屬于啟發(fā)式智能搜索優(yōu)化算法,具有全局搜索性能,保證全局最優(yōu),有效克服傳統(tǒng)優(yōu)化算法不收斂或局部最優(yōu)的缺點.但是模擬優(yōu)化算法未考慮參數(shù)取值、測量誤差等各種不確定性因素,得到參數(shù)的一個確定的估計值,容易丟掉真值.
近年來,不確定性分析方法在研究參數(shù)反演問題中應(yīng)用較多,主要包括廣義似然不確定性估計方法[4]、地質(zhì)統(tǒng)計學(xué)方法[5]、卡爾曼濾波方法[6]以及貝葉斯統(tǒng)計方法[7]等.貝葉斯統(tǒng)計方法應(yīng)用最為廣泛,常借助獨立抽樣的蒙特卡羅方法(MC)[8]進(jìn)行近似求解,其中馬爾科夫鏈蒙特卡羅方法(MCMC)[9-17]作為一種經(jīng)典抽樣方法得到廣泛應(yīng)用.單鏈MCMC算法[9-14]容易存在反演不收斂,或者反演結(jié)果局部最優(yōu)的問題[17].而多鏈MCMC算法[15-16]適用于參數(shù)維度高,有多個局部最優(yōu)值點,搜索量大的參數(shù)空間,能夠更好地解決馬爾科夫鏈局部收斂的問題[15-16].差分進(jìn)化Metropolis算法(DE- MC)[15]和差分進(jìn)化自適應(yīng)Metropolis算法(DREAM)[16]是常用多鏈MCMC算法,其中 DREAM算法是DE-MC算法的改進(jìn)版本,采用自適應(yīng)隨機子空間采樣技術(shù)及能夠自適應(yīng)調(diào)整的交叉概率,并且運用IQR統(tǒng)計方法[16]去除無用鏈,可顯著提高搜索效率和求解結(jié)果的精度.
在污染源及含水層參數(shù)反演過程中,往往出現(xiàn)監(jiān)測數(shù)據(jù)與反演參數(shù)關(guān)聯(lián)性較弱的問題,因此對監(jiān)測方案(包括監(jiān)測井的位置、數(shù)量及監(jiān)測頻率等)進(jìn)行優(yōu)化設(shè)計成為參數(shù)反演的關(guān)鍵.通常做法是定義某個目標(biāo)函數(shù),對監(jiān)測方案的信息含量進(jìn)行量化.常用的目標(biāo)函數(shù)有信噪比(SNR)[18],以及基于貝葉斯公式的相對熵[19-20].但信噪比僅考慮監(jiān)測誤差對監(jiān)測數(shù)據(jù)的干擾影響,而忽略了監(jiān)測數(shù)據(jù)對參數(shù)后驗分布的影響;相對熵未考慮參數(shù)先驗分布對后驗分布的影響.為此將信息熵概念[21-22]引入到地下水監(jiān)測方案優(yōu)化設(shè)計中.
另外在污染源及含水層參數(shù)反演過程中均需反復(fù)調(diào)用地下水水流及溶質(zhì)運移數(shù)值模擬模型,導(dǎo)致巨大的計算負(fù)荷.通常采用構(gòu)建數(shù)值模擬模型的替代模型的方法減少計算量.常用構(gòu)建替代模型方法有多項式回歸法[23],徑向基函數(shù)神經(jīng)網(wǎng)絡(luò)法[24], Kriging法[25-26]等.其中Kriging替代模型具有兼顧局部和全局的統(tǒng)計特性,適用于解決非線性程度較高的問題.
本文充分考慮模擬模型參數(shù)不確定性、替代模型誤差不確定性以及測量誤差的不確定性,運用地下水模擬軟件GMS(Groundwater Modeling System)建立二維非均質(zhì)各向同性含水層水流及溶質(zhì)運移數(shù)值模型.在確定初始監(jiān)測時刻、監(jiān)測間隔時間及監(jiān)測次數(shù)條件下,利用最優(yōu)拉丁超立方抽樣方法及Kriging法建立數(shù)值模擬模型的替代模型,采用基于信息熵最小累進(jìn)加井的方式制定監(jiān)測方案,并根據(jù)此優(yōu)化監(jiān)測方案監(jiān)測數(shù)據(jù),基于貝葉斯統(tǒng)計方法運用DREAM算法進(jìn)行污染源及含水層參數(shù)的同步反演.
貝葉斯公式如下:
本研究采用MCMC (Markov Chain Monte Carlo)算法對(5)式的后驗分布進(jìn)行求解.
式(11)的求解算法較為復(fù)雜,難以得出顯示表達(dá)式,本文運用文獻(xiàn)[20,22]中蒙特卡羅方法近似求解.
DREAM算法具體步驟可查閱文獻(xiàn)[16].采用Gelman-Rubin收斂診斷方法[27]對DREAM算法后50%采樣過程的收斂性進(jìn)行判斷,判斷指標(biāo)為:
如果馬爾科夫鏈滿足Gelman-Rubin收斂準(zhǔn)則,則計算終止,否則繼續(xù)進(jìn)化平行序列.
為了減少監(jiān)測方案優(yōu)化設(shè)計及參數(shù)反演過程中多次調(diào)用地下水溶質(zhì)運移數(shù)值模擬模型產(chǎn)生的計算量,運用Kriging法[25-26]構(gòu)造數(shù)值模擬模型的替代模型.
替代模型的建立與DREAM算法初始樣本的選取均需采取一定的抽樣方法抽取樣本.在替代模型構(gòu)建過程中,為保證替代模型能夠捕捉到研究對象函數(shù)的變化趨勢,保證抽取的樣本均勻分布在先驗分布的全空間,以最優(yōu)拉丁超立方抽樣方法[28]進(jìn)行樣本抽樣.而DREAM算法是多鏈MCMC算法,可采用拉丁超立方抽樣方法[25]在參數(shù)先驗分布范圍內(nèi)相對均勻的抽取隨機樣本作為初始樣本,降低了隨機選取樣本對反演結(jié)果產(chǎn)生的影響.
研究區(qū)內(nèi)共有10眼監(jiān)測井,地下水水質(zhì)背景值較好,含水層污染物的初始濃度為零,某日研究區(qū)下游發(fā)現(xiàn)污染物X,初步斷定污染源S在上游的某一區(qū)域范圍內(nèi),且在某時間段內(nèi)以注水井的形式(200m3/d)持續(xù)恒定地向含水層中排放污染物.研究區(qū)平面示意圖如圖1.
圖1 算例模型示意
表1 研究區(qū)域已知水文地質(zhì)參數(shù)
以西南角為坐標(biāo)原點建立坐標(biāo)系,根據(jù)研究區(qū)水文地質(zhì)條件,建立地下水水流數(shù)值模型:
建立的地下水水流及溶質(zhì)運移模型運用GMS軟件中的MODFLOW和MT3DMS模塊進(jìn)行計算求解.為保證每個網(wǎng)格中心對應(yīng)一個潛在污染源(污染源樣本)的位置,將研究區(qū)域剖分為100行250列的正方形有限差分網(wǎng)格,基本單元格邊長為4m.
表2 10眼監(jiān)測井的坐標(biāo)位置
表3 待求參數(shù)的先驗分布范圍
表4 從待求參數(shù)的先驗分布范圍中得到400組訓(xùn)練輸入樣本
(2)分別對10眼監(jiān)測井建立Kriging替代模型將表4中400組參數(shù)分別代入到GMS軟件中,得到10眼監(jiān)測井在[800d,1027d]內(nèi)每隔30d的污染物濃度計算值,作為Kriging替代模型輸出值.將400組輸入值與輸出值作為訓(xùn)練樣本代入MATLAB軟件中,運用DACE工具箱建立各備選監(jiān)測井的Kriging替代模型.
表5 從待求參數(shù)的先驗分布范圍中得到30組檢驗輸入樣本
圖2 Kriging替代模型與數(shù)值模擬模型輸出結(jié)果對比
表6 10眼備選監(jiān)測井Kriging替代模型的R2,RMSE及WRE
由表6數(shù)據(jù)可知, 10眼備選監(jiān)測井Kriging替代模型的平均決定系數(shù)2=0.9874,表明Kriging替代模型模擬精度較高.
首先優(yōu)化單井監(jiān)測方案,即求(11)式的最小值,問題可概化如下:
表7 從參數(shù)的先驗分布中得到20組“參數(shù)真值”
常用的參數(shù)反演效果評價指標(biāo)有相對均方根誤差和相對誤差均值[30],指標(biāo)具體公式如下:
②參數(shù)后驗均值估計MP與“參數(shù)真值”之間相對誤差均值
其中:表示第組“參數(shù)真值”,表示參數(shù)的第個分量.
表8 單井監(jiān)測方案,
由表8可知,6#備選監(jiān)測井的信息熵最小,故把其作為多井監(jiān)測方案選定的第1眼監(jiān)測井.后續(xù)選取多井監(jiān)測方案中的第2眼監(jiān)測井,步驟為:將選定6#井和其余9眼備選監(jiān)測井分別進(jìn)行組合,得到9種組合形式,并計算這9種組合監(jiān)測方案的信息熵.此問題可概化如下:
由表9可知,信息熵隨著監(jiān)測井?dāng)?shù)量的增加而減小,其中4眼井組合監(jiān)測方案信息熵、5眼井組合監(jiān)測方案信息熵以及10眼井組合監(jiān)測方案信息熵之間差距很小,故6眼井及以上組合監(jiān)測方案信息熵不再贅述.兼顧反演精度及監(jiān)測成本,并滿足每個參數(shù)分區(qū)至少有1眼監(jiān)測井的要求,本算例最優(yōu)監(jiān)測方案為5眼井組合監(jiān)測方案(6,5,1,2,8)(定義為方案1).為了驗證所選監(jiān)測方案參數(shù)反演效果,與10眼井組合監(jiān)測方案(1,2,3,4,5,6,7,8,9,10)(定義為方案2)的參數(shù)反演結(jié)果進(jìn)行對比分析.
表9 多井監(jiān)測方案優(yōu)化設(shè)計中各監(jiān)測方案的E(MP)
以表7中第1組參數(shù)真值為例,采用2.3節(jié)選定的2種監(jiān)測方案,運用DREAM算法(初始平行鏈數(shù)為11,反演穩(wěn)定后平行鏈數(shù)記為)進(jìn)行參數(shù)反演,其中每條馬爾科夫鏈長度是50000.
由表10可知,方案1中11個參數(shù)后驗均值偏離率的平均值為10.3%,方案2中11個參數(shù)后驗均值偏離率的平均值為9.1%,兩種監(jiān)測方案下參數(shù)后驗均值偏離率結(jié)果相近,即2種方案的反演效果類似.但是方案2用了10眼監(jiān)測井,方案1僅用了5眼監(jiān)測井,方案2所需監(jiān)測成本是方案1所需監(jiān)測成本的2倍.對于本文案例,基于監(jiān)測成本及反演精度雙重考慮,認(rèn)為方案1為最佳的監(jiān)測方案.
表10 2種監(jiān)測方案下模型參數(shù)后驗統(tǒng)計結(jié)果及收斂性判斷指標(biāo)
由表10可知,污染源位置參數(shù)S和S的后驗分布范圍較先驗分布范圍顯著減小(其中S的后驗分布范圍僅為先驗分布范圍的25%,S的后驗分布范圍為先驗分布范圍的10%),而其他參數(shù)的后驗分布范圍未見明顯減小,這主要是由于各參數(shù)對污染物濃度監(jiān)測值的靈敏度不同造成的.
為避免局部靈敏度分析方法沒有考慮參數(shù)之間相互作用對輸出結(jié)果影響的缺陷,運用全局靈敏度分析方法(Sobol’法)[31]及2.2節(jié)Kriging替代模型,得到參數(shù)對10眼備選監(jiān)測井10個監(jiān)測數(shù)據(jù)的1階靈敏度系數(shù),見表11.參照參數(shù)靈敏度分級[31](表12),參數(shù)S和S對于各監(jiān)測井監(jiān)測數(shù)據(jù)來講屬于中等靈敏參數(shù)或靈敏參數(shù),而其余參數(shù)對于各監(jiān)測井監(jiān)測數(shù)據(jù)來講多為不靈敏參數(shù).
本文算例采用理想地下水水流及溶質(zhì)運移模型進(jìn)行了監(jiān)測方案優(yōu)化設(shè)計及參數(shù)反演.但這些方法應(yīng)用到實際案例中時,首先需要保證建立的水流模型符合實際問題,需對水流模型各水文地質(zhì)參數(shù)進(jìn)行反演識別,然后再進(jìn)行監(jiān)測方案優(yōu)化設(shè)計及溶質(zhì)運移模型參數(shù)反演.
表11 參數(shù)a的1階靈敏度系數(shù)
表12 參數(shù)靈敏度分級
3.1 運用最優(yōu)拉丁超立方抽樣方法和Kriging法建立數(shù)值模擬模型的替代模型,其精度較高,能以較小的計算量得到和地下水?dāng)?shù)值模擬模型相近的輸入輸出關(guān)系,顯著降低監(jiān)測方案優(yōu)化設(shè)計與污染源及含水層參數(shù)識別過程中反復(fù)調(diào)用地下水溶質(zhì)運移數(shù)值模擬模型產(chǎn)生的計算負(fù)荷.
3.2 參數(shù)后驗均值偏離率有效避免了參數(shù)真值取值對反演結(jié)果評價的影響,適用性更廣.參數(shù)反演結(jié)果的后驗均值偏離率與監(jiān)測方案信息熵呈現(xiàn)較好的正相關(guān)關(guān)系.參數(shù)后驗均值偏離率越小,說明參數(shù)反演效果越好.信息熵是參數(shù)后驗分布不確定性的有效量度,信息熵越小,參數(shù)后驗范圍越小,基于貝葉斯公式與信息熵的監(jiān)測方案優(yōu)化設(shè)計方法是確定地下水污染監(jiān)測方案的有效方法.
3.3 多井監(jiān)測方案優(yōu)化設(shè)計過程中,不同數(shù)量監(jiān)測井組合產(chǎn)生的備選監(jiān)測方案數(shù)量較大,為尋求最優(yōu)監(jiān)測方案,需求解所有備選監(jiān)測方案的信息熵,導(dǎo)致計算量較大.而基于信息熵最小的累進(jìn)加井的多井監(jiān)測方案優(yōu)化設(shè)計方法,可有效降低優(yōu)化設(shè)計過程中求解信息熵的計算量, 監(jiān)測方案優(yōu)化效率得到顯著提高.
3.4 DREAM算法采用多條馬爾科夫鏈并行計算策略及自適應(yīng)隨機子空間采樣技術(shù),并且能夠自適應(yīng)調(diào)整交叉概率,可大大提高后驗概率分布的搜索效率,能夠有效解決局部收斂問題,而且能夠?qū)崿F(xiàn)多參數(shù)同步反演,計算效率顯著提高.
[1] Dougherty D E, Marryott R A. Optimal groundwater management: 1. Simulated annealing [J]. Water Resources Research, 1991,27(10): 2493-2508.
[2] Mahinthakumar G, Sayeed M. Hybrid genetic algorithm-local search methods for solving groundwater source identification inverse problems [J]. Journal of Water Resources Planning and Management, 2005,131(1):45-57.
[3] Glover F, Laguna M. Tabu Search [M]. Boston, MA: Kluwer Academic Publisher, 1997:125-151.
[4] Hassan A E, Bekhit H M, Chapman J B. Uncertainty assessment of a stochastic groundwater flow model using GLUE analysis [J]. Journal of Hydrology, 2008,362(1/2):89-109.
[5] Snodgrass M F, Kitanidis P K. A geostatistical approach to contaminant source identification [J]. Water Resources Research, 1997, 33(4):537-546.
[6] Rasmussen J, Madsen H, Jensen K H, et al. Data assimilation in integrated hydrological modeling using ensemble Kalman filtering: evaluating the effect of ensemble size and localization on filter performance [J]. Hydrology and Earth System Sciences, 2015,12(2): 2267-2304.
[7] Chen M, Izady A, Abdalla O A, et al. A surrogate-based sensitivity quanti?cation and Bayesian inversion of a regional groundwater ?ow model [J]. Journal of Hydrology, 2018,557:826-837.
[8] Roberts C P, Casella G. Monte Carlo statistical methods (Second edition) [M]. Berlin: Springer, 2004:79-122.
[9] Metropolis N, Rosenbluth A W, Rosenbluth M N, et al. Equation of state calculations by fast computing machines [J]. The Journal of Chemical Physics, 1953,21(6):1087-1092.
[10] Hastings W K. Monte Carlo sampling methods using Markov chains and their applications [J]. Biometrika, 1970,57(1):97-109.
[11] Haario H, Saksman E, Tamminen J. An adaptive Metropolis algorithm [J]. Bernoulli, 2001,7(2):223-242.
[12] Tierney L, Mira A. Some adaptive Monte Carlo methods for bayesian inference [J]. Statistics in Medicine, 1999,18:2507-2515.
[13] Mira A. Ordering and improving the performance of Monte CarloMarkov Chains [J]. Statistical Science, 2002,16:340-350.
[14] Haario H, Laine M, Mira A. DRAM: Efficient adaptive MCMC [J]. Statistics and Computing, 2006,16(4):339-354.
[15] Ter Braak C J F. A Markov chain Monte Carlo version of the genetic algorithm differential evolution: easy Bayesian computing for real parameter spaces [J]. Statistics and Computing, 2006,16(3):239-249.
[16] Vrugt J A, Ter Braak C J F, Diks C G H, et al. Accelerating Markov Chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling [J]. International Journal of Nonlinear Sciences and Numerical Simulation, 2009,10(3):273- 290.
[17] Vrugt J A. Markov chain Monte Carlo simulation using the DREAM software package: Theory, consepts, and Matlab implementation [J]. Enviromental Modeling & Software, 2016,75:273-316.
[18] Czanner G, Sarma S V, Eden U T, et al. A Signal-to-Noise Ratio Estimator for Generalized Linear Model Systems [J]. Lecture Notes in Engineering & Computer Science, 2008,2171:1063-1069.
[19] Lindley D V. On a measure of the information provided by an experiment [J]. The Annals of Mathematical Statistics, 1956,27(4): 986-1005.
[20] Huan X, Marzouk, Y M. Simulation-based optimal Bayesian experimental design for nonlinear systems [J]. Journal of Computational Physics, 2013,232(1):288-317.
[21] Shannon C E. A mathematical theory of communication [J]. The Bell System Technical Journal, 1948,27(3):379-423.
[22] 張雙圣,強 靜,劉漢湖,等.基于貝葉斯公式的地下水污染源識別 [J]. 中國環(huán)境科學(xué), 2019,39(4):1568-1578. Zhang S, Qiang J, Liu H, et al. Identification of groundwater pollution sources based on Bayes’ theorem [J]. China Environmental Science, 2019,39(3):1568-1578.
[23] Knill D L, Giunta A A, Baker C A, et al. Response surface models combining linear and Euler aerodynamics for supersonic transport design [J]. Journal of Aircraft, 1999,36(1):75-86.
[24] Li J, Chen Y, Pepper, D. Radial basis function method for 1-D and 2-D groundwater contaminant transport modeling [J]. Computational Mechanics, 2003,32(1):10-15.
[25] 安永凱,盧文喜,董海彪,等.基于克里格法的地下水流數(shù)值模擬模型的替代模型研究 [J]. 中國環(huán)境科學(xué), 2014,34(4):1073-1079. An Y, Lu W, Dong H, et al. Surrogate model of numerical simulation model of groundwater based on Kriging [J]. China Environmental Science, 2014,34(4):1073-1079.
[26] Lophaven S N, Nielsen H B, Sondergaard J. Dace: A MATLAB Kriging toolbox [R]. Kongens Lyngby: Technical University of Denmark, Technical Report No. IMM-TR-2002-12.
[27] Grlman A G, Rubin D B. Inference from iterative simulation using multiplesequences [J]. Statistical Science, 1992,7:457-472.
[28] Hickernell, F.A. A generalized discrepancy and quadrature error bound [J]. Mathematics of Computation of the American Mathematical Society, 1998,67(221):299-322.
[29] 初海波.基于自適應(yīng)替代模型的DNAPLs污染地下水修復(fù)優(yōu)化設(shè)計及其不確定性分析 [D]. 長春:吉林大學(xué), 2015.Chu H. The optimization design and uncertainty analysis of DNAPLs- contaminated groundwater remediation based on the adaptive surrogate model [D]. Changchun: Jilin University, 2015.
[30] 張江江.地下水污染源解析的貝葉斯監(jiān)測設(shè)計與參數(shù)反演方法 [D]. 杭州:浙江大學(xué), 2017. Zhang J. Bayesian monitoring design and parameter inversion for groundwater contaminant source identification [D]. Hangzhou: Zhejiang University, 2017.
[31] Lenhart L, Eckhardt K, Fohrer N, et al. Comparison of two different approaches of sensitivity analysis [J]. Physics and Chemistry of the Earth, 2002,27:645-654.
Synchronous inversion of groundwater pollution source and aquifer parameters based on Bayesian formula.
ZHANG Shuang-sheng1,3, LIU Han-hu1, QIANG Jing2*, LIU Xi-kun3, ZHU Xue-qiang1
(1.School of Environment Science and Spatial Informatics, China University of Mining and Technology, Xuzhou 221116, China;2.School of Mathematics, China University of Mining and Technology, Xuzhou 221116, China;3.Xuzhou City Water Resource Administrative Office, Xuzhou 221018, China)., 2019,39(7):2902~2912
Aiming at the optimization of monitoring schemes in the process of the identification of pollution source and the inversion of aquifer parameters in the heterogeneous underground aquifer, this paper proposes an optimization method for the multi-well monitoring schemes based on Bayesian formula and progressive addition of wells with minimum information entropy. Firstly, the two-dimensional heterogeneous isotropic subsurface groundwater flow and solute transport models under hypothetical case were constructed, and the numerical simulation models were solved by GMS software. The surrogate model of the numerical simulation model was established by the optimal Latin hypercube sampling method and Kriging method. Then Taking the minimum information entropy of the parameter posterior distribution as the objective function, the optimization design of multi-well monitoring schemes was carried out by means of progressive addition of wells. Finally, the differential evolution adaptive Metropolis algorithm was used to inverse the pollution source and aquifer parameters synchronously according to the optimized monitoring scheme. The case study results showed that: The 5combination monitoring scheme (6, 5, 1, 2, 8) under the condition of taking into account the inversion accuracy and monitoring cost and ensuring that there was at least one monitoring well in each parameter section was the optimal monitoring scheme. Compared with the 10combined monitoring scheme (1, 2, 3, 4, 5, 6, 7, 8, 9, 10) with the smallest information entropy, the 11parameters posterior mean deviation rate increased by 1.2%, but the monitoring cost was reduced by 50%.
monitoring well optimization;bayesian formula;information entropy;kriging surrogate model;differential evolution adaptive metropolis algorithm;parameter posterior mean deviation rate
X523
A
1000-6923(2019)07-2902-11
張雙圣(1983-),男,山東省昌邑市人,工程師,博士,主要從事地下水污染控制研究.發(fā)表論文20余篇.
2018-11-27
國家水體污染控制與治理科技重大專項基金資助項目(2015ZX07406005)
* 責(zé)任作者, 講師, 57591340@qq.com