周黎明,張 楊

(長江科學(xué)院 水利部巖土力學(xué)與工程重點(diǎn)實(shí)驗(yàn)室,武漢 430010)

0 引 言

長江源地區(qū)地處青藏高原腹地,素有“中華水塔”的美譽(yù),是中國重要的生態(tài)安全屏障,同時也是氣候變化的敏感響應(yīng)區(qū)和生態(tài)環(huán)境脆弱區(qū)。長江源地區(qū)平均海拔4 760 m,流域面積約13.82萬km2。長江源地區(qū)水系包括北源楚瑪爾河水系、正源沱沱河水系、南源當(dāng)曲水系以及干流通天河水系。近50 a來,在全球氣候變暖和人類活動的雙重影響下,長江源地區(qū)雪線上升、冰川退縮、水土流失、荒漠化和草地退化等問題日益凸顯,已對長江源區(qū)域乃至整個長江流域的生態(tài)安全和民生經(jīng)濟(jì)安全造成威脅。冰川作為響應(yīng)高原氣候變化的指示器,其厚度的精確探測和淡水資源儲量的準(zhǔn)確計(jì)算對于研究長江源地區(qū)氣候變化和推測長江徑流量具有重要意義。獲取精確的冰下地形對于推斷冰川的發(fā)育和運(yùn)動學(xué)過程具有重要作用。另外,探測長江源地區(qū)的永久凍土位置對于研究長江源地區(qū)水土保持狀況和淡水資源儲量具有重要價值,同時也為青藏高原地區(qū)公路、鐵路和水電工程的安全建設(shè)提供參考資料。

近年來,隨著探地雷達(dá)(Ground Penetrating Radar, GPR)設(shè)備精度和便攜性能的不斷提高,我國利用探地雷達(dá)在多個冰川展開了探測工作,其中多數(shù)為地形崎嶇的山地冰川,并取得了良好的效果。2008年,王璞玉等[1]對天山托木爾峰地區(qū)的青冰灘72號冰川進(jìn)行了GPR多剖面的厚度測量,并使用差分全球定位系統(tǒng)(Global Positioning System,GPS)對探測點(diǎn)位進(jìn)行精準(zhǔn)定位。2009年,吳利華等[2]對天山博格達(dá)峰地區(qū)的四工河4號冰川進(jìn)行了GPR測厚工作,并在地理信息系統(tǒng)(Geographic Information System,GIS)技術(shù)的支持下,采用Co-Kriging插值算法結(jié)合理想塑性體理論對冰川非測厚區(qū)域的厚度進(jìn)行了重建。2014年,王玉哲等[3]使用GPR對祁連山老虎溝12號冰川的中流線及多條橫剖面進(jìn)行了厚度測量,獲取了準(zhǔn)確的冰下地形信息。2018年,李亞楠等[4]結(jié)合雷達(dá)冰川測厚數(shù)據(jù)和遙感信息對東昆侖山煤礦冰川冰儲量進(jìn)行了估算并獲得冰床地形圖。2018年,靳勝強(qiáng)等[5]利用GPR對西藏阿里地區(qū)喀喇昆侖山脈南部的嘎尼冰川進(jìn)行了冰川厚度測量,并結(jié)合差分GPS數(shù)據(jù)和遙感影像數(shù)據(jù)繪制了冰下地形及冰厚分布圖。

土壤層、凍土活動層和多年凍土層的介電性質(zhì)存在明顯差異,為GPR研究凍融界面及多年凍土層上限提供有效地球物理前提。20世紀(jì)70年代,GPR技術(shù)在高寒地區(qū)的淺層地質(zhì)勘察中得以應(yīng)用,隨后在國內(nèi)外逐漸推廣[6]。杜二計(jì)等[7]在祁連山老虎溝和觀山河源頭等地開展了GPR凍土調(diào)查,獲得了多年凍土上限的準(zhǔn)確位置,并分析了凍土上限分布與海拔和氣候變化的關(guān)系。武小鵬等[8]將GPR應(yīng)用在青藏高原共玉高速公路多年凍土工程地質(zhì)勘察中,獲得了準(zhǔn)確的多年凍土上限位置以及沿公路設(shè)計(jì)線路的分布情況,為高原公路施工提供了參考依據(jù)。2018年,單波等[9]利用GPR技術(shù)為青海玉樹、新疆天山等高寒地區(qū)的輸電線路工程地質(zhì)勘察提供準(zhǔn)確的凍土物探解譯結(jié)果。

