張 俊 獨知行 杜 寧 張顯云
1 山東科技大學測繪科學與工程學院,青島市經濟技術開發區前灣港路579號,266590
2 貴州大學礦業學院,貴陽市吉林新村,550025
目前,約束反演和聯合反演是地球物理和大地測量反演研究的熱點[1-6]。由于反演分析方法本質上是由“現象”推求“原因”,其推求“原因”的能力取決于所建立的模型能多大程度上反映二者之間的真實聯系,因此,解的不適定性即成為反演分析方法的固有特性[7-8]。目前,解決不適定問題的主要手段是建立多種數據的聯合反演模型[1-12],并利用模型參數的某些先驗信息對參數可能的數值區間予以限制來實現的[1-2]。在這一過程中,由于對模型參數進行了一定的約束,從而使反演搜索區間大大縮小,在提高解算效率的同時,也改善了模型參數解的不適定問題。在模型具體解算時,不同數據對問題分析的貢獻是由各類數據在目標函數中的相對權比來決定的。文獻[1]指出,相對權比通過影響模型空間交集的大小,從而影響聯合反演模型空間的大小。相對權比最初基本上采用觀測值初始方差或簡單的試錯法確定,具有一定的主觀性[11-12]。文獻[13]提出將相對權比與模型參數一同作為待反演參數參與聯合反演模型函數極小化求解過程,從而避免主觀性;許才軍等[14-16]討論了聯合反演中各類觀測值相對權比的幾種確定方法,認為赫爾默特方差分量估計是最優方法。以上方法可歸納為4種:第一種,根據數據種類數的簡單平均分配;第二種,相對權比與模型參數一同反演;第三種,利用赫爾默特方差分量估計法原理,對不同數據的先驗精度進行調整,然后依據調整后的精度確定相對權比;第四種,根據先驗信息,基于某些經驗直接給出相對權比值。本文借鑒第二和第四種方法的思想,即采取相對權比與模型參數一起參與模型全局極小化反演過程來確定相對權比,但必須根據數據先驗信息,事先對相對權比給予一定約束,從而得到較為合理的反演分析結果。
大地測量聯合反演是指利用包括大地測量數據在內的兩種或兩種以上的觀測量(如測點位移和主應力資料等)求解未知參數的反演方法[1,6]。為正確反映不同觀測量對模型參數求解的貢獻度,在構造模型時,需要對各類數據賦予一定的權比,恰當的相對權比可使反演結果更合理、更真實。相對權比的數值大小反映了模型對不同觀測量的擬合程度。設有模型:

式中,假定有n類觀測數據,則有λ1+λ2+…+λn=1,且0≤λi≤1,其中λi為第i類觀測數據的相對權比值。m為參數空間,φi(m)為聯合反演模型中第i類目標的子函數,一般取φi(m)=(fi(m)-di)2,其中fi(m)為第i類數據與反演參數之間的正演函數,di為第i類觀測數據。當有可靠精度信息時,φi(m)可取加權形式,即取φi(m)=Pi(fi(m)-di)2。
假定各類數據已經過歸一化等數據預處理,從而建立模型(1),則根據多種數據進行聯合反演的問題,可歸結為求能使模型(1)的目標函數取得極小值的參數m的最優化問題。
建立模型(1)后,大地測量聯合各類數據(如地質、地震、地球物理數據)進行反演分析時,其任務即為選擇、利用某種全局最優化方法,來固定模型(1)的參數空間,并進行相關評價和解釋。本文僅通過數值實驗討論模型參數的獲取問題。后續評價與解釋一般需要結合具體實際才能進行,因此不作深入討論。
模型(1)中目標函數極小值的取得,明顯受相對權比值λ的影響。只有合理地確定各類數據的相對權比關系,才能獲得合理的參數反演解。文獻[13]提出,將模型(1)中的相對權比作為未知參數,與其他模型參數一同參與反演,取得了較好的應用效果。這種方法的優勢在于可避免人為確定相對權比的主觀性,但是,卻忽略了反演的本質是一種純粹的數學方法,有時得不到切合實際的反演結果。如果能夠根據數據本身的某些先驗信息,對相對權比值上下界附加某種約束,既可避免人為選擇的主觀性,又顧及了數據的原本特質。若對參數也有一定的了解,則對參數也附加一定的約束,可望得到較好的反演結果。為此,可在模型(1)的基礎上,對參數和相對權比都施加一定的約束,則整個反演模型即為模型(1)與如下約束模型的復合模型:

