本站小編為你精心準備了不等式約束的非線性共軛梯度法研究參考范文，愿這些范文能點燃您思維的火花，激發您的寫作靈感。歡迎深入閱讀并收藏。

《石油物探雜志》2016年第3期

摘要：

井中激電三維反演存在多解性且計算時間長，制約了其由理論走向實際應用。針對這一問題，采用非線性共軛梯度（NLCG）反演方法，在反演目標函數中直接施加約束條件壓制其解的非唯一性，也就是將介質電阻率取值范圍作為先驗信息和約束條件以懲罰函數法的方式引入到反演目標函數中。與常規三維電阻率反演目標函數相比，增加了不等式約束項的目標函數理論上可以壓制反演的多解性。理論模型的反演結果表明，基于不等式約束的井中激電三維反演方法有效地改善了反演結果的精度，以懲罰函數法施加不等式約束條件的方式現實可行。

關鍵詞：

井中激電；反演；非線性共軛梯度法；不等式約束；懲罰函數法

由SEIGEL［1］理論可知，時間域激發極化的正反演理論建立在電阻率法正反演基礎上，直流電阻率法正反演的發展與激發極化法正反演的發展有著緊密聯系，因此直流電阻率法的反演方法可以方便地應用到激發極化數據反演中。早期的直流電阻率法反演由于計算速度及內存的限制，大多數反演為一維、二維的最小二乘反演［2－4］；隨著計算機的快速發展及數值技術的進步，越來越多的研究針對三維反演［5－10］。反演中網格化的方式、偏導數矩陣的處理方法、反演方法技術的選擇是反演問題關鍵，直接決定著反演的效率和精度。SASAKI［11］利用互換定理計算偏導數矩陣對復雜地電模型進行了反演，但是計算時間依然較長；共軛梯度法能夠避免矩陣的直接求解，只需求解矩陣的向量積，大大節省了對計算時間和內存的需求，基于這個優勢，ZHANG等［12］使用線性共軛梯度法完成了直流電阻率法三維模型的反演，成為了直流電阻率法三維反演的主要方法；在國內，三維電阻率反演研究也取得了較多的研究成果［13－20］。但反演的計算速度、非唯一性以及穩定性仍然是三維反演的主要難題。SASAKI［11］提出的最小構造約束使相鄰網格之間電阻率差異極小，已經成為電阻率反演中常用的約束形式；黃俊革等［18］將體積因子作為光滑約束，使得三維電阻率反演的精度和深部網格的分辨率得到顯著提高；宛新林等［21］對粗糙度矩陣進行改進，使之適用于各種情況下光滑度矩陣的求取，將二階等間距光滑因子進行加權處理，改變光滑因子值大小以控制模型的擬合程度，并用最小二乘正交分解法求解反演方程，改善了三維反演的效果；劉斌等［22］引入自適應加權光滑約束控制深部分辨率，改善了深部網格的反演效果，隨后，討論了不等式約束條件和先驗結構約束條件的加入方式［23－24］，有效地去除了反演成像中的假異常和多余構造。非線性問題的線性化［25］往往會造成反問題的不適定性，因此迭代擬合法正逐漸發展為非線性反演方法，其中非線性共軛梯度法（NLCG）不僅減小了方法的不適定性，而且由于其將偏導數矩陣的計算融入梯度計算，因而提高了計算效率。RODI等［26］和NEWMAN等［27］利用非線性共軛梯度法分別解決了大地電磁二維和三維反演問題；胡祖志等［28］使用非線性共軛梯度法實現了大地電磁擬三維反演，計算過程中采用一維靈敏度矩陣來代替三維靈敏度矩陣，提高了反演速度；翁愛華等［29］和林昌洪等［30］分別使用非線性共軛梯度法解決了可控源頻率測深的三維反演問題；劉云鶴等［31－32］使用非線性共軛梯度法研究了三維頻率域航空電磁法反演；董浩等［33－34］利用非線性共軛梯度法對帶地形的三維大地電磁數據反演進行了研究。我們結合非線性共軛梯度反演方法以及ZHANG等［12］提出的偏導數矩陣計算方法，將介質電阻率取值范圍作為不等式約束以懲罰函數法的方式加入到反演目標函數中進行井中激電三維反演。利用理論模型驗證了基于不等式約束的三維電阻率反演方法的可行性與有效性。

