張 露，楊鵬輝，張德鑫，闞 劍

(安徽財經(jīng)大學(xué)統(tǒng)計與應(yīng)用數(shù)學(xué)學(xué)院，安徽 蚌埠 233030)

太陽影子定位技術(shù)研究

張 露，楊鵬輝，張德鑫，闞 劍

(安徽財經(jīng)大學(xué)統(tǒng)計與應(yīng)用數(shù)學(xué)學(xué)院，安徽 蚌埠 233030)

目的 針對太陽影子定位問題，建立一系列數(shù)學(xué)模型對視頻中的數(shù)據(jù)進行研究，并通過太陽影子定位技術(shù)確定視頻拍攝時間與地點。方法 以一段時長30 min拍攝時間與地點未知的直桿影長變化視頻為研究對象，根據(jù)全國大學(xué)生數(shù)學(xué)建模網(wǎng)提供的數(shù)據(jù)，運用最小二乘擬合、構(gòu)建方程組求解、畫圖分析以及圖片灰度處理、銳化處理等方法，分別構(gòu)建了影長變化模型、影子定位模型和影子定時模型等，并依據(jù)影子的坐標(biāo)點，結(jié)合MATLAB、Excel等軟件得出了影子長度的變化規(guī)律、目標(biāo)物所在地的經(jīng)緯度以及測量數(shù)據(jù)的時間。結(jié)果 經(jīng)畫圖分析得知2015年10月22日北京時間9∶00～15∶00之間天安門廣場(北緯39度54分26秒,東經(jīng)116度23分29秒)3 m高直桿的太陽影子先隨時間的推移而變短，在中午某一時刻達到最小值，后又隨時間的變化而變長；在只知直桿在水平地面上的太陽影子頂點坐標(biāo)的情況下，求得實例中直桿所在地點為(37.707 23 °N,71.788 99 °E)，測量日期為11月10號；在只獲取一段視頻的情況下，經(jīng)過對圖片截取、灰度化以及銳化等處理讀取出視頻中直桿影子頂端坐標(biāo)，并求得該視頻的拍攝地點與時間分別為(40.16 °N,122.70 °E)、2015年7月26日。結(jié)論 在視頻中物體所投射的影子坐標(biāo)可讀取的情況下，運用合理的數(shù)學(xué)模型，可以通過太陽影子的變化對視頻的拍攝時間及其地點進行精準(zhǔn)的判斷。

最小二乘擬合法；太陽高度角；MATLAB；Excel

在古代，人們就利用太陽投射的影子來測定時刻，并發(fā)明了聞名世界的計時儀器“日晷”。漢書天文志有“夏至至于東井，北近極，故晷短；立八尺之表，而晷景長尺五寸八分。冬至至于牽牛，遠極，故晷長；立八尺之表，而晷景長丈三尺一寸四分。春秋分日至婁、角，去極中，而晷中；立八尺之表，而晷景長七尺三寸六分。”可見古人對利用太陽影子進行定位進行過一定的研究。在現(xiàn)代，如何確定視頻的拍攝地點和拍攝日期則成為視頻數(shù)據(jù)分析的重要方面，而太陽影子定位技術(shù)就是通過分析物體的太陽影子變化，確定視頻拍攝地點和日期的一種重要方法。眾所周知，視頻數(shù)據(jù)的分析在現(xiàn)代生活中具有重要的意義，特別是對于刑事案件的偵查起著突破性的作用。國內(nèi)外學(xué)者基于日晷原理，結(jié)合現(xiàn)代數(shù)學(xué)軟件探索太陽影子定位技術(shù)。但這項技術(shù)還不夠成熟，有待學(xué)者們繼續(xù)深入研究，以便將其應(yīng)用于現(xiàn)代通信與監(jiān)測等領(lǐng)域。基于這樣的研究背景，運用赤緯角算法和時差算法等現(xiàn)代數(shù)學(xué)知識對太陽影子定位技術(shù)做了一定深度的研究。

1 數(shù)據(jù)來源與模型假設(shè)

數(shù)據(jù)主要來源于全國大學(xué)生數(shù)學(xué)建模比賽。為保證文章的嚴(yán)密性，現(xiàn)做出以下幾條假設(shè)：(1)設(shè)地球大氣層對太陽高度角沒有影響，即忽略大氣層對實際觀測的太陽高度角和理論值的不同影響;(2)設(shè)太陽、地球都是一個表面光滑的正球體；(3)太陽光為一束平行的光束；(4)目標(biāo)地直桿與水平地面呈垂直狀態(tài)；(5)時差算法以北京時間為準(zhǔn)。

