二維曲線計(jì)算機(jī)回歸方法在墑情傳感器律定中的應(yīng)用

尚德功

（河南省人民勝利渠管理局，河南 新鄉(xiāng) 453003）

1 引言

在實(shí)際工作中，經(jīng)常需要根據(jù)一些試驗(yàn)數(shù)據(jù)求出反映其變化規(guī)律的曲線，這些曲線可以是二維、三維等等。這里僅就二維曲線的計(jì)算機(jī)回歸進(jìn)行探討。這里的二維是指一個(gè)量隨另一個(gè)量的變化而變化的規(guī)律；這里所說的曲線是廣義上的曲線，包括直線和折線等。

如圖1所示小麥生長期與耗水強(qiáng)度的關(guān)系，這條曲線可根據(jù)觀測試驗(yàn)數(shù)據(jù)在坐標(biāo)紙上描點(diǎn)并根據(jù)其變化趨勢畫出平滑的曲線，還可以用插值法、及其它方法來粗略地求出其函數(shù)關(guān)系式。但如何利用計(jì)算機(jī)來精確地模擬出曲線，進(jìn)而準(zhǔn)確地求出其函數(shù)關(guān)系式呢？

圖1

對于某些變化規(guī)律呈線性或近似線性的函數(shù)關(guān)系，可用最小二乘法來作線性回歸，得出的是一次直線函數(shù)關(guān)系。對于某些參量之間呈現(xiàn)的不是線性關(guān)系，而接近初等幾何所包括的幾種曲線，如：拋物線、雙曲線、對數(shù)曲線，指數(shù)曲線、橢圓，圓弧等時(shí)，也可根據(jù)點(diǎn)狀圖呈現(xiàn)的規(guī)律，將曲線假設(shè)為相應(yīng)初等幾何曲線的標(biāo)準(zhǔn)方程，來求待定系數(shù)。

但是，在科研工作中，常見兩維參量之間的對應(yīng)關(guān)系不能用線性和上述初等幾何曲線來描述。對于這些函數(shù)關(guān)系又如何用函數(shù)表達(dá)式來描述呢？能否找出一種統(tǒng)一的方法，無論是線性還是某種初等幾何曲線，或是非初等幾何曲線，只要是二維參量函數(shù)關(guān)系都能根據(jù)其試驗(yàn)數(shù)據(jù)求出其對應(yīng)的函數(shù)表達(dá)式呢？進(jìn)而又能在計(jì)算機(jī)上得到統(tǒng)一解決呢？

2 二維參量曲線回歸的原理

數(shù)學(xué)中有這樣一個(gè)定理：若Y=f(X)在點(diǎn)X0的某一鄰域內(nèi)連續(xù)且有直到(n+1)階的連續(xù)導(dǎo)數(shù)，對X0某個(gè)鄰域內(nèi)的任一點(diǎn)X都有：

f(X)=a0+a1(X－X0)+a2(X－X0)2+…+an(X－X0)n+… X∈(X0－δ，X0+δ)

當(dāng)X0=0時(shí)，上式化為：

對任一 X∈[0，δ)，顯然 X∈(－δ，+δ)。 因此，該定理在[0，δ)上成立。

對于每一組試驗(yàn)數(shù)據(jù)（X，Y)，上式都是關(guān)于系數(shù) a0，a1，a2，…，an，…的方程。

即：Y=a0+a1X+a2X2+…+anXn+…

這里，a0，a1，a2，…，an，…皆為一次，問題又轉(zhuǎn)化為求線性方程待定系數(shù)問題，即可用最小二乘法來求解，并可用F—檢驗(yàn)來觀察其擬合情況，察看其顯著性水平。（有關(guān)用最小二乘法求解方程待定系數(shù)的問題本文不再贅述）

3 二維曲線回歸的方法及應(yīng)用軟件

將上述二維曲線回歸的方法編程在計(jì)算機(jī)上運(yùn)行，可以使問題輕松得到解決，可以用C語言來編寫，先將二維曲線在X=0點(diǎn)附近的展開式假設(shè)為：

f(X)=a0+a1X+a2X2+…+anXn

試驗(yàn)數(shù)據(jù)代入后，上式成為關(guān)于 a0，a1，a2，…，an的待定系數(shù)方程，用矩陣最小二乘法求解，得出回歸曲線，并顯示其相關(guān)系數(shù)R，（這里0≤R≤l，R值愈接近于0說明回歸效果越好；愈接近于1回歸效果越差）。可根據(jù)R值與0和1的接近程度來評(píng)價(jià)曲線的回歸教果，決定是否增加展開式項(xiàng)數(shù)n的值，作下一次回歸試驗(yàn)。若回歸效果較好，即可進(jìn)一步實(shí)施F一檢驗(yàn)。在檢驗(yàn)過程中，計(jì)算機(jī)先后要求輸入在0.1，0.05，0.01置信度下的F一檢驗(yàn)值。并與E=(N－K)U/KQ的值相比較，判斷在相應(yīng)顯著性水平下，曲線的回歸效果是否顯著。

