雅魯藏布江桑日—墨脫段巖石圈均衡研究以及虛擬水庫誘發地震危險性分析

郭凌冬， 付廣裕， 佘雅文,3， 王振宇， 楊君妍

1 中國地震局地震預測研究所(地震預測重點實驗室)， 北京 100036 2 中國地質大學(北京)地球物理與信息技術學院， 北京 100083 3 南京大學地球科學與工程學院， 南京 210023 4 中國地震局地球物理研究所， 北京 100081

0 引言

雅魯藏布大峽谷位于喜馬拉雅東構造結，該地水力資源豐沛，米林縣至墨脫縣之間河谷適宜建設大型梯級水庫(趙希濤和魏樂軍，2020)，但迄今為止尚未確定建設位置.上述梯級水庫一旦建成，總裝機容量以及年發電量預計均為世界第二，其巨大的蓄水量在地表形成的載荷勢必會對庫底及附近區域的巖石圈均衡造成一定影響.該區域構造背景復雜，特別是在1950年和2017年分別發生過察隅8.6級(Ben-Menahem et al., 1974)和米林6.9級(尹鳳嶺等, 2018)地震(圖1)，地震活動背景強烈.因此，研究雅魯藏布江上大型水庫蓄水導致的巖石圈均衡變化及其可能導致的水庫誘發地震問題，不僅可以對水庫選址提供一定參考，還對我國防震減災工作具有推動作用.

截至2000年，在全球范圍內確定由水庫誘發并且具有一定強度的地震超過90個(Gupta, 2002).迄今為止，由水庫蓄水誘發的最大地震是1967年印度戈伊納6.7級地震(Shashidhar et al., 2019; Gupta, 2002)；位于埃及尼羅河流域的阿斯旺水庫，在1981年誘發了5.6級地震(Gupta, 2002)；我國的新豐江水庫自1962年蓄水后，在水庫下游誘發了6.2級的地震，為我國最大的水庫誘發地震(丁原章,1989).水庫誘發地震的研究大部分是通過庫侖應力的變化、孔隙壓力的擴散以及斷層穩定性的動態變化等進行分析.Talwani等(1984)對斷層孔隙壓力進行計算，發現水庫蓄水導致的斷層孔隙壓力增加是水庫誘發地震的主要原因；Deng等(2010)通過計算水庫蓄水引起斷層庫倫應力的變化，得出紫坪鋪水庫蓄水誘發汶川地震可能性很低的結論；與之相反，Ge等(2009)認為紫坪鋪水庫蓄水使汶川地震的發生提前了幾十年，然而實際的地殼物性參數、斷層幾何形態以及震源深度等因素都具有不確定性，因此該文得出的結果尚未被學界普遍接受.與上述研究不同，本研究將從巖石圈均衡狀態的角度研究水庫建設位置和誘發地震之間的關系.

巖石圈均衡問題，研究人員利用重力/GNSS聯測數據得出了一定的研究成果.高尚華等(2016)利用重力/GNSS聯合觀測數據，提出了計算巖石圈均衡附加力的新方法.需要說明的是，高尚華等(2016)將表征巖石圈均衡狀態的物理量命名為“巖石圈垂向構造應力”，但由于該物理量是一個矢量而不是張量，使用“垂向構造應力”表述容易與巖石圈內部的構造應力場相混淆，因此本文將 “垂向構造應力”這一概念修改為“均衡附加力”.佘雅文等(2017)對巴顏喀拉塊體東緣的巖石圈均衡分布形態進行分析，發現均衡附加力為正值和負值區域分別容易發生淺源地震和深源地震，表明均衡附加力和地震之間存在著某種對應關系.均衡附加力可以表征巖石圈均衡調整的趨勢，大型水庫蓄水后在地表形成的巨大載荷必將對巖石圈均衡附加力產生影響，從而影響巖石圈均衡調整的趨勢.因此，通過均衡附加力的分布狀態尋找適宜于建設大型水庫的區域，對減少生命以及財產損失具有顯著的應用價值.