2 影子長度的變化規(guī)律

2.1 研究思路

建立影子長度變化的數(shù)學(xué)模型，并且分析影子長度關(guān)于各個參數(shù)的變化規(guī)律。在忽略大氣對光線折射的因素下，易知影子的長度與目標(biāo)所在地的經(jīng)度、緯度、太陽高度角等因素存在某種相關(guān)性。于是首先嘗試建立影子長度與經(jīng)緯度、太陽高度角等各因素之間的關(guān)系，分析影子長度關(guān)于各個參數(shù)的變化規(guī)律；然后基于各因素之間的變化規(guī)律，利用抽象幾何的畫圖研究方法，推導(dǎo)出太陽高度角與當(dāng)?shù)亟?jīng)緯度、太陽赤緯角及時角之間的函數(shù)變化關(guān)系，從而推導(dǎo)出太陽赤緯角以及時角的計算公式；最后利用時差算法及赤緯角算法對模型進行修正并推導(dǎo)出影長變化公式，帶入實例數(shù)據(jù)，運用MATLAB編程畫出實例中直桿的太陽影子長度的變化曲線，求其變化范圍。

2.2 研究方法

2.2.1 模型的建立

圖1 影長模型圖

若需求得直桿影長，需先求得角ε。又易知ε為太陽高度角，太陽高度角由目標(biāo)所在地的地理緯度γ、當(dāng)日的太陽赤緯ω以及當(dāng)時的太陽時角τ等因素共同決定[1]。各角度幾何圖示如圖2。

已知各個角度之間的關(guān)系滿足如下公式：

sinε=sinγsinω+cosγcosωcosτ

圖2 幾何角度圖

其中太陽赤緯ω是由觀測當(dāng)日日期在1年中的序號n所決定：

sinω=0.39795cos[0.98563(n-173)]

太陽時角τ則是由真太陽時T所決定：

τ=15×(T-12)

而真太陽時又是由北京時間T′及觀測地時差m兩者共同決定：

T=T′-m

這其中時差又是隨經(jīng)度φ變化的函數(shù)：

m=(120 °-φ)/15 °

建立以上關(guān)系式可以得出太陽高度角關(guān)于目標(biāo)所在地的地理緯度γ、觀測當(dāng)日日期在1年中的序號n、觀測點的北京時間T′以及觀測點經(jīng)度φ的變化關(guān)系式如下[2]：

其中

sinω=0.39795cos[0.98563(n-173)],m=(120°-φ)/15°

綜上可知，求1根直桿的影子長度，需了解該直桿長度h、所在地地理緯度γ、觀測當(dāng)日日期在一年中的序號n、觀測點的北京時間以及觀測點經(jīng)度φ。

圖3 直桿影子長度變化曲線圖

2.2.2 模型的求解

已知2015年10月22日北京時間9∶00～15∶00之間天安門廣場(北緯39度54分26秒,東經(jīng)116度23分29秒)有3 m高的直桿，即有h=3，γ=116.3914，n=295，φ=39.9072，T′=t(t∈[9,13])。將上述數(shù)據(jù)代入所求公式,運用MATLAB畫圖可得出該時間段觀測點直桿影子長度變化曲線如圖3。

觀察圖3中曲線可知，該時間段內(nèi)影子長度先隨時間的推移而變短，在中午某一時刻達到最小值，后又隨時間的推移而變長，整個時間段影長變化范圍為3.56～7.67 m。

3 基于影子變化的定位技術(shù)

3.1 研究思路