如果當(dāng)函數(shù)展開式項(xiàng)數(shù)n的值增加很大，以至于計(jì)算機(jī)內(nèi)存發(fā)生溢出，仍不能求出擬合效果良好的曲線時(shí)，可根據(jù)情況確定在非零點(diǎn)求出該曲線的展開式。這時(shí)只要根據(jù)程序運(yùn)行要求，結(jié)合實(shí)際情況，確定輸入X0的值，即可同上運(yùn)行，直到求出理想的曲線多項(xiàng)式為止。

4 在墑情傳感器曲線律定工作中的應(yīng)用

在科技攻關(guān)項(xiàng)目：翟坡試區(qū)作物水環(huán)境自動(dòng)監(jiān)測及預(yù)測預(yù)報(bào)研究課題中，購置的土壤墑情傳感器的頻率與土壤墑情變化曲線，是生產(chǎn)廠家在實(shí)驗(yàn)室里作的，與安裝現(xiàn)場實(shí)際情況差距很大。為保證所測墑情的準(zhǔn)確性，我們對購置的墑情傳感器逐個(gè)進(jìn)行了曲線律定。這些曲線就是用該軟件進(jìn)行回歸的。

在各傳感器的埋設(shè)現(xiàn)場取土樣烘干得到的土壤含水率數(shù)據(jù)，與相應(yīng)的傳感器頻率觀測數(shù)據(jù)資料如表1。

在應(yīng)用該軟件求其中某一個(gè)墑情傳感器的變化曲線時(shí)，首先將相應(yīng)的實(shí)測數(shù)據(jù)資料按順序輸入程序的數(shù)值語句中，根據(jù)實(shí)際情況選定適當(dāng)?shù)腦0點(diǎn)（一般情況下先選定X0=0），在X0點(diǎn)附近將函數(shù)展開成多頂式。然后，確定展開式的頂數(shù)，一般從一項(xiàng)開始試算，每次增加一項(xiàng)。試算期間可以只觀測相關(guān)系數(shù)R的值與0和1的接近程度，不作F－檢驗(yàn)。通過試算可以找出R值最接近于0的一次試算，確定其展開式的頂數(shù)n的值。再根據(jù)n的值進(jìn)一步運(yùn)行軟件求出曲線，并作F一檢驗(yàn)，若能通過F－檢驗(yàn)，則該曲線就是所需求的土壤墑情與傳感器頻率變化函數(shù)。若不能通過F－檢驗(yàn)，則需要更換展開點(diǎn)X0的值繼續(xù)試算。

用此方法和該軟件回歸曲線時(shí)，還應(yīng)注意的一個(gè)問題是：數(shù)據(jù)資料的選取應(yīng)稍超出實(shí)際應(yīng)用曲線定義域的范圍，防止回歸出來的曲線在定義域邊界附近出現(xiàn)變異現(xiàn)象，以保證所應(yīng)用曲線的正確性和完整性。

根據(jù)上述實(shí)際觀測資料描出數(shù)據(jù)點(diǎn)狀圖，依據(jù)點(diǎn)狀圖發(fā)展趨勢，補(bǔ)充邊界數(shù)據(jù)如表2。

表2 補(bǔ)充邊界數(shù)據(jù)表

將上述資料輸入該軟件的數(shù)據(jù)語句，運(yùn)行后分別得出回歸結(jié)果如下：

09#（33# 分站20cm埋深）

Y=－0.11+12.22(X-0)－1.48(X-0)2+7.18E-02(X-0)3－1.56E-03(X-0)4+1.26E-05(X-0)5

R=Qe/Syy=0.4792

Et=(n－k)U/KQe)=6.0863

當(dāng) α=0.1 時(shí),Fa(5,28)=2.06,Et＞Fa(5,28),OK!

當(dāng) α=0.05 時(shí),Fa(5,28)=2.56,Et＞Fa(5,28),OK!

當(dāng) α=0.01 時(shí),Fa(5,28)=3.75,Et＞Fa(5,28),OK!

11#（33#分站50cm埋深）

Y=－9.3948E-02+13.0948(X-0)－1.8511(X-0)2+0.11(X-0)3－2.9682E-03(X-0)4+2.9681E-05(X-0)5

R=Qe/Syy=0.4431

Et=(n－k)U/KQe)=5.7825

當(dāng) α=0.1 時(shí),Fa(5,23)=2.13,Et＞Fa(5,23),OK!

當(dāng) α=0.05 時(shí),Fa(5,23)=2.64,Et＞Fa(5,23),OK!

當(dāng) α=0.01 時(shí),Fa(5,23)=3.94,Et＞Fa(5,23),OK!