在國外,探地雷達(dá)已廣泛應(yīng)用于冰川和冰凍物質(zhì)的研究,包括冰川和凍土的含水量、冰川地下的排水路線、冰川和冰凍物質(zhì)的內(nèi)部結(jié)構(gòu)、凍土和冰川冰的破裂和變形、冰川中沉積物等[10]。Navarro等[11]基于探地雷達(dá)冰厚探測數(shù)據(jù),對斯瓦爾巴群島韋德爾-雅爾斯貝格冰川進(jìn)行了冰儲量估計(jì),研究結(jié)果表明對冰川、冷冰和溫帶冰層采用不同的速度進(jìn)行雷達(dá)成果時間-深度轉(zhuǎn)換與采用不同的波速估計(jì)冰儲量的誤差范圍是一致的。

冰川厚度探測方法還有合成孔徑雷達(dá)技術(shù)[12],凍土厚度探測技術(shù)有高密度電法[13]、航空電磁技術(shù)[14]、地震折射波法、地震映像法、瑞雷面波法[15]、瞬變電磁法[16]等。合成孔徑雷達(dá)技術(shù)能較好適用于冰川厚度探測,但不適用于凍土厚度探測?，F(xiàn)有凍土厚度探測方法存在如下弊端:航空電磁技術(shù)成本較高;不能單獨(dú)利用地震折射波數(shù)據(jù)檢測山地多年凍土;地震映像法抗干擾能力弱,勘探深度有限,對凍土上限確定的準(zhǔn)確度略低;瑞雷面波法地表介質(zhì)不均勻或激發(fā)條件差時,勘探深度減小;瞬變電磁法能較好適用于凍土厚度探測,但不適用于冰川厚度探測。綜合其他無損檢測技術(shù),GPR能用于冰川及凍土厚度探測,具有高效、高精度、能非常直觀地反映出地質(zhì)體的幾何位置及形態(tài)的優(yōu)點(diǎn)。

目前,針對長江源地區(qū)冰川和凍土的探測和研究僅限于遙感影像和地質(zhì)情況調(diào)查,對于長江源頭冰川厚度和長江源濕地多年凍土上限的GPR探測尚缺少研究。2022年和2023年7月,長江水利委員會長江科學(xué)院開展了自2012年以來的第12次、13次江源科考,并以江源“冰和碳”為這兩次科考的主要內(nèi)容,采用GPR技術(shù)進(jìn)行冰川厚度探測以及凍土上限調(diào)查。其中,冰川探測位置選取長江正源沱沱河的發(fā)源地格拉丹東主峰冰川,凍土探測位置選擇在查旦濕地。本文設(shè)計(jì)了多種典型冰川和凍土地質(zhì)體模型并進(jìn)行GPR波場數(shù)值模擬,研究了冰川和凍土探地雷達(dá)信號響應(yīng)特征。分析了在長江源地區(qū)獲取的實(shí)測冰川和凍土雷達(dá)數(shù)據(jù),結(jié)合地質(zhì)模型的GPR波場模擬結(jié)果,增加了實(shí)測數(shù)據(jù)的可靠性。對比分析了2022年和2023年冰川及凍土厚度的變化,為今后開展更大規(guī)模的長江源地區(qū)冰川凍土科考奠定了基礎(chǔ)。

1 研究地區(qū)及其物性條件

1.1 格拉丹東主峰

