三維數(shù)值流形方法研究及其在地學(xué)中的初步應(yīng)用

武艷強(qiáng)

（中國(guó)地震局地質(zhì)研究所，北京 100029）

（作者電子信箱，武艷強(qiáng)：chdqyw＠126．com）

伴隨地震孕育、發(fā)生與震后調(diào)整過(guò)程的應(yīng)力、應(yīng)變問(wèn)題在地震研究中占有重要地位。以GPS為代表的空間觀測(cè)技術(shù)為大尺度地殼運(yùn)動(dòng)監(jiān)測(cè)提供了高效、穩(wěn)定、精確的測(cè)量結(jié)果，并為研究上述應(yīng)力、應(yīng)變問(wèn)題提供了有效數(shù)據(jù)約束。GPS速度場(chǎng)和應(yīng)變率場(chǎng)結(jié)果既可用于變形分析，又可作為數(shù)值模擬的有效約束和檢驗(yàn)。為了獲取研究區(qū)域可靠的應(yīng)變率場(chǎng)結(jié)果，有必要開(kāi)展GPS應(yīng)變率場(chǎng)計(jì)算方法的對(duì)比研究。另一方面，由于在整個(gè)地震過(guò)程中應(yīng)力、應(yīng)變存在動(dòng)態(tài)演化現(xiàn)象，因此需要研究適合于動(dòng)態(tài)資料約束、大變形模擬、連續(xù)與非連續(xù)變形耦合的三維數(shù)值流形方法?；谏鲜龇治觯疚陌ㄈ缦?個(gè)方面研究?jī)?nèi)容。

（1）通過(guò)模擬與實(shí)際數(shù)據(jù)分析，討論了4種應(yīng)變率場(chǎng)計(jì)算方法的解算精度與抗差性等問(wèn)題，發(fā)展了能夠較客觀反映實(shí)際變形分布的最小二乘配置球面應(yīng)變率解算方法，并與其他幾種常用方法進(jìn)行了對(duì)比分析。

首先，對(duì)大空間尺度（75°～135°E，20°～50°N）模擬數(shù)據(jù)計(jì)算結(jié)果的分析表明，采用1°×1°采樣數(shù)據(jù)及其50%限定（對(duì)1°×1°數(shù)據(jù)進(jìn)行50%數(shù)據(jù)量的離散化采樣，并剔除2個(gè)5°×10°區(qū)域）采樣數(shù)據(jù)作為輸入的情況下，Delaunay三角形方法因噪聲對(duì)解算結(jié)果影響太大不可取，其他3種整體應(yīng)變率場(chǎng)解算結(jié)果具有一致性，但抗差性有所差別。通過(guò)分析理論結(jié)果與附加了不同誤差的計(jì)算結(jié)果的相關(guān)系數(shù)表明，抗差性由好到壞排列如下：最小二乘配置方法、球諧函數(shù)方法、多面函數(shù)方法和Delaunay三角形方法。從輸入數(shù)據(jù)的稀疏程度對(duì)不同應(yīng)變率計(jì)算方法影響程度看，在數(shù)據(jù)采樣率介于2°～1°網(wǎng)格之間時(shí)最小二乘配置方法受數(shù)據(jù)稀疏的影響最小，球諧函數(shù)和多面函數(shù)對(duì)輸入數(shù)據(jù)的密集性要求較高。

其次，中等空間尺度（90°～120°E，25°～40°N）模擬數(shù)據(jù)（1°～0．5°網(wǎng)格）結(jié)果表明，3種整體方法在測(cè)點(diǎn)分布足夠密的情況下均能滿足實(shí)際計(jì)算需求。此時(shí)，3種方法對(duì)附加的誤差敏感程度有一定差別但不顯著，對(duì)比而言最小二乘配置稍強(qiáng)于其他2種方法。通過(guò)對(duì)不同空間采樣數(shù)據(jù)的應(yīng)變率計(jì)算結(jié)果的分析表明，隨著輸入數(shù)據(jù)越來(lái)越稀疏，多面函數(shù)和球諧函數(shù)方法計(jì)算結(jié)果與理論值的相關(guān)性減弱的幅度要快于最小二乘配置方法，表明此2種方法對(duì)數(shù)據(jù)的分布密度要求較高。

