2022, Vol. 54
2. 中国科学院东北地理与农业生态研究所, 长春 130012;
3. 哈尔滨师范大学地理科学学院, 哈尔滨 150025;
4. 东北农业大学公共管理与法学院, 哈尔滨 150030
土壤有机质(soil organic matter，SOM)是土壤碳库的重要组成部分，它不仅能够为植物生长提供养分，而且能够改善土壤的物理条件和温度[1-2]，是土壤肥力分析与耕地质量评价的重要指标，对全球碳估算意义重大[3]。农牧交错带是防风固沙的天然屏障，是半干旱区向干旱区的过渡地带，是集农、林、牧于一体的最佳产业区。我国东北农牧交错带分布在内蒙古、黑龙江、吉林、辽宁4省，其生态环境脆弱，土地退化严重，准确地监测SOM对农业可持续发展具有重要意义。
传统SOM化学分析方法费时费力[4]。20世纪70年代起，遥感技术的快速发展为SOM反演提供了新思路。在350 ~ 2 500 nm光学遥感范围，SOM具有独特的光谱响应区域，其与土壤光谱反射率呈显著负相关关系[5]。但农牧交错带区域内SOM含量相对较低，这可能增加SOM遥感预测的不确定性。当前，针对低SOM含量反演主要采用高光谱技术。例如，王延仓等[6]采用室内高光谱数据，利用连续小波变换对含量为14.2 g/kg的北京东部潮土SOM进行反演，其验证集R2最高为0.69；李阳等[7]对采集于南疆地区的荒漠土样采用3种不同建模方法和6种数据处理方式构建SOM高光谱反演模型，其最优模型R2可达0.76。这些研究主要采用室内高光谱数据，无法反映连续区域内的野外监测结果[5]。光学遥感影像数据源众多，特别是多光谱影像，其成像周期短、空间分辨率高，更容易获取，探索其适用于农牧交错带低SOM含量反演，具有重要意义。
目前已有大量研究使用多种方法用于SOM制图，主要有克里格插值模型(Kriging)、多元线性逐步回归模型(MLSR)、随机森林模型(RF)、BP神经网络模型(BPNN)等[4, 8-11]。不同方法具有不同的优势与不足，Kriging方法需要较多样品来揭示SOM空间相关性实现区域预测；MLSR模型通过筛选最优输入量，可以解决输入量间的多元共线问题[12]。BPNN和RF模型已经广泛用于SOM反演，李耀翔等[13]基于近红外光谱与BPNN模型对森林土壤有机碳进行了预测；刘焕军等[14]引入多时相信息，建立基于BP神经网络的SOM遥感反演模型。张欢[15]利用随机RF预测了华东沿海滩涂围垦区土壤有机碳分布；Nabiollahi等[16]利用RF模型对不同土地利用背景下的土壤有机碳进行了预测。
本研究以典型农牧交错带内蒙古自治区赤峰市敖汉旗耕地为研究区，结合Landsat 8 OLI多光谱遥感数据与野外调查采样，提取了波段反射率、波段反射率对数、亮度指数与相关地形因子，分别采用多元线性逐步回归(MLSR)、随机森林(RF)和BP神经网络(BPNN)模型，分别构建农牧交错带SOM多光谱反演模型，为农牧交错带低SOM区域的多光谱遥感反演提供借鉴。
1 材料与方法 1.1 研究区概况敖汉旗位于内蒙古自治区赤峰市东南部(图 1)，地理位置41°69′~ 43°03’N，119°53′ ~ 120°89′E，地处燕山山脉东段努鲁尔虎山北麓、科尔沁沙地南缘，属于典型的半干旱区农牧交错带，温带大陆性季风气候，年降水量374 ~ 524 mm，年平均气温6 ℃，面积约8 316.2 km2，土壤以褐土、栗钙土和风沙土为主。地形和SOM空间差异较大，地势南高北低，平均高程517 m，SOM总体含量较低。耕地以旱地为主，主要种植作物为玉米和谷子。敖汉旗大部分秸秆可以实现机器还田。
|
(该图基于国家测绘地理信息局标准地图服务网站下载的审图号为GS(2016)2923号的标准地图制作，底图无修改) 图 1 研究区及采样点分布图 Fig. 1 Location of study area and sampling points |
2018年5月14日，沿公路两侧采集102个裸土样本，采样深度为0 ~ 20 cm，每个样点设置3次重复，将3个土样混合，作为一个土壤样品，包含的土壤类型主要有褐土、栗钙土、风沙土等。在室内对土壤样品进行研磨、风干和过2 mm筛，采用重铬酸钾容量法测定SOM含量[17]，结果如表 1。2018年研究区SOM均值为10.20 g/kg，处于较低水平，轻度右偏离，基本符合正态分布，变异系数离散程度较高。本文随机选取SOM样本，按照2:1比例建模，其中，建模集68个样本，预测集34个样本(表 1)。建模集与验证集的样本含量分布比较均匀，均值相近，因此，建模集与验证集样本合理，具有代表性。
|
表 1 采样点有机质数据描述统计(n = 102) Table 1 Description statistics of SOM content |
从美国地质调查局(USGS，https://glovis.usgs.gov/)中获取2018年3月25日的Landsat 8 OLI影像，影像在研究区无云覆盖，由于敖汉旗气候干旱，2018年3月1日—4月30日平均降水量为0.30 mm(http://data.cma.cn/)，2018年3月降水数据由澳大利亚插值软件auspline得到(图 2)，空间分辨率为0.1° × 0.1°。且本研究选取的影像处于裸土期，植被覆盖与降水对影像的影响可以忽略不计。SOM随时间短期变化不大，故2018年3月25日的影像可以用于研究区SOM遥感反演研究。影像预处理过程包括采用ENVI 5.3进行辐射定标、大气校正(大气校正后的结果为扩大了10 000倍的反射率)等，经过穗帽变换获取亮度指数。从美国国家航空航天局(https://search.asf.alaska.edu/#/)获取ALOS 12.5m DEM，利用ArcGIS 10.6提取高程(E)、坡度(S)、坡向(A)、坡向变率(sope of aspect，SOA)、坡度变率(slope of slope，SOS)。敖汉旗耕地范围数据来源于中国科学院资源环境科学数据中心数据注册与出版系统(http://www.resdc.cn/DOI)。
|
图 2 2018年3月降水量分布图 Fig. 2 Precipitation in March, 2018 |
本文使用Landsat 8 OLI B1 ~ B7波段反射率、波段反射率对数、亮度指数和地形因子(高程、坡度、坡向、坡度变率、坡向变率)共计20个因子作为建模输入变量(表 2)，按照2:1的比例随机选取68个建模集，根据RF重要性，首先选择4个重要性较大的影像输入变量，然后增加2个地形输入变量以便探究地形因子对模型精度的影响，并分别构建MLSR、RF、BPNN模型。
|
|
表 2 输入变量列表 Table 2 List of input variables |
MLSR通过优化自变量的组合来解释与因变量的线性关系，根据前馈方式输入和删除自变量，确定了输入量满足选择标准时应包括在回归模型中的最佳子集[18]。该模型还避免了自变量的多重共线性。选取SOM的实测值作为因变量，设置变量方差的显著性水平0.05作为进入与剔除的标准，选择B6、E作为自变量，利用SPSS Statistics 22中多元线性逐步回归分析建立预测模型，SOM与自变量的显著性均小于0.001。
1.3.2 RF模型Breiman[19]将多个决策树模型与套袋算法结合提出了RF模型。本文使用R语言平台开发的Random forest包构建模型，主要涉及3个参数：决策树的数量(mtree)，分割节点的随机变量的数量(mtry)，自变量对SOM影响的重要性(importance)[10]。从训练样本中采用随机有放回的抽样方式构建每一个树模型，使用袋外数据误差确定最佳回归树的数量和最优分裂节点数，最终的预测值是所有树模型预测值的均值，本文mtree设置为500，mtry为6[20]。RF模型可处理高维数据，训练速度快，调解参量少，异常值和噪声具有较高的容忍度，且不会过度拟合。
1.3.3 BPNN模型BPNN是一种按误差逆传播算法训练的多层前馈网络，可以完成复杂的非线性计算。它由输入层、隐藏层和输出层组成。在前向传播过程中，输入信号通过隐藏层从输入层处理并传输到输出层。一层神经元的状态只影响下一层的神经元。如果在输出层中无法获得期望的输出，则误差值随后通过网络向后传播，并对每层中的权重进行更改[21]。利用Matlab 2014b建立BP神经网络SOM预测模型，采用3层BP神经网络，选用反向传播算法进行训练，隐含层为双曲正切S型函数，3层输出层节点均为16，输出层为线性传递函数，1个输出层节点，学习效率为0.01，迭代次数为100次。
1.4 精度评价随机选取68个训练样本构建SOM预测模型，使用决定系数(R2)和均方根误差(RMSE)对预测结果进行精度评价[22]。R2越高、RMSE越低，模型预测精度越高。
| $ {R^2} = \frac{{\sum\limits_{i = 1}^n {{{({x_1} - \bar x)}^2}} }}{{\sum\limits_{i = 1}^n {{{({x_2} - \bar x)}^2}} }} $ | (1) |
| $ {\text{RMSE}} = \sqrt {\frac{{\sum\limits_{i = 1}^n {{{({x_1} - {x_2})}^2}} }}{n}} $ | (2) |
式中：x1表示SOM实测值，x2表示SOM预测值，
水分对土壤反射率具有很强的干扰作用，持续的干旱天气，降低了水分对土壤反射率的影响。去除SOM样点的极大值和极小值，选取有代表性的3个样点。从Landsat 8 OLI影像提取不同SOM含量的土壤反射光谱曲线(图 3)，从整体上看，随着SOM含量增加，光谱反射率整体降低，B1 ~ B3波段数值均较小，说明SOM吸收作用强烈，而反射率差异较大的波段可以作为特征波段，B4 ~ B7波段差异明显，因此，Landsat 8 OLI的4 ~ 7波段能较好地反映SOM信息，是敖汉旗SOM遥感反演的特征波段。
|
图 3 Landsat8 OLI不同SOM含量土壤反射率光谱曲线 Fig. 3 Soil spectral reflectance curves with different SOM contents of Landsat 8 OLI |
图 4A为研究区高程图，南部丘陵区高程较高，为1 200 m左右；低高程地区为300 m左右，主要分布在北部平原、沙漠区，南北部高程差约700 m (图 4B)。随着高程的升高，SOM含量逐渐升高后不变(表 3)。图 5为研究区坡度图与不同坡度下的SOM含量图，根据不同坡度将敖汉旗耕地分为0° ~ 3°、3° ~ 4°、4° ~ 5°、5° ~ 6°、> 6°五级(图 5A)。坡度在0° ~ 3°的SOM含量最低，4° ~ 5°的SOM含量最高(图 5B)。结合研究区坡度分析(图 5A)，坡度在0° ~ 3°的地区主要分布在沙漠区和山脚下，因此该区域SOM含量相对较低；而从4° ~ 5°开始，随着坡度的增加，受土壤侵蚀影响，SOM含量降低。
|
图 4 研究区高程分析 Fig. 4 Analysis of DEM of study area |
|
|
表 3 不同高程下的SOM平均值 Table 3 Mean values of SOM in different elevations |
|
图 5 研究区坡度分析 Fig. 5 Analysis of slope of study area |
本文使用Landsat 8 OLI B1 ~ B7波段反射率、波段反射率对数、亮度指数和地形因子作为建模输入量，按照与SOM的重要性排序，结果如表 4。
|
|
表 4 最佳输入量重要性排序 Table 4 Order of importance of optimal input |
选取Landsat 8 OLI lgB4、B5、B6、BI与地形因子(E、SOA)作为输入量，利用MLSR、RF、BPNN构建SOM多光谱反演模型，通过R2、RMSE两个指标评价SOM反演模型的精度(表 5)。
|
|
表 5 MLSR、RF、BPNN模型验证结果 Table 5 Validation results of MLSR model, RF model and BPNN model |
未加入地形因子之前，所有模型中BPNN模型的R2最高，为0.60，RMSE最小，为1.80 g/kg，说明BPNN模型的精度最好。加入地形因子E、SOA后，3种预测模型精度均有所提高，其中BPNN模型精度提高最大，R2为0.82，提高了0.22，RMSE为1.40 g/kg，下降了0.40 g/kg，R2与文献[23]的研究结果相近(R2=0.86)，而RMSE明显优于文献[23] (RMSE=2.29 g/kg)，原因是RMSE与样本数成反比，本文较文献[23]样本数较多，所以RMSE较低。本研究构建的模型能够较好地预测农牧交错带SOM的变化。
2.4.2 有机质反演结果lgB4、B5、B6、BI与E、SOA共同作为输入变量时，BPNN模型反演精度最高(R2：0.82，RMSE：1.40 g/kg)，因此，选择该模型进行SOM制图(图 6)。敖汉旗地区的SOM含量基本低于20 g/kg，北部及其中部地区的SOM含量基本低于10 g/kg，东部地区的SOM含量较高，多集中于15 ~ 20 g/kg。SOM表现出南高北低的特点，其中，高值区均分布在南部高高程地区，低值区分布在北部风沙土区，这与李龙等[24]空间插值的研究结果一致。
|
图 6 2018年SOM遥感反演图 Fig. 6 SOM content inversion map in 2018 |
研究聚焦于农牧交错带区域SOM多光谱反演，以Landsat 8 OLI为数据源，构建适于该地区的SOM反演模型。结果表明，lgB4、B5、B6和BI可作为农牧交错带SOM反演的主要输入变量。SOM对可见光–近红外–短波红外(400 ~ 2 500 nm)比较敏感，已经被广泛用于土壤理化参数反演。王延仓等[25]利用反射光谱与TM5模拟光谱反演北方潮土的SOM，发现室内实测原始光谱与SOM的最大相关系数为0.608 (433 nm)；经吸收峰深度变换后，最大相关系数为0.748(791 nm)。王海峰等[26]探究了以色列南部不同土壤类型下的SOM反演模型，发现砂质土SOM的光谱响应区间为820 ~ 860 nm和940 ~ 970 nm，黏壤土SOM光谱响应区间为730 ~ 790 nm和800 ~ 820 nm。SOM与土壤颜色有关，BI可以表征遥感影像中的土壤颜色和反射率大小，SOM含量越低，反射率越大，BI值越高，图像越亮[27]，对研究区低SOM具有很好的适用性。BPNN和RF模型的反演精度高于MLSR模型(表 5)。这是因为SOM与土壤光谱之间存在着某些线性或非线性关系，而BPNN模型和RF模型具有强大的非线性拟合能力和数据挖掘能力，能更好地模拟SOM与遥感光谱信息及地形因子之间的复杂多元非线性关系。通过反复训练的BPNN模型具有更强的预测效果，可以实现区域SOM更高精度的反演。
敖汉旗SOM空间分布呈南高北低的趋势，与地形和降水有关。南部海拔高的地区有机质含量高，这与Teng等[28]的研究结果一致。原因是南部丘陵地区植被茂盛，SOM容易积累。北部沙漠区气候干旱，土壤质地差，SOM很难积累。但本文中，高程大于700 m后，坡度的增加加剧了土壤侵蚀，水土流失会带走一部分SOM，SOM含量不再持续升高。降水与高程基本呈现一致趋势，降水可以影响区域土壤湿润程度，湿润程度越高，SOM含量越高[14]。而坡度主要影响土壤侵蚀程度，有研究表明[29]，随着坡度的增加，SOM减少。但研究区的SOM随坡度的增大出现先增加后减少的趋势。因为0° ~ 3°主要分布在北部沙漠区和南部耕地(图 5A)，沙漠区环境条件差，坡度不起主导作用，南部耕地和3° ~ 4°部分主要分布在山脚下，SOM随雨水冲刷等作用在山脚沉积，故0° ~ 4°研究区有机质随着坡度的增大出现一定程度的增加。从4° ~ 5°开始，随着坡度的增加，SOM含量降低(图 5B)。坡度大的地区，土壤侵蚀较严重，水土流失加快，造成SOM的损失，SOM含量较低。
以往研究多基于波段反射率构建SOM反演模型，如陈思明等[5]对比全波段与显著波段，构建了基于不同机器学习方法的SOM反演模型，最高R2可达0.74；王琼等[30]通过比较不同波段反射率、波段反射率变换形式得到的回归模型R2最高可达0.72。本研究在对前人SOM反演研究的基础上进行了改进，构建了针对农牧交错带耕地的SOM反演模型。虽然地形因子E和SOA的加入在敖汉旗对SOM预测精度有明显提高，但应用于其他区域是否可行，还需要进一步验证。未来的研究应主要集中于：①多种地形因子在农牧交错带对SOM遥感反演精度提高的作用；②运用最优模型分析农牧交错带SOM的时空变化特征。
4 结论受裸土窗口期较短的影响，针对农牧交错带区域SOM遥感反演研究相对较少。本文以典型农牧交错带内蒙古自治区赤峰市敖汉旗为例，分别利用MLSR、RF、BPNN模型进行SOM反演，BPNN模型的精度最高(R2：0.82，RMSE：1.40 g/kg)，表明BPNN模型在SOM预测中更具有优势，可以得到较高的预测精度。地形因子对SOM反演精度提升明显，增加了SOM预测模型的精度。研究表明，lgB4、B5、B6、BI、E、SOA适用于作为农牧交错带低SOM反演的输入变量。研究成果可以为农牧交错带SOM遥感制图提供借鉴。
| [1] |
Wang J, He T, Lv C, et al. Mapping soil organic matter based on land degradation spectral response units using Hyperion images[J]. International Journal of Applied Earth Observation and Geoinformation, 2010, 12: S171-S180 DOI:10.1016/j.jag.2010.01.002 (  0) |
| [2] |
Wood S A, Tirfessa D, Baudron F. Soil organic matter underlies crop nutritional quality and productivity in smallholder agriculture[J]. Agriculture, Ecosystems&Environment, 2018, 266: 100-108 ( 0) |
| [3] |
Gu X H, Wang Y C, Sun Q, et al. Hyperspectral inversion of soil organic matter content in cultivated land based on wavelet transform[J]. Computers and Electronics in Agriculture, 2019, 167: 105053 DOI:10.1016/j.compag.2019.105053 ( 0) |
| [4] |
杨顺华, 张海涛, 郭龙, 等. 基于回归和地理加权回归Kriging的土壤有机质空间插值[J]. 应用生态学报, 2015, 26(6): 1649-1656 ( 0) |
| [5] |
陈思明, 毛艳玲, 邹小兴, 等. 基于不同建模方法的湿地土壤有机质含量多光谱反演[J]. 土壤通报, 2018, 49(1): 16-22 ( 0) |
| [6] |
王延仓, 张兰, 王欢, 等. 连续小波变换定量反演土壤有机质含量[J]. 光谱学与光谱分析, 2018, 38(11): 3521-3527 ( 0) |
| [7] |
李阳, 刘新路, 彭杰, 等. 基于可见光近红外光谱的南疆荒漠土壤有机质反演研究[J]. 土壤通报, 2018, 49(4): 767-772 ( 0) |
| [8] |
Dou X, Wang X, Liu H J, et al. Prediction of soil organic matter using multi-temporal satellite images in the Songnen Plain, China[J]. Geoderma, 2019, 356: 113896 DOI:10.1016/j.geoderma.2019.113896 ( 0) |
| [9] |
Alloa M, Todoroffa P, Jameuxa M, et al. Alain albrechtc science-soil science and hydrology; researchers at French Agricultural Research Center for International Development (CIRAD) release new data on soil science and hydrology (prediction of tropical volcanic soil organic carbon stocks by visible-near-and mid-infrared spectroscopy)[J]. Agriculture Week, 2020, 189 ( 0) |
| [10] |
刘焕军, 张美薇, 杨昊轩, 等. 多光谱遥感结合随机森林算法反演耕作土壤有机质含量[J]. 农业工程学报, 2020, 36(10): 134-140 DOI:10.11975/j.issn.1002-6819.2020.10.016 ( 0) |
| [11] |
Wang X P, Zhang F, Kung H T, et al. New methods for improving the remote sensing estimation of soil organic matter content (SOMC) in the Ebinur Lake Wetland National Nature Reserve (ELWNNR) in northwest China[J]. Remote Sensing of Environment, 2018, 218: 104-118 DOI:10.1016/j.rse.2018.09.020 ( 0) |
| [12] |
侯艳军, 塔西甫拉提·特依拜, 买买提·沙吾提, 等. 荒漠土壤有机质含量高光谱估算模型[J]. 农业工程学报, 2014, 30(16): 113-120 DOI:10.3969/j.issn.1002-6819.2014.16.016 ( 0) |
| [13] |
李耀翔, 汪洪涛, 耿志伟, 等. 基于近红外光谱及BP神经网络分析法预测森林土壤有机碳含量[J]. 西部林业科学, 2014, 43(3): 1-6 DOI:10.3969/j.issn.1672-8246.2014.03.002 ( 0) |
| [14] |
刘焕军, 潘越, 窦欣, 等. 黑土区田块尺度土壤有机质含量遥感反演模型[J]. 农业工程学报, 2018, 34(1): 127-133 ( 0) |
| [15] |
张欢. 华东沿海滩涂围垦区土壤有机碳动态及其模型预测[D]. 南京: 南京大学, 2017.
( 0) |
| [16] |
Nabiollahi K, Eskandari S, Taghizadeh-Mehrjardi R, et al. Assessing soil organic carbon stocks under land-use change scenarios using random forest models[J]. Carbon Management, 2019, 10(1): 63-77 DOI:10.1080/17583004.2018.1553434 ( 0) |
| [17] |
吴才武, 夏建新, 段峥嵘. 土壤有机质测定方法述评与展望[J]. 土壤, 2015, 47(3): 453-460 ( 0) |
| [18] |
Zare S, Fallah Shamsi S R, Abtahi S A. Weakly-coupled geo-statistical mapping of soil salinity to Stepwise Multiple Linear Regression of MODIS spectral image products[J]. Journal of African Earth Sciences, 2019, 152: 101-114 ( 0) |
| [19] |
Breiman L. Random forests[J]. Machine Learning, 2001, 45(1): 5-32 DOI:10.1023/A:1010933404324 ( 0) |
| [20] |
Liaw A, Wiener M. Classification and regression by randomforest[J]. The Newsletter of the R Project, 2002(2): 18-22 ( 0) |
| [21] |
Xiao H F, Tian Y L. Prediction of mine coal layer spontaneous combustion danger based on genetic algorithm and BP neural networks[J]. Procedia Engineering, 2011, 26: 139-146 DOI:10.1016/j.proeng.2011.11.2151 ( 0) |
| [22] |
Ćalasan M, Abdel Aleem S H E, Zobaa A F. On the root mean square error (RMSE) calculation for parameter estimation of photovoltaic models: A novel exact analytical solution based on Lambert W function[J]. Energy Conversion and Management, 2020, 210: 112716 ( 0) |
| [23] |
马驰. 基于Landsat8吉林中北部地区土壤有机质定量反演研究[J]. 干旱区资源与环境, 2017, 31(2): 167-172 ( 0) |
| [24] |
李龙. 赤峰市敖汉旗土壤有机碳时空变异研究[D]. 呼和浩特: 内蒙古农业大学, 2017.
( 0) |
| [25] |
王延仓, 顾晓鹤, 朱金山, 等. 利用反射光谱及模拟多光谱数据定量反演北方潮土有机质含量[J]. 光谱学与光谱分析, 2014, 34(1): 201-206 DOI:10.3964/j.issn.1000-0593(2014)01-0201-06 ( 0) |
| [26] |
王海峰, 张智韬, Arnon Karnieli, 等. 基于灰度关联-岭回归的荒漠土壤有机质含量高光谱估算[J]. 农业工程学报, 2018, 34(14): 124-131 ( 0) |
| [27] |
Gholizadeh A, Žižala D, Saberioon M, et al. Soil organic carbon and texture retrieving and mapping using proximal, airborne and Sentinel-2 spectral imaging[J]. Remote Sensing of Environment, 2018, 218: 89-103 ( 0) |
| [28] |
Teng M J, Zeng L X, Xiao W F, et al. Spatial variability of soil organic carbon in Three Gorges Reservoir area, China[J]. Science of the Total Environment, 2017, 599/600: 1308-1316 ( 0) |
| [29] |
张枝枝, 张福平, 燕玉超, 等. 渭河两岸缓冲带的土壤有机质含量分布特征及其影响因子[J]. 土壤, 2017, 49(2): 393-399 ( 0) |
| [30] |
王琼, 陈兵, 王方永, 等. 基于HJ卫星的棉田土壤有机质空间分布格局反演[J]. 农业工程学报, 2016, 32(01): 174-180 ( 0) |
2. Northeast Institute of Geography and Agroecology, Chinese Academy of Sciences, Changchun 130012, China;
3. School of Geographical Sciences, Harbin Normal University, Harbin 150025, China;
4. School of Public Administration and Law, Northeast Agricultural University, Harbin 150030, China