根據(jù)某固定直桿在水平地面上的太陽影子的變化，建立數(shù)學(xué)模型,確定直桿所處的地點。首先要獲取影子的頂點坐標(biāo)。對于在水平地面影子頂點坐標(biāo)的獲取，是以直桿底端為原點、水平地面為xy平面建立的坐標(biāo)系，雖然坐標(biāo)軸方向的選取影響了影子的坐標(biāo)，但是對于1天內(nèi)影子頂點的運動軌跡以及影長的變化規(guī)律不會產(chǎn)生影響，并且影長的最短處對應(yīng)的當(dāng)?shù)貢r間應(yīng)為正午12點。以時間變化為橫坐標(biāo)，影長變化為縱坐標(biāo)建立直角坐標(biāo)系，在坐標(biāo)系內(nèi)運用影子定點坐標(biāo)擬合出影長的變化曲線圖，求出影長的最低點；再根據(jù)測量時的背景時間以及測量的間隔時間，推算出測量地的正午12點和北京時間的時間差，從而根據(jù)時差和經(jīng)度的關(guān)系推算出測量地的經(jīng)度。由于直桿長度與位置無關(guān)，所以可以通過與太陽高度角的比值忽略掉直桿長度的影響，進而就可代入數(shù)據(jù)求出直桿所在的緯度，確定測量點的位置。

3.2 研究方法

3.2.1 對經(jīng)度的求解

I 數(shù)據(jù)處理

表1 不同時刻影子長度對應(yīng)表

II 最小二乘擬合

以時間(T)為橫坐標(biāo)，影子長度(L)為縱坐標(biāo)建立影子長度與時間變化的坐標(biāo)系，在建立這個坐標(biāo)系之前，需要對橫坐標(biāo)即時間變化做一定的處理，把時間的分鐘數(shù)同除以60換算成小時制，然后結(jié)合影子長度建立影子長度變化的坐標(biāo)系。在坐標(biāo)系內(nèi)將測量所得的數(shù)據(jù)運用MATLAB進行擬合，得出時間與影長的函數(shù)y=ax2+bx+c(a≠0，并且利用最小二乘法判斷其擬合是否合理。

3.2.2 對緯度的求解

在上述太陽高度角、緯度、日期和太陽赤緯的關(guān)系中，已知日期、太陽赤緯是與真太陽時(當(dāng)?shù)貢r間)構(gòu)成某種函數(shù)關(guān)系。若以當(dāng)?shù)貢r間作為自變量，則只有太陽高度角和緯度是需要求解的。當(dāng)建立太陽方位角公式時，太陽方位角就可以通過影子頂點的縱坐標(biāo)與影長的比值解出，所得結(jié)果即為太陽高度角與緯度之間的關(guān)系。

已知太陽方位角公式為：

建立公式sinε=sinγsinω+cosγcosωcosτ，化簡可得到如下方程：

代入數(shù)值即可求出當(dāng)?shù)鼐暥圈肹4]。

4 基于影子變化的定時技術(shù)

4.1 研究思路

根據(jù)某固定直桿在水平地面上的太陽影子頂點坐標(biāo)數(shù)據(jù)，建立數(shù)學(xué)模型，確定直桿所處的地點和日期。首先根據(jù)太陽高度角和太陽方位角的公式，推導(dǎo)出日期和緯度的聯(lián)合公式，將所測得的影子坐標(biāo)進行處理，得出太陽方位角的余弦值，再將不同的影子長度和時間所對應(yīng)的坐標(biāo)帶入方程，運用MATLAB程序得出緯度值和日期，從而得出可能的測量地點和日期。

4.2 研究方法

4.2.1 對時間的求解

在地球自轉(zhuǎn)和公轉(zhuǎn)的過程中，某一地點的正午太陽高度角既和所在地的緯度有關(guān)，又與測量的日期以及在該日期測量的時間有關(guān)。同一緯度不同日期的正午太陽高度角不同，同一日期不同緯度的正午太陽高度角也不同，因此在給定影長時，測量地的位置以及測量的日期都是可以確定的。

在如下公式中：

sinε=sinγsinω+cosγcosωcosτ

對太陽高度角求其反函數(shù)，即可得出太陽高度角ε的表達式為：

ε=arcsin(sinγsinω+cosγcosωcosτ)

由圖1可看出，太陽高度角ε也可表示為：

因此同一地點不同時刻的太陽高度角的比值方程為：

此時，同一地點不同時刻的太陽高度角ε正切值的比值與影長相關(guān)，與直桿的高度無關(guān)。

將ε的表達式帶入太陽高度角的比值方程，即可得出緯度日期模型：[5]

模型中的緯度、測量日期、太陽赤緯共同決定了測量地的太陽高度角的數(shù)值，將實際測量的影子長度和時間數(shù)據(jù)帶入模型，運用MATLAB程序可解出測量地的緯度和測量日期。