最后，通過(guò)對(duì)中國(guó)大陸1999—2004年期間應(yīng)變率計(jì)算結(jié)果的分析，發(fā)現(xiàn)最小二乘配置方法的穩(wěn)定性最高，即使僅用50%的數(shù)據(jù)作為輸入，得到的應(yīng)變率結(jié)果與全部數(shù)據(jù)輸入基本一致，誤差水平也沒(méi)有明顯的變化。球諧函數(shù)方法受點(diǎn)位分布稀疏程度的影響較大，表現(xiàn)為邊緣效應(yīng)的量值有所增大，影響范圍也有所擴(kuò)大。多面函數(shù)方法受測(cè)點(diǎn)稀疏程度的影響較大，表現(xiàn)為計(jì)算的不穩(wěn)定性。另外，多面函數(shù)方法和球諧函數(shù)方法受點(diǎn)位分布的幾何形狀影響較大，比如中國(guó)東北地區(qū)應(yīng)變率計(jì)算結(jié)果的失真。

因此總體而言，從抗差性、邊緣效應(yīng)、誤差分布、穩(wěn)定性角度來(lái)看，最小二乘配置方法最佳。究其原因在于，最小二乘配置的協(xié)方差函數(shù)是對(duì)實(shí)際觀測(cè)數(shù)據(jù)統(tǒng)計(jì)計(jì)算得到的，能夠反映數(shù)據(jù)的真實(shí)分布特征；而球諧函數(shù)的展開(kāi)階數(shù)，多面函數(shù)的光滑因子、核函數(shù)、平差結(jié)點(diǎn)等的選擇都是通過(guò)反復(fù)試算得到的，很難做到對(duì)輸入數(shù)據(jù)的最優(yōu)描述。在輸入數(shù)據(jù)分布密度滿足的情況下，最小二乘配置方法的計(jì)算過(guò)程無(wú)需人工進(jìn)行參數(shù)選擇，即使不同人員進(jìn)行解算也能得到一致的結(jié)果。

（2）通過(guò)單純形積分特性研究、具體的矩陣元素表達(dá)式推導(dǎo)、模擬結(jié)果的理論測(cè)試等過(guò)程，從多角度對(duì)三維數(shù)值流形方法開(kāi)展了研究。

首先，通過(guò)對(duì)單純形積分公式的分析，給出了二維和三維單純形積分的C＋＋算法實(shí)現(xiàn)，并結(jié)合具體實(shí)例描述了單純形積分的求解過(guò)程。在對(duì)中空規(guī)則形狀積分區(qū)域的面積（體積）和重心的計(jì)算值與理論值對(duì)比分析的基礎(chǔ)上，對(duì)單純形積分的精度進(jìn)行了討論。通過(guò)對(duì)比不同邊長(zhǎng)比（10－6～106）情況下計(jì)算值與理論值的差異，討論了積分區(qū)域圖形條件對(duì)積分結(jié)果的影響；通過(guò)對(duì)中空非規(guī)則形狀積分區(qū)域計(jì)算結(jié)果的分析，討論了單純形積分的普適性特征；最后，分析了高階多項(xiàng)式在中空非規(guī)則形狀積分區(qū)域的積分精度。結(jié)果表明，積分項(xiàng)為低價(jià)或高價(jià)多項(xiàng)式情況下，計(jì)算值與理論值的相對(duì)誤差均約為10－15～10－14，圖形條件對(duì)積分結(jié)果的影響極小且不存在系統(tǒng)性。綜合分析認(rèn)為，計(jì)算結(jié)果與理論結(jié)果的差異是由計(jì)算機(jī)的計(jì)算誤差造成的，單純形積分結(jié)果具有高度的穩(wěn)定性和精確性。該積分算法為本文三維數(shù)值流形方法的研究奠定了基礎(chǔ)。