在長江科學(xué)院歷年的江源科考成果中,已有包括對冰川地區(qū)地理特征、地質(zhì)特征、冰面高程變化監(jiān)測等多方面的考察研究基礎(chǔ)。2022年和2023年科考,在以往的基礎(chǔ)上,登上冰川表面對冰川厚度進(jìn)行精準(zhǔn)探測,以掌握長江源地區(qū)冰川的分布狀況。探測位置選擇在長江正源沱沱河的發(fā)源地格拉丹東雪山,位于青海省西南部青藏邊境唐古拉山中段,由一大片南北長達(dá)50余千米、東西寬約30余千米、共50余條山岳冰川群所組成。主峰格拉丹東峰位于91°E、33.5°N,海拔6 621 m,為唐古拉山脈的主峰(圖1)。格拉丹東雪山冰川的伸縮變化明顯地揭示了當(dāng)?shù)啬酥燎嗖馗咴腿驓夂虻淖冞w,極具科研價值。

圖1 格拉丹東主峰冰川地理位置及冰川衛(wèi)星照片

1.2 查旦濕地

2022年和2023年江源科考,還進(jìn)行了長江源地區(qū)凍土調(diào)查研究工作。調(diào)查地點(diǎn)選擇在地下水含量豐富的查旦濕地,也是長江南源當(dāng)曲的發(fā)源地。為了驗(yàn)證凍土探測的準(zhǔn)確性,2022年和2023年科考凍土探測測線選擇在查旦鄉(xiāng)一個常年設(shè)置的地下水監(jiān)測井附近,如圖2所示。

圖2 凍土探測測線位置及地下水監(jiān)測井

1.3 長江源地區(qū)常見介質(zhì)物性條件

長江源地區(qū)常見介質(zhì)的電性特征見表1[17]。

表1 長江源地區(qū)常見介質(zhì)的電性特征[17]

根據(jù)電磁波在介質(zhì)中的波速和旅行時間可以計(jì)算界面深度為

h=vt/2 。

(1)

式中:h表示反射界面深度;t表示傳播時間;v表示電磁波在介質(zhì)中的傳播速度,即

v=c/(εr)1/2。

(2)

式中:c為電磁波在真空中的傳播速度,c=0.3 m/ns;εr為介質(zhì)的相對介電常數(shù)。根據(jù)所探測區(qū)域介質(zhì)的不同選擇適當(dāng)?shù)碾娦蕴卣鲄?shù)(表1)可以獲取目標(biāo)體準(zhǔn)確的位置信息。

電磁波在傳播過程中,遇到不同的阻抗界面(介質(zhì)分界面)時將產(chǎn)生反射波和透射波,其反射與透射遵循反射與透射定律,反射波能量大小取決于反射系數(shù)。反射系數(shù)r的數(shù)學(xué)表達(dá)式為

(3)

式中ε1、ε2為反射界面兩側(cè)的相對介電常數(shù)。

2 冰川和凍土GPR波場模擬

2.1 模型建立

本文采用GprMax[18]軟件,基于時域有限差分法,設(shè)置滿足穩(wěn)定性條件的適當(dāng)參數(shù)并結(jié)合PML[19]衰減邊界條件,模擬GPR電磁波在各向同性無限半空間地質(zhì)模型中的傳播過程。本文設(shè)計(jì)地質(zhì)模型尺寸均為10 m×20 m,網(wǎng)格間距為0.01 m,GPR在地表采用剖面法觀測,天線間距1.0 m、探測點(diǎn)距0.2 m,相關(guān)介質(zhì)參數(shù)見表1。圖3中0～2 m深度的界面為直達(dá)波。

圖3 冰川模型及GPR正演記錄

圖3(a)為2層水平層狀冰川模型,其中沿深度方向0.0～9.0 m為凍結(jié)冰層、9.0～20.0 m為冰下基巖,且冰下基巖地形呈水平狀。圖3(b)為2層向右傾斜狀冰川模型,其中冰巖界面的深度從左到右,從9.0 m變化到11.0 m。圖3(c)為3層水平層狀冰川模型,其中沿深度方向0.0～5.0 m為含水的融化冰層、5.0～9.0 m為凍結(jié)冰層、9.0～20 m為冰下基巖。冰下基巖地形呈水平狀;凍結(jié)冰層的相對介電常數(shù)ε3=3;融化冰層的相對介電常數(shù)ε4=80;基巖的相對介電常數(shù)ε5=8。

2.2 GPR波場模擬結(jié)果

