王澤民 張曲辰 金爽 艾松濤
(武漢大學中國南極測繪研究中心,湖北 武漢 430079)
冰川槽谷是由冰川侵蝕作用過的谷地, 這種地貌的主要成因是冰川侵蝕作用, 冰川槽谷的侵蝕作用主要分為冰川下蝕和側蝕拓寬, 冰川侵蝕作用的不同也會影響冰川槽谷形態發育的差異,由冰川槽谷的橫剖面輪廓判別槽谷形態發育的程度, 進而判斷冰川侵蝕的作用模式是一種研究的新思路。根據槽谷外部輪廓的形狀差異, 冰川槽谷分為V型谷和U型谷[1], 從V型谷到U型谷的演化過程稱為冰川槽谷的完整發育過程[2]。
1959年Svensson[3]建立了冪函數模型, 該模型可以使用冪函數方程對冰川槽谷的外部形態進行數值模擬進而得到冰川槽谷的定量化描述。相較于其他模型, 冪函數模型在比較不同槽谷橫剖面及研究槽谷形態發展中更具優勢[2,4]。但是這個模型存在系統性的缺陷, 冪函數模型參數相同的槽谷剖面, 形態上可能存在的巨大差異,Graf[5]引入的槽谷寬深比(FR)限定了槽谷的寬度與深度比例, 完善了最初的模型。1984年Wheeler[6]對位于阿倫島的Glen Rosa冰川槽谷的兩段分別進行了擬合, 得出了用冪函數描述冰川槽谷橫剖面時, 槽谷的不同坡段冪函數模型參數存在差異的結論。
在冰川槽谷定量研究的發展歷程中, 冪函數模型法自提出以來已被國內外學者廣泛應用、研究。此方法的優勢在于簡化了方程的計算, 除此之外還有懸鏈線模型法及二次多項式方法也可用于冰川槽谷發育的定量研究中。三種常見的方法均可應用于大多數具備V型至U型演化特點的冰川槽谷, 但因懸鏈線模型法方程式復雜不便計算,二次多項式方法無法適用于對稱性差的冰川槽谷的原因, 冪函數模型法成為冰川槽谷形態研究的主流方法。國內學者[7-8]對冪函數模型法已有廣泛應用, 該方法成功應用于天山烏魯木齊河源1號冰川及白馬雪山冰川的研究。
國內外學者[9-12]對Austre Lovénbreen冰川和Pedersenbreen冰川的研究主要集中于冰川物質平衡、冰川地形測量、冰川表面積雪覆蓋率及冰流速研究中, 因 Austre Lovénbreen 冰川和Pedersenbreen冰川的底部溫冰層含水量較大難以直接觀測冰下地形[13], 冰面之下的冰川槽谷發育的形態特征缺乏定量性的分析研究。本次北極山地冰川的數值模擬實驗對冰川槽谷的形態發育進行了基于冪函數模型的定量分析, 計算了冰川槽谷各個橫剖面與冪函數模型擬合的各項參數, 對比分析了兩個冰川槽谷發育過程存在的差異并探究了差異的成因。
北極斯瓦爾巴地區是目前國際冰川研究的熱點地區, 該群島總面積約59250 km2, 其中超過57%的面積被2100多個冰川覆蓋。該區域的大多數冰川都是多溫冰川, 而本文選取的研究對象Austre Lovénbreen冰川和Pedersenbreen冰川是兩個典型的多溫山谷冰川[13](圖1)。