其次，在真實(shí)坐標(biāo)下推導(dǎo)了三維數(shù)值流形方法的所有矩陣的元素表達(dá)式，具體包括單元?jiǎng)偠染仃?、初?yīng)力矩陣、點(diǎn)荷載矩陣、體荷載矩陣、慣性矩陣、速度矩陣、點(diǎn)位移矩陣、接觸矩陣、摩擦矩陣。將所有矩陣裝配起來(lái)即可實(shí)現(xiàn)連續(xù)、非連續(xù)三維數(shù)值流形方法模擬。由于數(shù)值流形方法不做等參變換，因此該方法的點(diǎn)荷載、點(diǎn)位移等約束可以施加到單元內(nèi)任意一點(diǎn)。同時(shí)，在三維數(shù)值流形方法實(shí)現(xiàn)過(guò)程中，添加了多點(diǎn)位移組合約束算法，可實(shí)現(xiàn)相對(duì)運(yùn)動(dòng)約束，其求解策略包括最小二乘方法和強(qiáng)制約束等。

最后，通過(guò)連續(xù)變形和非連續(xù)變形具體計(jì)算實(shí)例分析，驗(yàn)證了固定點(diǎn)矩陣、點(diǎn)荷載矩陣、體荷載矩陣、慣性矩陣、速度矩陣、接觸矩陣、摩擦矩陣等的正確性。16單元懸臂梁模擬結(jié)果與理論結(jié)果的對(duì)比分析表明程序具有高效性，同時(shí)還驗(yàn)證了應(yīng)力逐步累積模塊的有效性；多接觸面剪切破壞實(shí)例模擬了彈性回跳現(xiàn)象，驗(yàn)證了開(kāi)閉迭代算法的有效性，證明現(xiàn)有程序可以在極限平衡條件下開(kāi)展靜－動(dòng)單次轉(zhuǎn)化模擬。

通過(guò)對(duì)單純形積分算法、矩陣元素表達(dá)式推導(dǎo)、關(guān)鍵算法描述、實(shí)例測(cè)試等研究，對(duì)三維數(shù)值流形方法的算法實(shí)現(xiàn)及模擬效果進(jìn)行了討論，驗(yàn)證了程序的有效性和模擬結(jié)果的高精度特性。

（3）汶川地震的發(fā)生與青藏高原的運(yùn)動(dòng)與變形密切相關(guān)，由于巴顏喀拉地塊東向運(yùn)動(dòng)受阻，致使龍門(mén)山斷裂帶積累了大量的應(yīng)變能。因此，開(kāi)展震前區(qū)域變形特征研究有利于更加清楚地認(rèn)識(shí)汶川地震孕育過(guò)程。該部分從主應(yīng)力與主應(yīng)變率方向差異、不同空間尺度的GPS應(yīng)變率場(chǎng)分布特征和多步三維數(shù)值流形方法模擬等3個(gè)方面對(duì)汶川地震前地殼變形特征開(kāi)展了針對(duì)性研究。

首先，利用最小二乘配置方法對(duì)中國(guó)大陸及周邊World Stress Map（WSM）計(jì)劃應(yīng)力方向數(shù)據(jù)和GPS速度場(chǎng)數(shù)據(jù)進(jìn)行了處理，經(jīng)過(guò)數(shù)據(jù)預(yù)處理、經(jīng)驗(yàn)協(xié)方差參數(shù)估計(jì)和擬合推估等過(guò)程，得到了中國(guó)大陸范圍內(nèi)應(yīng)力方向和應(yīng)變率場(chǎng)結(jié)果。數(shù)據(jù)處理過(guò)程中，對(duì)應(yīng)力方向數(shù)據(jù)和GPS速度場(chǎng)擬合殘差的無(wú)偏性和正態(tài)性進(jìn)行了檢驗(yàn)，并對(duì)應(yīng)力方向和主壓應(yīng)變率網(wǎng)格結(jié)果的誤差分布特征及其成因進(jìn)行了分析。在充分考慮數(shù)據(jù)結(jié)果精度的基礎(chǔ)上，討論了中國(guó)大陸主壓應(yīng)力和主壓應(yīng)變率方向特征的異同，結(jié)果表明二者分布總體一致。在青藏地塊東部存在顯著差異，應(yīng)力方向表現(xiàn)為近NE向而GPS主壓應(yīng)變率表現(xiàn)為近EW向，該特征表明青藏地塊東部地區(qū)近年來(lái)東西向擠壓應(yīng)變積累顯著；另一顯著差異區(qū)域位于西域地塊西部，應(yīng)力方向表現(xiàn)為自西向東由NW向NE的轉(zhuǎn)變過(guò)程，GPS主壓應(yīng)變率以近NS向?yàn)橹鳌?/p>