1.基于不等式約束的三維電阻率反演方法

井中激電反演包括電阻率反演和極化率反演兩個階段，極化率反演建立在電阻率反演基礎之上，因此激發極化反演首先需要對電阻率進行反演計算，完成電阻率反演后再進行極化率反演。

1．1目標函數依據正則化理論，最小二乘意義下的電阻率反演的目標函數為：Φ（m）＝Φd（m）＋λΦm（m）＝‖Wd［dobs－F（m）］‖22＋λ‖Wmm‖22（1）式中：F（m）為正演響應函數；m（mi，i＝1，2，…，m）為模型參數；dobs（dj，j＝1，2，…，n）為觀測數據；Wd為N×N階加權矩陣；Wd＝diag（1／σ1，1／σ2，…，1／σn），σi為第i個數據的標準差；λ為正則化參數；Wm為M×M階光滑度矩陣。設反演的觀測參數為單極－單極電位φobs，三維反演參數為各網格單元的電導率值σ，由于這兩個參數變化范圍很大，為了反演的穩定性，通常采用對數來標定反演數據及模型參數，即d＝lnφobs，m＝lnσ。

1．2不等式約束三維電阻率反演中的解通常存在不適定性問題，例如在源點附近會出現奇異點，或者出現負值（采用取對數的方式消除）以及很大的電阻率異常值，為了降低這種不適定性，往往會在程序中人為設定某些參數的上限值，甚至下限值：ρmini≤mi≤ρmaxi（2）式中：mi為第i個網格單元的電阻率；ρmini和ρmaxi分別為第i個網格單元的電阻率下限和上限。這種經驗性的在程序中加入上、下限約束的過程和反演計算本身沒有任何關系，這種方式也不適于復雜的實際反演情況，所以我們以懲罰函數的形式將不等式約束引入到目標函數中很好地解決了三維電阻率反演的不適定問題。將不等式約束（（2）式）加入到常規的電阻率反演目標函數，構造出帶有不等式約束的反演目標函數（不等式約束施加的方式詳見附錄A）：Z（m）＝Φd（m）＋λΦm（m）＋μP（m）＝［d－F（m）］TC－1d［d－F（m）］＋λmTCTmCmm＋μ［min（m－mmin，0）Tmin（m－mmin，0）＋min（mmax－m，0）Tmin（mmax－m，0）］（3）式中：mmin，mmax為每個單元網格電阻率的上、下限。公式（3）即為本文帶不等式約束條件的三維電阻率反演的目標函數，與常規三維電阻率反演目標函數相比，該方法增加了不等式約束項，但表達式簡單、容易理解、易于實現。理論上講，不等式約束對抑制多解性有重要作用，能有效地改善反演效果。

1．3非線性共軛梯度反演非線性共軛梯度反演需要解決的關鍵問題有3個：①目標函數梯度的計算；②最優步長的計算；③預條件因子的選取與計算。1．3．1梯度的求解從NLCG算法的計算過程來看，需要計算目標函數的梯度，目標函數（（3）式）的梯度表達式為：g＝－2JTC－1d［d－F（m）］＋2λCTmCmm＋2μ［min（m－mmin，0）－min（mmax－m，0）］（4）從（4）式中可以看出，梯度g的計算關鍵在于雅克比矩陣J的計算，而且只需要計算雅克比矩陣的轉置與任意一維向量的乘積JTy，因此無需考慮雅克比矩陣的存儲問題，并且每次JTy的計算均可在每次反演計算中的正演計算后一起得到，從而大大加快了反演的計算速度。ZHANG等［12］對JTy計算進行了詳細的分析。

1．3．2步長的確定在極小化目標函數Φ（m（i）＋αu（i））時，需要確定搜索步長α，因為只包含一個變量α，所以可以用全局最優化的一維線性搜索方法，如精確線性搜索法、黃金分割法、插值法和非精確線性搜索法等求取。而通常精確線性搜索花費的時間較長，因此一般采用非精確的一維線性搜索。我們采用充分下降條件和二次插值法來求取反演步長，具體的搜索方式見參考文獻［31］。非精確線性搜索中的Wofe－Powell準則為：Φ（mk＋αkpk）≤Φ（mk）＋c1αkΔΦ（mk）TpkΔΦ（mk＋1＋αkpk）Tpk≥c2ΔΦ（mk）Tp烅烄烆k（5）式中：Φ為正演算子；mk為第k次迭代的模型參數向量；αk為迭代步長；pk為搜索方向；c1，c2為常數，滿足0＜c1＜c2＜1；k為迭代次數。在實際應用中，一般取c1＝0．0001，c2＝0．1000。