本研究通過重力/GNSS觀測數據，對雅魯藏布江林芝段附近的均衡附加力進行計算和分析，從均衡附加力的角度給出適合建設大型水庫的位置.首先根據地表重力/GNSS聯測數據計算雅魯藏布江河谷的重力異常場，然后在地殼均一和分層模型下分別計算該河谷的莫霍面深度和均衡面深度，并進一步計算相應模型下均衡附加力的大小和方向.然后對比兩套計算結果，探討均一模型帶來的計算誤差.最后通過均衡附加力的分布形態圈定適宜建設大型水庫的位置，為我國在雅魯藏布江上水庫大壩的選址工作提供參考.

1 雅魯藏布江河谷的重力/GNSS觀測以及數據處理

Fu和She(2017)在雅魯藏布江大峽谷附近建成了一個由107個測站組成的重力觀測網，并進行了一期流動重力觀測.在上述觀測網的基礎上，She和Fu(2020)于2019年進行了擴充觀測，形成一個由408個測站組成的觀測網(圖1b)，為本研究奠定了數據基礎.整個觀測網中，相鄰兩個測站之間的平均間距約為5 km，對區域重力異常場有較強的約束作用.本文主要研究雅魯藏布江桑日—墨脫段的重力異常，所以只采用河谷附近的重力/GNSS觀測數據，其分布如圖1b中綠色三角形所示.需要特別指出的是，雅魯藏布大峽谷地勢險峻、交通不便，難以進行重力觀測，因此大峽谷地區存在數據空白(圖1b).圖1a中，藍色五角星表示地殼運動觀測網絡CMONOC(Crustal Movement Observation Network of China)在喜馬拉雅東構造地區的重力觀測站(祝意青等, 2012).通過上述四個CMONOC重力觀測站，可將喜馬拉雅東構造重力觀測網連接到拉薩絕對重力觀測站，獲得絕對重力控制.

重力數據的觀測過程遵循Fu和She(2017)的思路，使用兩臺精度為0.02 mGal的Burris重力儀(B105和B115)進行同步觀測.為了減少觀測誤差，觀測方式采取了往返測量，即a→b→c→…→c→b→a.為了保證數據的精度，對交叉點以及測線末端點進行了多次測量.使用LGADJ軟件(劉冬至等, 1991)對觀測數據進行平差處理，發現重力數據的總體精度為0.02mGal.同時，通過三個Leica GX1230雙頻GNSS接收器獲取重力觀測點的高程信息，每個觀測點的觀測時間為40 min，采樣間隔為30 s，垂向觀測誤差在分米量級以內.

通過如下改正實現雅魯藏布江河谷布格重力異常的計算：(1)固體潮改正，(2)正常場改正，(3)高度改正，(4)中間層改正，(5)地形改正.其中，中間層改正使用的地殼密度為2.67 g·cm-3.基于ASTER GDEM 2009地形模型(Abrams et al., 2010)，以觀測點為中心將區域劃分為近場、中場和遠場，采用有限元方法對重力數據進行地形改正；并根據Kane(1962)和Nagy(1966)給出的公式，計算觀測點周邊各個區域地形起伏對該點的影響.對于近場區域(0～2 km)，使用間距為1′×1′的網格地形數據進行改正.類似的，對于中場區域(2～20 km)，使用間距為5′×5′的網格；對于遠場區域(20～167 km)，使用間距為10′×10′的網格.雅魯藏布江大峽谷上游地形起伏較小的區域，地形改正量約為20mGal；雅魯藏布江大峽谷下游，墨脫段地形起伏較大，地形改正量達到60mGal.同時，使用上述重力改正方法的1至3步，得到雅魯藏布江河谷的自由空氣重力異常.結果顯示，布格和自由空氣重力異常的均方根誤差均為0.2 mGal左右.

2 雅魯藏布江桑日—墨脫段重力異常剖面

