最小二乘法在地質(zhì)測(cè)量中的應(yīng)用

吳小偉

（西安石油大學(xué)地球科學(xué)與工程學(xué)院，西安710065）

1 引 言

礦區(qū)的地質(zhì)測(cè)量工作為礦區(qū)的前期規(guī)劃和實(shí)際開(kāi)采提供有力的數(shù)據(jù)支撐。目前常見(jiàn)的GPS 測(cè)量技術(shù)、數(shù)字?jǐn)z影測(cè)量技術(shù)等被廣泛應(yīng)用于礦區(qū)的地質(zhì)測(cè)量工作中。此類(lèi)技術(shù)可結(jié)合具體地形地貌來(lái)計(jì)算開(kāi)采空間、開(kāi)采時(shí)間、開(kāi)采設(shè)備等，也可通過(guò)數(shù)字影像測(cè)出坍塌區(qū)的具體形態(tài)和特征，如體積、深度、形態(tài)等[1-3]。由于地形地貌差異和數(shù)據(jù)儀器測(cè)量精度的影響，會(huì)存在一定的誤差[4]，以較大誤差值預(yù)測(cè)實(shí)際的礦產(chǎn)量、開(kāi)采率等數(shù)值，將對(duì)生產(chǎn)實(shí)踐造成不利影響。最小二乘法作為最佳的觀測(cè)組合工具之一[5]，具有精確性高、誤差小的特點(diǎn)，廣泛應(yīng)用于諸多領(lǐng)域中。在現(xiàn)階段大數(shù)據(jù)迅猛發(fā)展的環(huán)境下，以相同地形、地理?xiàng)l件為基準(zhǔn)，可對(duì)開(kāi)采礦區(qū)的數(shù)據(jù)進(jìn)行精準(zhǔn)記錄和分析。采用最小二乘法算法，把開(kāi)采數(shù)據(jù)隨深度的改變，以體積為參考對(duì)實(shí)際開(kāi)采量數(shù)據(jù)進(jìn)行建模。根據(jù)最小二乘法令平方和最小化的方式減小誤差，通過(guò)相關(guān)系數(shù)，利用MATLAB 軟件編程擬合成具體的函數(shù)。后期在具體的評(píng)估和預(yù)測(cè)中，可以從相同的地形和地理位置下從數(shù)據(jù)庫(kù)中匹配曲線函數(shù)，計(jì)算新礦區(qū)的開(kāi)采量和固定體積開(kāi)采率等數(shù)值。

2 算法應(yīng)用原理

最小二乘法是高斯于1809 年所著的《天體運(yùn)動(dòng)論》中提出的。它是數(shù)學(xué)中用于數(shù)據(jù)匹配的一種最佳方法[6]，運(yùn)用逼近原則，尋求數(shù)據(jù)，讓實(shí)測(cè)數(shù)據(jù)與估計(jì)數(shù)據(jù)平方和最小，即距離和最小，來(lái)尋測(cè)和匹配最佳的函數(shù)數(shù)據(jù)，從而擬合成函數(shù)曲線。

2.1 算法應(yīng)用的必要性

礦區(qū)第一深度空間定義為由地表到深度500m處；第二深度空間為500m～2000m；第三深度空間為2000m 以下[7]。以具體深度空間為界限，對(duì)每一深度空間的礦區(qū)開(kāi)采數(shù)據(jù)進(jìn)行收集、整理、分析。由于地形地貌的不同，隨著深度的增加，體積的變化，開(kāi)采總量整體呈上升趨勢(shì)。同一定義空間下，由于地表平面的不規(guī)則形狀，選取固定深度會(huì)造成體積的差異，從而造成固定體積中開(kāi)采率的差異，所以開(kāi)采總量的改變不能用具體的函數(shù)表示。統(tǒng)計(jì)體積V 和開(kāi)采總量F，用散點(diǎn)圖標(biāo)記，利用最小二乘法算法，對(duì)散點(diǎn)圖采取擬合方式，建模構(gòu)造函數(shù)曲線。同時(shí)，針對(duì)地形地貌的改變，平面圖形的建立也可采用最小二乘法算法減小誤差，進(jìn)行建模，模擬函數(shù)曲線。

2.2 最小二乘法原理

在對(duì)地表平面面積的運(yùn)算中，長(zhǎng)度和寬度的改變，以及礦產(chǎn)資源中固定體積V 與礦產(chǎn)量F 的改變，皆符合建模構(gòu)造函數(shù)曲線。結(jié)合擬合次數(shù)的方法[8-9]，對(duì)散點(diǎn)圖進(jìn)行數(shù)值分析，并通過(guò)MATLAB 軟件數(shù)據(jù)建模，尋求最佳的函數(shù)圖像。理論推導(dǎo)過(guò)程為：對(duì)給定數(shù)(xi,yi) （i=0,1,2,...,m），共 m+1 個(gè)數(shù)據(jù)坐標(biāo)點(diǎn)，選取多項(xiàng)式P(x)，滿足：

