張世輝韓德偉何 歡（燕山大學信息科學與工程學院 秦皇島 066004）

②（河北省計算機虛擬技術與系統集成重點實驗室 秦皇島 066004）

利用視覺目標遮擋和輪廓信息確定下一最佳觀測方位

張世輝*①②韓德偉①何 歡①①（燕山大學信息科學與工程學院 秦皇島 066004）

②（河北省計算機虛擬技術與系統集成重點實驗室 秦皇島 066004）

下一最佳觀測方位的確定是視覺領域一個比較困難的問題。該文提出一種基于視覺目標深度圖像利用遮擋和輪廓信息確定下一最佳觀測方位的方法。該方法首先對當前觀測方位下獲取的視覺目標深度圖像進行遮擋檢測。其次根據深度圖像遮擋檢測結果和視覺目標輪廓構建未知區(qū)域，并采用類三角剖分方式對各未知區(qū)域進行建模。然后根據建模所得的各小三角形的中點、法向量、面積等信息構造目標函數。最后通過對目標函數的優(yōu)化求解得到下一最佳觀測方位。實驗結果表明所提方法可行且有效。

深度圖像；遮擋；輪廓；未知區(qū)域；類三角剖分；下一最佳觀測方位

1 引言

下一最佳觀測方位（Next Best View， NBV）的確定是指根據攝像機（傳感器）在當前觀測方位下獲取到的視覺目標信息確定出攝像機的下一觀測方位，使得攝像機在所求方位下能夠最大量地獲取到當前觀測方位下未采集到的信息。下一最佳觀測方位的確定一直是數據采集、3維重建、機器人導航、自動裝配、目標識別等領域重要且困難的研究課題之一，受到學者們的廣泛關注［1-4］。

目前，下一最佳觀測方位確定算法所基于的圖像信息類型主要有兩種：圖像亮度信息［5，6］和圖像深度信息［7-11］。同亮度圖像相比，深度圖像更利于獲得場景的3維信息，因此，現有的下一最佳觀測方位確定算法通常都基于深度圖像實現。文獻［7］采用八叉樹模型描述視覺目標，并根據每個節(jié)點的實際情況對其進行標記，進而確定下一最佳觀測方位。該方法是求解下一最佳觀測方位問題方面最早的方法，但該方法由于未考慮遮擋因素，故當視覺目標出現遮擋現象時，無法保證所得結果的正確性和有效性。文獻［8］通過對視覺目標遮擋區(qū)域進行多邊形建模，利用遮擋信息確定下一最佳觀測方位，但該方法需要特定的設備，適用范圍有限。文獻［9］提出一種基于整體觀測策略確定下一最佳觀測方位的方法，該方法要求對視覺目標有一定的了解且該方法確定出的攝像機位置被限定在一個固定的表面，方法的通用性受到了限制。文獻［10］提出一種基于 B樣條計算信息熵的視點規(guī)劃方法，該方法將攝像機獲得模型最大信息熵的方位作為下一最佳觀測方位，但該方法沒有充分利用視覺目標背面的信息，并且計算出的攝像機位置同樣被限定在一個固定的表面。文獻［11］提出了一種基于射線追蹤的下一最佳觀測方位確定方法，該方法的思想是從若干候選觀測方位中確定出下一最佳觀測方位，導致所求結果被限定在事先確定的候選觀測方位中。文獻［12］利用當前觀測方位下未觀測到區(qū)域的置信模型來估計每個可能的觀測點的信息增益，提出一種靈活的概率框架來估計下一最佳觀測方位，但該方法涉及到的可觀測性參數需要人為設定。

針對已有方法存在的未考慮遮擋因素、攝像機位置受限等方面的不足，本文提出一種基于深度圖像利用遮擋和輪廓信息確定下一最佳觀測方位的方法。

2 方法概述

2.1相關概念

當攝像機觀測視覺目標時，視覺目標的一部分表面會遮擋另外一部分表面。遮擋邊界是當前觀測方位下視覺目標上遮擋部分與被遮擋部分的分界線，該分界線由遮擋部分的像素點組成，且每個遮擋邊界點都在各自的八鄰域點內對應一個深度差最大點，這個深度差最大點稱為該遮擋邊界點的下鄰接邊界點，下鄰接邊界點組合形成下鄰接邊界。圖1為深度攝像機觀測視覺目標并獲取視覺目標深度圖像的示意圖。其中，紅色線表示當前觀測方位下獲取到的視覺目標的輪廓，黃色線表示遮擋邊界，藍色線表示下鄰接邊界。