式(2a)中,hj(m)表達了對參數的某種約束。h(m)可以是線性也可以是非線性的,當有多個約束時,hj(m)中下標j表達了對參數的多種約束。式(2b)是對相對權比的約束,ci表示第i類數據與模型參數一起作為待反演參數時的下界值。
聯合反演模型相對權比的約束反演值就是在同時滿足模型(1)和模型(2)的條件下所得到的相對權比的反演值,而最終聯合反演模型參數m值就是在相應的相對權比值下的反演值。具體實施過程如下。
第一步,建立模型(1)。模型(1)的建立,關鍵在于能否正確建立各類數據與待反演參數的數學物理描述。一般來講,可直接利用現有的大地測量、地質、地球物理中已建立的成熟模型,如利用基于力學模式的彈性力學模型,附加一定邊值條件,結合有限元方法,來反演物質彈性模量、泊松比、邊界力等。
第二步,對原始數據進行分析,獲取數據和待反演參數的先驗信息。如針對觀測數據,可考慮從數據內外符合精度、是否包含系統誤差等數據質量信息著手,分析數據相對質量好壞,從而為第三步確定相對權比約束提供依據;針對模型參數,可根據已有實驗值或經驗,確定模型(2)中參數或參數函數的約束端點值a、b,如地殼彈性模量或密度等數值一般都有相對可靠的范圍,利用這些信息很容易建立約束模型(式(2a))。
第三步,根據第二步分析結果或已有經驗,對相對權比附加適當約束,即確定模型(2b)中的c值以建立模型(2b),從而對相對權比進行有效的下界約束。
第四步,利用某種全局最優化方法,如模擬退火或改進的模擬退火法、遺傳算法和區間算法等作為反演分析的數學方法,確定模型參數初始值X(0),并設定解的搜索方式,在模型(2)的約束下更新模型參數X(0)和相對權比值λ0,得到相對權比和參數更新值λ(k)和X(k)。
第五步,比較第四步中X(k)與X(0)作用于模型(1)的目標函數值Φ(λ,X),當相鄰兩次目標函數值差值的范數小于事先設定的閾值ε時,停止迭代過程;否則,重復第四、五步,不斷更新X(0)和λ0并得到更新值X(k)和λ(k),直到滿足|Φ(λ(k),X(k))-Φ(λ(0),X(0))|≤ε時為止。迭代終止后,取Xk作為最終反演模型參數解,λ(k)為對應最終相對權比的反演值。
本文算例源于文獻[1]。設某一大地測量反演問題建立的線性反演模型為V=Bx-l,在該模型中,x為待反演的模型參數(某滑坡的幾何變化信息),包括4個分量:x1(長度)、x2(寬度)、x3(深度)和x4(傾斜)。所附加的線性不等式約束為Gx≤W。同時,由監測得到的參數變化范圍是-0.2≤xi≤2(i=1,2,3,4),反演參數的真值為=[-0.07-0.2 0.20 0.37]T,具體參見文獻[1]。這里對原算例進行如下改化:令T=,并將L分為3組,即令L,相應地T=[BTGT],則原問題可簡單地看作由3類數據建立的不同線性反演模型構成的聯合反演問題。鑒于實際中模型和觀測數據都不可避免地含有誤差,因此在3組觀測數據及相應系數矩陣中分別模擬了大小不同的正態隨機誤差(表1)。