圖3(d)為圖3(a)所示模型的GPR正演記錄,其中在深度方向9.0 m處出現(xiàn)明顯反射信號水平狀同相軸。由表1和式(3)可知,由于凍結(jié)的冰層與冰下基巖相對介電常數(shù)差異較大,冰、巖分界面為電磁波強(qiáng)反射界面,在探地雷達(dá)正演記錄中,反射波同相軸出現(xiàn)的位置及形狀與圖3(a)所示模型介質(zhì)分界面一致。

圖3(e)為圖3(b)所示模型的GPR正演記錄,其中在深度方向7.0 m處開始出現(xiàn)明顯反射信號向右傾斜狀同相軸,反射波同相軸出現(xiàn)的位置及形狀與圖3(c)中所示模型介質(zhì)分界面基本一致。圖3(c)和圖3(e)所示傾斜狀冰下地形為冰川槽谷常見地形,多呈“V”形或“U”形。

圖3(f)為圖3(c)所示模型的GPR正演記錄,其中在深度方向5.0 m和15.0 m處各出現(xiàn)強(qiáng)反射信號水平狀同相軸,冰層(融化)和冰層(凍結(jié))界面能較好對應(yīng)理論模型深度,而冰層(凍結(jié))和基巖界面與理論模型深度偏差較大,其原因?yàn)楸鶎?融化)的電磁波傳播速度較小,而進(jìn)入冰層(凍結(jié))波速較大,導(dǎo)致計(jì)算結(jié)果深度偏大。

圖4(a)為2層水平層狀凍土模型,其中沿深度方向0.0～2.0 m為含水土壤層(即季節(jié)性融土層)、2.0～20 m為凍土層(多年)。

圖4 凍土模型及GPR正演記錄

圖4(b)為圖4(a)所示模型的探地雷達(dá)正演模擬記錄,其中在深度方向2.0 m處出現(xiàn)明顯的水平狀反射同相軸。由表1和式(3)可知,由于融土層與凍土層之間存在介電性質(zhì)差異,因此在探地雷達(dá)正演記錄中,反射波同相軸出現(xiàn)的位置及形狀與圖4(a)中所示模型介質(zhì)分界面一致。在凍土模型中,融土層與凍土層之間的電磁波反射界面振幅強(qiáng)弱還與兩層介質(zhì)中含水率有關(guān)。

3 長江源地區(qū)冰川和凍土探測

3.1 探測儀器設(shè)備及基本數(shù)據(jù)處理流程

冰川、凍土科考所使用設(shè)備為加拿大sensor &software公司制造的EKKO PRO型探地雷達(dá)儀,天線主頻為100 MHz。GPR數(shù)據(jù)處理[20]要經(jīng)過零點(diǎn)校正、背景噪聲消除、直流濾波、Laplace濾波、歸一化、調(diào)節(jié)增益、時深轉(zhuǎn)換等處理步驟,得到最終突出地下介質(zhì)特征的波形信息。其中,零點(diǎn)校正主要消除天線離地和天線延時的影響,背景噪聲消除主要去除背景中含有的電磁噪聲信號;直流濾波主要切除電磁波中含有的直流成分;Laplace帶通濾波主要是消除地表環(huán)境產(chǎn)生的噪聲干擾;歸一化是為了在信號強(qiáng)度大的地方(通常為近地表地區(qū))采用小增益,在信號衰減的地方(通常為深部)采用高增益;調(diào)節(jié)增益主要補(bǔ)償電磁波在地下介質(zhì)中傳播時的衰減能量;時深轉(zhuǎn)換是將時間剖面轉(zhuǎn)換為深度剖面。GPR的地質(zhì)解釋是在數(shù)據(jù)處理后得到的雷達(dá)剖面上,根據(jù)反射波組的波場特征,通過同相軸追蹤,確定反射波組的地質(zhì)意義。

3.2 格拉丹東主峰冰川探測

格拉丹東主峰冰川探測所設(shè)探地雷達(dá)參數(shù)如表2所示。綜合考慮了采集的原始數(shù)據(jù)的質(zhì)量和人員設(shè)備安全等因素,設(shè)計(jì)了與格拉丹東主峰冰川中流線夾角呈61.52°的橫剖面測線一條(圖5),測線參數(shù)如表3所示。