分析圖1可知，本文把由遮擋邊界及其對應的下鄰接邊界組成的區(qū)域及當前觀測方位下獲取到的視覺目標輪廓所圍成的視覺目標的背面區(qū)域統稱為未知區(qū)域。對未知區(qū)域的觀測可以得到視覺目標更多的信息。為了獲取視覺目標更多的信息，需要對視覺目標進行下一次觀測。在所有下一觀測方位中能夠最多地獲取到視覺目標未知信息的觀測方位稱為下一最佳觀測方位。

2.2 總體思想

本文方法總體思想如下。首先，獲取視覺目標的一幅深度圖像，同時記錄下攝像機的內外參數；其次，根據深度圖像中的遮擋和輪廓信息構建未知區(qū)域；再次，對未知區(qū)域進行類三角剖分建模，并計算出各小三角形的中點、法向量和面積；然后，基于所有小三角形信息構造出以攝像機位置為自變量，以該觀測方位下可以被實際觀測到的小三角形面積為因變量的函數，將下一最佳觀測方位求解問題轉化為函數優(yōu)化問題；最后，通過優(yōu)化算法求解目標函數確定出下一最佳觀測方位。

圖1 攝像機觀測并獲取視覺目標深度圖像示意圖

3 下一最佳觀測方位確定方法

3.1 構建未知區(qū)域并計算各未知區(qū)域的中點

首先利用文獻［13］中的遮擋邊界檢測方法檢測出深度圖像中的遮擋邊界，然后把每條遮擋邊界與其對應的下鄰接邊界組成一個未知區(qū)域。此外，通過分析深度圖像中的像素點可知，八鄰域點內有且不全為背景點或數據點的數據點為輪廓點，根據輪廓點的這個特征提取出視覺目標的輪廓。

本文定義未知區(qū)域的中點為圍成該區(qū)域的所有像素點的3維中心點，即未知區(qū)域中點的計算公式為

3.2 未知區(qū)域建模

為了實現對遮擋邊界產生的未知區(qū)域的建模，還需確定出未知區(qū)域對應的遮擋邊界的端點。通過分析深度圖像中遮擋邊界點之間的位置關系可知，如果遮擋邊界點P'的八鄰域內只有一個遮擋邊界點或只有兩個相鄰的遮擋邊界點，那么遮擋邊界點P'一定為該遮擋邊界的端點。除以上兩種情況外遮擋邊界點P'均不為端點。判斷遮擋邊界點P'的八鄰域內的兩個遮擋邊界點是否相鄰的公式為

圖2展示了某一未知區(qū)域按上述方法建模后的效果。其中，黃色點代表遮擋邊界上的遮擋邊界點，藍色點代表與之對應的下鄰接邊界點，綠色點代表該區(qū)域的中點。

3.3 計算小三角形的中點、法向量和面積

未知區(qū)域建模后為了求解下一最佳觀測方位，需要確定各小三角形的中點、法向量和面積。小三角形的中點為構成小三角形的3個頂點的3維中心點，計算公式為

下面以圖3所示的某一未知區(qū)域中的一個小三角形ABC為例說明小三角形法向量的計算方法。

圖3中A， B為相鄰的遮擋邊界點或下鄰接邊界點，C為該未知區(qū)域的中點，T為小三角形的中點。構成小三角形的3個頂點的3維坐標兩兩相減可以得到6個共面的向量。根據向量叉乘的定義可以推出，上述6個向量中任意兩個不共線向量的叉積即為小三角形的法向量。為了保證計算出的小三角形實際可觀測面積的有效性，本文選取指向視覺目標外部的法向量作為小三角形的法向量。下面以6個向量中點A到點B的向量a和點A到點C的向量c為例詳細介紹小三角形法向量的計算方式。

為了得到指向視覺目標外部的法向量就要確定兩個向量的叉乘順序，即 a ×c還是 c ×a。下面以為例給出其計算公式：