1．3．3預條件因子的選擇預條件因子能夠改善系數矩陣的特征值分布，減少其條件數以提高反演的穩定性，在NLCG法中，NEWMAN等［27］認為近似Hessian矩陣的逆作為預條件因子使得搜索方向具有與牛頓方向類似的性質，因此預條件因子M－1越近似Hes－sian矩陣的逆，越能加速收斂。預條件矩陣可以在每步反演迭代中固定或者更新變化，NEW－MAN等［27］認為Hessian矩陣不是一個常數，因此預條件矩陣也應該時刻更新變化；而計算全Hessian矩陣或者近似Hessian矩陣代價較大，因此只需要計算它的對角元素來作為預條件因子，NEWMAN等［27］選擇擬牛頓法中的BFGS算法產生預條件因子，只需要計算并存儲Hessian矩陣的對角元素，其對角元素可以依靠目標函數的梯度和上一次迭代的結果給出［27，33］。我們選取的預條件因子為：Mk＋1＝Mk＋gkgTkgTkuk＋ykyTkαkyTkuk（6）式中：gk＝ΔΦ（mk），yk＝ΔΦ（mk＋1）－ΔΦ（mk）。通常第一次迭代時，M1設為單位矩陣I。按公式（6）更新預條件因子能夠確保Mk＋1為正定且NLCG法中的搜索方向是下降方向。

2.井中激電的三維反演

根據體激發極化理論，假設σi（i＝1，2，…，M）與ηi（i＝1，2，…M）分別表示地下介質的電導率和極化率，由等效電阻率法可知，激發極化總場φ（σ＊）由等效電導率σ＊＝（1－η）σ決定，且視極化率可表示為：ηs＝φ（σ＊）－φ（σ）φ（σ＊）（7）由（7）式可推導出視極化率ηs與地下介質真極化率ηi之間的近似線性關系：ηs＝－∑Mi＝1lnφlnσiηi＝－Jη（8）其中，J剛好是電阻率三維反演中偏導數矩陣的對數形式，因此根據正則化理論，得到最小二乘意義下極化率反演的目標函數為：Φ（m）＝Φd（m）＋λΦm（m）＝‖Wd（ηs＋Jη）‖22＋λ‖Wmη‖22（9）對公式（9）求極小，即令Φη＝0，有：g（η）＝Φη＝2JTWTdWd（ηs＋Jη）＋2λWTmWmη＝0（10）（JTWTdWdJ＋λWTmWm）η＝－JTWTdWdηs（11）可以看出，（11）式滿足Ax＝b的形式，因此（11）式可以直接利用線性共軛梯度法求解。（11）式中各項參數和電阻率反演中的各項參數相同。

3.數值模型反演算例

為了驗證方法的可靠性和計算效率，我們給出了幾個理論模型來進行反演計算。三維非線性共軛梯度反演算法具體流程如圖1所示。需要特別說明的是，本文在目標函數中施加不等式約束均是在電阻率三維反演中實現的，極化率的反演直接利用線性共軛梯度法反演完成。流程圖中各參數含義如下：初始預處理矩陣H0一般取為單位矩陣；收斂系數ε1為正則化因子的閾值；收斂系數ε2為前、后兩次相對擬合差Rms的絕對值，即目標函數值的下降量小于某個閾值時，說明目標函數值不再下降，達到極小值，因為理論數據也不能使數據擬合差降低到0；收斂系數ε3為目標函數的相對擬合差閾值；收斂系數ε4為目標函數的梯度值的閾值，如果目標函數的梯度很小，目標函數值也達到極小值。相對擬合差Rms定義如下：Rms＝∑Mj＝1dobsj－d0jdobsj2槡M（12）式中：M是觀測數據的個數；dobsj和d0j分別為合成的觀測數據與正演理論值。值得注意的是，在大多數情況下用公式（12）的值來表示三維電阻率反演精度［35］，而屈有恒［36］認為其值的大小只對反演的過程有意義，不能用來衡量反演精度，任何一個給定的初始模型的反演精度的估算都是相對的，因此對于給定的初始模型，相對擬合差只要能下降到初始模型誤差的1／20～1／10即可認為達到了較高的精度。