圖1 研究區域。a)北極黃河站、A冰川、P冰川位置示意圖; b)斯瓦爾巴群島示意圖及研究區域所處的位置Fig.1.The study area.a) location map of the Yellow River Station, Austre Lovénbreen and Pedersenbreen; b) location of study region in Svalbard
2004年, 我國首次北極黃河站科學考察后, 將Austre Lovénbreen冰川(下文中簡稱為A冰川)、Pedersenbreen冰川(下文中簡稱為P冰川)列為長期監測冰川。2005年, 科考隊員在兩條冰川上埋設標桿,其中A冰川表面埋設17根, P冰川表面埋設5根, 用于監測冰川物質平衡的變化及高精度GPS測量[14]。2009年, 我國科考隊員利用單頻GPS接收機和探冰雷達(Ground Penetrating Radar, GPR)在兩條冰川表面采集了密集的GPS及GPR數據, 開展了冰面地形和冰下地形的研究[15-16]。
本次實驗所用數據是我國北極科考隊員采集于2009年4月的GPS數據及GPR數據。在數據采集的過程中, 將考察設備固定在自制的木質雪橇上, 前雪橇組安裝GPR主機、接收天線和GPS接收機(圖2), 后雪橇組安裝GPR發射天線, 用雪地摩托車牽引自制的雪橇在冰面上同步采集GPS數據和GPR數據[15,17]。其中GPS數據是由加拿大NovAtel公司生產的SMATY-V1型一體化GPS接收機采集, GPR數據則由加拿大Sensor&Software公司生產的pulseEKKO PRO型探地雷達采集。與此同時使用雙頻GPS對測樁進行多次測量, 獲得精度優于1 cm的冰面流速[12,18]。

圖2 雪地摩托牽引雪橇組工作現場: 木質雪橇搭載的GPR和GPS測量系統Fig.2.Field work by snowmobile and sleds: GPR & GPS surveying instruments were loaded in wooden sleds
在采集過程中, 獲得A冰川表面23500個測點的數據, 測線里程達到了83.2 km; P冰川表面采集了17200個測點的數據[16], 測線里程達到了47.4 km, 具體采集區域及測線分布見圖3。

圖3 北極A冰川、P冰川處GPS/GPR測線分布圖[19]Fig.3.Acquired GPS/GPR points in surveying lines in Austre Lovénbreen and Pedersenbreen, Arctic[19]
在上述實驗所用的GPS數據及GPR數據采集完畢后, 利用TGO軟件和EKKO Deluxe軟件分別對GPS數據和GPR數據進行預處理, 之后利用自主開發的GPRead軟件, 實現GPS數據讀取、平差處理、殘差平滑和精度評定, 以及對GPR數據的深度提取[15,17]。經殘差平滑后, A冰川表面的GPS測線交叉點高程最大互差為1.43 m, P冰川表面的GPS測線交叉點高程最大互差為-0.71 m,GPS交叉點高程平均誤差達到分米級。對實時動態測量數據(RTK)的自然鄰近插值獲取了冰川表面數字高程模型(DEM)。
使用探底雷達對冰層厚度進行測量, 本次數據采集所使用的探地雷達分為兩種頻率類型, 在冰厚小于100 m的絕大多數區域使用頻率為100 MHz的雷達采集冰厚數據; 在冰厚大于100 m的局部地區使用頻率為5 MHz的低頻雷達,避免了冰層底部處于壓融狀態的溫冰無法被高頻信號穿透的問題[13,17]。冰面DEM與探地雷達測量得到的冰厚相減得到冰川底部DEM, 冰川表面、底部的DEM空間分辨率均為10 m[15,17], 數據處理結果如圖4所示。

圖4 經GPS、GPR數據處理的冰川表面、底部地形DEM。a) A冰川表面DEM; b) A冰川冰下DEM; c) P冰川表面DEM; d) P冰川冰下DEMFig.4.Surface DEM and bedrock DEM of glaciers.a) surface DEM of Austre Lovénbreen; b) bedrock DEM of Austre Lovénbreen; c) surface DEM of Pedersenbreen; d) bedrock DEM of Pedersenbreen
冰川槽谷在地貌學的分類中以拋物線形態為主, 本文利用計算拋物線的研究方法[3], 構建冰川槽谷形態的冪函數模型, 假設冰川槽谷的橫截面近似為y=axb, 其中x,y分別為原點到槽谷橫剖面上點的水平距離和垂直距離, 通過確定系數a及指數b可以判斷冰川槽谷的形態和發育情況。
為了方便計算, 對冪函數模型兩邊取對數,將函數y=axb轉化為Y=A+bX, 其中Y=lny,A=lna,X=lnx。之后采用最小二乘法確定線性關系中的常數A與b, 在殘差平方和出現最小值時常數A、b為數值模擬的最優解:

