管守奎 瞿 偉 蔣 軍

1 武漢大學(xué)測繪學(xué)院，武漢市珞瑜路129號，430079

2 長安大學(xué)地質(zhì)工程與測繪學(xué)院，西安市雁塔路126號，710054

3 武漢大學(xué)GNSS中心，武漢市珞瑜路129號，430079

應(yīng)變場是描述區(qū)域變形的一項重要指標。各國學(xué)者基于GPS觀測資料對應(yīng)變場開展研究，取得許多成果［1－5］。但是當研究區(qū)域僅有少量觀測數(shù)據(jù)時，對應(yīng)變場的精確求解就會變得困難，而最小二乘配置能夠較好地解決觀測數(shù)據(jù)較少的問題。本文提出使用剛體旋轉(zhuǎn)模型（RRM），采用最小二乘配置解法獲得較為充分的虛擬觀測速度值，之后結(jié)合REHSM 方法，解決觀測數(shù)據(jù)量較少時應(yīng)變場的解算問題。

1 最小二乘配置與REHSM 解算應(yīng)變場

1.1 REHSM 解算應(yīng)變場

歐拉定理是現(xiàn)代板塊運動定量描述的基本定理。如果把板塊看成剛性，地球近似為球體，球心看成強制在地球表面運動的剛體板塊運動的固定點，則地殼運動滿足歐拉定理。設(shè)一點的經(jīng)緯度為（λ，φ），則有：

式中，r為地球半徑，Ve、Vn為站心參考系的經(jīng)向和緯向速度，ωe、ωen、ωn為板塊的旋轉(zhuǎn)參數(shù)分量［6］。該式稱為塊體的剛體旋轉(zhuǎn)模型（RRM）。

一個塊體在周圍塊體作用下發(fā)生整體旋轉(zhuǎn)的同時，內(nèi)部也將發(fā)生應(yīng)變。塊體內(nèi)部的應(yīng)變可分解為兩部分［7］：一是塊體的純應(yīng)變，二是由應(yīng)變引起的塊體旋轉(zhuǎn)。表示為：

x、y分別為觀測點至研究區(qū)域中心經(jīng)線與緯線間的平行圈弧長與子午線弧長，北東為正。εe、εen、εn分別是塊體的東西向線應(yīng)變、東西向和南北向之間的剪應(yīng)變、南北向線應(yīng)變。該式稱為塊體的整體旋轉(zhuǎn)與均勻應(yīng)變模型（REHSM）。

由上式可知，對于一個觀測點，ε矩陣不能直接求得。本文對整個研究區(qū)進行柵格劃分，以每個柵格4個頂點的速率信息求出的應(yīng)變參數(shù)來代替柵格中心的應(yīng)變參數(shù)。求得應(yīng)變參數(shù)后，即可進一步計算應(yīng)變場的特征值［8］。最大主應(yīng)變率：

最大剪應(yīng)變率：

面膨脹率：

1.2 最小二乘配置原理

最小二乘平差將所有待估參數(shù)當作非隨機量，或不考慮參數(shù)的隨機性質(zhì)。最小二乘配置除包括經(jīng)典最小二乘的非隨機參數(shù)外，還包含隨機參數(shù)，按照極大驗后估計、最小方差估計或廣義最小二乘原理求定參數(shù)的最佳估值［9］。函數(shù)模型為：

式中，Y′是非觀測向量，在式（6）中其系數(shù)為0，并未包含在l內(nèi)，而是依賴于協(xié)方差相聯(lián)系，在對濾波模型作參數(shù)估計的同時對Y′進行推估。

配置模型的最優(yōu)線性無偏估計為：

在最小二乘配置理論中，關(guān)鍵是求取已知點信號的方差－協(xié)方差矩陣、已知點信號與未知點信號的協(xié)方差矩陣。這兩個矩陣都是由同一個協(xié)方差函數(shù)計算得到的，協(xié)方差函數(shù)選擇恰當與否直接影響估計值的精度。本文分別采用以下6個協(xié)方差擬合函數(shù)進行計算。

1）多項式函數(shù)：

2）高斯函數(shù)：

3）似高斯函數(shù)：

4）希爾沃寧函數(shù)：

5）空間調(diào)和函數(shù)：