基于圖1所示的重力觀測結果，通過“陸態網絡”的相對重力觀測數據以及拉薩觀測站的絕對重力觀測值，得出雅魯藏布江桑日—墨脫段河谷的自由空氣和布格重力異常場(圖2).延續Fu等(2014)和Fu和She(2017)的表達方式，相對可靠的重力異常值在以觀測點為中心、10 km為半徑所連成的緩沖區內用彩色顯示.圖2a所示的自由空氣重力異常表達的是真實地球相對于參考橢球體的密度變化，與地球表面的地形起伏具有一定相關性.沿著河谷將重力觀測站所在位置展開成一條剖面，可清晰地展示自由空氣重力異常與地形的相關性.由圖2c可以看出，除了米林以西部分地區以外，地形對自由空氣重力異常具有較大的影響.與自由空氣重力異常相比，布格重力異常消除了地形的影響，反映了地球內部的密度變化.因此，布格重力異常與地表地形的相關性較弱(圖2d).圖2c和2d中，沿著河谷的觀測站存在高程起伏現象，而不是完全遵循沿著河谷前進高程逐步降低的規律.其中林芝以東部分觀測站的起伏高達百米量級，該起伏反應了觀測站處的地形的變化，也大體反映了觀測站與附近雅魯藏布江河面高程的差異，與分米級以內的觀測誤差關系極小.

圖1 雅魯藏布江河谷的重力觀測以及區域構造簡圖 (a) 研究區域的地理位置. 藍色五角星表示陸態網絡重力觀測站，紅框為本研究區域； (b) 研究區域構造簡圖.綠色三角形表示雅魯藏布江河谷附近重力觀測點,棕色三角形表示非河谷段觀測點.紅色圓點表示1950年察隅8.6級地震和2017年林芝6.9級地震.藍色線條表示雅魯藏布江.紅色虛線和深紅色的線條分別表示縫合帶和斷裂帶,數據來源于宋鍵等(2011)、丁林和鐘大賚(2013)和Chang等(2015). IYS雅魯藏布江縫合帶; JLF嘉黎斷裂帶; BNS班公—怒江縫合帶; STDS藏南拆離系; V1喜馬拉雅塊體; V2-1 拉薩塊體; V2-2那曲塊體; V2-3墨脫察隅塊體; V3羌塘塊體.Fig.1 The gravity network along the Yarlung Zangbo River and the local tectonic structures (a) Location of the study area. The five-pointed stars denote the CMONOC stations. Red box denotes the study area； (b) Tectonic sketch of the study area. The green and brown triangles denote the gravity stations. The red dots denote the 1950 M8.6 Zayu earthquake and 2017 M6.9 Nyingchi earthquake, respectively. The blue line denotes the Yarlung Zangbo River. The red and crimson lines denote suture zone and fault zone, respectively(after Song et al., 2011; Ding and Zhong, 2013; Chang et al., 2015). IYS denotes the Yarlung Zangbo suture zone; JLF denotes the Jiali Fault; BNS denotes the Bangong-Nujiang suture; STDS denotes the Southern Tibet detachment system; V1 denotes the Himalayan block; V2-1 denotes the Lhasa block; V2-2 denotes the Nagqu block; V2-3 denotes the Mêdog-Zayu block; V3 denotes the Qiangtang block.

據圖2b可知，雅魯藏布江河谷桑日至墨脫段的

圖2 雅魯藏布江桑日—墨脫段的重力異常 (a) 自由空氣重力異常； (b) 布格重力異常； (c) 自由空氣重力異常剖面； (d) 布格重力異常剖面.Fig.2 Gravity anomaly of the Yarlung Zangbo Valley from Sangri to Mêdog (a) The free-air gravity anomaly; (b) The Bouguer gravity anomaly; (c) The profile of free-air gravity anomaly; (d) The profile of Bouguer gravity anomaly.

布格重力異常均為負值，在雅魯藏布大峽谷上游從西到東表現為由小到大再到小的變化趨勢，大峽谷下游墨脫段從東北到西南表現為增大趨勢.上游的布格重力異常在西端最小達到-500 mGal，在米林西側的高值區達到-330 mGal；下游墨脫段最大布格重力異常值達到-350 mGal.雅魯藏布大峽谷上游河谷海拔變化基本平穩(圖2d)，因此米林西側的布格重力異常高值不是由均衡補償引起的.前人研究表明，米林西側既無高密度礦藏(Hu and Yao, 2018)，也無地殼流(Zhao et al., 2013).因此，該地的布格重力異常高值主要反映了莫霍面的相對隆升.