表1 分組后的原始數據及模擬誤差Tab.1 Grouped original data and simulation error
考慮到參數滿足約束條件-0.2≤xi≤2(i=1,2,3,4)(相當于模型(2a)中,h(x)=x,a=-0.2,b=2.0情形),采用下列方案進行反演計算,各方案反演結果列于表2(其中為參數反演值,‖Δx‖ =‖x-‖ 為參數反演值與真值之差的范數)。
方案1 3類數據按各自模擬方差定權,相對權比均取為1/3進行聯合反演;
方案2 3 類數據按各自模擬方差定權,式(2)中相對權比下界約束值分別取c1=0.1,c2=1/3,c3=0.1,將相對權比與模型參數一并進行聯合反演;
方案3 3類數據按各自模擬方差定權,式(2)中相對權比下界值分別取c1=0.1,c2=1/3,c3=0.2,將相對權比與模型參數一并進行聯合反演;
方案4 3類數據按各自模擬方差定權,式(2)中相對權比下界值分別取c1=0.2,c2=1/3,c3=0.1,將相對權比與模型參數一并進行聯合反演;
方案5 3類數據按各自模擬方差定權,相對權比與模型參數一并進行聯合反演,但不對相對權比作任何約束,即式(2)中相對權比下界值分別取c1=c2=c3=0。
對表2各方案反演結果進行對比可以看出,根據模擬誤差精度對相對權比進行不同約束反演,得到不同的反演結果。在對相對權比進行約束時,考慮到第二組數據模擬精度最高、第三組次之、第一組最差,在3類數據聯合反演時,充分利用先驗精度信息,考慮3組數據的貢獻度。比如對于第二組數據,由于其精度最高,可考慮使其相對權比具有不小于1/3的貢獻度,而對其余兩組數據進行不同下界閾值的約束實驗。結果顯示,當對相對權比按精度進行符合實際的約束反演時,反演結果有顯著改善(方案3),這是因為方案3對相對權比約束的下界值很好地顧及了3組數據的精度情況,正確反映了各類數據在聯合反演中應有的貢獻度。當相對權比約束不符合實際時,反演結果反而不及將相對權比簡單平均分配的反演結果好。如方案2認為第一組數據和第三組數據至少應該有相同的貢獻,而方案4認為第一組數據應有高于第三組數據的貢獻。顯然,根據模擬誤差先驗信息,這兩種約束都不合理。當對相對權比不作任何約束并將相對權比與待反演參數一起反演時,結果最差(方案5)。并且從結果看,幾乎完全擬合了第一組數據,而第一組數據恰恰是精度最差的。這說明,僅將相對權比簡單地與其他參數一起作為反演參數參與反演計算是一種純粹的數學方法,沒有顧及數據本身的實際情況。