Z＊＝Hi＋Hj，Hi、Hj為數(shù)據(jù)點的海拔高，考慮到d＝0的特殊情況，取Z＊＝｜Hi－Hj｜，但并不影響函數(shù)的對稱、正定和協(xié)調(diào)性。該函數(shù)可以作為經(jīng)驗協(xié)方差函數(shù)：

式中，D為協(xié)方差，d為兩測站點間的距離，D（0）為所有觀測點方差的平均值，k為待定系數(shù)。通過不同的協(xié)方差擬合函數(shù)比較，獲得最優(yōu)的速度場擬合信息，為應(yīng)變場的解算提供數(shù)據(jù)支持［10］。

2 計算實例

2.1 速度場構(gòu)建

將位于運城盆地的基于“數(shù)字地震網(wǎng)絡(luò)工程”的19個監(jiān)測點分為兩組，13個為已測點，6個為檢核點。通過6種協(xié)方差擬合函數(shù)，比較獲得最優(yōu)的擬合方式。根據(jù)表1，采用希爾沃寧函數(shù)為協(xié)方差擬合函數(shù)時，E方向上為1.046 4mm／a，N方向上為2.009 6mm／a，兩個方向的殘差最大值均小于其他協(xié)方差擬合函數(shù)的結(jié)果。采用希爾沃寧函數(shù)為速率場構(gòu)建時的協(xié)方差擬合函數(shù)，對所有的觀測速度進行擬合，實測與擬合的速度場如圖1，具體速度殘差見表2。

由圖1可直觀看出，19個監(jiān)測點的實測結(jié)果與擬合結(jié)果大小與方向較為符合。由表2 可看出，E方向上的殘差絕對值最大值為0.775 6 mm／a，N方向為1.218 7mm／a，擬合結(jié)果良好。

表1 速率檢核統(tǒng)計Tab.1 The statistics of velocity check

圖1 實測速率與擬合速率比較Fig.1 Compared the measured velocity and fitted velocity

表2 速率擬合結(jié)果Tab.2 The result of velocity fitting

協(xié)方差擬合函數(shù)確定后，根據(jù)最小二乘配置方法構(gòu)建的速率場信息如圖2。由圖2 可知，速度場擬合結(jié)果比較緊密。研究區(qū)域SN 方向的速率較EW 方向小，整體為水平偏南走向，同武艷強［3］、李延興［6］的結(jié)果基本一致。

圖2 擬合速率場Fig.2 Fitted velocity field

2.2 應(yīng)變場構(gòu)建

有了充足的速度信息，結(jié)合REHSM 與柵格化的解算方法，即可得到研究區(qū)域應(yīng)變場特征參數(shù)的等值線圖。

由圖3～5 可以看出，高值區(qū)域主要分布在WS、EN 區(qū)域，低值區(qū)域主要分布在WN－ES 區(qū)域。高值區(qū)域為張性或應(yīng)變率較為強烈的地區(qū)，數(shù)量級最大達1×10－8，反映這些區(qū)域是應(yīng)變能高積聚區(qū)，發(fā)生中強地震的可能性較大。結(jié)合運城盆地的地質(zhì)背景可知，WS、EN 區(qū)域為中條山斷裂帶，驗證了應(yīng)變參數(shù)求解的有效性。圖6中黑色箭頭為張性劇烈程度增加走向。可以看出，中條山斷裂兩側(cè)向中條山區(qū)域兩側(cè)膨脹加劇，中條山兩側(cè)均受擠壓，在碰撞處向WS、EN 兩個方向釋放能量，導(dǎo)致WS、EN 即中條山斷裂走向區(qū)域張性提升。膨脹區(qū)域在膨脹交匯處，因各自膨脹方向均不相同，地殼運動活躍。據(jù)山西省地震局通報，2010－12－14 11：45在山西省運城市聞喜縣、夏縣交界發(fā)生3.3級地震，震中位于35.3°N、111.1°E。此區(qū)域正是圖6中的橢圓區(qū)域，故而驗證了該區(qū)域地殼運動的活躍性。

圖3 面膨脹率等值線圖Fig.3 The isogram of superficial expansion rate

圖4 最大剪應(yīng)變等值線圖Fig.4 The isogram of shear strain rate

圖5 最大主應(yīng)變等值線圖Fig.5 The isogram of maximum principal strain rate

3 結(jié) 語