滕吉文等(2006)通過重力異常場對喜馬拉雅東構造結的地殼結構進行研究，發現在亞歐板塊與印度板塊的擠壓碰撞下，東構造結的高密度物質向北移動.Jiménez-Munt等(2008)依據重力異常等數據，計算了西藏地區的巖石圈厚度.Bai等(2013)和Jin等(1996)基于西藏地區重力異常數據，發現印度板塊俯沖到了歐亞板塊下方.Jordan和Watts(2005)通過重力異常等數據確定了印度板塊與歐亞板塊碰撞區巖石圈的有效彈性層厚度.He等(2014)利用布格重力異常資料對藏北的均衡異常進行了研究，約束了藏北歷史演化的地球動力學模型.上述研究均覆蓋了此次的研究區域，但本研究依據密集的高精度地表觀測數據展開，因此具備更高的可信度.

3 雅魯藏布江桑日—墨脫段巖石圈密度結構

依據圖2所示的布格重力異常數據和高程數據，可進一步反演區域地殼厚度分布和上地幔總體結構等深部信息.考慮到重力反演的多解性，利用布格重力異常數據反演巖石圈的密度構造時，需要建立相對可靠的巖石圈密度初始模型.本文以Hu等(2018)基于瑞利波相速度得到的喜馬拉雅東構造結地殼速度結構為基礎，使用速度和密度之間關系的經驗公式(Wang et al., 2018)：

(1)

本文以圖1所示的重力測線西端點為起點，對雅魯藏布江河谷展開成的剖面進行密度結構反演.采用NGA(Northwest Geophysical Associates)開發的GMSYS(Gravity/Magnetic Modeling Software)軟件反演二維地殼密度結構，具體步驟如下：(1)基于地表觀測獲得河谷布格重力異常數據(圖3a)；(2)根據地殼速度結構(Hu et al., 2018)、莫霍面深度(Xu et al., 2013)和沉積層數據(Wang et al., 2014)建立初始地殼密度模型；(3)使用GMSYS正演初始密度模型對應的布格重力異常場；(4)利用非線性參數的最小二乘估計算法(Marquardt, 1963)反演地殼密度結構.經過反復擬合，得到最優的地殼密度結構(圖3b).反演結果表明，雅魯藏布江桑日—墨脫段河谷下部的地殼結構有明顯的分層特征，上地殼密度的變化范圍是2.6 g·cm-3至2.8 g·cm-3，下地殼密度為2.85 g·cm-3，地幔密度為3.2 g·cm-3.上述密度結構均是通過Hu和Yao(2018)給出的區域地殼速度結構推算得來，界面的差異主要是通過重力數據反演得到.

圖3 雅魯藏布江桑日—墨脫段分層地殼模型 (a) 觀測與模擬布格重力異常; (b) 地殼密度結構,紅色三角形為各個城市在測線上的投影位置.Fig.3 Layered density model along the Yarlung Zangbo River from Sangri to Mêdog (a) The observed and simulated Bouguer gravity anomalies; (b) The inversed crustal density structure. Red triangles denote the location of cities.

正演的布格重力異常結果與觀測結果之間擬合程度反映了地殼結構的可靠性(Deng et al., 2014; Zhang et al., 2017).本文正演值與觀測值之間的擬合結果如圖3a所示，雅魯藏布大峽谷上游正演值與觀測值的均方根誤差為1.33 mGal，下游的均方根誤差為3.98 mGal，整條河谷的均方根誤差為1.79 mGal.反演得到的莫霍面深度與前人的研究結果基本一致(Wang et al., 2019; Xu et al., 2013; Peng et al., 2016; She and Fu, 2020)，與Wang等(2019)反演的地殼密度結構也具有很好的一致性.