函數(shù)P(x)即稱(chēng)為擬合函數(shù)，或最小二乘法解，設(shè)：

式中，a0,a1,a2,...,an為待求的未知數(shù)，n 為多項(xiàng)式的最高次冪，因此該問(wèn)題化為求I=I(a0,a1,a2,...,an)的極值問(wèn)題。由多元函數(shù)求極值的必要條件：

可以得到：

這是一個(gè)關(guān)于a0,a1,a2,...,an的線性方程組，用矩陣表示如下：

因此，只要給出坐標(biāo)點(diǎn)(xi,yi)及其個(gè)數(shù)，再給出所要擬合的參數(shù)n，則可求出未知數(shù)矩陣[a0,a1,a2,...,an]。

3 地質(zhì)測(cè)量實(shí)施應(yīng)用

3.1 地質(zhì)測(cè)量設(shè)計(jì)思路

地質(zhì)測(cè)量的整體實(shí)施步驟及部分環(huán)節(jié)的設(shè)計(jì)理念如圖1 所示。根據(jù)地表形狀測(cè)繪出圖形，對(duì)圖形求解面積。由于平面面積的不規(guī)則性，按積分求解平面面積。為提高精確性，選取固定深度求解體積，統(tǒng)計(jì)礦產(chǎn)開(kāi)采量，利用所記錄體積內(nèi)的開(kāi)采量，運(yùn)用二乘法對(duì)數(shù)據(jù)建模，運(yùn)用MATLAB 軟件編程，擬合成具體的函數(shù)曲線。

圖1 地質(zhì)測(cè)量流程圖

3.2 平面面積計(jì)算理論

通過(guò)地表平面形狀，按比例縮小，測(cè)繪出具體的平面圖紙。利用平面直角坐標(biāo)系計(jì)算圖紙面積，按比例放大，計(jì)算地表平面實(shí)際面積。平面面積計(jì)算可采取過(guò)蒙特卡洛方法，然而由于地表平面的不規(guī)則形狀，類(lèi)似幾何概型面積比值求解的思路不適用于地質(zhì)礦產(chǎn)精準(zhǔn)的測(cè)量。在此種情況下，對(duì)平面面積的求解采用以X 軸為寬，Y 軸為長(zhǎng)建立平面直角坐標(biāo)系，為減小誤差寬度改變、遵循等差數(shù)列及公差較小的增長(zhǎng)方式，以平面圖形長(zhǎng)度的邊緣建立平面直角坐標(biāo)系（x,y）。利用MATLAB 軟件對(duì)數(shù)據(jù)建模處理，擬合成y 與x（即平面長(zhǎng)度與寬度）的具體函數(shù)表達(dá)式，利用積分思想求解面積。

不規(guī)則形狀的平面圖形積分法如圖2 所示。在區(qū)間[a,b]內(nèi)分割成 n 份（n →∞），每份矩形面積的加和為最后整個(gè)區(qū)間的面積，即：

由于f(x)數(shù)據(jù)建模是根據(jù)冪指數(shù)擬合具體的函數(shù)圖像，所以f(x)在區(qū)間內(nèi)是可積分的。把該積分的極限值定為該區(qū)間的定積分，即：

以此求出測(cè)繪圖形的具體面積，根據(jù)縮放比例求出地表平面的實(shí)際面積。

圖2 不規(guī)則形狀的平面圖形積分法

3.3 模擬面積仿真實(shí)驗(yàn)

對(duì)數(shù)據(jù)的面積求解，模擬出擬合函數(shù)曲線。遵循上述平面面積求解理論，根據(jù)長(zhǎng)度和寬度確定7組（x,y）坐標(biāo)，具體為：（x1,y1）＝（-3,4）；（x2,y2）＝（-2,6）；（x3,y3）＝（-1,7）；（x4,y4）＝（0,9）；（x5,y5）＝（1,12）；（x6,y6）＝（2,34）；（x7,y7）＝（3,31）。運(yùn)用數(shù)據(jù)建模的思想，結(jié)合MATLAB 軟件進(jìn)行編程設(shè)計(jì)，對(duì)數(shù)據(jù)做多項(xiàng)式最小二乘擬合。利用已知的坐標(biāo)呈現(xiàn)散點(diǎn)圖的形式，確定模擬參數(shù)的冪指數(shù)，根據(jù)數(shù)據(jù)確定矩陣，研究?jī)纭⑿?、列之間的關(guān)系式。根據(jù)方程確定系數(shù)，制作出具體的擬合函數(shù)曲線。實(shí)驗(yàn)結(jié)果如圖3、圖4 所示，以多次擬合結(jié)果的兩種情況為例做出對(duì)比分析，圖3 是六次冪的擬合曲線，圖4 是五次冪的擬合曲線。基于最小二乘法圖像擬合原理，使散點(diǎn)圖中坐標(biāo)點(diǎn)盡量多地落在回歸方程上，可果斷舍去五次冪的擬合圖像，選取六次冪擬合圖像。