式中，i， j， k分別為直角坐標系下與坐標軸相對應的3個單位向量分別為向量a和c的坐標表示。通過分析遮擋邊界和下鄰接邊界的性質可知：由位置關系不同的遮擋邊界與其下鄰接邊界所圍成的未知區(qū)域內的小三角形的法向量計算方式是不同的，但同一個區(qū)域內的小三角形的法向量的計算方式是相同的。確定兩個向量的叉乘順序的方法描述如下：

圖2 未知區(qū)域建模示意圖

圖3 小三角形的中點及其法向量示意圖

由圖4可以看出，為了保證所求得的小三角形法向量指向視覺目標的外部，圖4（a）， 4（b）， 4（c）， 4（d）4種情形下對應的小三角形法向量的計算方式應為點A'到點B'的向量叉乘上點A'到點C'的向量，對應為 a ×c。圖4（e）， 4（f）， 4（g）， 4（h）4種情形下對應的小三角形法向量的計算方式應為點A'到點 C'的向量叉乘上點A'到點B'的向量，對應為 c ×a。因此判斷某一未知區(qū)域內小三角形法向量的計算方式問題就可轉換成判斷該未知區(qū)域屬于圖4中哪種情形的問題。

然后，出于在一條遮擋邊界的兩個端點橫坐標相等時斜率k依然有意義以及簡化計算的考慮，根據上述分析把斜率k的計算方式改進為

式中，e為一個大于 0的小數?？紤]到極端情況，也就是當一條遮擋邊界只有兩個像素點組成時，如果e大于或等于 0.5就有可能造成遮擋邊界與其下鄰接邊界的位置關系的誤判，因此本文限定0 ＜ e＜ 0.5。這樣就可以把情形圖4（a）， 4（c）， 4（e），4（g）分別歸并到情形圖4（b）， 4（d）， 4（f）， 4（h）。點與直線位置關系的確定公式為

最后，通過判斷某一未知區(qū)域計算出的k和D的值來確定該未知區(qū)域內的小三角形法向量的計算方式。也就是說，當 k＜ 0且 D＞ 0時，即情形圖4（b），對應的法向量的計算方式為 a ×c；當 k＞ 0且D＜ 0時，即情形圖 4（d），對應的法向量的計算方式為 a ×c；當 k＜ 0且 D＜ 0時，即情形圖4（f），對應的法向量的計算方式為 c ×a；當 k＞ 0且 D＞0時，即情形圖 4（h），對應的法向量的計算方式為c ×a?？紤]到輪廓點所圍成的未知區(qū)域只有一種情形且比較簡單，又由于輪廓所代表的未知區(qū)域為視覺目標的背面，為了得到指向視覺目標外部的法向量，計算該區(qū)域內的小三角形的法向量時兩個向量的叉乘順序應為 a ×c。

圖4 遮擋邊界與其下鄰接邊界的位置關系示意圖

得到小三角形的法向量后就可以求出小三角形的面積。第i個小三角形的面積 si的計算公式為

3.4 下一最佳觀測方位問題轉化

由此可見只要計算出攝像機位置即可完成對下一最佳觀測方位問題的求解?；谏鲜龇治觯蠼庀乱蛔罴延^測方位下攝像機位置的公式為

式中，w為所有小三角形的個數，u為常數，u值越大函數越精確，但考慮到時空開銷，經多組實驗對比，本文選取 u= 300。 si為第i個小三角形的面積。 mi為第i個小三角形的法向量與向量夾角的余弦值。是一個激活函數：當 mi為正時，約等于1，即小三角形能被攝像機觀測到時小三角形參與計算，此時 mi與面積 si的乘積表示小三角形實際被攝像機觀測到的面積；當 mi為負時，約等于0，表示小三角形不能被攝像機觀測到，此時小三角形不參與計算。L為當前觀測方位下攝像機的觀測距離。

經上述處理后，下一最佳觀測方位求解問題就轉化為式（9）的函數優(yōu)化問題，通過對函數的優(yōu)化求解即可得出下一最佳觀測方位。

3.5 下一最佳觀測方位函數優(yōu)化求解

根據3.4節(jié)的描述可知，式（9）的最優(yōu)解即為下一最佳觀測方位下攝像機的位置。由于物體的全局幾何信息是未知的，所以本文要解決的是一個非凸問題，即式（9）是一個非凸函數。本文采用一種當前公認的處理中小規(guī)模非線性規(guī)劃問題的最好算法中的序列二次規(guī)劃［14］算法對式（9）進行求解。