3．1算例一：2個低阻立方體模型首先建立三維反演模型，然后利用常規反演目標函數（公式（1））對電阻率進行非線性共軛梯度反演計算，再進行極化率的反演。模型如圖2所示，在電阻率為ρ0＝100Ω•m，極化率為η0＝0的均勻半空間中存在2個大小為20m×20m×20m（邊長a＝20m），電阻率ρ1＝ρ2＝10Ω•m，極化率η1＝η2＝0．1的低阻體，頂面埋深均為h＝5m，底面埋深25m，異常體邊界離x軸與y軸距離相等且距離d＝5m，紅色圓點代表井口所在位置，井口坐標為（50m，50m，0）。分別在井中不同深度布設多個源：A0（z＝0），A1（z＝－10m），A2（z＝－20m），A3（z＝－30m），A4（z＝－40m），A5（z＝－50m），在這些源處固定A極供電；觀測電極距2m，在地面51條測線上逐點移動M極，通過三維有限元正演程序共得到15606個一次場電位數據和視極化率的“實測數據”。為了提高反演速度，將整個區域分為邊界區域和目標區域，目標區域選為100m×100m×100m，擴展區域分別向3個方向延展200m。反演初始模型選取電阻率ρ＝100Ω•m，極化率為η＝0的均勻半空間。反演參數選擇如下：正則化因子初始值l0＝0．05，減小系數q＝0．1，給定反演終止收斂系數ε2＝ε3＝1．0×10－5。圖3為模型的電阻率與極化率反演結果。圖4為y＝40m和y＝60m處的電阻率反演切片（xoz垂直斷面）。圖5為y＝40m和y＝60m處的極化率反演切片（xoz垂直斷面）。低阻異常體反演的電阻率數值基本接近真實值，高阻體差異較大；高阻體的極化率數值接近真實值；反演結果基本圈定了高、低阻異常體的水平中心位置和邊界，反演的高低阻異常體的中心位置都有不同程度的向上偏移現象。為了對比施加約束條件前后的反演效果，我們應用基于不等式約束的反演目標函數（（13）式）對電阻率進行非線性共軛梯度反演計算。首先，利用有限元法對模型進行正演模擬，井中點電源供電點一共11個，點源間隔為5m，觀測電極距2m，在地面51條測線上逐點移動M極，通過三維有限元正演程序共得到28611個視電阻率的“實測數據”。以點電源深度為縱軸，得到視電阻率斷面圖（圖6），圖6中虛線方框分別為低阻異常體在xoz平面上的投影。由視電阻率斷面圖可以看到觀測數據中的視電阻率基本在120Ω•m以內，因此將整個模型中網格單元的電阻率上限設為150Ω•m（有一定冗余），下限設為0。即mmin＝0，mmax＝150Ω•m。為了加快反演計算速度，提取A0（z＝0），A1（z＝－10m），A2（z＝－20m），A3（z＝－30m），A4（z＝－40m），A5（z＝－50m）6個供電點處15606個一次場電位數據進行反演。計算區域與模型大小的設置、反演參數與單個低阻立方體模型一樣，其中懲罰因子μ取為常數1，反演結果見圖7至圖9。由未施加不等式約束條件下的反演結果（圖3至圖5）可看出，圖4電阻率反演切片圖中源點附近出現了較大的電阻率值（大于200Ω•m），導致整個三維反演成像效果較差，尤其對2個低阻異常體的成像精度較低。而在加入不等式約束條件后的反演結果（圖7至圖9）中，從圖7中可以看出兩個低阻異常體的成像效果大大提高，在規模、形狀、電阻率值等特征方面與原模型基本一致。

