摘 要:研究地下水中溶質(zhì)的遷移規(guī)律可以為地下水環(huán)境質(zhì)量預(yù)報(bào)、控制和管理服務(wù)。首次采用數(shù)值方法探討與位置坐標(biāo)相關(guān)的溶質(zhì)運(yùn)移參數(shù)模型,采用即時(shí)時(shí)間域顯式差分法格式雙精度Fortran語(yǔ)言編程計(jì)算了變?nèi)苜|(zhì)運(yùn)移參數(shù)模型的一維溶質(zhì)運(yùn)移問題,并進(jìn)行了算法可靠性驗(yàn)證;然后對(duì)本文首次提出的3種與位置坐標(biāo)相關(guān)的變?nèi)苜|(zhì)運(yùn)移參數(shù)模型進(jìn)行了參數(shù)分析,獲得了這3種模型的初步認(rèn)識(shí)。還將時(shí)間域顯式差分方程計(jì)算結(jié)果與解析解和實(shí)驗(yàn)結(jié)果進(jìn)行了驗(yàn)證對(duì)比分析,說明算法可靠,顯式法列式簡(jiǎn)單,更適于復(fù)雜模型計(jì)算。本文提出的算法和處理技巧可以應(yīng)用到許多與地質(zhì)相關(guān)的溶質(zhì)運(yùn)移問題的研究,通過建立不同模式的溶質(zhì)運(yùn)移參數(shù)模型,從而實(shí)現(xiàn)復(fù)雜條件下溶質(zhì)運(yùn)移過程的預(yù)測(cè)和模擬。
關(guān)鍵詞:一維溶質(zhì)遷移;與深度相關(guān);變?nèi)苜|(zhì)運(yùn)移參數(shù);建模與模擬;即時(shí)顯式差分法
中圖分類號(hào):O357.3
文獻(xiàn)標(biāo)志碼:A
由于工農(nóng)業(yè)生產(chǎn)、雨水徑流沖刷和人為因素等作用,造成水流中溶質(zhì)瞬間注入或連續(xù)注入從而引起河流或地下水的水質(zhì)改變。研究溶質(zhì)瞬時(shí)注入或連續(xù)注入情況下不同時(shí)刻溶質(zhì)濃度在河流或地下水流中的分布狀況、及不同位置時(shí)溶質(zhì)濃度隨時(shí)間的變化情況,可通過數(shù)學(xué)模型來模擬或預(yù)報(bào)水質(zhì)在時(shí)間或空間的變化,達(dá)到為水環(huán)境質(zhì)量預(yù)報(bào)、控制和水資源管理服務(wù)的目的[1-22]。
數(shù)值模擬水流溶質(zhì)運(yùn)移的主要方法有有限元法、邊界單元法、有限差分法和其他數(shù)值法等。由于進(jìn)行水質(zhì)模擬過程中涉及的參數(shù)在不同的時(shí)間和空間上變化較大,在水質(zhì)模擬過程中,會(huì)造成數(shù)值彌散、數(shù)值波動(dòng)和計(jì)算量巨大等問題。國(guó)內(nèi)主要的研究有徐玉佩[2-3]用Galerkin法推導(dǎo)了二維地下水運(yùn)動(dòng)的有限元方程和溶質(zhì)運(yùn)移有限元方程,并編制了計(jì)算程序應(yīng)用于明溝排水和田面保持淹灌水層情況下的沖洗種稻改良鹽堿地的分析,還與現(xiàn)有解析解和實(shí)驗(yàn)資料進(jìn)行了對(duì)照。徐玉佩還用水動(dòng)力彌散理論進(jìn)行了地下水中水鹽動(dòng)態(tài)的研究,可直接用于鹽堿地的預(yù)防和改良。黃康樂[4]則提出了一種求解二維飽和-非飽和溶質(zhì)運(yùn)移問題的交替方向特征有限單元法(ADCG),在克服數(shù)值彌散、數(shù)值波動(dòng)和提高計(jì)算速度方面具有顯著的優(yōu)越性。張效先等[5]提出了一種計(jì)算地下水中溶質(zhì)運(yùn)移的數(shù)值法,將多維問題分解成幾個(gè)一維問題, 減少了計(jì)算工作量, 并消除了數(shù)值波動(dòng)現(xiàn)象,具有較高的精度。馮紹元[6]提出了特征有限元法的基本理論和實(shí)施過程,其計(jì)算實(shí)例表明能消除大部分的數(shù)值彌散而又不增加解的振蕩性。姚磊華[7]首先基于分步求解的思想, 用廣義迎風(fēng)對(duì)偶單元均衡法求解對(duì)流定解問題, 對(duì)擴(kuò)散定解問題采用一般的Galerkin有限元法求解, 不僅避免了數(shù)值彌散和過量問題, 而且避免了求節(jié)點(diǎn)速度這一步, 簡(jiǎn)化了運(yùn)算步驟。邱克儉等[8]結(jié)合拉格朗日觀點(diǎn)的運(yùn)動(dòng)坐標(biāo)與歐拉觀點(diǎn)的固定坐標(biāo),分別求解對(duì)流-彌散方程的對(duì)流效應(yīng)和彌散效應(yīng), 并通過與理論解和試驗(yàn)結(jié)果做比較,驗(yàn)證了模型的可靠性。任理等[9]采用混合拉普拉斯變換有限單元法求解對(duì)流占優(yōu)的地下水溶質(zhì)運(yùn)移問題,能有效地消除數(shù)值擴(kuò)散和過量現(xiàn)象,特別適于大區(qū)域地下水污染的長(zhǎng)期預(yù)報(bào)。李煥榮等[10-11]利用廣義差分法建立了一維非飽和水流問題的守恒形式的數(shù)值模型, 具有計(jì)算量小和穩(wěn)定性好的特點(diǎn), 其提出的非粘性土壤水中溶質(zhì)運(yùn)移問題的守恒混合元格式,精度高且數(shù)值穩(wěn)定。段德宏等[14]以黃河蘭州段為例,采用一維穩(wěn)態(tài)河流水質(zhì)模型,對(duì)岸邊19個(gè)污水排放口及共約32.9 km長(zhǎng)的全河段污染狀況成功地進(jìn)行了模擬。徐文彬等[15]根據(jù)數(shù)值隨時(shí)間步長(zhǎng)變化的特點(diǎn),提出了針對(duì)數(shù)值波動(dòng)、彌散和過量現(xiàn)象的新改進(jìn)方案,并用實(shí)例驗(yàn)證了方法的可行性。何丕文[16]利用一維水質(zhì)模型,采用有限差分法求解河流中污染物濃度,以淮河淮南段為例進(jìn)行了驗(yàn)證。馬東豪等[17]研究了土壤溶質(zhì)遷移的兩流區(qū)與兩區(qū)模型,結(jié)果表明與兩區(qū)模型相比,兩流區(qū)模型可以更好地描述優(yōu)先流情況下的溶質(zhì)穿透曲線。對(duì)于無優(yōu)先流情況, 連續(xù)流輸入情況下兩個(gè)模型的參數(shù)一致性較好;而在脈沖輸入時(shí)一致性較差。國(guó)外在溶質(zhì)運(yùn)移模型及其數(shù)值計(jì)算方法上的研究也很多,如JAISWAL等[18]運(yùn)用Chebyshev配點(diǎn)法求解一維溶質(zhì)運(yùn)移問題的數(shù)值解。SILAVWE等[19]討論了某些數(shù)值方法用于一維溶質(zhì)運(yùn)移模型的參數(shù)估計(jì),通過將數(shù)值解與試驗(yàn)結(jié)果比較,反演溶質(zhì)運(yùn)移模型參數(shù)。JAISWAL等[20]采用算子矩陣法對(duì)多孔介質(zhì)的非線性微分方程進(jìn)行了數(shù)值分析,該法應(yīng)用切比雪夫多項(xiàng)式并結(jié)合切比雪夫?qū)?shù)運(yùn)算矩陣和譜配置法對(duì)解進(jìn)行逼近。其優(yōu)點(diǎn)是可以將這些問題轉(zhuǎn)化為易于求解的代數(shù)方程組。CARR[21]利用拉普拉斯變換導(dǎo)出了層狀多孔介質(zhì)溶質(zhì)運(yùn)移的新半解析解,該法在相鄰層的界面引入表示運(yùn)移量的未知函數(shù),將多層問題置于拉普拉斯域上分層解決,然后再將其數(shù)值反演回時(shí)域。ALLWRIGHT等[22] 對(duì)地下水溶質(zhì)運(yùn)移方程在對(duì)流明顯強(qiáng)于擴(kuò)散時(shí)出現(xiàn)的數(shù)值波動(dòng)提出了兩種新的數(shù)值近似計(jì)算方法,即只適用于對(duì)流項(xiàng)的迎風(fēng)Crank-Nicolson格式和加權(quán)迎風(fēng)-下風(fēng)格式,并對(duì)這些新提出的格式進(jìn)行了數(shù)值穩(wěn)定性分析和比較。發(fā)現(xiàn)如果Crank-Nicolson格式只用于對(duì)流項(xiàng)時(shí)很合宜。此外,顯式加權(quán)迎風(fēng)-下風(fēng)有限差分格式是對(duì)傳統(tǒng)顯式一階迎風(fēng)格式的改進(jìn),而隱式加權(quán)一階迎風(fēng)-下風(fēng)有限差分格式在給定適當(dāng)?shù)臋?quán)重因子時(shí)無條件穩(wěn)定。
本文主要貢獻(xiàn)在于對(duì)一維溶質(zhì)對(duì)流彌散基本方程,首先,建立了一維溶質(zhì)運(yùn)移3種變?nèi)苜|(zhì)運(yùn)移參數(shù)模型,即線性變化模型/多項(xiàng)式變化模型、指數(shù)變化模型和冪變化模型;其次,采用顯式時(shí)域差分法成功地進(jìn)行了溶質(zhì)濃度時(shí)空分布的數(shù)值計(jì)算,有效克服了很多研究人員遭遇的數(shù)值波動(dòng)、數(shù)值彌散和計(jì)算量大等難題,并與解析解和實(shí)驗(yàn)結(jié)果做了對(duì)比,驗(yàn)證了該數(shù)值方法的正確性;最后,對(duì)3種溶質(zhì)運(yùn)移參數(shù)模型進(jìn)行了參數(shù)分析,獲得了3種模型的初步認(rèn)識(shí)。
1 基本原理
研究溶質(zhì)在水流中的遷移規(guī)律,試驗(yàn)測(cè)定可采用電解質(zhì)脈沖示蹤法等來進(jìn)行。電解質(zhì)脈沖在水流中的運(yùn)移規(guī)律符合溶質(zhì)對(duì)流彌散基本方程,根據(jù)夏衛(wèi)生等[12-13]的研究,其溶質(zhì)運(yùn)移參數(shù)隨水流位置而變化,在大規(guī)模的地下水流問題中,溶質(zhì)運(yùn)移參數(shù)亦會(huì)隨空間位置而變化。對(duì)其溶質(zhì)運(yùn)移規(guī)律的模擬可采用近似解析解和數(shù)值解。
3 計(jì)算結(jié)果與討論
根據(jù)上面的理論基礎(chǔ),利用雙精度FORTRAN語(yǔ)言編寫程序,分別采用方程(19)所代表的數(shù)值方法和方程(18)所代表的解析方法求得溶質(zhì)濃度分布,以檢驗(yàn)數(shù)值方法的可靠性。并對(duì)不同溶質(zhì)運(yùn)移參數(shù)變化模式下瞬時(shí)注入溶質(zhì)后水流中溶質(zhì)濃度的時(shí)空分布進(jìn)行了數(shù)值模擬和參數(shù)分析。
3.1 一維溶質(zhì)運(yùn)移參數(shù)線性變化時(shí)方程數(shù)值解與解析解比較
由夏衛(wèi)生等[12-13]的試驗(yàn)結(jié)果得知,坡面薄層水流的溶質(zhì)擴(kuò)散系數(shù)及水流速度均隨位置而變化。本計(jì)算的變?nèi)苜|(zhì)運(yùn)移參數(shù)模型采用D(x)=D0(1+bx/l),u(x)=u0(1+bux/l)。根據(jù)夏衛(wèi)生等[12-13]的試驗(yàn)數(shù)據(jù),取D0=5.1×10-2m2/s,b=-0.55,u0=1.02 m/s和bu=0.042。l=3.5 m,t0=0.6 s, 注入溶液濃度c1=0.63,水流起始溶質(zhì)濃度和平衡濃度c0=0.20,Δt=0.04 sec和Δx=0.025 m,分別求得了水流方向上不同位置處溶質(zhì)濃度的數(shù)值解和近似解析解。
圖1給出了數(shù)值解、近似解析解與文獻(xiàn)[12]實(shí)測(cè)結(jié)果的比較,從圖中可以看出兩種模擬結(jié)果同夏衛(wèi)生等[12-13]的結(jié)果一致,測(cè)量曲線和模擬曲線的溶質(zhì)濃度最大值均呈逐漸減小的趨勢(shì);在距離較
大處其模擬結(jié)果和實(shí)測(cè)結(jié)果吻合較好。近似解析解和數(shù)值解能較好地吻合,特別是距離越大,二者吻合得越好。這說明顯式差分法求解溶質(zhì)運(yùn)移方程的精度已經(jīng)能夠達(dá)到要求。
表1給出了一維常溶質(zhì)運(yùn)移參數(shù)模型在脈沖注入后溶質(zhì)濃度隨深度、時(shí)間變化的數(shù)值解式(20)與解析解式(18)的比較,其中,常溶質(zhì)運(yùn)移參數(shù)D0=5.1×10-2m2/s,u0=1.02 m/s,其余計(jì)算參數(shù)同圖1取值,從表中可以看出二者吻合非常好,說明了數(shù)值方法是可靠的。
由圖1和表1可知,常溶質(zhì)運(yùn)移參數(shù)一維溶質(zhì)運(yùn)移基本方程的解析解與變?nèi)苜|(zhì)運(yùn)移參數(shù)一維溶質(zhì)運(yùn)移基本方程的數(shù)值解基本一致。因此,求解溶質(zhì)運(yùn)移基本方程可以采用本文數(shù)值方法。
3.2 一維溶質(zhì)運(yùn)移參數(shù)多項(xiàng)式變化時(shí)的數(shù)值結(jié)果
給定D0=0.03 m2/s,u0=1.00 m/s,計(jì)算水流位置l=3.5 m,脈沖注入溶質(zhì)濃度時(shí)間t0=0.6 s,注入溶液濃度c1=1.00,水流起始溶質(zhì)濃度和平衡濃度c0=0.20,計(jì)算步長(zhǎng)Δt=0.04 sec和Δx=0.025 m,當(dāng)式(8)和(9)中模型參數(shù)c=0.1,b=±0.6,±0.4,±0.2,對(duì)cu,bu和c,b同步一樣取完全相同值時(shí),其數(shù)值解法求出的溶質(zhì)濃度沿水流位置的分布規(guī)律如圖2所示,圖2中c和b同時(shí)代表cu,bu和c,b的取值。從圖2可以看出,c,cu的零值附近微小變化對(duì)溶質(zhì)濃度分布影響較計(jì)算方法帶來的差別要小。c,cu值一定時(shí),隨著b,bu從-0.6到0.6的增加,某時(shí)刻溶質(zhì)濃度位置上的分布,其峰值和峰寬增加,波峰位置滯后。時(shí)間越長(zhǎng),這種濃度分布的波峰變化規(guī)律越明顯。在1 s時(shí)間以內(nèi),本文對(duì)b,bu 和c,cu 的不同取值對(duì)溶質(zhì)濃度分布影響不大。
3.3 一維溶質(zhì)運(yùn)移參數(shù)冪函數(shù)變化時(shí)的數(shù)值結(jié)果
溶質(zhì)運(yùn)移參數(shù)采用式(10)和(11)的冪函數(shù)變化模型,給定模型參數(shù)b=2且c=±0.9,±0.6,±0.3下,cu,bu和c,b取值完全相同,其余計(jì)算參數(shù)D0,u0,l,t0, c1, c0,Δt和Δx同以上多項(xiàng)式模式中參數(shù)的取值,一維溶質(zhì)運(yùn)移濃度分布數(shù)值解如圖3所示。圖3表明,參數(shù)中指數(shù)b/c,bu/cu較大時(shí),c,cu值的大小對(duì)溶質(zhì)濃度分布的影響較小,時(shí)間越短,這種影響越小。
3.4 一維溶質(zhì)運(yùn)移參數(shù)指數(shù)變化時(shí)的數(shù)值結(jié)果
圖4所示為溶質(zhì)運(yùn)移參數(shù)模型采用式(12)和(13)模型的數(shù)值計(jì)算結(jié)果,計(jì)算中取模型參數(shù)b=0.125且c=±0.9,±0.6,±0.3下,cu,bu和c,b同步取完全相同的值,其余計(jì)算參數(shù)D0,u0,l,t0,c1,c0,Δt和Δx同以上多項(xiàng)式模式中參數(shù)的取值。從圖4看出模型參數(shù)的變化對(duì)濃度分布形態(tài)的影響較小,對(duì)于本例,當(dāng)b,bu一定時(shí),c,cu的增大會(huì)導(dǎo)致濃度分布峰值位置滯后,濃度峰值增加。c,cu一定時(shí),b,bu在一定范圍內(nèi)的變化對(duì)溶質(zhì)濃度分布的影響很小。
4 結(jié)論
本文采用顯式時(shí)間差分格式計(jì)算變?nèi)苜|(zhì)運(yùn)移參數(shù)模型一維水流溶質(zhì)運(yùn)移過程,計(jì)算結(jié)果通過解析解進(jìn)行驗(yàn)證,并和試驗(yàn)結(jié)果進(jìn)行對(duì)比分析,說明算法可靠,顯式法列式簡(jiǎn)單,更適于復(fù)雜模型計(jì)算[1];模型越復(fù)雜,本法的計(jì)算效率越高。本文首次進(jìn)行了變?nèi)苜|(zhì)運(yùn)移參數(shù)模型溶質(zhì)運(yùn)移問題的數(shù)值計(jì)算,并建立不同的溶質(zhì)運(yùn)移參數(shù)模型,對(duì)該模型進(jìn)行了參數(shù)分析。提出的系列算法和處理技巧可以應(yīng)用到溶質(zhì)注入水流后,變?nèi)苜|(zhì)運(yùn)移參數(shù)模型的溶質(zhì)遷移規(guī)律研究,通過建立不同模式下溶質(zhì)運(yùn)移參數(shù)模型,從而實(shí)現(xiàn)復(fù)雜條件下溶質(zhì)運(yùn)移過程的預(yù)測(cè)和模擬。
參考文獻(xiàn):
[1] 邵明安, 王全九, 黃明彬.土壤物理學(xué)[M]. 北京: 高等教育出版社, 2006.
[2] 徐玉佩. 地下水中溶質(zhì)運(yùn)移數(shù)值計(jì)算(有限元法)及其在沖洗種稻改良鹽堿地中的應(yīng)用[J], 水利學(xué)報(bào), 1983,14(4):1-15.
[3] 徐玉佩. 地下水中溶質(zhì)運(yùn)移數(shù)值計(jì)算的沿流線分段追蹤法[J], 水利學(xué)報(bào), 1986, 17(12):38-47.
[4] 黃康樂. 求解二維飽和-非飽和溶質(zhì)運(yùn)移問題的交替方向特征有限單元法[J].水利學(xué)報(bào), 1988,19(7):1-7.
[5] 張效先, 張志剛, 李法虎. 計(jì)算地下水中溶質(zhì)運(yùn)移的時(shí)間分裂Eulerian-Lagrangian有限差分法[J].水利學(xué)報(bào),1991,22(8):35-45.
[6] 馮紹元. 解地下水中溶質(zhì)運(yùn)移問題的特征有限元法[J].武漢水利電力學(xué)院學(xué)報(bào), 1991,24(4):410-417.
[7] 姚磊華. 三維溶質(zhì)運(yùn)移問題的分步廣義迎風(fēng)解法[J].長(zhǎng)春地質(zhì)學(xué)院學(xué)報(bào),1997,27(4):424-428.
[8] 邱克儉, 戚隆溪, 秦寧. 飽和非飽和土壤中溶質(zhì)運(yùn)移的數(shù)值模擬[J].力學(xué)學(xué)報(bào), 1993,25(1):40-47.
[9] 任理,李春友,李韻珠. 粘性土壤溶質(zhì)運(yùn)移新模型的應(yīng)用[J].水科學(xué)進(jìn)展,1997,8(4):321-328.
[10]李煥榮, 羅振東, 謝正輝, 等. 非飽和土壤水流問題的廣義差分法及其數(shù)值模擬[J].計(jì)算數(shù)學(xué),2006,28(3):321-336.
[11]李煥榮, 羅振東. 非粘性土壤中溶質(zhì)運(yùn)移問題的守恒混合有限元法及其數(shù)值模擬[J].計(jì)算數(shù)學(xué), 2010,32(2): 183-194.
[12]夏衛(wèi)生, 雷廷武, 趙軍,等. 薄層水流速度測(cè)量系統(tǒng)的研究[J]. 水科學(xué)進(jìn)展, 2003, 14(6): 781-784.
[13]夏衛(wèi)生, 雷廷武, 張晴雯,等. 坡面薄層水流中電解質(zhì)脈沖遷移模型[J].水利學(xué)報(bào),2003,34 (11),90-95.
[14]段德宏, 王根霞, 王萍. 一維穩(wěn)態(tài)河流水質(zhì)的計(jì)算機(jī)模擬[J]. 水資源與水工程學(xué)報(bào), 2005, 16(3): 54-56.
[15]徐文彬, 徐維生, 方濤. 溶質(zhì)運(yùn)移數(shù)值模擬彌散問題再探[J]. 災(zāi)害與防治工程, 2007, 1(6): 56-61.
[16]何丕文. 有限差分法在河流水質(zhì)預(yù)測(cè)中的應(yīng)用[J]. 長(zhǎng)江大學(xué)學(xué)報(bào)(自然科學(xué)版),2006,3(1): 38-40.
[17]馬東豪,王全九. 土壤溶質(zhì)遷移的兩區(qū)模型與兩流區(qū)模型對(duì)比分析[J]. 水利學(xué)報(bào),2004,35(6):92-97.
[18]JAISWAL S,CHOPRA M,SENG-HUAT O,et al. Numerical solution of one-dimensional finite solute transport system with first type source boundary condition[J].International Journal of Applied and Computational Mathematics,2016,3(4):1-11.
[19]SILAVWE D D,BRINK I C,WALLIS S G. Assessment of some numerical methods for estimating the parameters of the one-dimensional advection-dispersion model[J]. Acta Geophysica,2019,67(3):999-1016.
[20]JAISWAL S,CHOPRA M. Numerical solution of non-linear partial differential equation for porous media using operational matrices[J]. Mathematics and Computers in Simulation,2019,160:138-154.
[21]CARR E J. New semi-analytical solutions for advection-dispersion equations in multilayer porous media[J].Transport in Porous Media,2020,135(1):39-58.
[22]ALLWRIGHT A, ATANGANA A. Augmented upwind numerical schemes for the groundwater transport advection-dispersion equation with local operators[J].International Journal for Numerical Methods in Fluids, 2018,87(9):437-462.
(責(zé)任編輯:曾 晶)
A Numerical Study on Depth-Relevant Solute
Advection-dispersion Model
CAI Songbai1, SHAO Ming’an*2, QI Libin3
(1.College of Resources and Environmental Science, Hunan Normal University, Changsha 410081, China;2.Institute of Soil and Water Conservation, Chinese Academy of Science, Yangling 712100, China;3.Department of Biochemical Engineering, Sanmenxia Polytechnic, Sanmenxia 472000, China)
Abstract:
Research on modelling of solute transport behaviour in groundwater could serve for environmental quality prediction, controlling and management of groundwater. This work has for the first time studied depth-relevant solute transport parameter model by adopting numerical method. By coding a double precision FORTRAN program, this work has calculated the solution of one-dimensional solute transport model of changing solute transport parameter by instant explicit time domain difference method. Then parameter analysis has been done for the 3 changing solute transport models of groundwater and some conceptual knowledge of the model has been obtained.The results of explicit difference equation in time domain are verified and compared with the analytical solution and experimental results. The results show that the algorithm is reliable, the explicit formula is simple, and it is more suitable for complex model calculation. Finally, the 3 changing solute transport models and its numerical analysis methods presented can be applied for many fields such as 2-D and/or 3-D solute transport problem for modelling and prediction of complicated geological condition of solute transport in groundwater.
Key words:
one dimension solute transport; depth-relevant; changing solute transport coefficient; modelling and simulation; instant explicit difference method
收稿日期:2020-11-11
基金項(xiàng)目:國(guó)家自然科學(xué)基金資助項(xiàng)目(40371060)
作者簡(jiǎn)介:蔡松柏(1962—),男,副教授,博士,博士后,研究方向:土壤物理學(xué),E-mail:songbai_cai@sina.com.
通訊作者:邵明安,E-mail:mashao@ms.iswc.ac.cn.