其次，利用前文GPS應(yīng)變率計(jì)算方法對(duì)比分析擇優(yōu)推薦的最小二乘配置球面應(yīng)變率解算方法，基于1999—2007年中國(guó)大陸西部約700個(gè)GPS測(cè)點(diǎn)的速度場(chǎng)結(jié)果，分析了汶川8．0級(jí)地震前的區(qū)域應(yīng)變率場(chǎng)分布特征。GPS主應(yīng)變率場(chǎng)結(jié)果顯示青藏高原主要表現(xiàn)為南北—北北東向擠壓、東西—北西西向拉伸的特性。青藏中東部地區(qū)呈現(xiàn)由南到北、自西向東主壓應(yīng)變率方向逐步向東偏轉(zhuǎn)的有序分布特征。青藏塊體東西向GPS應(yīng)變率場(chǎng)結(jié)果表明，塊體西部（92°E以西）以東西向拉張變形為主；青藏塊體東部（92°E以東、100°E以西）以東西向擠壓為主，該結(jié)果反映了青藏高原物質(zhì)東向流動(dòng)受到了華北—華南地塊阻擋的動(dòng)力學(xué)特征。因此，應(yīng)變能在塊體邊界處積累，有利于汶川地震的發(fā)生。雖然汶川8．0級(jí)地震震源區(qū)不處于剪應(yīng)變率的高值區(qū)，但是處于面應(yīng)變率擠壓區(qū)和東西向擠壓應(yīng)變率高值區(qū)；發(fā)震斷裂在震前存在擠壓兼右旋剪切變形背景，且震源區(qū)西側(cè)的擠壓變形幅度明顯大于東側(cè)；在斷裂帶尺度上，應(yīng)變積累較為緩慢，存在變形趨于極限現(xiàn)象。

最后，通過(guò)多步三維數(shù)值流形模擬方法研究了汶川地震前擠壓應(yīng)力特征和變形特征。結(jié)果表明，最小二乘配置GPS速度場(chǎng)、應(yīng)變率場(chǎng)結(jié)果和三維數(shù)值流形方法模擬結(jié)果一致性較高，驗(yàn)證了三維數(shù)值流形方法模擬結(jié)果的可靠性。GPS主應(yīng)變率場(chǎng)和三維NMM模擬結(jié)果均顯示在1999—2007年間龍門(mén)山斷裂帶南段的應(yīng)變積累速率大于中北段，中北段的擠壓應(yīng)變積累已趨于極限，南段依然可以繼續(xù)進(jìn)行應(yīng)變積累。

汶川地震前地殼變形特征表明，孕震后期龍門(mén)山斷裂帶及其附近區(qū)域的彈性變形已接近極限，與鄰近區(qū)域相比處于顯著弱變形狀態(tài)。龍門(mén)山斷裂帶南段的應(yīng)變積累速率大于中北段的結(jié)果表明，龍門(mén)山斷裂帶中北段的閉鎖程度要大于南段，更有利于地震破裂發(fā)生。因此，在研究GPS應(yīng)變率場(chǎng)分布特征時(shí)，應(yīng)該結(jié)合構(gòu)造變形的背景從不同空間尺度進(jìn)行分析。

總體而言，本文側(cè)重于方法研究，其中GPS應(yīng)變率場(chǎng)既可用于地殼變形分析又可用于檢驗(yàn)?zāi)M結(jié)果的正確性；三維數(shù)值流形方法的研究側(cè)重公式推導(dǎo)、結(jié)果測(cè)試等，以保證算法的正確性為首要任務(wù)。文章最后給出的三維數(shù)值流形模擬實(shí)例較為簡(jiǎn)單，在地學(xué)領(lǐng)域進(jìn)行更廣泛的應(yīng)用尚需開(kāi)展更加深入的研究工作。