圖3b顯示，雅魯藏布大峽谷上游的上地殼與下地殼密度偏大，中間地殼的密度偏小，下游墨脫段的地殼密度由上到下呈現增加的趨勢.雅魯藏布大峽谷上游莫霍面從西到東的變化趨勢為“深-淺-深”，在米林西側莫霍面深度達到最小值，下游的莫霍面從東北至西南逐漸升高.雅魯藏布大峽谷上游的地殼在20～40 km深處存在一個較大的低密度層(圖3b)，該低密度層的空間分布與地震波層析成像的低速層以及大地電磁反演的低電阻率帶相一致(Hu and Yao, 2018; Lin et al., 2017; Peng et al., 2016).在含水量、高溫和化學成分的綜合作用下，該地區的中地殼地震波速和電阻率出現異常(Hu and Yao, 2018)，同時也出現地殼低密度異常的現象.

上述地殼結構表明，莫霍面在米林西側存在明顯的上凸.該地區位于構造單元交界處，米林西側的雅魯藏布江縫合帶(IYS)基本呈東西走向，米林東側急轉至東北西南走向，藏南拆離系(STDS)與雅魯藏布江縫合帶(IYS)在米林處交會(圖1).沿著雅魯藏布大峽谷上游的重力測線在米林處與雅魯藏布江縫合帶(IYS)交會，從拉薩塊體穿過縫合帶到達喜馬拉雅塊體.布格重力異常在米林附近呈現高值，因此反演得到的莫霍面呈現出隆起的趨勢.米林西側的雅魯藏布江斷裂是右旋走滑斷裂(宋鍵等，2011)，米林東側的東久—米林斷裂為左旋走滑斷裂(Zhang et al., 2004)，藏南拆離系(STDS)為一系列的低角度正斷層(Burchfiel et al., 1992).Chang等(2015)發現喜馬拉雅東構造結的巖石圈形變具有垂直連貫性，地下塊體的運動與地表運動具有相同的趨勢.通過米林附近斷裂帶性質可以推斷地下塊體的運動趨勢，認為米林南部地體分別向東西移動，使得地殼破裂變得易碎并形成了地幔上涌通道，導致此處莫霍面抬升.該解釋也得到了前人研究的證實，Peng等(2016)對喜馬拉雅東構造結的地震P波速度異常進行研究，發現米林附近的地幔物質存在上涌特征.Dong等(2016)對該區域的大地電磁進行了研究，該研究認為在印度板塊的俯沖作用下深部熱地幔物質通過薄弱區域上涌，流體狀物質的上涌和橫向擴散降低了局部黏度，逐漸將地殼從薄弱區域推離，于是地殼呈現出“東-西”伸展的趨勢，該過程類似于洋中脊的擴張.綜上所述，米林西側莫霍面的上凸可歸因于地幔物質侵入縫合帶.

4 雅魯藏布江桑日—墨脫段均衡附加力分析

在地質構造時間尺度下，密度較小但相對堅固的地殼漂浮在密度較大但相對柔軟的地幔之上.當地殼無外力作用而處于均衡狀態時，均衡面和莫霍面處于同一位置(Fu et al., 2014).反之，當地殼受到垂向外力作用時，莫霍面與均衡面發生分離，地殼處于非均衡狀態.此時莫霍面和均衡面之間的剩余物質因殼幔之間的密度差而產生附加浮力，需要巖石圈撓曲變形產生的彈性力來維持平衡.該彈性力即為均衡附加力，原則上與附加浮力大小相等且方向相反，因此可以通過附加浮力來計算巖石圈承載的均衡附加力.均衡附加力可以定量地表征巖石圈均衡狀態，從而有助于確定適合建設水電站的位置.

4.1 均一地殼模型

