尚德功
(河南省人民勝利渠管理局,河南 新鄉 453003)
在實際工作中,經常需要根據一些試驗數據求出反映其變化規律的曲線,這些曲線可以是二維、三維等等。這里僅就二維曲線的計算機回歸進行探討。這里的二維是指一個量隨另一個量的變化而變化的規律;這里所說的曲線是廣義上的曲線,包括直線和折線等。
如圖1所示小麥生長期與耗水強度的關系,這條曲線可根據觀測試驗數據在坐標紙上描點并根據其變化趨勢畫出平滑的曲線,還可以用插值法、及其它方法來粗略地求出其函數關系式。但如何利用計算機來精確地模擬出曲線,進而準確地求出其函數關系式呢?

圖1
對于某些變化規律呈線性或近似線性的函數關系,可用最小二乘法來作線性回歸,得出的是一次直線函數關系。對于某些參量之間呈現的不是線性關系,而接近初等幾何所包括的幾種曲線,如:拋物線、雙曲線、對數曲線,指數曲線、橢圓,圓弧等時,也可根據點狀圖呈現的規律,將曲線假設為相應初等幾何曲線的標準方程,來求待定系數。
但是,在科研工作中,常見兩維參量之間的對應關系不能用線性和上述初等幾何曲線來描述。對于這些函數關系又如何用函數表達式來描述呢?能否找出一種統一的方法,無論是線性還是某種初等幾何曲線,或是非初等幾何曲線,只要是二維參量函數關系都能根據其試驗數據求出其對應的函數表達式呢?進而又能在計算機上得到統一解決呢?
數學中有這樣一個定理:若Y=f(X)在點X0的某一鄰域內連續且有直到(n+1)階的連續導數,對X0某個鄰域內的任一點X都有:
f(X)=a0+a1(X-X0)+a2(X-X0)2+…+an(X-X0)n+… X∈(X0-δ,X0+δ)
當X0=0時,上式化為:

對任一 X∈[0,δ),顯然 X∈(-δ,+δ)。 因此,該定理在[0,δ)上成立。
對于每一組試驗數據(X,Y),上式都是關于系數 a0,a1,a2,…,an,…的方程。
即:Y=a0+a1X+a2X2+…+anXn+…
這里,a0,a1,a2,…,an,…皆為一次,問題又轉化為求線性方程待定系數問題,即可用最小二乘法來求解,并可用F—檢驗來觀察其擬合情況,察看其顯著性水平。(有關用最小二乘法求解方程待定系數的問題本文不再贅述)
將上述二維曲線回歸的方法編程在計算機上運行,可以使問題輕松得到解決,可以用C語言來編寫,先將二維曲線在X=0點附近的展開式假設為:
f(X)=a0+a1X+a2X2+…+anXn
試驗數據代入后,上式成為關于 a0,a1,a2,…,an的待定系數方程,用矩陣最小二乘法求解,得出回歸曲線,并顯示其相關系數R,(這里0≤R≤l,R值愈接近于0說明回歸效果越好;愈接近于1回歸效果越差)。可根據R值與0和1的接近程度來評價曲線的回歸教果,決定是否增加展開式項數n的值,作下一次回歸試驗。若回歸效果較好,即可進一步實施F一檢驗。在檢驗過程中,計算機先后要求輸入在0.1,0.05,0.01置信度下的F一檢驗值。并與E=(N-K)U/KQ的值相比較,判斷在相應顯著性水平下,曲線的回歸效果是否顯著。
如果當函數展開式項數n的值增加很大,以至于計算機內存發生溢出,仍不能求出擬合效果良好的曲線時,可根據情況確定在非零點求出該曲線的展開式。這時只要根據程序運行要求,結合實際情況,確定輸入X0的值,即可同上運行,直到求出理想的曲線多項式為止。
在科技攻關項目:翟坡試區作物水環境自動監測及預測預報研究課題中,購置的土壤墑情傳感器的頻率與土壤墑情變化曲線,是生產廠家在實驗室里作的,與安裝現場實際情況差距很大。為保證所測墑情的準確性,我們對購置的墑情傳感器逐個進行了曲線律定。這些曲線就是用該軟件進行回歸的。
在各傳感器的埋設現場取土樣烘干得到的土壤含水率數據,與相應的傳感器頻率觀測數據資料如表1。
在應用該軟件求其中某一個墑情傳感器的變化曲線時,首先將相應的實測數據資料按順序輸入程序的數值語句中,根據實際情況選定適當的X0點(一般情況下先選定X0=0),在X0點附近將函數展開成多頂式。然后,確定展開式的頂數,一般從一項開始試算,每次增加一項。試算期間可以只觀測相關系數R的值與0和1的接近程度,不作F-檢驗。通過試算可以找出R值最接近于0的一次試算,確定其展開式的頂數n的值。再根據n的值進一步運行軟件求出曲線,并作F一檢驗,若能通過F-檢驗,則該曲線就是所需求的土壤墑情與傳感器頻率變化函數。若不能通過F-檢驗,則需要更換展開點X0的值繼續試算。
用此方法和該軟件回歸曲線時,還應注意的一個問題是:數據資料的選取應稍超出實際應用曲線定義域的范圍,防止回歸出來的曲線在定義域邊界附近出現變異現象,以保證所應用曲線的正確性和完整性。
根據上述實際觀測資料描出數據點狀圖,依據點狀圖發展趨勢,補充邊界數據如表2。