式中,n為觀測值樣本數。
在冰川槽谷的發育過程中, 冰川規模和侵蝕能力決定了槽谷的形態[20-21]。因冰川侵蝕作用在不同位置的差異較大, 僅用一組冪函數模型參數A、b來描述冰川槽谷的橫剖面是不準確的, 所以需要在不同海拔位置提取多個橫剖面, 利用冪函數模型開展分析。冪函數模型需要以各橫剖面最低點為原點, 在原點兩側分別進行擬合, 以獲取對應坡段的冪函數模型參數[6]。
在y=axb中,a值描述谷底的寬闊程度,b值反映槽谷的發育情況及陡峭程度[21]。在數值范圍上b值具有很大的變化范圍, 通常范圍為1≤b≤5, 當b值靠近1時, 冰川槽谷接近于V型谷;b值大于2時, 冰川槽谷更接近于U型谷, 一個冰川槽谷的完整發育過程就是V型谷向U型谷的演化過程[2]。
冪函數模型y=axb作為冰川槽谷定量描述的數學模型具有片面性, 這是由冪函數曲線的性質所決定的, 針對同一冪函數模型, 可完成對形態不同(主要指寬度、深度不同)的冰川槽谷曲線的擬合。因此, 為了更加精確地用數值模型表示冰川槽谷的形態, 需引入其他參數描述冰川槽谷橫截面的幾何形狀, 以避免兩段不同的冰川槽谷回歸曲線相同。
槽谷寬深比保證了基于冪函數模型的冰川槽谷的唯一性, 其本質是使用正比例函數對符合冪函數模型的冰川槽谷形態進行約束, 只有槽谷寬深比與冪函數一起使用才能描述和限定冰川槽谷的形態。槽谷寬深比的引入克服了冪函數模型的局限性, 更加準確地表示了冰川槽谷的橫截面幾何形狀, 并忽視量綱的影響, 精確了先前的冪函數模型[5]。槽谷寬深比(FR)被定義為槽谷深度(D)和谷頂寬度(WT)之比:

FR保證了基于冪函數模型的冰川槽谷的唯一性, 其本質是使用正比例函數對符合冪函數模型的冰川槽谷形態進行約束, 只有FR與冪函數一起使用才能描述和限定冰川槽谷的形態。

式中,ymax為冰川槽谷的槽谷深度D,xmax為谷頂寬度WT的一半。
式中正比例函數用于約束冰川槽谷的形態,FR的幾何意義是正比例函數斜率的一半, 圖5為FR對槽谷形態限制示意圖。