莫霍面是近地表最大的一個密度不連續面，是介于地殼和地幔之間的一個化學界面，也是地震波速產生顯著跳躍的界面(Fu et al., 2014).莫霍面深度的起伏變化可以通過重力數據反演得到(王謙身等, 2009).為了簡化計算，布格重力異常可以簡化為由布格板(密度層)厚度變化引起的(Bagherbandi, 2012)，而地殼與地幔之間的密度差異為密度層的密度.因此，可以用密度層的變化近似表示莫霍面的變化.假設雅魯藏布江河谷下方的地殼均一且密度為2.67 g·cm-3，該密度相當于大地水準面以上物質的平均密度，則莫霍面深度附近密度層的變化可以通過布格重力異常轉化得到.如果有一個范圍比較大且厚度為h的密度層出現在地表以下深度為H的地方，以Δρ表示此密度層的密度，則該密度層在地表所引起的重力變化可以由如下公式表示(郭俊義, 2001)：

(2)

(3)

同樣，將重力觀測站所在位置展開成一條剖面(圖4)，大峽谷上游從桑日西端至林芝東端，下游從墨脫東北至墨脫西南.圖4b中藍色曲線表示均一模型下的莫霍面深度，大峽谷上游地殼厚度(莫霍面深度)從西到東展現出由厚到薄再到厚的整體分布態勢，下游從東北到西南表現為由厚到薄的分布態勢.在整條測線上，桑日東側莫霍面深度達到最深為68 km，最淺處在米林西側為62 km.

目前存在多種維持巖石圈均衡狀態的密度補償模型，其中最著名的是Airy均衡模型.該模型假設地殼和地幔的密度均為常數，莫霍面的深度與地表載荷成正比(Zhang et al., 2014).基于Airy均衡模型，當巖石圈處于均衡狀態時，莫霍面深度和地表高程之間存在如下關系(郭俊義, 2001)：

(4)

其中，r是山根的厚度，H0是大地水準面以下地殼的平均厚度，ρt是大地水準面以上物質的平均密度(2.67 g·cm-3)，Δρ是殼幔密度差異(0.63 g·cm-3).基于Xu等(2013)的接收函數研究結果，喜馬拉雅東構造結平均地殼厚度H0約為67 km.圖4b中，紅色曲線表示均衡面(均衡條件下的莫霍面)深度，可以看出均衡面深度受地形影響很大；其中大峽谷上游的均衡面深度基本在62～68 km之間變化，下游的均衡面基本在52～57 km之間變化.由于觀測站的高程有變化，導致大峽谷上游最東端均衡面有小范圍變深的異常.重力觀測站總體上沿著河谷布設，但河谷與旁邊道路的高差較大，導致圖2c所示的地形起伏變化.圖2c表示的是重力觀測站的高程而不是河面的高程，因此出現圖2c中下游高程比上游高程高的情況.在Airy均衡模型下地形與均衡面具有很大相關性(郭俊義, 2001)，因此林芝以東均衡面出現了異常變化.

圖4 雅魯藏布江桑日—墨脫段的布格重力異常、均衡面以及莫霍面深度 (a) 實測布格重力異常,紅色三角形表示各個城市在測線上投影的位置； (b) 均一地殼模型下莫霍面和均衡面深度, 藍色線條表示莫霍面深度,紅色線條表示均衡面深度； (c) 分層地殼模型下莫霍面和均衡面深度,圖例同(b).Fig.4 The BGAs, isostatic surface and Moho depths along the Yarlung Zangbo River between Sangri and Mêdog (a) The observed Bouguer gravity anomalies. Red triangles denote the locations of the cities; (b) The depths of the Moho and the isostatic surface under the homogeneous crust model. Blue and red lines denote the Moho and the isostatic surface depth, respectively; (c) The depths of the Moho and the isostatic surface under the layered crust model. Legends are the same as (b).