3．2算例二：單個低阻立方體模型模型如圖10所示，在極化率η0＝0，電阻率ρ0＝100Ω•m的均勻半空間中存在1個大小10m×10m×10m（a＝10m），電阻率ρ1＝10Ω•m，極化率η1＝0．1的低阻低極化體，頂面埋深h＝45m，底面埋深55m，立方體兩邊距離x軸與y軸相等，d1＝d2＝5m，紅色圓點代表井口，井口坐標為（50m，50m，0）。4個方位電極點離井口20m；網格剖分單元總數為289280，節點總數為49797。采用二級裝置進行地－井五方位觀測，分別在地面井口方位A0（50m，50m，0），主方位A1（30m，50m，0），A2（70m，50m，0），輔助方位A3（50m，30m，0），A4（50m，70m，0）處固定A極供電進行五方位觀測，井中觀測電極距2m，在井中逐點移動M極，通過三維有限元正演程序共得到255個一次場電位數據和視極化率的“實測數據”。最常用的地井激化法工作方式是地井五方位觀測，由于地井方式單次觀測獲取的觀測數據量非常少，所以三維反演有很強的多解性，很難實現精確的成像反演，較常規的反演方法是通過正演擬合的方式針對某些特定的參數進行最優化反演；而井地方式中，由于測點在地面，可以得到更多的觀測數據，獲取更多的地下信息，能夠更好地進行三維成像反演。地井方式利用單次觀測的數據進行三維電阻率和極化率的反演結果多解性較為嚴重，本文嘗試采取兩種途徑來減小傳統反演方式帶來的多解性問題：1）將地井五方位觀測的數據進行同時反演，增加觀測數據的個數；2）增加更多的先驗信息約束。整個區域分為邊界區域和目標區域，目標區域選為100m×100m×100m，擴展區域分別向3個方向延展200m。反演初始模型選取電阻率ρ＝100Ω•m，極化率η＝0的均勻半空間。反演參數選擇：正則化因子初始值λ0＝0．05，減小系數q＝0．1，給定反演終止收斂系數ε2＝ε3＝1．0×10－5。首先進行地井五方位觀測數據的同時反演，不施加約束條件，利用非線性共軛梯度法對模型進行電阻率反演計算，圖11與圖12為地井五方位觀測數據的同時反演電阻率結果。從圖12中可以看出，基本能反映有低阻異常體的存在，異常體的中心位置基本吻合，無多余構造產生，但是異常體的邊界范圍模糊，延伸較遠；由于異常體埋深較大，在地面供電的情況下，地下電性異常體產生的電位異常幅值很小，因此反演出來的電阻率值與真實電阻率值相差較大。結合五方位觀測數據，將不等式約束引入到反演目標函數中。觀測數據中的視電阻率都在110Ω•m以內，因此將整個模型中網格單元的電阻率上限設為150Ω•m，下限設為0；在使用電阻率探測時，經常會配合其它地球物理探測方法（如淺層地震勘探法、地質雷達探測法），這些方法對異常體界面識別和定位具有較好的效果，由此獲得的異常體的位置與形態的信息是一種極為重要的先驗構造信息。因此將這種構造位置的先驗信息作為局部約束加入到三維反演目標函數中。假設由其它地球物理探測方法可知異常體的構造位置與形態信息，由探區礦石標本物性參數可得：55m≤x≤65m，55m≤y≤65m，－55m≤z≤－45m處已知異常體的電阻率值在50Ω•m以內。因此將這處網格單元的電阻率上限設為50Ω•m，下限設為0。懲罰因子μ取為常數1。圖13至圖15為施加不等式約束后的電阻率與極化率的反演結果，從反演結果可以看出，異常體水平和縱向分辨率較五方位同時反演明顯提高，在位置、規模、形狀等特征方面均與原模型基本一致，驗證了局部不等式約束能夠得到更為精確的三維反演成像結果。

4.結論

1）我們以懲罰函數法的方式將介質電阻率取值范圍的不等式約束施加到三維井中激電反演目標函數中，改善了反演的效果，因此本文提出的以懲罰函數法施加不等式約束條件的方式是可行的，為約束條件的施加方式提供了一條新的途徑；2）在地井觀測模式中，我們嘗試采取兩種途徑來減小傳統反演方法帶來的多解性問題，理論模型反演試算表明，利用構造位置的先驗信息作為局部約束加入到三維反演目標函數中得到的反演結果更為精確，反演的多解性顯著降低；3）懲罰因子的取值是本文的難點，本文懲罰因子的取值采用的是實驗法，過大或者過小的懲罰因子都會影響反演的效果，有時甚至會出現不收斂的情況，因此給出一個合理的懲罰因子非常重要；4）本文的非線性共軛梯度反演算法計算依然很耗時，在反演多個點源多次觀測數據的情況下顯得更為突出，下一步需要研究并行反演算法進行計算，實現反演算法的快速計算。