圖5 槽谷寬深比對冰川槽谷形態的限制[5]。紅線為受寬深比影響的正比例函數曲線;藍線為與冪函數模型曲線擬合的冰川槽谷Fig.5.The constraint of Form Ratio (FR) to morphology of glacial valley[5].The red line is constraint function curve affected by FR; the blue line is curve of glacial valley which fitted to the power law model
使用2009年科考隊員所采集的A冰川、P冰川的GPS數據和GPR數據, 分別制作顯示冰川槽谷發育形態的剖面圖[15,17]。在A冰川主流線上選取9個橫剖面, 另外在東部支流截取4個橫剖面, 橫剖面截取見圖6a所示; P冰川則只在主流線截取10個橫剖面, 橫剖面截取見圖7a所示。在剖面的截取過程中沿主流線做大量剖面, 根據剖面形態特征手動提取符合冪函數曲線的剖面,再由冪函數坐標原點附近的沉積情況為判據, 將剖面中冰下地形曲線受沉積物影響嚴重的對象剔除, 避免沉積物覆蓋的影響, 因此在參數的擬合過程中保證了數據的精確性。
截取所得的A冰川、P冰川的各橫剖面數據折線圖分別見圖6b、圖7b, 之后將處理過的GPS數據和GPR數據與取對數的冪函數模型Y=A+bX擬合, 經求解后獲得各橫剖面的b值, 如表1~3所示, A冰川的b平均值范圍為0.7554~1.7289, P冰川的b平均值范圍為1.0250~1.8738, 由于兩冰川槽谷b值范圍均小于2, 可得出兩者更接近于V型谷而非U型谷。通過分析b值的分布情況, 發現在兩條冰川中部雪線附近, 即冰流速最快的區域均出現b值的極大值; 另外, 通過高精度GPS點監測獲取P冰川及A冰川主流線上的流速, 發現P冰川整體冰流速是A冰川整體冰流速的2倍以上[12,18], 且P冰川相較于A冰川在下游b值更大, 證明P冰川槽谷發育情況優于A冰川, 即冰川槽谷發育程度與冰流速之間存在著正相關的關系。此外, 在同一條冰川內部, 隨著冰流速增大,b值也增大, 表明在冰流速越大的位置發育成型為U型谷的可能性越高。A冰川、P冰川沿主流線上各位置流速數據如圖8所示, 該流速數據分別為2009—2015年A冰川、P冰川的平均速度, 圖8a中各標號是A冰川沿主流線自下游至上游布設的6根監測標桿的編號, 圖8b中P1—P5是P冰川沿主流線自下游至上游布設的5根監測標桿的編號, 監測標桿具體位置可參考圖6a、圖7a主流線上黑點位置。

表1 A冰川干流橫剖面冪函數模型參數Table 1.Parameters of power law model in central line of Austre Lovénbreen

圖6 A冰川槽谷剖面圖。a)A冰川的槽谷橫剖面位置示意圖, 黑點為冰流速監測的觀測標桿布設位置; b)A冰川的橫剖面圖, 其中藍色為冰面曲線, 綠色為槽谷形狀Fig.6.Cross-sections in Austre Lovénbreen.a) the location map of cross-section in Austre Lovénbreen, the black sign is the location of monitoring stakes; b) profiles in Austre Lovénbreen, the blue line shows ice surface and the green line shows the shape of glacial valley

圖7 P冰川槽谷剖面圖。a) P冰川的槽谷橫剖面位置示意圖; b) P冰川的橫剖面圖, 其中藍色為冰面曲線, 綠色為槽谷形狀Fig.7.Cross-sections in Pedersenbreen.a) the location map of cross-section in Pedersenbreen; b) profiles in Pedersenbreen,the blue line shows ice surface and the green line shows the shape of glacial valley

圖8 A冰川與P冰川流速對比。a)A冰川主流線速度分布; b)P冰川主流線速度分布Fig.8.The comparison of ice flow velocities between Austre Lovénbreen and Pedersenbreen.a) the velocities alone the central flow line in Austre Lovénbreen; b) the velocities alone the central flow line in Pedersenbreen
FR是冰川槽谷發育形態研究中的一個重要影響因素, 可以通過FR與冪函數指數b的相關性判斷該冰川槽谷發育過程中所受冰川侵蝕的具體作用形式[22-25]。當FR值與b值呈正相關關系時,表明冰川槽谷的形成過程中下蝕作用為冰川侵蝕的主要形式;FR值與b值呈負相關關系時, 表明冰川槽谷的形成過程中側蝕作用為冰川侵蝕的主要形式。
使用A冰川、P冰川各橫剖面b值和FR值作散點圖, 由圖9a可知A冰川b-FR關系總體呈正相關, 冰川槽谷受冰川侵蝕作用主要以下蝕作用為主, 表現為典型的山地冰川模式; 由圖9b可知P冰川中b-FR關系總體呈現一定的負相關關系, 側蝕作用強于下蝕作用; P冰川下游的槽谷發育優于A冰川, 形態更接近于U型谷。

