磁法數(shù)據(jù)多目標(biāo)粒子群反演方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明具體涉及一種用于磁法數(shù)據(jù)的多目標(biāo)粒子群反演方法。
【背景技術(shù)】
[0002] 地球物理反演存在不穩(wěn)定性和多解性，吉洪諾夫正則化方法[1'2]已被廣泛接受并 應(yīng)用于解決此類(lèi)病態(tài)反問(wèn)題。正則化反演的目標(biāo)函數(shù)由數(shù)據(jù)擬合、模型約束和正則化因子 構(gòu)成，引入模型約束和正則化因子可穩(wěn)定反演過(guò)程，降低反演的多解性。其中，正則化因子 是數(shù)據(jù)擬合和模型約束之間的折衷系數(shù)，其取值直接影響到反演結(jié)果的好壞，諸多學(xué)者對(duì) 正則化參數(shù)選取問(wèn)題進(jìn)行了深入的研究。如何選取最合適的正則化因子，目前仍然是正則 化反演研究的熱點(diǎn)與難點(diǎn)之一。
[0003] 從反演目標(biāo)函數(shù)的求解方法來(lái)看，地球物理反演方法分為兩類(lèi)：局部搜索法和全 局搜索法。傳統(tǒng)的局部搜索法以目標(biāo)函數(shù)梯度信息為指導(dǎo)，線(xiàn)性迭代求解。該類(lèi)方法穩(wěn)定 性好、計(jì)算效率高，但過(guò)于依賴(lài)初始模型的選擇，不易加入約束條件。蒙特卡羅、模擬退火、 遺傳算法、人工神經(jīng)網(wǎng)絡(luò)、粒子群優(yōu)化、蟻群優(yōu)化、差分進(jìn)化等經(jīng)典的全局搜索法則不需要 目標(biāo)函數(shù)梯度信息，采用啟發(fā)式隨機(jī)搜索策略直接搜索整個(gè)解空間，具有不依賴(lài)初始模型、 不易陷入局部最優(yōu)解等優(yōu)點(diǎn)。然而，目前該類(lèi)方法更多地考慮數(shù)據(jù)擬合，而與正則化反演的 結(jié)合較少。由于全局搜索法比局部搜索法更耗時(shí)，其正則化因子的選取也更為困難。如何 既解決正則化因子選取難題，又克服初始模型依賴(lài)，是地球物理反演領(lǐng)域亟待解決的問(wèn)題。

【發(fā)明內(nèi)容】

[0004] 有鑒于此，有必要研究一種既能解決正則化因子選取難題，又克服初始模型依賴(lài) 的地球物理反演方法。
[0005] -種磁法數(shù)據(jù)多目標(biāo)粒子群反演算法，其包括如下步驟：
[0006] S1、在磁化率上下界范圍內(nèi)，隨機(jī)初始化各粒子的初始位置，各粒子初始速度為 〇 ；
[0007] S2、將每個(gè)粒子的位置看作模型向量，計(jì)算所述模型向量的理論響應(yīng)，進(jìn)而計(jì)算得 到每個(gè)模型向量的數(shù)據(jù)擬合和模型約束；
[0008] S3、初始化每個(gè)粒子的個(gè)體最優(yōu)，選擇與非受控目標(biāo)函數(shù)向量對(duì)應(yīng)的粒子構(gòu) 成初始多目標(biāo)反演解集；
[0009] S4、從初始多目標(biāo)反演解集中選擇一個(gè)全集最優(yōu)LEADER ;
[0010] S5、根據(jù)全集最優(yōu)LEADER和個(gè)體最優(yōu)更新每個(gè)粒子的速度和位置；
[0011] S6、根據(jù)更新后的速度和位置計(jì)算粒子的模型向量的理論響應(yīng)，進(jìn)而得到新的模 型向量的數(shù)據(jù)擬合和模型約束；
[0012] S7、根據(jù)新的數(shù)據(jù)擬合和模型約束對(duì)每個(gè)粒子的個(gè)體最優(yōu)進(jìn)行更新，
[0013] S8、相應(yīng)更新初始多目標(biāo)反演解集；
[0014] S9、如果所述磁法數(shù)據(jù)多目標(biāo)粒子群反演算法未收斂且未達(dá)到最大迭代次數(shù)，返 回步驟S4,進(jìn)行下一輪迭代；
[0015] S10、直至所述磁法數(shù)據(jù)多目標(biāo)粒子群反演算法收斂且達(dá)到最大迭代次數(shù)，則最終 得到最優(yōu)多目標(biāo)反演解集，同時(shí)輸出最優(yōu)多目標(biāo)反演解集在目標(biāo)函數(shù)空間的映射Pareto 前沿；
[0016] S11、根據(jù)最優(yōu)多目標(biāo)反演解集在目標(biāo)函數(shù)空間的映射Pareto前沿，權(quán)衡數(shù)據(jù)擬 合和模型約束的相對(duì)重要程度，自多目標(biāo)反演解集中選出一個(gè)最優(yōu)解向量，從而得到最優(yōu) 反演結(jié)果。
[0017] 本發(fā)明所述的磁法數(shù)據(jù)多目標(biāo)粒子群反演算法，通過(guò)將正則化反演轉(zhuǎn)換成多目標(biāo) 最優(yōu)化問(wèn)題，先采用全局最優(yōu)化方法同時(shí)求數(shù)據(jù)擬合和模型約束的多目標(biāo)反演解集，再權(quán) 衡兩者的相對(duì)重要程度，最后從反演解集中優(yōu)選出最終反演結(jié)果，從而起到正則化因子的 作用，能同時(shí)解決正則化因子選取困難和初始模型依賴(lài)問(wèn)題。與傳統(tǒng)正則化反演方法只得 到一個(gè)反演結(jié)果不同，該算法能盡可能多地保留可行解，得到一個(gè)反演解集。地球物理研究 人員通過(guò)分析該解集，既可更加深入的理解反演的過(guò)程，又可靈活地從數(shù)據(jù)擬合和模型約 束兩方面進(jìn)行權(quán)衡與選擇，得到優(yōu)于正則化反演的結(jié)果。
【附圖說(shuō)明】
[0018] 圖1為二維磁法正演示意圖；
[0019] 圖2為多目標(biāo)反演示意圖；
[0020] 圖3為本發(fā)明所述的磁法數(shù)據(jù)多目標(biāo)粒子群反演方法的流程圖；
[0021 ] 圖4為圖3中步驟S4的子流程圖；
[0022] 圖5為圖3中步驟S8的子流程圖；
[0023] 圖6為相鄰兩棱柱體模型剖面圖；
[0024] 圖7為相鄰兩個(gè)棱柱體模型無(wú)噪聲數(shù)據(jù)正則化反演結(jié)果；
[0025] 圖8是相鄰兩棱柱體模型無(wú)噪聲多目標(biāo)粒子群反演算法中反演解集的Pareto前 沿；
[0026] 圖9是相鄰兩棱柱體模型無(wú)噪聲多目標(biāo)粒子群反演算法中最優(yōu)粒子的收斂曲線(xiàn)；
[0027] 圖10是相鄰兩棱柱體模型無(wú)噪聲多目標(biāo)粒子群反演算法的注重?cái)?shù)據(jù)擬合的反演 結(jié)果；
[0028] 圖11是相鄰兩棱柱體模型無(wú)噪聲多目標(biāo)粒子群反演算法的注重模型約束的反演 結(jié)果；
[0029] 圖12是相鄰兩棱柱體模型無(wú)噪聲多目標(biāo)粒子群反演算法的最終反演結(jié)果。
[0030] 具體實(shí)施方法
[0031] 為了使本發(fā)明的目的、技術(shù)方案及優(yōu)點(diǎn)更加清楚明白，以下結(jié)合附圖及實(shí)施例，對(duì) 本發(fā)明進(jìn)行進(jìn)一步詳細(xì)說(shuō)明，應(yīng)當(dāng)理解，此處所述的具體實(shí)施案例僅僅用于解釋本發(fā)明，并 不用于限定本發(fā)明。
[0032] 首先簡(jiǎn)要回顧經(jīng)典的正則化反演方法，然后闡述本發(fā)明提出的磁法數(shù)據(jù)多目標(biāo)粒 子群反演算法。
[0033] (1)現(xiàn)有技術(shù)的磁法數(shù)據(jù)正則化反演方法
[0034] 磁法數(shù)據(jù)二維物性正則化反演涉及到規(guī)則網(wǎng)格單元的正演計(jì)算，因此需對(duì)正演計(jì) 算進(jìn)行討論。將地下半空間剖分為規(guī)則排列的二維棱柱體單元，任意一個(gè)棱柱體單元如圖 1所示。當(dāng)磁化方向與地磁場(chǎng)方向一致時(shí)，磁異常計(jì)算公式如下：
[0035] F = 2kFe[sin (21) sin β In (r2r3/r4rl) + (cos2I + sin2 β - sin2I) (Φ1-Φ2-Φ3+Φ4)]. (I)
[0036] 其中，F(xiàn)是磁異常，F(xiàn)f3是地磁場(chǎng)強(qiáng)度，k是磁化率，I是地磁傾角，β是棱柱體走向 與磁北的夾角，^，Φi分別是圖中所示的距離和夾角。
[0037] 磁異常向量F可視為離散網(wǎng)格單元磁化率向量m與常量核矩陣A的乘積，因此， (1)式的矩陣形式如下：
[0038] F = Am. (2)
[0039] 其中，m = Qc1,…h(huán)，…kn)，Ici是第i個(gè)棱柱體的磁化率。
[0040] 反演是找到一個(gè)模型m，使其理論磁異常與實(shí)測(cè)磁異常擬合得最好，因此，數(shù)據(jù)擬 合目標(biāo)函數(shù)定義如下 [2]:
[0041] f! (m) = I I d-Am | 1min, (3)
[0042] 其中，d是實(shí)測(cè)磁異常，Am是模型m的理論磁異常。
[0043] 為穩(wěn)定反演過(guò)程，克服多解性，正則化反演引入模型約束f2(m)和正則化因子λ， 構(gòu)成正則化反演目標(biāo)函數(shù)如下[2]:
[0044] Ρλ (m) = f1 (m) + λ f2 (m) min. (4)
[0045] 最常見(jiàn)的模型約束是模型范數(shù)最?。?br>[0046] f2(m) = I |m| |2, (5)
[0047] 此時(shí)，正則化反演的矩陣形式目標(biāo)函數(shù)為：
[0048] Pλ (m) = Π d-Am | |2+ λ | | m | |2 = (d-Am) T (d-Am) + λ mTm. (6)
[0049] 其正則化解為[1°]:
[0050] m = [ΑΤΑ+λ I] 1AtCL (7)
[0051] 公式（4) (6)表明，正則化反演的實(shí)質(zhì)是將數(shù)據(jù)擬合和模型約束兩個(gè)泛函最優(yōu)化 問(wèn)題轉(zhuǎn)換為單目標(biāo)最優(yōu)化問(wèn)題，其中，正則化因子決定數(shù)據(jù)擬合和模型約束的相對(duì)權(quán)重，其 取值直接影響到反演結(jié)果的好壞。
[0052] (2)本發(fā)明所述磁法數(shù)據(jù)多目標(biāo)粒子群反演算法
[0053] 本發(fā)明所述磁法數(shù)據(jù)多目標(biāo)粒子群反演算法包括如下步驟：
[0054] S1、在磁化率上下界范圍內(nèi)，隨機(jī)初始化各粒子的初始位置，各粒子初始速度為 〇 ；
[0055] S2、將每個(gè)粒子的位置看作模型向量，根據(jù)磁異常計(jì)算公式計(jì)算所述模型向量的 理論響應(yīng)，所述磁異常計(jì)算公式如下：
[0056] F = 2kFe[sin (21) sin β In (r2r3/r4rl) + (cos2I + sin2 β - sin2I) (Φ1-Φ2-Φ3+Φ4)]. (I)
[0057] 其中，F(xiàn)是磁異常向量，F(xiàn)f3是地磁場(chǎng)強(qiáng)度，k是磁化率，I是地磁傾角，β是棱柱體 走向與磁北的夾角， ri，O1分別是圖中所示的距離和夾角；
[0058] 并根據(jù)矩陣形式公式計(jì)算得到每個(gè)模型向量的數(shù)據(jù)擬合和模型約束，所述矩陣形 式公式為：
[0059] F = Am. (2)
[0060] 其中，A為常量核矩陣；m為離散網(wǎng)格單元磁化率向量；m= d · · · h，· · · kn)，1^是 第i個(gè)棱柱體的磁化率；
[0061] S3、初始化每個(gè)粒子的個(gè)體最優(yōu)Phe<，選擇與非受控目標(biāo)函數(shù)向量對(duì)應(yīng)的粒子構(gòu) 成初始多目標(biāo)反演解集；
[0062] S4、從初始多目標(biāo)反演解集中選擇一個(gè)全集最優(yōu)LEADER ;
[0063] S5、根據(jù)全集最優(yōu)LEADER和個(gè)體最優(yōu)pliest吏新每個(gè)粒子的速度和位置；
[0064] S6、根據(jù)更新后的速度和位置，按照公式公式（1) (2)計(jì)算粒子的模型向量的理論 響應(yīng)，進(jìn)而按照數(shù)據(jù)擬合目標(biāo)函數(shù)定義公式及最常見(jiàn)的模型約束是模型范數(shù)最小公式計(jì)算 得到新的模型向量的數(shù)據(jù)擬合和模型約束；
[0065] 所述數(shù)據(jù)擬合目標(biāo)函數(shù)定義公式如下[2]:
[0066] f! (m) = I I d-Am | 1min, (3)
[0067] 其中，d是實(shí)測(cè)磁異常，Am是模型m的理論磁異常。
[0068] 所述最常見(jiàn)的模型約束是模型范數(shù)最小公式：
[0069] f2(m) = I |m| I2, (5)；
[0070] S7、根據(jù)新的數(shù)據(jù)擬合和模型約束對(duì)每個(gè)粒子的個(gè)體最優(yōu)進(jìn)行更新，
[0071] S8、相應(yīng)更新初始多目標(biāo)反演解集；
[0072] S9、如果所述磁法數(shù)據(jù)多目標(biāo)粒子群反演算法未收斂且未達(dá)到最大迭代次數(shù)，返 回步驟S4,進(jìn)行下一輪迭代；
[0073] S10、直至所述磁法數(shù)據(jù)多目標(biāo)粒子群反演算法收斂且達(dá)到最大迭代次數(shù)，則最終 得到最優(yōu)多目標(biāo)反演解集，同時(shí)輸出最優(yōu)多目標(biāo)反演解集在目標(biāo)函數(shù)空間的映射Pareto 前沿；
[0074] S11、根據(jù)最優(yōu)多目標(biāo)反演解集在目標(biāo)函數(shù)空間的映射Pareto前沿，權(quán)衡數(shù)據(jù)擬 合和模型約束的相對(duì)重要程度，自多目標(biāo)反演解集中選出一個(gè)最優(yōu)解向量，從而得到最優(yōu) 反演結(jié)果。
[0075] 其中，所述步驟S4包括以下子步驟：
[0076] S41、將數(shù)據(jù)誤差函數(shù)和模型約束函數(shù)張成的空間劃分為IOX 10均勻網(wǎng)格；
[0077] S42、將REP中個(gè)粒子按其目標(biāo)函數(shù)值對(duì)應(yīng)到不同網(wǎng)格；
[0078] S43、給每個(gè)非空的網(wǎng)格賦予一個(gè)權(quán)值，該權(quán)值等于該網(wǎng)格包含的粒子數(shù)目的倒 數(shù)；
[0079] S44、以輪盤(pán)賭的方式從非空的網(wǎng)格中選擇一個(gè)包含最少粒子的網(wǎng)格；
[0080] S45、從該網(wǎng)格中隨機(jī)選擇一個(gè)粒子作為L(zhǎng)EADER。
[0081] 其中，所述步驟S8包括以下子步驟：
[0082] S81、檢查所有粒子之間的控制和非控制關(guān)系