表2 補充邊界數據表
將上述資料輸入該軟件的數據語句,運行后分別得出回歸結果如下:
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
當 α=0.1 時,Fa(5,28)=2.06,Et>Fa(5,28),OK!
當 α=0.05 時,Fa(5,28)=2.56,Et>Fa(5,28),OK!
當 α=0.01 時,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
當 α=0.1 時,Fa(5,23)=2.13,Et>Fa(5,23),OK!
當 α=0.05 時,Fa(5,23)=2.64,Et>Fa(5,23),OK!
當 α=0.01 時,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
當 α=0.1 時,Fa(5,23)=2.13,Et>Fa(5,23),OK!
當 α=0.05 時,Fa(5,23)=2.64,Et>Fa(5,23),OK!
當 α=0.01 時,Fa(5,23)=3.94,Et>Fa(5,23),OK!
由上述各自的F-撿驗可知:上述三條曲線均在α=0.01顯著性水平下回歸效果顯著,墑情傳感器能夠反映土壤墑情的變化規律,是要求的函數表達式。
為驗證所求出的土壤墑情與傳感器頻率變化曲線,我們繼續進行取土樣烘干試驗,對烘干得出的土壤含水率值與相應的變化函數求出的土壤含水率值進行了比較。如表3。
33#分站20cm埋深09#傳感器頻率與土壤含水率變化函數
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 計算值與墑情烘干值比較表
從表3可以看出:雖然回歸曲線函數計算值與土樣烘干值絕對誤差小于0.5的僅占45.45%,但是差值小于l的卻占81.82%,對±壤墑情來說,這是可以接受的。因此,該曲線函數能反映土壤墑情傳感器頻率值與土壤含水率之間的變化規律。
33#分站50cm理深11#傳感器頻率與土壤含水率變化函數
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 計算值與墑情烘干值比較表
從表4可以看出:回歸曲線函數計算值與土樣烘干值絕對誤差小于0.5的占63.64%,絕對差值小于l的占81.82%,說明用該軟件回歸出來的曲線函數能較好地反映土壤墑情傳感器頻率值與土壤含水率之間的變化規律。
33#分站80cm埋深08#傳感器頻率與土壤含水率變化函數
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 計算值與墑情烘干值比較表
從表5可以看出:回歸曲線函數計算值與土樣烘干值絕對誤差小于1的占72.73%,說明用該軟件回歸出來的曲線函數能較好地反映土壤墑情傳感器頻率值與土壤含水率之間的變化規律。
通過上述驗證,說明用二維參量曲線回歸方法和軟件得出的曲線函數能較好地反映傳感器頻率與土壤墑情的變化規律。
用該軟件作曲線回歸時,只需將原始數據資料中的個別明顯非點刪除,即可輸入計算機進行回歸運算,避免了人為干預的因素,曲線的走向完全根據數據分布規律而定。不象插值法和其它方法那樣,先要人為地刪除許多所謂不規律的點,再根據觀察保留有規律的曲線附近的點,這樣難免存在觀察者厚此薄彼的現象,以致影響理論曲線函數與客觀規律的吻合程度。
[1]馬玉芳;陳建華;郝揚滿.基于C語言對三次樣條函數的求解及程序.價值工程.2011.11
[2]蘇俊杰.基于矩陣樣條函數的二階矩陣微分方程數值解法研究.合肥工業大學.2010.03
[3]張曉丹;邵帥;劉欽圣.基于樣條函數的光滑支持向量機模型.北京科技大學學報.2012.7
[4]邱慧敏;楊濟安.樣條函數與樣條小波.重慶郵電學院學報(自然科學版).2000-03