本文方法根據攝像機在當前觀測方位下獲取到的視覺目標的輪廓及遮擋情況的不同利用式（10）計算得出初始迭代位置。通過分析小三角形的位置、大小和法向量定義初始迭代位置為

4 實驗及分析

4.1 實驗方案和數據集

為了驗證所提方法的可行性和有效性，本文做了3組實驗：（1）所提方法各階段結果的可視化實驗；（2）不同方法計算攝像機觀測點時對最終結果影響的對比實驗；（3）不同下一最佳觀測方位確定方法之間的對比實驗。實驗硬件環(huán)境為 Intel（R）Core（TM）2 Duo CPU E7500、內存2.0 G。算法采用C++編程實現。實驗分為仿真實驗和真實實驗，仿真實驗所用 Bunny， Duck， Dragon， Mole和Rocker 3維模型來自深度圖像領域較著名的Stuttgart Range Image Database［15］，真實實驗所用Bonsai， Box和Kettle是現實生活中的3個實物。仿真實驗中，通過OpenGL模擬攝像機觀測3維模型并獲取其深度圖像，設置投影矩陣的參數為（60，1，200，600），攝像機觀測距離為 300 mm。真實實驗中，用Kinect相機采集深度圖像，設定攝像機觀測距離為1.2 m，其余設置與仿真實驗一致。

4.2 實驗結果及分析

圖5 本文方法各階段結果的可視化圖

4.2.1 所提方法的可視化實驗為了更加形象地說明本文方法的求解過程，圖5給出了本文方法各階段結果的可視化圖。

由圖5展示的實驗結果可以看出，對于遮擋情況不太顯著的Duck和Kettle而言，由于其輪廓所圍成的未知區(qū)域內小三角形的面積在建模所得總的小三角形的面積中所占比重較大，故所求下一最佳觀測方位更加偏向于直接觀測視覺目標的背面；對于遮擋情況相對顯著的 Bunny， Dragon， Mole，Rocker， Bonsai和Box而言，由于其由遮擋邊界及下鄰接邊界建模所得三角形的面積在建模所得總的小三角形的面積中也占有較大比重，故這6個視覺目標所求下一最佳觀測方位兼顧遮擋和輪廓信息。由此可見，基于視覺目標遮擋及輪廓進行建模，并利用建模所得三角形信息進行下一最佳觀測方位的求解方案是可行的，且無論是數據集中的模型還是現實生活中的實物，本文方法所求下一最佳觀測方位都與人類視覺的觀測習慣相符。

4.2.2 不同觀測點計算方法對下一最佳觀測方位影響的對比實驗

攝像機觀測點的選取將對最終計算出的下一最佳觀測方位下獲取的信息產生較大的影響。為了驗證本文攝像機觀測點計算方法的有效性。我們對基于兩種不同的觀測點計算方法的實驗結果進行了比較。第1種為本文3.4節(jié)所提的攝像機觀測點計算方法，在此不再贅述。第2種攝像機觀測點的計算方法為

式（11）表明該計算方法下攝像機觀測點為以各小三角形面積為權重的小三角形中點坐標的平均值。

表1給出了攝像機在兩種不同觀測點計算方法所求解出的下一最佳觀測方位下獲取的視覺目標表面點的量化結果。其中 Nn代表攝像機在下一最佳觀測方位下獲取到的視覺目標表面點的個數， No代表被重復觀測到的視覺目標表面點的個數， Nnew表示實際新增點的個數。由表1中可以看出，采用本文所提第1種觀測點計算方法時的結果明顯優(yōu)于采用第2種觀測點計算方法時的結果。

4.2.3 不同下一最佳觀測方位方法間的對比實驗

為了更好地評估本文方法的效果，我們將本文方法與同樣基于深度圖像且考慮了遮擋信息的文獻［9］，文獻［10］中的方法做了對比實驗。圖6給出了部分具有代表性的視覺目標在不同下一最佳觀測方位求解方法下對應的實驗結果。

由圖6可以看出，文獻［9］方法計算出的下一最佳觀測方位側重觀測當前觀測方位下視覺目標的背面。文獻［10］方法計算出的下一最佳觀測方位側重觀測當前觀測方位下信息增益較大點的相鄰未知區(qū)域。本文方法通過充分利用遮擋、輪廓的信息，計算出了一個更加合理的下一觀測方位。