表2 冰川探測探地雷達(dá)參數(shù)設(shè)置

表3 冰川探測測線參數(shù)

圖5 冰川探測測線起點(diǎn)位置

格拉丹東主峰冰川GPR探測現(xiàn)場工作由多人完成,分別負(fù)責(zé)儀器搬運(yùn)、點(diǎn)位標(biāo)定、距離測量、雷達(dá)操作、現(xiàn)場記錄和安全保障(圖6)。

圖6 冰川探測工作現(xiàn)場

2022年GPR探測數(shù)據(jù)處理與解譯如圖7所示。由圖7可看出:從冰面深度d=0.0 m至冰下深度d=3.0 m處雷達(dá)信號表現(xiàn)為振幅較強(qiáng),頻帶較低,同相軸連續(xù)性好,在d=3.0 m處存在振幅突變界面,結(jié)合第2節(jié)中圖3(f)所示冰川模型GPR波場模擬記錄,推斷d=0.0～3.0 m的范圍內(nèi)為融化冰層;沿地表測線方向從x=0.0～9.5 m對應(yīng)在深度d=11.0～12.5 m處存在振幅突變界面,振幅明顯增強(qiáng),在該界面上方無明顯雷達(dá)反射信號,符合冰介質(zhì)對電磁波信號弱衰減作用的特征,該界面下方存在復(fù)雜反射信號及多組線性同相軸(多次波),符合冰下基巖的GPR信號響應(yīng)特征,推斷從融化冰層底界面至該強(qiáng)振幅界面中間為凍結(jié)冰層,該強(qiáng)振幅分界面為冰、巖分界面。

圖7 2022年冰川探測探地雷達(dá)數(shù)據(jù)解譯

從圖7解譯結(jié)果中讀取的冰厚信息由測線起點(diǎn)(遠(yuǎn)離冰川中流線一側(cè))的11.0 m逐漸變化為測線終點(diǎn)(靠近冰川中流線一側(cè))的12.5 m。冰下基巖坡度由式(4)計(jì)算得到。

(4)

2022年探測所得冰下基巖坡度為9.08°,地形走勢呈向右傾斜狀,結(jié)合圖3(c)所示冰川模型和圖3(d)所示GPR波場模擬記錄,此次探測數(shù)據(jù)中解譯的冰下地形符合冰川槽谷地形特征。冰川槽谷是冰川長期作用地表的結(jié)果,其形態(tài)受冰川蠕動、作用時間和基巖巖性等因素的影響。不同的槽谷形態(tài)對冰川運(yùn)動速度的影響很大。通常需根據(jù)冰川槽谷形態(tài)(如“U”形、“V”形) 對側(cè)向應(yīng)力進(jìn)行參數(shù)化,從而建立預(yù)測冰川運(yùn)動過程的數(shù)值模型。由于2022年科考受地形因素限制,考慮到冰上作業(yè)的安全問題,橫剖面測線未穿越冰川中流線,尚不能完全掌握格拉丹東主峰冰川冰下地形的整體情況。

2023年GPR探測數(shù)據(jù)處理與解譯如圖8所示,其中:對應(yīng)在深度d=6.0～7.0 m處存在振幅突變界面,振幅明顯增強(qiáng),該界面下方存在復(fù)雜反射信號及多組線性同相軸(多次波),符合冰下基巖的探地雷達(dá)信號響應(yīng)特征,推斷該強(qiáng)振幅分界面為冰、巖分界面。冰川下部基巖不完整,導(dǎo)致基巖裂縫中有冰層,深度d=7.0～8.0 m范圍內(nèi)可能為冰層,深度d=8.0 m以下才為完整基巖層,導(dǎo)致在深度d=8.0 m的位置也存在振幅突變界面。冰川探測的厚度有待進(jìn)一步采用實(shí)測資料加以佐證。

圖8 2023年冰川探測探地雷達(dá)數(shù)據(jù)解譯

3.3 查旦濕地凍土探測