表2 各方案的反演結果Tab.2 The inversion results of each scheme
實驗表明,將相對權比與待反演參數一起作為未知參數參與反演的方法是可行的,但必須根據數據質量情況或地球物理的某些先驗信息,對相對權比進行適當的約束??紤]到通常事先難以對各類數據有精確的了解,這種約束一般也不宜太強。如針對第二組數據,因其精度最高,僅根據模擬誤差精度情況對其相對權比下限進行不低于平均貢獻度的約束,而并未對其上限作進一步約束。事實上,當對數據缺乏足夠了解時,不適當的約束反而會導致結果扭曲。因此,當缺乏數據先驗信息時,聯合反演分析中相對權比如何確定仍然是值得進一步研究的問題。
[1]王樂洋,許才軍.等式約束反演與聯合反演的對比研究[J].大地測量與地球動力學,2009,29(1):74-78(Wang Leyang,Xu Caijun.Comparative Research on Equality Constraint Inversion and Joint Inversion[J].Journal of Geodesy and Geodynmics,2009,29(1):74-78)
[2]王樂洋,朱建軍.附不等式約束的大地測量反演[J].大地測量與地球動力學,2008,28(1):109-114(Wang Leyang,Zhu Jianjun.Geodetic Inversion Theory Constrained with Inequalty[J].Journal of Geodesy and Geodynmics,2008,28(1):109-114)
[3]Zeng Y H,Chen C H.Fault Rupture Process of the 20 September 1999Chi-Chi,Taiwan Earthquake[J].Bulletin of the Seismological Society of America,2009,99(5):1 088-1 098
[4]Ma K,Mori J,Jong L S,et al.Spatial and Temporal Distribution of Slip for the 1999 Chi-Chi,Taiwan Earthquake[J].Bulletin of the Seismological Society of America,2001,91(5):1 069-1 087
[5]Wald D J,Heaton T H.Spatial and Temporal Distribution of Slip for the 1992Landers,California Earthquake[J].Bulletin of the Seismological Society of America,1994,84(3):668-691
[6]Wright T J,Parsons B E,Jackson J A,et al.Source Parameters of the 1 October 1995 Dinar(Turkey)Earthquake from SAR Interferometer and Seismic Body Wave Modeling[J].Earth and Planetary Science Letters,1999,172:23-37
[7]趙少榮,於宗儔,陶本藻,等.動態大地測量反演理論及其若干應用[J].測繪學報,1993,22(4):241-248(Zhao Shaorong,Yu Zongchou,Tao Benzao,et al.Theory of Dynamic Geodetic Inversion and Its Applications[J].Acta Geodaetic et Cartographica Sinica,1993,22(4):241-248)
[8]獨知行,盧秀山.基于力學模式的大地測量反演理論及應用[M].北京:地震出版社,2003(Du Zhixing,Lu Xiushan.Theory and Application of Geodesy Inversion Based on Mechanical Models[M].Beijing:Earthquake Press,2003)
[9]獨知行,盧秀山,溫興水,等.中國大陸主要塊體現今運動狀態聯合反演研究[J].山東科技大學學報:自然科學版,2006,25(1):17-20(Du Zhixing,Lu Xiushan,Wen Xingshui,et al.Joint Inversion on Present Motion State of the Main-Plates in Chinese Continent[J].Journal of Shandong University of Science and Technology:Natural Science,2006,25(1):17-20)
[10]Hartzell S H,Heaton T H.Inversion of Strong Ground Motion and Teleseismic Waveform Data for the Fault Rupture History of the 1979Imperial Valley California,Earthquake[J].Bulletin of the Seismological Society of America,1983,73(6):1 553-1 583
[11]Yagi Y,Kikuchi M.Source Rupture Process of the Kocaeli,Turkey,Earthquake of August 17,1999,Obtained by Joint Inversion of Near-Field Data and Teleseismic Data[J].Geophys Res Lett,2000,27:1 969-1 972
[12]Delouis B,Giardini D,Lundgren P,et al.Joint Inversion of InSAR,GPS,Teleseismic,and Strong-Motion Data for the Spatial and Temporal Distribution of Earthquake Slip:Application to the 1999Izmit Mainshock[J].Bulletin of the Seismological Society of America,2002,92(1):278-299
[13]獨知行,歐吉坤,靳奉祥,等.聯合反演模型中相對權比的優化反演[J].測繪學報,2003,32(1):15-19(Du Zhixing,Ou Jikun,Jin Fengxiang,et al.Optimization Inversion of the Relative Weight Ratioλin the Joint Inversion Models[J].Acta Geodaetic et Cartographica Sinica,2003,32(1):15-19)
[14]段虎榮,張永志.斷層滑動特征與多種反演數據的相對權比關系研究[J].大地測量與地球動力學,2009,29(6):18-21(Duan Hurong,Zhang Yongzhi.On Relative Weight Ratio Relation Between Characteristics of Fault Slip and Various Type of Inversion Data[J].Journal of Geodesy and Geodynamics,2009,29(6):18-21)
[15]王樂洋,許才軍,張朝玉.一種確定聯合反演中相對權比的兩步法[J].測繪學報,2012,41(1):19-24(Wang Leyang,Xu Caijun,Zhang Chaoyu,et al.A Two-Step Method to Determine Relative Weight Ratio Factors in Joint Inversion[J].Acta Geodaetic et Cartographica Sinica,2012,41(1):19-24)
[16]Xu Caijun,Ding Kaihua,Cai Jianqing,et al.Methods of Determining Weight Scaling Factors for Geodetic-Geophysical Joint Inversion[J].Journal of Geodynamics,2009,47:39-46