為了進一步衡量不同方法的性能，我們對不同方法的結果進行了量化評估。表2給出了不同方法所求下一最佳觀測方位下的 Nn，No，Nnew，重合率和新增率 Rnew。

表2 不同方法在所求下一最佳觀測方位下實驗結果的量化評估

分析表2可知，就文獻［9］方法而言，因為其主要觀測當前觀測方位下視覺目標的背面，因此在其所求得的下一最佳觀測方位下獲得的表面點重合率較低，新增點個數與當前觀測方位下獲取到的圖像中視覺目標輪廓的大小有直接關系。例如 Bunny，Duck和Mole，因為在當前觀測方位下獲取的視覺目標的輪廓相對較大，所以其獲得的新增點數較多；而對于Dragon和Rocker，因為在當前觀測方位下獲取的視覺目標的輪廓相對較小，所以其獲得的新增點數較少。由此可見，文獻［9］方法具有較大的局限性。就文獻［10］方法而言，其根據已觀測點中信息增益較大點及其法向量計算下一最佳觀測方位。該方法求得的最佳觀測這些已觀測信息增益較大點的觀測方位即為下一最佳觀測方位，這就導致了該方法求得的下一最佳觀測方位與當前觀測方位之間的角度相差不大，因此在該所求觀測方位下獲取到的表面點與當前觀測方位下獲取到的表面點重合率較高，新增點個數相對較少。對于某些特殊的視覺目標例如 Dragon而言，因為其表面比較復雜，所以重合率相對較低。而對于其他視覺目標而言，其重合率都較高。本文方法通過建立綜合考慮視覺目標遮擋及輪廓信息的求解模型，采用優(yōu)化算法對模型求解獲得下一最佳觀測方位。由表2的量化評估結果可以看出，在最能反映方法效果的新增點個數指標上，本文方法結果明顯優(yōu)于文獻［9］和文獻［10］方法的結果，更好地達到了下一最佳觀測方位問題的求解要求。

綜合圖6和表2的結果可知，在相同的實驗條件（如觀測距離）下，本文方法較其他方法獲得的下一最佳觀測方位更為合理。

5 結束語

本文提出一種基于深度圖像利用遮擋和輪廓信息求解下一最佳觀測方位的方法。該方法無需限定攝像機的觀測位置，具有較好的通用性。所提方法的貢獻在于：（1）提出一種綜合考慮遮擋和物體輪廓信息的視覺目標未知區(qū)域建模方法，使得計算出的下一最佳觀測方位更加合理，克服了只利用遮擋或視覺目標背面信息求解下一最佳觀測方位方法的局限性；（2）根據建模所得小三角形的中點、法向量和面積等信息構造了目標函數，使下一最佳觀測方位確定問題轉化為優(yōu)化問題；（3）給出一種根據當前觀測方位下所獲得的視覺目標信息確定優(yōu)化函數初始迭代位置的計算方法，該計算方法對視覺目標具有自適應性，有助于獲得較好的下一最佳觀測方位求解結果。

圖6 不同方法在所求下一最佳觀測方位下的視覺目標深度圖像

［1］ Bottino A and Laurentini A. What is next an interactive next best view approach［J］. Pattern Recognition， 2006， 39（1）: 126-132.

［2］ Pintilie G D and Stuerzlinger W. An evaluation of interactive and automated next best view methods in 3D scanning［J］. Computer-Aided Design and Applications， 2013， 10（2）: 279-291.

［3］ 曾明勇， 吳澤民， 田暢， 等. 基于外觀統計特征融合的人體目標再識別［J］. 電子與信息學報， 2014， 36（8）: 1844-1851. Zeng Ming-yong， Wu Ze-min， Tian Chang， et al.. Fusing appearance statistical features for person re-identification［J］. Journal of Electronics & Information Technology， 2014， 36（8）: 1844-1851.

［4］ 李烈辰， 李道京. 基于壓縮感知的連續(xù)場景稀疏陣列 SAR三維成像［J］. 電子與信息學報， 2014， 36（9）: 2166-2172. Li Lie-Chen and Li Dao-jing. Sparse array SAR 3D imaging for continuous scene based on compressed sensing［J］. Journal of Electronics & Information Technology， 2014， 36（9）: 2166-2172.