4.2.2 舉例說明

有一直竿垂直于地面，以直竿底端為原點、水平地面為xy平面建立坐標(biāo)系，其影子坐標(biāo)變化如下表2：

表2 影子坐標(biāo)變化表

圖4 影子落點坐標(biāo)圖

以直桿底端為原點，東西方向為橫軸，南北方向為縱軸，建立影子落點的坐標(biāo)系，將上表中的數(shù)據(jù)在坐標(biāo)系中畫出，如圖4所示：

根據(jù)影子落點坐標(biāo)圖，對應(yīng)表2中測量每個點的北京時間，首先對北京時間進行處理，把時間的分鐘數(shù)換算成小時制，如表3。

表3 處理時間表

表4 不同時刻影子長度對應(yīng)表

將不同時刻影子長度對應(yīng)表中的數(shù)據(jù)，運用MATLAB程序進行曲線擬合，擬合曲線如圖5所示：

圖5 影子長度變化圖

得出影子長度和時間的二次方程：

L=0.09812T2-2.985T=23.319

誤差的平方和為4.402588e-0.07，由于誤差和極小于1，故此函數(shù)擬合程度較高。

在此方程中，影長l的最低點對應(yīng)的時間即為測量地的正午12點所對應(yīng)的北京時間：

此時測量地與北京時間的時差為：

則當(dāng)測量地時間為正午12點時，對應(yīng)的北京時間為15.214 07，由于Δt>0，故測量地在北京時間所在時區(qū)的經(jīng)度120 °E的西邊，由于地理上經(jīng)度相隔1 °時，時差為4min，故測量地經(jīng)度由公式：

可求得：

φ=71.78899°E

因此測量地所在的經(jīng)度為71.78899 °。

由表2數(shù)據(jù)中，選取前2個點對應(yīng)的影長和當(dāng)?shù)貢r間，帶入緯度日期模型，因為：

τ=15×(T-12)

借助MATLAB程序，運行后結(jié)果為：

x=37.7072，y=162.5577

x即所求緯度，y是太陽赤緯。

由公式：

sinω=0.39795cos[0.98563(n-173)]

可反推出測量日期在1年中的序號n的表達式為：

將程序運行出的y(即ω)值帶入反推出的公式，計算出n≈314，即11月10日。

結(jié)合經(jīng)度模型和緯度日期模型，案例中數(shù)據(jù)求解結(jié)果為：該地點地理位置約為：37.70723 °N,71.78899 °E，測量日期為11月10日。

5 基于影長變化對視頻拍攝時間地點的判斷

5.1 研究思路

在上述問題中已經(jīng)解決了在目標(biāo)物影子坐標(biāo)已知的情況下，根據(jù)太陽影子求目標(biāo)物所在地及拍攝時間等問題。下面探討如何根據(jù)太陽影子直接從一段視頻中讀取其拍攝時間及地點等信息。利用太陽影子定位技術(shù)，首先要從視頻中讀取影子的信息，即影子隨時間變化其坐標(biāo)點的變化。然后利用太陽高度角和太陽方位角的公式，推導(dǎo)出日期和緯度的聯(lián)合公式，將所測得的影子坐標(biāo)進行處理，得出太陽方位角的余弦值，再將不同影子長度和時間所對應(yīng)的坐標(biāo)帶入方程，運用MATLAB程序得出緯度值和日期，從而得出可能的測量地點和日期。在此問題中對影子坐標(biāo)的讀取最為關(guān)鍵[6]。

5.2 研究方法

5.2.1 視頻數(shù)據(jù)的讀取方法

首先要將視頻數(shù)據(jù)轉(zhuǎn)換成圖片數(shù)據(jù)，即等時間段對視頻進行截圖處理，但憑主觀感覺進行截取誤差太大，故需要運用MATLAB對視頻數(shù)據(jù)進行讀取處理。由于視頻是由一幀一幀的圖片組合而成，所以要想讀取視頻中的圖片數(shù)據(jù)需先弄清楚這段視頻由多少幀圖片組合而成。

把所需讀取的視頻導(dǎo)入MATLAB中，讀取視頻信息如圖6。

圖6 視頻讀取信息

從視頻信息中可知該視頻時長為2 440s，該視頻的幀率為25Hz，即每秒25幀圖片，所以該視頻總共有61 000幀。若每兩分鐘截取一張圖片，總共可得21張圖片。[7]