08#(33#分站80cm埋深)

Y=4.4937E-02+12.1704(X-0)-1.5222(X-0)2+7.9693E-02(X-0)3-1.873E-03(X-0)4+1.6311E-05(X-0)5

R=Qe/Syy=0.4886

Et=(n－k)U/KQe)=4.8149

當(dāng) α=0.1 時(shí),Fa(5,23)=2.13,Et＞Fa(5,23),OK!

當(dāng) α=0.05 時(shí),Fa(5,23)=2.64,Et＞Fa(5,23),OK!

當(dāng) α=0.01 時(shí),Fa(5,23)=3.94,Et＞Fa(5,23),OK!

由上述各自的F－撿驗(yàn)可知：上述三條曲線均在α=0.01顯著性水平下回歸效果顯著，墑情傳感器能夠反映土壤墑情的變化規(guī)律，是要求的函數(shù)表達(dá)式。

5 對二維參量曲線回歸方法及軟件的評(píng)價(jià)

為驗(yàn)證所求出的土壤墑情與傳感器頻率變化曲線，我們繼續(xù)進(jìn)行取土樣烘干試驗(yàn)，對烘干得出的土壤含水率值與相應(yīng)的變化函數(shù)求出的土壤含水率值進(jìn)行了比較。如表3。

33#分站20cm埋深09#傳感器頻率與土壤含水率變化函數(shù)

Y=－0.11+12.22(X-0)－1.48(X-0)2+7.18E-02(X-0)3－1.56E-03(X-0)4+1.26E-05(X-0)5

表3 計(jì)算值與墑情烘干值比較表

從表3可以看出：雖然回歸曲線函數(shù)計(jì)算值與土樣烘干值絕對誤差小于0.5的僅占45.45%，但是差值小于l的卻占81.82%，對±壤墑情來說，這是可以接受的。因此，該曲線函數(shù)能反映土壤墑情傳感器頻率值與土壤含水率之間的變化規(guī)律。

33#分站50cm理深11#傳感器頻率與土壤含水率變化函數(shù)

Y=－9.3948E-02+13.0948(X-0)－1.8511(X-0)2+0.11(X-0)3－2.9682E-03(X-0)4+2.9681E-05(X-0)5

表4 計(jì)算值與墑情烘干值比較表

從表4可以看出：回歸曲線函數(shù)計(jì)算值與土樣烘干值絕對誤差小于0.5的占63.64%，絕對差值小于l的占81.82%，說明用該軟件回歸出來的曲線函數(shù)能較好地反映土壤墑情傳感器頻率值與土壤含水率之間的變化規(guī)律。

33#分站80cm埋深08#傳感器頻率與土壤含水率變化函數(shù)

Y=4.4937E-02+12.1704(X-0)-1.5222(X-0)2+7.9693E-02(X-0)3-1.873E-03(X-0)4+1.6311E-05(X-0)5

表5 計(jì)算值與墑情烘干值比較表

從表5可以看出：回歸曲線函數(shù)計(jì)算值與土樣烘干值絕對誤差小于1的占72.73%，說明用該軟件回歸出來的曲線函數(shù)能較好地反映土壤墑情傳感器頻率值與土壤含水率之間的變化規(guī)律。

通過上述驗(yàn)證，說明用二維參量曲線回歸方法和軟件得出的曲線函數(shù)能較好地反映傳感器頻率與土壤墑情的變化規(guī)律。

用該軟件作曲線回歸時(shí)，只需將原始數(shù)據(jù)資料中的個(gè)別明顯非點(diǎn)刪除，即可輸入計(jì)算機(jī)進(jìn)行回歸運(yùn)算，避免了人為干預(yù)的因素，曲線的走向完全根據(jù)數(shù)據(jù)分布規(guī)律而定。不象插值法和其它方法那樣，先要人為地刪除許多所謂不規(guī)律的點(diǎn)，再根據(jù)觀察保留有規(guī)律的曲線附近的點(diǎn)，這樣難免存在觀察者厚此薄彼的現(xiàn)象，以致影響理論曲線函數(shù)與客觀規(guī)律的吻合程度。

[1]馬玉芳;陳建華;郝?lián)P滿.基于C語言對三次樣條函數(shù)的求解及程序.價(jià)值工程.2011.11

[2]蘇俊杰.基于矩陣樣條函數(shù)的二階矩陣微分方程數(shù)值解法研究.合肥工業(yè)大學(xué).2010.03

[3]張曉丹;邵帥;劉欽圣.基于樣條函數(shù)的光滑支持向量機(jī)模型.北京科技大學(xué)學(xué)報(bào).2012.7

[4]邱慧敏;楊濟(jì)安.樣條函數(shù)與樣條小波.重慶郵電學(xué)院學(xué)報(bào)(自然科學(xué)版).2000-03