圖9 各橫剖面b-FR關系。a)A冰川各橫剖面b-FR關系散點圖; b)P冰川各橫剖面b-FR關系散點圖Fig.9.Relation between b and FR.a) relation between b and FR of Austre Lovénbreen; b) relation between b and FR of Pedersenbreen
在用冪函數模型表示冰川槽谷的過程中, 參數a、b是描述冰川槽谷發育形態的重要參數, 故a、b的關系對冰川槽谷形態特征具有一定的研究意義[21]。由于a值數值過小, 甚至會小于b值數個數量級,a、b之間線性關系不顯, 故研究a值的對數值A與b值的關系。
通常情況冰川槽谷的冪函數模型中冪函數系數a值范圍為0<a<1, 可得出a的對數化的結果A為負數, 取兩冰川各橫剖面所得數據|A|、b分別作各自的回歸曲線, 如圖10所示。

圖10 |A|與b關系。a) A冰川; b)P冰川Fig.10.Relation between |A| and b.a) Austre Lovénbreen; b) Pedersenbreen
A冰川各橫剖面|A|與b擬合的回歸方程為y=5.3061x-3.9089, 相關系數r=0.951, P冰川各橫剖面|A|與b擬合的回歸方程為y=4.9584x- 3.1028,相關系數r=0.963。根據所得數據可知兩條冰川的|A|與b表現出了很強的線性相關,|A|隨b值增大而增大, 兩條冰川的趨勢相一致; 比較分析兩條冰川的曲線斜率, 發現兩者曲線斜率相近, 表明兩條冰川槽谷形態特征為同一類型, 均為V型谷; A冰川的擬合曲線可以近似看做是P冰川擬合曲線向上平移的結果, 如果取相同b值的情況下, A冰川的|A|值將大于P冰川, 則A冰川的a值更小, 谷底相較P冰川更寬闊。

表2 A冰川東部支流橫剖面冪函數模型參數Table 2.Parameters of power law model in the eastern of Austre Lovénbreen
本文使用2009年北極科考隊員所采集的Austre Lovénbreen冰川、Pedersenbreen冰川的GPS數據及GPR數據進行處理后, 采用冪函數模型, 引用槽谷寬深比等概念對兩冰川的冰川槽谷的發育形態進行系列數值模擬工作, 研究分析后結論如下。
1.A冰川、P冰川各橫剖面所得b值反映冰川槽谷的發育情況, 兩者b值范圍均小于2, 其中A冰川b值范圍為0.7554~1.72895, P冰川b值范圍為1.0250~1.8738, 地貌形態都表現為V型谷, P冰川槽谷相較A冰川槽谷更加接近于U型谷。同一條冰川內,b值與相應位置的冰流速呈正相關關系。
2.引入槽谷寬深比FR后發現, A冰川槽谷剖面b-FR關系總體呈正相關, 由此判斷A冰川槽谷為下蝕作用為主的山地冰川模式。P冰川槽谷剖面的b-FR關系呈現一定的負相關, 由此判斷P冰川側蝕作用強于下蝕作用。
3.取對數后的冪函數模型系數|A|與b呈線性相關關系, 且兩冰川槽谷|A|-b回歸方程的曲線斜率數值差異不大, 再次證明兩條冰川槽谷同為V型谷。從回歸曲線截距判斷, A冰川槽谷|A|-b回歸曲線截距更大, 表明A冰川谷底更加寬闊。
致謝感謝中國北極黃河站考察隊員在本研究野外調查工作中的幫助和支持。