據圖4b可知，通過布格重力異常計算得到的莫霍面深度和基于Airy模型計算得到的均衡面深度之間存在差異.如果莫霍面和均衡面的深度十分接近或者相等，則該地處于均衡狀態，即地殼的自重等于其所受浮力(Forsyth, 1985).但是，若兩者的深度差異較大，則該地區處于不均衡狀態.二者差異越大巖石圈越不均衡.根據前人的研究(Fu and She, 2017; 佘雅文等, 2017; 高尚華等, 2016)，均衡問題可以轉化為密度界面擾動補償的加載問題.當均衡面和莫霍面深度存在差異時，兩者之間的剩余物質就會產生附加浮力，附加浮力需要巖石圈承載的彈性力去平衡，該彈性力即為均衡附加力.該附加浮力與巖石圈受到的均衡附加力大小相等、方向相反.考慮均衡面和莫霍面深度的不同、殼幔密度的差異，基于均一地殼模型的均衡附加力可表示為(高尚華等, 2016)：

Pvertical=(MGra-MTopo)(ρM-ρC)g

=(MGra-MTopo)Δρg,

(5)

其中，g為重力加速度，ρM為地幔密度，ρC為地殼密度，Δρ為剩余密度(0.63 g·cm-3).通過公式(5)可求出雅魯藏布江河谷桑日—墨脫段的巖石圈均衡附加力,其中正值表示該處受到的均衡附加力向上，負值表示均衡附加力向下.

通過圖5中紅色曲線可知，雅魯藏布大峽谷上游的均衡附加力在-20～20 MPa之間變化.從最西端至測線220 km的位置，均衡附加力基本為正值，220 km以東基本為負值.大峽谷下游的均衡附加力在50～70 MPa之間變化，遠大于上游的均衡附加力.測線180～220 km之間河谷的均衡附加力向上且數值較小，甚至出現了0 MPa的測點，在地殼均一模型下這段河谷相比較于其他河谷更接近均衡狀態.本研究與Fu和She(2017)的計算結果基本一致.大峽谷上游東側均衡附加力計算結果出現急劇減小的異常點，亦為觀測點與河平面海拔出現顯著差異所致.Airy均衡模型下均衡面的深度與地形高度相關(郭俊義, 2001)，在求解均衡附加力時將這一信息帶入了最終結果.

圖5 雅魯藏布江桑日—墨脫段河谷的巖石圈均衡附加力 紅色線條表示基于均一地殼模型的巖石圈均衡附加力,藍色線條表示基于分層地殼模型的均衡附加力, 紅色三角形表示各個城市在剖面上的投影位置.Fig.5 The isostatic additional force of lithosphere along the Yarlung Zangbo River between Sangri and Mêdog The red line represents theisostatic additional force under the homogeneous crustal model, the blue line represents the one under the layered crustal model. Red triangles denote the locations of the cities.

4.2 分層地殼模型

基于均一地殼模型的均衡附加力計算相對簡單，但是分層地殼模型更符合實際，因此本文基于分層地殼模型計算了均衡附加力.與均一地殼模型類似，當某處的均衡面與莫霍面深度相同，則該處巖石圈處于均衡狀態.當莫霍面與均衡面的深度之間存在偏差時，兩者之間的剩余物質就會產生附加浮力，需要巖石圈承載的彈性力去平衡.某處莫霍面與均衡面的偏差越大，巖石圈就越不均衡.基于Airy均衡模型理論，根據地殼分層的密度結構，可以得到在地殼分層情況下的均衡面深度(王振宇等, 2018)：

+Hn-1ρn-HρM]/(ρn-ρM),

(6)

分層地殼模型下，雅魯藏布江河谷的巖石圈均衡附加力剖面如圖5中藍線所示，正值表示均衡附加力向上，負值表示均衡附加力向下.通過對比均一模型與分層模型下求得的均衡附加力，可以發現兩者的變化趨勢基本一致.從測線最西端至測線230 km處的均衡附加力的方向向上，米林段的均衡附加力向下，下游墨脫段的計算結果為正值；從林芝到大峽谷之間，均衡附加力在兩種模型下的計算結果在方向上不一致，但總體變化趨勢基本相同.上述結果與He等(2017)對西藏地區巖石圈均衡狀態的研究結果具有一致性，即朗縣以東附近的巖石圈基本處于均衡狀態，其他地區為非均衡狀態(圖5).