參考文獻：

［4］蘇朱劉，胡文寶，嚴良?。娮杪屎蜆O化率測深法的正演修正法反演［J］．石油物探，2005，44（2）：194－198

［10］張樹彬，孫東，郭秀娟．航空電磁法三維電阻率反演成像技術［J］．石油物探，2006，45（1）：98－101

［13］吳小平，汪彤彤．電阻率三維反演穩定性和可靠性研究［J］．煤田地質與勘探，2002，30（1）：57－60

［14］吳小平，徐果明．利用共軛梯度法的電阻率三維反演研究［J］．地球物理學報，2000，43（3）：420－427

［15］吳小平，徐果明．電阻率三維反演中偏導數矩陣的求取與分析［J］．石油地球物理勘探，1999，34（4）：363－372

［16］吳小平，汪彤彤．利用共軛梯度方法的電阻率三維反演若干問題研究［J］．地震地質，2001，23（2）：321－327

［17］吳小平．利用共軛梯度方法的激發極化三維快速反演［J］．煤田地質與勘探，2004，32（5）：62－64

［18］黃俊革．三維電阻率／極化率有限元正演模擬與反演成像［D］．長沙：中南大學，2003

［19］黃俊革，阮百堯，鮑光淑．基于有限單元法的三維地電斷面電阻率反演［J］．中南大學學報（自然科學版），2004，35（2）：295－299

［20］強建科．起伏地形三維電阻率正演模擬與反演成像研究［D］．武漢：中國地質大學（武漢），2006

［21］宛新林，席道瑛，高爾根，等．用改進的光滑約束最小二乘正交分解法實現電阻率三維反演［J］．地球物理學報，2005，48（2）：439－444

［22］劉斌，李術才，聶利超，等．基于自適應加權光滑約束與PCG算法的三維電阻率探測反演成像［J］．巖土工程學報，2012，34（9）：1646－1653

［23］劉斌，李術才，李樹忱，等．基于不等式約束的最小二乘法三維電阻率反演及其算法優化［J］．地球物理學報，2012，55（1）：260－268

［24］劉斌，聶利超，李術才，等．三維電阻率空間結構約束反演成像方法［J］．巖石力學與工程學報，2012，31（11）：2258－2268

［25］吳小平，汪彤彤．電阻率三維反演方法研究進展［J］．地球物理學進展，2002，17（1）：156－162

［28］胡祖志，胡祥云，何展翔．大地電磁非線性共軛梯度擬三維反演［J］．地球物理學報，2006，49（4）：1226－1234

［29］翁愛華，劉云鶴，賈定宇，等．地面可控源頻率測深三維非線性共軛梯度反演［J］．地球物理學報，2012，55（10）：3506－3515

［30］林昌洪，譚捍東，舒晴，等．可控源音頻大地電磁三維共軛梯度反演研究［J］．地球物理學報，2012，55（11）：3829－3838

［31］劉云鶴，殷長春．三維頻率域航空電磁反演研究［J］．地球物理學報，2013，56（12）：4278－4287

［32］劉云鶴．三維可控源電磁法非線性共軛梯度反演研究［D］．吉林長春：吉林大學，2011

［33］董浩，魏文博，葉高峰，等．基于有限差分正演的帶地形三維大地電磁反演方法［J］．地球物理學報，2014，57（3）：939－952

［34］董浩．基于有限差分法正演的大地電磁測深帶地形三維反演研究［D］．北京：中國地質大學（北京），2013

［35］胡朝俊．三維直流電法共軛梯度反演算法研究［D］．北京：中國地質大學（北京），2006

［36］屈有恒．井地電阻率法及雙頻激電三維數值模擬與反演研究［D］．北京：中國地質大學（北京）

作者：王智 潘和平 吳愛平 李剛 方思南 單位：長江大學電子信息學院 中國地質大學（武漢）地球物理與空間信息學院