本文使用的地殼形變參數(shù)求解方法，能夠有效解決監(jiān)測點數(shù)據(jù)較少時應(yīng)變場的解算問題。通過塊體的剛體旋轉(zhuǎn)模型，采用最小二乘配置以不同的協(xié)方差擬合函數(shù)獲得較優(yōu)的速率場，并同之前的全國性速度場結(jié)果進行比較。再將研究區(qū)域進行柵格化，通過塊體的均勻旋轉(zhuǎn)與應(yīng)變模型，解決區(qū)域地殼應(yīng)變場的求解問題。

圖6 面膨脹率分析Fig.6 The superficial expansion rate analysis

［1］Dumak B S，Kotlia B S，Kumar K，et al.Crustal Deformation Revealed by GPS in Kumaun Himalaya India［J］.Journal of Mountain Science，2014（1）：41－50

［2］江在森，馬宗晉，張希，等.GPS初步結(jié)果揭示的中國大陸水平應(yīng)變場與構(gòu)造變形［J］.地球物理學(xué)報，2003（3）：352－358（Jiang Zaisen，Ma Zongjin，Zhang Xi，et al.Horizontal Strain Field and Tectonic Deformation of China Mainland Revealed by Preliminary GPS Result［J］.Chinese Journal of Geophysics，2003（3）：352－358）

［3］武艷強，江在森，楊國華，等.利用最小二乘配置在球面上整體解算GPS應(yīng)變場的方法及應(yīng)用［J］.地球物理學(xué)報，2009（7）：1 707－1 714（Wu Yanqiang，Jiang Zaisen，Yang Guohua，et al.The Application and Method of GPS Strain Calculation in Whole Mode Using Least Square Collocation in Sphere Surface［J］.Chinese Journal of Geophysics，2009（7）：1 707－1 714）

［4］黃立人，王敏.中國大陸構(gòu)造塊體的現(xiàn)今活動和變形［J］.地震地質(zhì)，2003（1）：23－32（Huang Liren，Wang Min.Present－Day Activity and Deformation of Tectonic Blocks in Chinese Continent［J］.Seismology and Geology，2003（1）：23－32）

［5］瞿偉，張勤，王慶良，等.渭河盆地現(xiàn)今地殼水平形變特征及區(qū)域構(gòu)造活動性［J］.武漢大學(xué)學(xué)報：信息科學(xué)版，2011（7）：830－834（Qu Wei，Zhang Qin，Wang Qingliang，et al.Research on Present Crustal Horizontal Deformation Feature of Weihe Basin and Its Tectonic Activity［J］.Geomatics and Information Science of Wuhan University，2011（7）：830－834）

［6］李延興，李智，張靜華，等.中國大陸及周邊地區(qū)的水平應(yīng)變場［J］.地球物理學(xué)報，2004，47（2）：222－231（Li Yanxing，Li Zhi，Zhang Jinghua，et al.Horizontal Strain Field in the Chinese Mainland and Its Surrounding Areas［J］.Chinese Journal of Geophysics，2004，47（2）：222－231）

［7］許才軍，張朝玉.地殼形變測量與數(shù)據(jù)處理［M］.武漢：武漢大學(xué)出版社，2009（Xu Caijun，Zhang Chaoyu.Crustal Deformation Measurement and Data Processing［M］.Wuhan：Wuhan University Press，2009）

［8］瞿偉.基于GPS速度場分析研究青藏塊體東北緣的形變特征［D］.西安：長安大學(xué)，2008（Qu Wei.The Deformation Characteristic of the Northeastern Margin of Qinghai－Tibet Plateau Constrained by GPS Observations［D］.Xi’an：Chang’an University，2008）

［9］崔希璋，陶本藻.廣義測量平差［M］.北京：測繪出版社，1992（Cui Xizhang，Tao Benzao.Generalized Survey Adjustment［M］.Beijing：Surveying and Mapping Press，1992）

［10］李沖，李建成，瞿偉.基于移動原理的擬合推估模型建立區(qū)域地殼運動速率場［J］.大地測量與地球動力學(xué)，2012（5）：33－36（Li Chong，Li Jiancheng，Qu Wei.Establishing Regional Crustal Movement Velocity Field with Collocation Based Displacement Principle［J］.Journal of Geodesy and Geodynamics，2012（5）：33－36）