5.2.2 影子頂點坐標(biāo)的讀取方法

1張圖片是由無數(shù)個像素點組成，可以利用MATLAB讀取各像素點坐標(biāo)，從而推導(dǎo)出直桿影子頂點坐標(biāo)，為降低計算量，先對圖片進行一定的處理。根據(jù)主觀判斷選取直桿底部一點讀取該像素點坐標(biāo)為(866,882)，然后再讀取21張圖片影子頂點的像素點坐標(biāo)，發(fā)現(xiàn)其變化范圍為x∈(867,1 700),y∈(780,930)，故截取圖片影子所在的一部分進行灰度處理，截取圖片如圖7。

圖7 影子截圖

圖8 灰度處理效果圖

為了使影子圖像更加清晰，需要利用圖像銳化技術(shù)，使直桿影子的邊緣變得清晰，從而使得影子頂點坐標(biāo)的讀取更加精準(zhǔn)。圖像銳化處理的目的是使圖像的邊緣、輪廓線以及圖像的細節(jié)變的更加清晰，但平滑的圖像變得模糊的根本原因是圖像受到了平均或積分運算的影響，因此需要對其進行逆運算，使圖像變得更為清晰[9]。利用Photo-Shop的圖像銳化處理，銳化處理效果如圖9所示。

圖9 銳化處理效果圖

經(jīng)過上述處理，影子頂點已經(jīng)變得清晰可見，只需讀取出影子頂點坐標(biāo)即可。但讀取出的坐標(biāo)僅為像素坐標(biāo)，還需把其轉(zhuǎn)化為實際的坐標(biāo)才能進行影子定位技術(shù)運算。實際坐標(biāo)可由攝像機坐標(biāo)系與圖像坐標(biāo)系共同推導(dǎo)得出。設(shè)在圖像坐標(biāo)系中，以圖片中直桿底端為原點，以毫米為單位的影子頂點坐標(biāo)為(α,β)，以像素為單位的影子頂點坐標(biāo)為(x,y)，且影子頂點像素點在x與y軸方向的物理尺寸為Δx,Δy；在攝像機坐標(biāo)系中，以攝像機為坐標(biāo)原點的影子頂點坐標(biāo)為(φ,φ,γ)，且焦距f近似為攝像機與圖像坐標(biāo)系原點之間的直線距離；在實物坐標(biāo)系中，以實際的直桿底端為坐標(biāo)原點，影子頂點坐標(biāo)為(x,y,z)，則三個坐標(biāo)系下的影子頂點坐標(biāo)之間的轉(zhuǎn)化關(guān)系為：

(x,y,z,1)T=w(φ,φ,γ,1)T(w為4×4階的矩陣，是對攝影機坐標(biāo)系下的坐標(biāo)進行旋轉(zhuǎn)變換所得)

運用MATLAB進行上述坐標(biāo)變換，最終得到影子頂點坐標(biāo)，如表5所列。

表5 影子頂點坐標(biāo)

根據(jù)上面敘述的方法，用影子頂點坐標(biāo)求得該視頻的拍攝地點與時間為(40.16 °N,122.70 °E)、2015年7月26日。

6 討論與結(jié)論

物體的影長受到物體自身的體積、物體所在地的經(jīng)度和緯度、觀測日期及時間點等因素的影響。要通過物體影子的長度及影子坐標(biāo)來確定物體的經(jīng)緯度，需要根據(jù)觀測日期、太陽高度角以及方位角的公式來建立模型，確定物體的位置；在失去觀測日期等條件時，進一步優(yōu)化模型，只需利用目標(biāo)物影子頂點坐標(biāo)變化也可求出拍攝時間與地點；而當(dāng)影子頂點坐標(biāo)這一唯一條件也失去時，可以利用圖片截取、圖片灰度處理、圖片銳化處理以及坐標(biāo)系轉(zhuǎn)換等方法重現(xiàn)影子頂點坐標(biāo)，從而求取視頻的拍攝時間與地點。經(jīng)過上述分析與討論可知在只獲取一段視頻，而其他任何條件未知的情況下，通過一系列的圖片處理與數(shù)學(xué)演算，求取視頻的拍攝時間與地點是可行的。

[1]劉偉峰,謝永杰,陳若望,等.天頂亮度與太陽高度角關(guān)系的觀測[J].光電工程,2012,(07):49-54.