圖3 六次冪擬合曲線

圖4 五次冪擬合曲線

用MATLAB 求出具體的系數(shù) a(1), a(2)，...，a(7)數(shù)值，具體計(jì)算語(yǔ)句為：f(x)=(9.0000)+(-2.2737e-013)*x+(-0.8056)*x.^2+(-2.7500)*x.^3+(1.4444)*x.^4+(0.2500)*x.^5+(0.1389 )*x.^6。根據(jù)函數(shù)解析式，利用積分求解平面面積。以下代碼為MATLAB 軟件編寫(xiě)圖像擬合過(guò)程的程序設(shè)計(jì)代碼：

4 礦產(chǎn)開(kāi)采量測(cè)量

求得地表平面的面積，以深度的增長(zhǎng)遵循等差數(shù)列公差（d＜1m）進(jìn)行增長(zhǎng)，通過(guò) V=S·h 求得體積，并在開(kāi)采過(guò)的礦產(chǎn)中統(tǒng)計(jì)該體積中的礦產(chǎn)量F，以體積作為X 坐標(biāo)，實(shí)際的礦產(chǎn)量做為縱坐標(biāo)建立多組數(shù)據(jù)（V,F）。運(yùn)用最小二乘法算法，同上述地表面積求解過(guò)程的數(shù)據(jù)擬合方法類(lèi)似，進(jìn)行數(shù)據(jù)建模。采用冪指數(shù)形式擬合成具體曲線函數(shù)。根據(jù)不同冪指數(shù)的形式做多次擬合，選取散點(diǎn)圖盡量多的在曲線上選擇最高冪指數(shù)的擬合曲線。也可在MATLAB中進(jìn)行編程設(shè)置，輸入數(shù)據(jù)后一次輸出多個(gè)最高冪指數(shù)擬合圖像。觀察坐標(biāo)點(diǎn)在擬合曲線的多少，選擇擬合曲線函數(shù)，具體求解擬合參數(shù)的系數(shù)。為提高精確度、減小誤差，可在Excel 文件中建立多組數(shù)據(jù)，在MATLAB 中通過(guò)讀入的方式錄入數(shù)據(jù)。確定擬合曲線后，可根據(jù)曲線函數(shù)計(jì)算固定深度、固定面積礦產(chǎn)開(kāi)采量、開(kāi)采率、固定深度總開(kāi)采量等數(shù)據(jù)。

經(jīng)以上步驟，在大數(shù)據(jù)的支持下，可對(duì)不同地理位置和環(huán)境下的礦產(chǎn)數(shù)據(jù)進(jìn)行數(shù)據(jù)建模，擬合成不同的曲線函數(shù)。相同條件下，遇到未開(kāi)采的礦區(qū)時(shí)，可以從已建成的擬合曲線函數(shù)數(shù)據(jù)庫(kù)中，匹配函數(shù)曲線與解析式，通過(guò)深度和地形面積數(shù)據(jù)預(yù)測(cè)未開(kāi)采礦區(qū)的開(kāi)采量及固定體積下的開(kāi)采數(shù)據(jù)等。

開(kāi)采過(guò)程中，通過(guò)匹配擬合函數(shù)，可對(duì)某一體積下的開(kāi)采率進(jìn)行比對(duì)，研究開(kāi)采狀態(tài)和未知狀況，為合理、科學(xué)開(kāi)采提供參考及技術(shù)支持。

5 結(jié)束語(yǔ)

最小二乘法作為匹配數(shù)據(jù)的最優(yōu)算法之一，廣泛的應(yīng)用到各個(gè)領(lǐng)域，本研究將其應(yīng)用到地質(zhì)測(cè)量中。在大數(shù)據(jù)下，為減小誤差、遵循實(shí)用性的準(zhǔn)則，針對(duì)開(kāi)采過(guò)的礦區(qū)，運(yùn)用最小二乘法法，以MATLAB軟件編程，對(duì)地表平面的面積、開(kāi)采量等數(shù)據(jù)建模處理，以冪指數(shù)擬合方式確定具體擬合圖像，構(gòu)造具體曲線函數(shù)解析式，建立多種地理環(huán)境下擬合函數(shù)數(shù)據(jù)庫(kù)。地表面積求解中擬合函數(shù)圖像與實(shí)際圖形結(jié)合，運(yùn)用積分巧妙解決了不規(guī)則圖形面積求解問(wèn)題。數(shù)據(jù)建模擬合的具體函數(shù)對(duì)今后新礦區(qū)的開(kāi)采提供有利的、精確的數(shù)據(jù)支撐。數(shù)據(jù)模擬實(shí)驗(yàn)對(duì)數(shù)據(jù)的建模處理提供了思路和方法，對(duì)理論數(shù)據(jù)和實(shí)際數(shù)據(jù)逼近提供了有利的保證，同時(shí)也驗(yàn)證此算法在地址礦產(chǎn)測(cè)量中的實(shí)用性和可行性。