［5］ Trummer M， Munkelt C， and Denzler J. Online next-best-view planning for accuracy optimization using an extended E-criterion［C］. Proceedings of the 20th International Conference on Pattern Recognition， Istanbul，Turkey， 2010: 1642-1645.

［6］ Haner S and Heyden A. Covariance propagation and next best view planning for 3D reconstruction［C］. Proceedings of the 12th European Conference on Computer Vision， Florence，Italy， 2012: 545-556.

［7］ Connolly C. The determination of next best views［C］. Proceedings of the IEEE International Conference on Robotics and Automation， Missouri， USA， 1985: 432-435.

［8］ Maver J and Bajcsy R. Occlusion as a guide for planning the next view［J］. IEEE Transactions on Pattern Analysis and Machine Intelligence， 1993， 15（5）: 417-432.

［9］ Banta J E， Wong L M， Dumont C， et al.. A next-best-view system for autonomous 3-D object reconstruction［J］. IEEE Transactions on Systems， Man， and Cybernetics， Part A: Systems and Humans， 2000， 30（5）: 589-598.

［10］ Li Y F and Liu Z G. Information entropy-based viewpoint planning for 3-D object reconstruction［J］. IEEE Transactions on Robotics， 2005， 21（3）: 324-337.

［11］ Vasquez-Gomez J I， Sucar L E， and Murrieta-Cid R. Hierarchical ray tracing for fast volumetric next-best-view planning［C］. Proceedings of the International Conference on Computer and Robot Vision， Regina， SK， Canada， 2013: 181-187.

［12］ Potthast C and Sukhatme G. A probabilistic framework for next best view estimation in a cluttered environment［J］. Journal of Visual Communication and Image Representation，2014， 25（1）: 148-164.

［13］ 張世輝， 劉建新， 孔令富. 基于深度圖像利用隨機森林實現遮擋檢測［J］. 光學學報， 2014， 34（9）: 0915003: 1-12. Zhang Shi-hui， Liu Jian-xin， and Kong Ling-fu. Using random forest for occlusion detection based on depth image［J］. Acta Optica Sinica， 2014， 34（9）: 0915003: 1-12.

［14］ Curtis F and Overton M. A sequential quadratic programming algorithm for nonconvex， nonsmooth constrained optimization［J］. SIAM Journal on Optimization，2012， 22（2）: 474-500.

［15］ Hetzel G， Leibe B， Levi P， et al.. 3D object recognition from range images using local feature histograms［C］. Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition， Piscataway， USA， 2001: II394-II399.

張世輝： 男，1973 年生，博士，教授，博士生導師，研究方向為視覺信息處理、模式識別.

韓德偉： 男，1989 年生，碩士生，研究方向為遮擋檢測與規(guī)避.

何 歡： 男，1989 年生，碩士生，研究方向為計算機視覺.

Determining Next Best View Using Occlusion and Contour Information of Visual Object

Zhang Shi-hui①②Han De-wei①He Huan①
①（School of Information Science and Engineering， Yanshan University， Qinhuangdao 066004， China）②（The Key Laboratory for Computer Virtual Technology and System Integration of Hebei Province，Qinhuangdao 066004， China）

Determining camera’s next best view is a difficult issue in visual field. A next best view approach based on depth image of visual object is proposed by using occlusion and contour information in this paper. Firstly， the occlusion detection is accomplished for the depth image of visual object in current view. Secondly， the unknown regions are constructed according to the occlusion detection result of the depth image and the contour of the visual object， and then the unknown regions are modeled with triangulation-like. Thirdly， the midpoint， normal vector and area of each small triangle and other information are utilized to establish the objective function. Finally， the next best view is obtained by optimizing objective function. Experimental results demonstrate that the approach is feasible and effective.

Depth image； Occlusion； Contour； Unknown regions； Triangulation-like； Next best view

s: The National Natural Science Foundation of China （61379065）； The Natural Science Foundation of Hebei Province （F2014203119）

TP391

A

1009-5896（2015）12-2921-08

10.11999/JEIT150190

2015-02-02；改回日期：2015-08-19；網絡出版：2015-11-01

*通信作者：張世輝 sshhzz@ysu.edu.cn

國家自然科學基金（61379065）和河北省自然科學基金（F2014203119）