[2]栗琳,胡勇,鞏彩蘭,等.太陽高度角對圖像能量的影響及其校正[J].大氣與環(huán)境光學(xué)學(xué)報,2013,(01):11-17.

[3]王國安,米鴻濤,鄧天宏,等.太陽高度角和日出日落時刻太陽方位角一年變化范圍的計算[J].氣象與環(huán)境科學(xué),2007,(01):161-164.

[4]李先華,黃雪樵,池天河,等.衛(wèi)片像元太陽高度角和方位角的計算原理與方法[J].測繪學(xué)報,1993,(02):149-154.

[5]楊長青,李華,徐寶芳.“正午太陽高度角模型”的設(shè)計、制作與運用[J].中學(xué)地理教學(xué)參考,2011,(08):26-28.

[6]楊桂元.數(shù)學(xué)建模方法[M].合肥：中國科學(xué)技術(shù)出版社，2014：1-58.

[7]任利平.視頻中關(guān)鍵幀提取技術(shù)的研究[D].蘭州：蘭州大學(xué),2011.

[8]李貞培,李平,郭新宇,等.三種基于GDI+的圖像灰度化實現(xiàn)方法[J].計算機技術(shù)與發(fā)展,2009,19(07):73-75+79.

[9]靳佳澍.一種針對彩色二維碼圖像的二值化方法[J].科技與企業(yè),2016,(04):83-84.

[10]潘國榮,周躍寅.兩種坐標(biāo)系轉(zhuǎn)換計算方法的比較[J].大地測量與地球動力學(xué),2011,(03):58-62.

[責(zé)任編輯：劉守義 英文編輯：劉彥哲]

Sun Shadow Positioning Technology

ZHANG Lu,YANG Peng-hui,ZHANG De-xin,KAN Jian

(School of Statistics and Applied Math，Anhui University of Finance and Economics,Bengbu,Anhui 233030,China)

Objective The paper,according to the sun shadow positioning issue,establishes a series of mathematical models to study the data of the video,and determines the shooting time and site of the video through the sun shadow positioning technology.Methods Taking the video of 30-minute-long and unknown shooting site reflecting the length variation of the straight rod’s shadow as the research object,this paper respectively constructed the shadow’s length variation model,the shadow positioning model,the shadow timing model and other models by applying the least square fitting,the equation set,the graph drawing analysis,and the image processing of the gray scale and sharpening,based on the data provided by the contemporary undergraduate mathematical website in modeling.The paper also concluded the law of the shadow’s length variation,the longitude and latitude of the site in the object,and the time of the data measurement by MATLAB,Excel and other software,based on the coordinates of the shadow.Results Through the graph drawing analysis,it was found that on October 22,2015,Beijing time between 9∶00-15∶00,the sun shadow of the 3-meter-high straight rod firstly got short with time passing by,reached the minimum in a certain moment of the noon,and got long as the time varied.Under the circumstance that only the coordinate of the shadow’s vertex on the rod horizontal ground was known,it was calculated that in instance 1,the site of the straight rod was in the coordinate and the date of the measurement was November 10th.Under the circumstance there was only one piece of the video,the coordinate of the shadow’s vertex of the straight rod,the shooting site ofand the shooting date of October 22,2015.Conclusion In the case that the coordinate of the shadow projected by the object can be read;the shooting time and site of the video can be judged concisely through analyzing the changes of the sun shadow by applying the reasonable mathematical models.

least square fitting;solar elevation angle;MATLAB;Excel

國家自然科學(xué)基金項目(11301001)；安徽財經(jīng)大學(xué)教研項目(acjyzd201429)

張露(1995-)，女，安徽六安人，安徽財經(jīng)大學(xué)統(tǒng)計與應(yīng)用數(shù)學(xué)學(xué)院在讀學(xué)生，研究方向：數(shù)學(xué)與應(yīng)用數(shù)學(xué)。

楊鵬輝(1981-)，女，安徽淮南人，安徽財經(jīng)大學(xué)統(tǒng)計與應(yīng)用數(shù)學(xué)學(xué)院講師，碩士，研究方向：應(yīng)用數(shù)學(xué)與數(shù)學(xué)建模。

O 242

A

10.3969/j.issn.1673-1492.2016.11.007

來稿日期：2016-04-18