相對于均一地殼模型下的計算結果(圖4b)，分層模型下莫霍面和均衡面深度之間的差異更大(圖4c)，該差異與地殼分層結構相關.兩種模型下均衡附加力的計算結果差異較大，最大差異達到20 MPa.總體來說，基于分層模型的計算結果與前人結論更加接近(Peng et al., 2016; Wang et al., 2019).因此，若希望得到精度較高的巖石圈均衡結果，應使用分層地殼模型.

4.3 雅魯藏布江水庫誘發地震風險分析

水庫誘發地震與大型水庫的蓄排水有關，是一種發生在庫區及其附近的特殊地震活動(周斌等，2010).水庫蓄水會在地表形成載荷，影響周圍巖石圈均衡，因此可能誘發地震.

水庫蓄水在地表形成的載荷對地殼具有向下的加載作用，通過對巖石圈均衡附加力的方向與水庫加載的方向進行比較，可以確定河谷的哪些地方適合建設大型水庫.如圖6a，若水庫建設在巖石圈均衡附加力向上的區域，大量水體的自重會與部分均衡附加力相互抵消，使區域巖石圈更加均衡，從而不容易誘發地震；相反，如果在均衡附加力向下的區域建設水庫(圖6b)，大量水體的自重會使巖石圈變得更加不均衡，從而容易誘發地震.因此，巖石圈均衡附加力向上的區域適宜建設大型水庫.

圖6 水庫蓄水前后巖石圈均衡附加力變化示意圖 (a) 均衡附加力向上時水庫蓄水后均衡附加力的變化； (b) 均衡附加力向下時水庫蓄水后均衡附加力的變化.Fig.6 The variation of the isostatic additional force before and after impoundment (a) The variation of the isostatic additional force after impoundment when the isostatic additional force is upward. (b) The variation of the isostatic additional force after impoundment when the isostatic additional force is downward.

圖5表示巖石圈均衡附加力的大小與方向，在雅魯藏布江米林段的均衡附加力為負值.假如在米林附近建設大型水庫，水庫蓄水在地表形成的載荷會使巖石圈更加不均衡，容易誘發觸發地震.因此，從巖石圈均衡的角度來講，米林附近不適合建設大型水庫.相反，在米林上游約100 km的朗縣附近，巖石圈均衡附加力為正值(圖5)，適宜建設大型水庫，該處水庫蓄水會使巖石圈更趨于均衡，不易觸發地震.因此，若在雅魯藏布江上建設大型水庫，建議建在朗縣附近，以減輕水庫誘發地震的風險.

5 結論

本文首先基于雅魯藏布江河谷的高密度地表重力觀測數據，經過一系列改正，得到桑日至墨脫段的重力異常場.然后分別基于均一和分層地殼模型，通過布格重力異常數據反演河谷底部的莫霍面深度，根據Airy均衡模型計算出河谷區域的均衡面深度.隨后依據莫霍面與均衡面之間的差異計算巖石圈所承載的均衡附加力，并據此對雅魯藏布江未來大型水庫建設伴隨的誘發地震風險進行了討論.得到如下主要結論：

(1)桑日至墨脫之間河谷的布格重力異常在-500 mGal與-330 mGal之間變化.沿著河流的方向，大峽谷上游的布格重力異常呈現出由小到大再到小的變化趨勢，下游則表現出逐漸增大的趨勢.

(2)大峽谷上游中地殼存在低密度層，河谷下游的地殼密度從淺到深逐漸增大.大峽谷上游的莫霍面深度自西向東呈現“深-淺-深”的變化趨勢，最淺處位于米林西側，深度約為52 km.大峽谷下游莫霍面深度自東北至西南逐漸變淺.

(3)基于分層地殼模型的巖石圈均衡附加力計算結果表明，米林附近的均衡附加力為負值，其他地區均為正值.下游均衡附加力自東北至西南呈現出逐漸變小的趨勢，最小值為39 MPa.

(4)米林附近的巖石圈均衡附加力為負值，不適宜建設大型水庫.從巖石圈均衡的角度來看，若將大型水庫建設在朗縣附近，將顯著降低誘發地震風險.