查旦濕地凍土探測所設(shè)GPR參數(shù)如表4所示。 探測測線長度9.0 m。 2022年探測數(shù)據(jù)處理與解譯如圖9所示。 為突出界面分層效果采用波形圖顯示, 其中從土壤表面d=0.0 m至地下d=4.7 m處雷達(dá)信號表現(xiàn)為振幅較強(qiáng), 頻帶較低, 同相軸連續(xù)性較好, 局部波形細(xì)密, 在d=4.7 m處存在振幅突變界面, 在d>4.7 m處電磁波頻率升高, 振幅變小, 且局部波形細(xì)密。 結(jié)合第2節(jié)中圖4(a)所示凍土模型和圖4(b)所示GPR波場模擬記錄, 推斷在d=0.0～4.7 m的范圍內(nèi)為土壤層和季節(jié)性融土層; 在d=4.7 m為凍土上限位置。 凍土探測現(xiàn)場地下水監(jiān)測井中可以清晰看到凍土層表面結(jié)冰, 經(jīng)深度測量得到結(jié)冰位置位于d=4.73 m(此處未考慮氣溫對地下水監(jiān)測井中的凍土層的影響)。 由此驗(yàn)證了科考采用探地雷達(dá)對長江源地區(qū)凍土調(diào)查的有效性。

表4 凍土探測探地雷達(dá)參數(shù)設(shè)置

圖9 2022年凍土GPR數(shù)據(jù)解譯

2023年探測結(jié)果如圖10所示, 從土壤表面d=0.0 m至地下d=3.6 m處雷達(dá)信號表現(xiàn)為振幅較強(qiáng), 頻帶較低, 同相軸連續(xù)性較好;在d=3.6 m處存在振幅突變界面, 在d>3.6 m處電磁波振幅突然減弱, 幾乎消失。 推斷在d=0.0～3.6 m的范圍內(nèi)為土壤層和季節(jié)性融土層; 在d=3.6 m為凍土厚度上限位置。

圖10 2023年凍土探測探地雷達(dá)數(shù)據(jù)解譯

4 結(jié)論和展望

本文基于時域有限差分方法和完全匹配層的吸收邊界條件,通過設(shè)置多種冰川和凍土理論模型,對冰川和凍土進(jìn)行探地雷達(dá)正演模擬及江源科考探測分析,采用GprMax軟件進(jìn)行數(shù)值模擬,得到以下結(jié)論:

(1)不同冰川、凍土模型所對應(yīng)的GPR信號特征同相軸形態(tài)與介質(zhì)分布形態(tài)一致,且反射信號振幅由不同介質(zhì)之間的介電性質(zhì)差異決定,其中融冰層與凍結(jié)冰層之間為較弱電磁波反射界面,凍結(jié)冰層與冰下基巖為強(qiáng)電磁波反射界面,融土層和凍土層之間的電磁波反射界面強(qiáng)度還與兩層之間的含水率有關(guān)。

(2)探測結(jié)果表明,格拉丹東雪山主峰冰川厚度和查旦濕地凍土厚度上限均有不同程度降低,冰川厚度和凍土厚度上限觀測是一個常年累月的結(jié)果,后續(xù)還需要持續(xù)進(jìn)行觀測,積累更多數(shù)據(jù),分析變化趨勢,以估算探測區(qū)域內(nèi)冰儲量,研究氣候變化對冰川的影響。

本文依托2022年和2023年長江源科考對格拉丹東主峰冰川探測和查旦濕地凍土探測,探測解譯結(jié)果均與理論模型研究結(jié)果符合,獲得了冰川測線處準(zhǔn)確冰厚信息、冰下地形信息和凍土上限的準(zhǔn)確位置,驗(yàn)證了GPR方法在長江源地區(qū)進(jìn)行冰川和凍土精準(zhǔn)探測方面的可行性。在未來的江源科考工作中,將進(jìn)一步提高GPR在長江源冰川探測和凍土調(diào)查中的應(yīng)用范圍,獲取精確的三維冰川和凍土資料,估算冰儲量和地下淡水資源儲量,并引進(jìn)航空探地雷達(dá)對受地形因素制約的冰川和凍土地區(qū)進(jìn)行探測試驗(yàn)。